5 * Created by Pat Schloss on 7/12/10.
6 * Copyright 2010 Schloss Lab. All rights reserved.
10 #include "qualityscores.h"
12 /**************************************************************************************************/
14 QualityScores::QualityScores(){
16 m = MothurOut::getInstance();
22 m->errorOut(e, "QualityScores", "QualityScores");
27 /**************************************************************************************************/
29 QualityScores::QualityScores(ifstream& qFile){
32 m = MothurOut::getInstance();
41 m->mothurOut("Error reading quality file, name blank at position, " + toString(qFile.tellg()));
42 m->mothurOutEndLine();
45 seqName = seqName.substr(1);
47 cout << seqName << endl;
48 string qScoreString = m->getline(qFile);
49 cout << qScoreString << endl;
50 istringstream qScoreStringStream(qScoreString);
52 while(!qScoreStringStream.eof()){
53 if (m->control_pressed) { break; }
54 qScoreStringStream >> score;
55 qScores.push_back(score);
56 cout << score << '\t' << count << endl;
67 if (seqName.length() != 0) {
68 seqName = seqName.substr(1);
69 while (!qFile.eof()) {
72 if (c == 10 || c == 13){ break; }
79 char letter=qFile.get();
80 if((letter == '>')){ qFile.putback(letter); break; }
81 else if (isprint(letter)) { scores += letter; }
89 //convert scores string to qScores
90 istringstream qScoreStringStream(scores);
93 while(!qScoreStringStream.eof()){
95 if (m->control_pressed) { break; }
97 qScoreStringStream >> score;
98 qScores.push_back(score);
102 seqLength = qScores.size();
106 catch(exception& e) {
107 m->errorOut(e, "QualityScores", "QualityScores");
113 /**************************************************************************************************/
115 string QualityScores::getName(){
120 catch(exception& e) {
121 m->errorOut(e, "QualityScores", "getName");
126 /**************************************************************************************************/
128 void QualityScores::printQScores(ofstream& qFile){
131 double aveQScore = calculateAverage();
133 qFile << '>' << seqName << '\t' << aveQScore << endl;
135 for(int i=0;i<seqLength;i++){
136 qFile << qScores[i] << ' ';
141 catch(exception& e) {
142 m->errorOut(e, "QualityScores", "printQScores");
147 /**************************************************************************************************/
149 void QualityScores::trimQScores(int start, int end){
154 hold = vector<int>(qScores.begin()+start, qScores.end());
158 if(qScores.size() > end){
159 hold = vector<int>(qScores.begin(), qScores.begin()+end);
164 seqLength = qScores.size();
166 catch(exception& e) {
167 m->errorOut(e, "QualityScores", "trimQScores");
172 /**************************************************************************************************/
174 void QualityScores::flipQScores(){
177 vector<int> temp = qScores;
178 for(int i=0;i<seqLength;i++){
179 qScores[seqLength - i - 1] = temp[i];
183 catch(exception& e) {
184 m->errorOut(e, "QualityScores", "flipQScores");
189 /**************************************************************************************************/
191 bool QualityScores::stripQualThreshold(Sequence& sequence, double qThreshold){
193 string rawSequence = sequence.getUnaligned();
194 int seqLength = sequence.getNumBases();
196 if(seqName != sequence.getName()){
197 m->mothurOut("sequence name mismatch btwn fasta: " + sequence.getName() + " and qual file: " + seqName);
198 m->mothurOutEndLine();
202 for(int i=0;i<seqLength;i++){
204 if(qScores[i] < qThreshold){
210 if (end == (seqLength-1)) { end = seqLength; }
212 sequence.setUnaligned(rawSequence.substr(0,end));
213 trimQScores(-1, end);
217 catch(exception& e) {
218 m->errorOut(e, "QualityScores", "flipQScores");
224 /**************************************************************************************************/
226 bool QualityScores::stripQualRollingAverage(Sequence& sequence, double qThreshold){
228 string rawSequence = sequence.getUnaligned();
229 int seqLength = sequence.getNumBases();
231 if(seqName != sequence.getName()){
232 m->mothurOut("sequence name mismatch btwn fasta: " + sequence.getName() + " and qual file: " + seqName);
233 m->mothurOutEndLine();
237 double rollingSum = 0.0000;
239 for(int i=0;i<seqLength;i++){
241 rollingSum += (double)qScores[i];
243 if(rollingSum / (double)(i+1) < qThreshold){
249 if(end == -1){ end = seqLength; }
252 sequence.setUnaligned(rawSequence.substr(0,end));
253 trimQScores(-1, end);
258 catch(exception& e) {
259 m->errorOut(e, "QualityScores", "flipQScores");
265 /**************************************************************************************************/
267 bool QualityScores::stripQualWindowAverage(Sequence& sequence, int stepSize, int windowSize, double qThreshold){
269 string rawSequence = sequence.getUnaligned();
270 int seqLength = sequence.getNumBases();
272 if(seqName != sequence.getName()){
273 m->mothurOut("sequence name mismatch between fasta: " + sequence.getName() + " and qual file: " + seqName);
274 m->mothurOutEndLine();
277 int end = windowSize;
280 if(seqLength < windowSize) { return 0; }
282 while(start < seqLength){
283 double windowSum = 0.0000;
285 for(int i=start;i<end;i++){
286 windowSum += qScores[i];
288 double windowAverage = windowSum / (double)(end-start);
290 if(windowAverage < qThreshold){
291 end = end - stepSize;
295 end = start + windowSize;
296 if(end >= seqLength){ end = seqLength - 1; }
299 if(end == -1){ end = seqLength; }
301 sequence.setUnaligned(rawSequence.substr(0,end));
302 trimQScores(-1, end);
306 catch(exception& e) {
307 m->errorOut(e, "QualityScores", "flipQScores");
313 /**************************************************************************************************/
315 double QualityScores::calculateAverage(){
317 double aveQScore = 0.0000;
319 for(int i=0;i<seqLength;i++){
320 aveQScore += (double) qScores[i];
322 aveQScore /= (double) seqLength;
327 /**************************************************************************************************/
329 bool QualityScores::cullQualAverage(Sequence& sequence, double qAverage){
331 string rawSequence = sequence.getUnaligned();
332 bool success = 0; //guilty until proven innocent
334 if(seqName != sequence.getName()) {
335 m->mothurOut("sequence name mismatch btwn fasta: " + sequence.getName() + " and qual file: " + seqName);
336 m->mothurOutEndLine();
339 double aveQScore = calculateAverage();
341 if(aveQScore >= qAverage) { success = 1; }
342 else { success = 0; }
346 catch(exception& e) {
347 m->errorOut(e, "QualityScores", "cullQualAverage");
352 /**************************************************************************************************/
354 void QualityScores::updateQScoreErrorMap(map<char, vector<int> >& qualErrorMap, string errorSeq, int start, int stop, int weight){
357 int seqLength = errorSeq.size();
359 int qIndex = start - 1;
361 for(int i=0;i<seqLength;i++){
363 if(errorSeq[i] == 'm') { qualErrorMap['m'][qScores[qIndex]] += weight; }
364 else if(errorSeq[i] == 's') { qualErrorMap['s'][qScores[qIndex]] += weight; }
365 else if(errorSeq[i] == 'i') { qualErrorMap['i'][qScores[qIndex]] += weight; }
366 else if(errorSeq[i] == 'a') { qualErrorMap['a'][qScores[qIndex]] += weight; }
367 else if(errorSeq[i] == 'd') { /* there are no qScores for deletions */ }
369 if(errorSeq[i] != 'd') { qIndex++; }
370 if(qIndex > stop){ break; }
373 catch(exception& e) {
374 m->errorOut(e, "QualityScores", "updateQScoreErrorMap");
379 /**************************************************************************************************/
381 void QualityScores::updateForwardMap(vector<vector<int> >& forwardMap, int start, int stop, int weight){
385 for(int i=start-1;i<stop;i++){
386 forwardMap[index++][qScores[i]] += weight;
390 catch(exception& e) {
391 m->errorOut(e, "QualityScores", "updateForwardMap");
396 /**************************************************************************************************/
398 void QualityScores::updateReverseMap(vector<vector<int> >& reverseMap, int start, int stop, int weight){
402 for(int i=stop-1;i>=start;i--){
403 reverseMap[index++][qScores[i]] += weight;
407 catch(exception& e) {
408 m->errorOut(e, "QualityScores", "updateForwardMap");
413 /**************************************************************************************************/