Imported Debian patch 0.12.7-3
[bowtie.git] / pat.cpp
1 #include <cmath>
2 #include <iostream>
3 #include <string>
4 #include <stdexcept>
5 #include <seqan/sequence.h>
6 #include <seqan/file.h>
7
8 #include "pat.h"
9 #include "filebuf.h"
10
11 using namespace std;
12 using namespace seqan;
13
14 /**
15  * Parse a single quality string from fb and store qualities in r.
16  * Assume the next character obtained via fb.get() is the first
17  * character of the quality string.  When returning, the next
18  * character returned by fb.peek() or fb.get() should be the first
19  * character of the following line.
20  */
21 int parseQuals(ReadBuf& r,
22                FileBuf& fb,
23                int readLen,
24                int trim3,
25                int trim5,
26                bool intQuals,
27                bool phred64,
28                bool solexa64)
29 {
30         int qualsRead = 0;
31         int c = 0;
32         assert(fb.peek() != '\n' && fb.peek() != '\r');
33         _setBegin (r.qual, (char*)r.qualBuf);
34         _setLength(r.qual, 0);
35         if (intQuals) {
36                 while (c != '\r' && c != '\n' && c != -1) {
37                         bool neg = false;
38                         int num = 0;
39                         while(!isspace(c = fb.peek()) && !fb.eof()) {
40                                 if(c == '-') {
41                                         neg = true;
42                                         assert_eq(num, 0);
43                                 } else {
44                                         if(!isdigit(c)) {
45                                                 char buf[2048];
46                                                 cerr << "Warning: could not parse quality line:" << endl;
47                                                 fb.getPastNewline();
48                                                 cerr << fb.copyLastN(buf);
49                                                 buf[2047] = '\0';
50                                                 cerr << buf;
51                                                 throw 1;
52                                         }
53                                         assert(isdigit(c));
54                                         num *= 10;
55                                         num += (c - '0');
56                                 }
57                                 fb.get();
58                         }
59                         if(neg) num = 0;
60                         // Phred-33 ASCII encode it and add it to the back of the
61                         // quality string
62                         r.qualBuf[qualsRead++] = ('!' + num);
63                         // Skip over next stretch of whitespace
64                         c = fb.peek();
65                         while(c != '\r' && c != '\n' && isspace(c) && !fb.eof()) {
66                                 fb.get();
67                                 c = fb.peek();
68                         }
69                 }
70         } else {
71                 while (c != '\r' && c != '\n' && c != -1) {
72                         c = fb.get();
73                         r.qualBuf[qualsRead++] = charToPhred33(c, solexa64, phred64);
74                         c = fb.peek();
75                         while(c != '\r' && c != '\n' && isspace(c) && !fb.eof()) {
76                                 fb.get();
77                                 c = fb.peek();
78                         }
79                 }
80         }
81         if (qualsRead < readLen-1 ||
82             (qualsRead < readLen && !r.color))
83         {
84                 tooFewQualities(r.name);
85         }
86         qualsRead -= trim3;
87         if(qualsRead <= 0) return 0;
88         int trimmedReadLen = readLen-trim3-trim5;
89         if(trimmedReadLen < 0) trimmedReadLen = 0;
90         if(qualsRead > trimmedReadLen) {
91                 // Shift everybody left
92                 for(int i = 0; i < readLen; i++) {
93                         r.qualBuf[i] = r.qualBuf[i+qualsRead-trimmedReadLen];
94                 }
95         }
96         _setLength(r.qual, trimmedReadLen);
97         while(fb.peek() == '\n' || fb.peek() == '\r') fb.get();
98         return qualsRead;
99 }
100
101 void wrongQualityFormat(const String<char>& read_name) {
102         cerr << "Encountered a space parsing the quality string for read " << read_name << endl
103              << "If this is a FASTQ file with integer (non-ASCII-encoded) qualities, please" << endl
104              << "re-run Bowtie with the --integer-quals option.  If this is a FASTQ file with" << endl
105              << "alternate basecall information, please re-run Bowtie with the --fuzzy option." << endl;
106         throw 1;
107 }
108
109 void tooFewQualities(const String<char>& read_name) {
110         cerr << "Too few quality values for read: " << read_name << endl
111                  << "\tare you sure this is a FASTQ-int file?" << endl;
112         throw 1;
113 }
114
115 void tooManyQualities(const String<char>& read_name) {
116         cerr << "Reads file contained a pattern with more than 1024 quality values." << endl
117                  << "Please truncate reads and quality values and and re-run Bowtie" << endl;
118         throw 1;
119 }
120
121 void tooManySeqChars(const String<char>& read_name) {
122         cerr << "Reads file contained a pattern with more than 1024 sequence characters." << endl
123                  << "Please truncate reads and quality values and and re-run Bowtie." << endl
124                  << "Offending read: " << read_name << endl;
125         throw 1;
126 }
127
128 /**
129  * C++ version char* style "itoa":
130  */
131 char* itoa10(int value, char* result) {
132         // Check that base is valid
133         char* out = result;
134         int quotient = value;
135         do {
136                 *out = "0123456789"[ std::abs( quotient % 10 ) ];
137                 ++out;
138                 quotient /= 10;
139         } while ( quotient );
140
141         // Only apply negative sign for base 10
142         if (value < 0) *out++ = '-';
143         std::reverse( result, out );
144
145         *out = 0; // terminator
146         return out;
147 }