← Back to branch summary

~ubuntu-branches/ubuntu/jaunty/amap-align/jaunty

« back to all changes in this revision

Viewing changes to align/Sequence.h

  • Committer: Bazaar Package Importer
  • Author(s): Charles Plessy, Charles Plessy
  • Date: 2008-04-10 09:17:01 UTC
  • mfrom: (1.1.1 upstream) (2.1.1 gutsy)
  • Revision ID: james.westby@ubuntu.com-20080410091701-ioqvofeuaify075g
Tags: 2.2-1
http://bugs.debian.org/468060
[ Charles Plessy ]
* New upstream release:
  - Corrected DNA alignment parameters file. Improved debug output. (2.01)
  - Includes new output for the AMAP Display Java based GUI. (2.1)
* debian/amap.1, debian/amap.1.xml: updated the manpage.
* debian/watch: changed URI to monitor.
* debian/patches, debian/rules, debian/control: switched to quilt.
* debian/patches/fix-gcc-4.3.diff: updated for version 2.2 (Closes: #468060).
* debian/control:
  - Moved the 'Conflicts:' field from the source to the binary section.

Show diffs side-by-side

added added

removed removed

Lines of Context:
align/Sequence.h
 
1
/////////////////////////////////////////////////////////////////
 
2
// Sequence.h
 
3
//
 
4
// Class for reading/manipulating single sequence character data.
 
5
/////////////////////////////////////////////////////////////////
 
6
 
 
7
#ifndef SEQUENCE_H
 
8
#define SEQUENCE_H
 
9
 
 
10
#include <string>
 
11
#include <fstream>
 
12
#include <iostream>
 
13
#include <cctype>
 
14
#include <cstdlib>
 
15
#include "SafeVector.h"
 
16
#include "FileBuffer.h"
 
17
 
 
18
/////////////////////////////////////////////////////////////////
 
19
// Sequence
 
20
//
 
21
// Class for storing sequence information.
 
22
/////////////////////////////////////////////////////////////////
 
23
 
 
24
class Sequence {
 
25
 
 
26
  bool isValid;                // a boolean indicating whether the sequence data is valid or not
 
27
  string header;               // string containing the comment line of the FASTA file
 
28
  SafeVector<char> *data;      // pointer to character data
 
29
  int length;                  // length of the sequence
 
30
  int sequenceLabel;           // integer sequence label, typically to indicate the ordering of sequences
 
31
                               //   in a Multi-FASTA file
 
32
  int inputLabel;              // position of sequence in original input
 
33
 
 
34
  /////////////////////////////////////////////////////////////////
 
35
  // Sequence::Sequence()
 
36
  //
 
37
  // Default constructor.  Does nothing.
 
38
  /////////////////////////////////////////////////////////////////
 
39
 
 
40
  Sequence () : isValid (false), header (""), data (NULL), length (0), sequenceLabel (0), inputLabel (0) {}
 
41
 
 
42
 public:
 
43
 
 
44
  /////////////////////////////////////////////////////////////////
 
45
  // Sequence::Sequence()
 
46
  //
 
47
  // Constructor.  Reads the sequence from a FileBuffer.
 
48
  /////////////////////////////////////////////////////////////////
 
49
 
 
50
  Sequence (FileBuffer &infile, bool stripGaps = false) : isValid (false), header ("~"), data (NULL), length(0), sequenceLabel (0), inputLabel (0) {
 
51
 
 
52
    // read until the first non-blank line
 
53
    while (!infile.eof()){
 
54
      infile.GetLine (header);
 
55
      if (header.length() != 0) break;
 
56
    }
 
57
 
 
58
    // check to make sure that it is a correct header line
 
59
    if (header[0] == '>'){
 
60
 
 
61
      // if so, remove the leading ">"
 
62
      header = header.substr (1);
 
63
 
 
64
      // remove any leading or trailing white space in the header comment
 
65
      while (header.length() > 0 && isspace (header[0])) header = header.substr (1);
 
66
      while (header.length() > 0 && isspace (header[header.length() - 1])) header = header.substr(0, header.length() - 1);
 
67
 
 
68
      // get ready to read the data[] array; note that data[0] is always '@'
 
69
      char ch;
 
70
      data = new SafeVector<char>; assert (data);
 
71
      data->push_back ('@');
 
72
 
 
73
      // get a character from the file
 
74
      while (infile.Get(ch)){
 
75
 
 
76
        // if we've reached a new comment line, put the character back and stop
 
77
        if (ch == '>'){ infile.UnGet(); break; }
 
78
 
 
79
        // skip whitespace
 
80
        if (isspace (ch)) continue;
 
81
 
 
82
        // substitute gap character
 
83
        if (ch == '.') ch = '-';
 
84
        if (stripGaps && ch == '-') continue;
 
85
 
 
86
        // check for known characters
 
87
        if (!((ch >= 'A' && ch <= 'Z') || (ch >= 'a' && ch <= 'z') || ch == '*' || ch == '-')){
 
88
          cerr << "ERROR: Unknown character encountered: " << ch << endl;
 
89
          exit (1);
 
90
        }
 
91
 
 
92
        // everything's ok so far, so just store this character.
 
93
        data->push_back(ch);
 
94
        ++length;
 
95
      }
 
96
 
 
97
      // sequence must contain data in order to be valid
 
98
      isValid = length > 0;
 
99
      if (!isValid){
 
100
        delete data;
 
101
        data = NULL;
 
102
      }
 
103
    }
 
104
  }
 
105
 
 
106
  
 
107
  /////////////////////////////////////////////////////////////////
 
108
  // Sequence::Sequence()
 
109
  //
 
110
  // Constructor.  Builds a sequence from existing data.  Note
 
111
  // that the data must use one-based indexing where data[0] should
 
112
  // be set to '@'.
 
113
  /////////////////////////////////////////////////////////////////
 
114
 
 
115
  Sequence (SafeVector<char> *data, string header, int length, int sequenceLabel, int inputLabel) :
 
116
    isValid (data != NULL), header(header), data(data), length (length), sequenceLabel (sequenceLabel), inputLabel (inputLabel) {
 
117
      assert (data);
 
118
      assert ((*data)[0] == '@');
 
119
  }
 
120
 
 
121
  /////////////////////////////////////////////////////////////////
 
122
  // Sequence::Sequence()
 
123
  //
 
124
  // Destructor.  Release allocated memory.
 
125
  /////////////////////////////////////////////////////////////////
 
126
 
 
127
  ~Sequence (){
 
128
    if (data){
 
129
      assert (isValid);
 
130
      delete data;
 
131
      data = NULL;
 
132
      isValid = false;
 
133
    }
 
134
  }
 
135
 
 
136
  /////////////////////////////////////////////////////////////////
 
137
  // Sequence::GetHeader()
 
138
  //
 
139
  // Return the string comment associated with this sequence.
 
140
  /////////////////////////////////////////////////////////////////
 
141
 
 
142
  string GetHeader () const {
 
143
    return header;
 
144
  }
 
145
 
 
146
  /////////////////////////////////////////////////////////////////
 
147
  // Sequence::GetName()
 
148
  //
 
149
  // Return the first word of the string comment associated with this sequence.
 
150
  /////////////////////////////////////////////////////////////////
 
151
 
 
152
  string GetName () const {
 
153
    char name[1024];
 
154
    sscanf (header.c_str(), "%s", name);
 
155
    return string(name);
 
156
  }
 
157
 
 
158
  /////////////////////////////////////////////////////////////////
 
159
  // Sequence::GetDataPtr()
 
160
  //
 
161
  // Return the iterator to data associated with this sequence.
 
162
  /////////////////////////////////////////////////////////////////
 
163
 
 
164
  SafeVector<char>::iterator GetDataPtr(){
 
165
    assert (isValid);
 
166
    assert (data);
 
167
    return data->begin();
 
168
  }
 
169
 
 
170
  /////////////////////////////////////////////////////////////////
 
171
  // Sequence::GetPosition()
 
172
  //
 
173
  // Return the character at position i.  Recall that the character
 
174
  // data is stored with one-based indexing.
 
175
  /////////////////////////////////////////////////////////////////
 
176
 
 
177
  char GetPosition (int i) const {
 
178
    assert (isValid);
 
179
    assert (data);
 
180
    assert (i >= 1 && i <= length);
 
181
    return (*data)[i];
 
182
  }
 
183
 
 
184
  /////////////////////////////////////////////////////////////////
 
185
  // Sequence::SetLabel()
 
186
  //
 
187
  // Sets the sequence label to i.
 
188
  /////////////////////////////////////////////////////////////////
 
189
 
 
190
  void SetLabel (int i){
 
191
    assert (isValid);
 
192
    sequenceLabel = i;
 
193
    inputLabel = i;
 
194
  }
 
195
 
 
196
  /////////////////////////////////////////////////////////////////
 
197
  // Sequence::SetSortLabel()
 
198
  //
 
199
  // Sets the sequence sorting label to i.
 
200
  /////////////////////////////////////////////////////////////////
 
201
 
 
202
  void SetSortLabel (int i){
 
203
    assert (isValid);
 
204
    sequenceLabel = i;
 
205
  }
 
206
 
 
207
  /////////////////////////////////////////////////////////////////
 
208
  // Sequence::GetLabel()
 
209
  //
 
210
  // Retrieves the input label.
 
211
  /////////////////////////////////////////////////////////////////
 
212
 
 
213
  int GetLabel () const {
 
214
    assert (isValid);
 
215
    return inputLabel;
 
216
  }
 
217
 
 
218
  /////////////////////////////////////////////////////////////////
 
219
  // Sequence::GetSortLabel()
 
220
  //
 
221
  // Retrieves the sorting label.
 
222
  /////////////////////////////////////////////////////////////////
 
223
 
 
224
  int GetSortLabel () const {
 
225
    assert (isValid);
 
226
    return sequenceLabel;
 
227
  }
 
228
 
 
229
  /////////////////////////////////////////////////////////////////
 
230
  // Sequence::Fail()
 
231
  //
 
232
  // Checks to see if the sequence successfully loaded.
 
233
  /////////////////////////////////////////////////////////////////
 
234
 
 
235
  bool Fail () const {
 
236
    return !isValid;
 
237
  }
 
238
 
 
239
  /////////////////////////////////////////////////////////////////
 
240
  // Sequence::Length()
 
241
  //
 
242
  // Returns the length of the sequence.
 
243
  /////////////////////////////////////////////////////////////////
 
244
 
 
245
  int GetLength () const {
 
246
    assert (isValid);
 
247
    assert (data);
 
248
    return length;
 
249
  }
 
250
 
 
251
  /////////////////////////////////////////////////////////////////
 
252
  // Sequence::WriteMFA()
 
253
  //
 
254
  // Writes the sequence to outfile in MFA format.  Uses numColumns
 
255
  // columns per line.  If useIndex is set to false, then the
 
256
  // header is printed as normal, but if useIndex is true, then
 
257
  // ">S###" is printed where ### represents the sequence label.
 
258
  /////////////////////////////////////////////////////////////////
 
259
 
 
260
  void WriteMFA (ostream &outfile, int numColumns, bool useIndex = false) const {
 
261
    assert (isValid);
 
262
    assert (data);
 
263
    assert (!outfile.fail());
 
264
 
 
265
    // print out heading
 
266
    if (useIndex)
 
267
      outfile << ">S" << GetLabel() << endl;
 
268
    else
 
269
      outfile << ">" << header << endl;
 
270
 
 
271
    // print out character data
 
272
    int ct = 1;
 
273
    for (; ct <= length; ct++){
 
274
      outfile << (*data)[ct];
 
275
      if (ct % numColumns == 0) outfile << endl;
 
276
    }
 
277
    if ((ct-1) % numColumns != 0) outfile << endl;
 
278
  }
 
279
 
 
280
  /////////////////////////////////////////////////////////////////
 
281
  // Sequence::Clone()
 
282
  //
 
283
  // Returns a new deep copy of the seqeuence.
 
284
  /////////////////////////////////////////////////////////////////
 
285
 
 
286
  Sequence *Clone () const {
 
287
    Sequence *ret = new Sequence();
 
288
    assert (ret);
 
289
 
 
290
    ret->isValid = isValid;
 
291
    ret->header = header;
 
292
    ret->data = new SafeVector<char>; assert (ret->data);
 
293
    *(ret->data) = *data;
 
294
    ret->length = length;
 
295
    ret->sequenceLabel = sequenceLabel;
 
296
    ret->inputLabel = inputLabel;
 
297
 
 
298
    return ret;
 
299
  }
 
300
 
 
301
  /////////////////////////////////////////////////////////////////
 
302
  // Sequence::GetRange()
 
303
  //
 
304
  // Returns a new sequence object consisting of a range of
 
305
  // characters from the current seuquence.
 
306
  /////////////////////////////////////////////////////////////////
 
307
 
 
308
  Sequence *GetRange (int start, int end) const {
 
309
    Sequence *ret = new Sequence();
 
310
    assert (ret);
 
311
 
 
312
    assert (start >= 1 && start <= length);
 
313
    assert (end >= 1 && end <= length);
 
314
    assert (start <= end);
 
315
 
 
316
    ret->isValid = isValid;
 
317
    ret->header = header;
 
318
    ret->data = new SafeVector<char>; assert (ret->data);
 
319
    ret->data->push_back ('@');
 
320
    for (int i = start; i <= end; i++)
 
321
      ret->data->push_back ((*data)[i]);
 
322
    ret->length = end - start + 1;
 
323
    ret->sequenceLabel = sequenceLabel;
 
324
    ret->inputLabel = inputLabel;
 
325
 
 
326
    return ret;
 
327
  }
 
328
 
 
329
  /////////////////////////////////////////////////////////////////
 
330
  // Sequence::AddGaps()
 
331
  //
 
332
  // Given an SafeVector<char> containing the skeleton for an
 
333
  // alignment and the identity of the current character, this
 
334
  // routine will create a new sequence with all necesssary gaps added.
 
335
  // For instance,
 
336
  //    alignment = "XXXBBYYYBBYYXX"
 
337
  //    id = 'X'
 
338
  // will perform the transformation
 
339
  //    "ATGCAGTCA" --> "ATGCC---GT--CA"
 
340
  //                    (XXXBBYYYBBYYXX)
 
341
  /////////////////////////////////////////////////////////////////
 
342
 
 
343
  Sequence *AddGaps (SafeVector<char> *alignment, char id){
 
344
    Sequence *ret = new Sequence();
 
345
    assert (ret);
 
346
 
 
347
    ret->isValid = isValid;
 
348
    ret->header = header;
 
349
    ret->data = new SafeVector<char>; assert (ret->data);
 
350
    ret->length = (int) alignment->size();
 
351
    ret->sequenceLabel = sequenceLabel;
 
352
    ret->inputLabel = inputLabel;
 
353
    ret->data->push_back ('@');
 
354
 
 
355
    SafeVector<char>::iterator dataIter = data->begin() + 1;
 
356
    for (SafeVector<char>::iterator iter = alignment->begin(); iter != alignment->end(); ++iter){
 
357
      if (*iter == 'B' || *iter == id){
 
358
        ret->data->push_back (*dataIter);
 
359
        ++dataIter;
 
360
      }
 
361
      else
 
362
        ret->data->push_back ('-');
 
363
    }
 
364
 
 
365
    return ret;
 
366
  }
 
367
 
 
368
  /////////////////////////////////////////////////////////////////
 
369
  // Sequence::GetString()
 
370
  //
 
371
  // Returns the sequence as a string with gaps removed.
 
372
  /////////////////////////////////////////////////////////////////
 
373
 
 
374
  string GetString (){
 
375
    string s = "";
 
376
    for (int i = 1; i <= length; i++){
 
377
      if ((*data)[i] != '-') s += (*data)[i];
 
378
    }
 
379
    return s;
 
380
  }
 
381
 
 
382
 
 
383
  /////////////////////////////////////////////////////////////////
 
384
  // Sequence::GetMapping()
 
385
  //
 
386
  // Returns a SafeVector<int> containing the indices of every
 
387
  // character in the sequence.  For instance, if the data is
 
388
  // "ATGCC---GT--CA", the method returns {1,2,3,4,5,9,10,13,14}.
 
389
  /////////////////////////////////////////////////////////////////
 
390
 
 
391
  SafeVector<int> *GetMapping () const {
 
392
    SafeVector<int> *ret = new SafeVector<int>(1, 0);
 
393
    for (int i = 1; i <= length; i++){
 
394
      if ((*data)[i] != '-') ret->push_back (i);
 
395
    }
 
396
    return ret;
 
397
  }
 
398
 
 
399
  /////////////////////////////////////////////////////////////////
 
400
  // Sequence::Highlight()
 
401
  //
 
402
  // Changes all positions with score >= cutoff to upper case and
 
403
  // all positions with score < cutoff to lower case.
 
404
  /////////////////////////////////////////////////////////////////
 
405
 
 
406
  void Highlight (const SafeVector<float> &scores, const float cutoff){
 
407
    for (int i = 1; i <= length; i++){
 
408
      if (scores[i-1] >= cutoff)
 
409
        (*data)[i] = toupper ((*data)[i]);
 
410
      else
 
411
        (*data)[i] = tolower ((*data)[i]);
 
412
    }
 
413
  }
 
414
};
 
415
 
 
416
#endif