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();
21 m->errorOut(e, "QualityScores", "QualityScores");
26 /**************************************************************************************************/
28 QualityScores::QualityScores(ifstream& qFile, int l){
31 m = MothurOut::getInstance();
38 //m->getline(qFile, line);
39 //istringstream nameStream(line);
42 while (!qFile.eof()) { char c = qFile.get(); if (c == 10 || c == 13 || c == -1){ break; } } // get rest of line
44 if (seqName == "") { m->mothurOut("Error reading quality file, name blank at position, " + toString(qFile.tellg())); m->mothurOutEndLine(); }
46 seqName = seqName.substr(1);
49 //m->getline(qFile, line);
50 //istringstream qualStream(line);
53 // qualStream >> score;
54 // qScores.push_back(score);
58 //seqLength = qScores.size();
66 //cout << name << endl;
67 if (name.length() != 0) {
68 saveName = name.substr(1);
71 if (c == 10 || c == 13){ break; }
78 char letter= in.get();
79 if(letter == '>'){ in.putback(letter); break; }
80 else{ scores += letter; }
83 //istringstream iss (scores,istringstream::in);
85 //int count = 0; int tempScore;
86 //while (iss) { iss >> tempScore; count++; }
87 //cout << saveName << '\t' << count << endl;
94 for(int i=0;i<seqLength;i++){
96 qScores.push_back(score);
101 catch(exception& e) {
102 m->errorOut(e, "QualityScores", "QualityScores");
107 /**************************************************************************************************/
109 /**************************************************************************************************/
111 string QualityScores::getName(){
116 catch(exception& e) {
117 m->errorOut(e, "QualityScores", "getName");
122 /**************************************************************************************************/
124 void QualityScores::printQScores(ofstream& qFile){
127 double aveQScore = calculateAverage();
129 qFile << '>' << seqName << '\t' << aveQScore << endl;
131 for(int i=0;i<seqLength;i++){
132 qFile << qScores[i] << ' ';
137 catch(exception& e) {
138 m->errorOut(e, "QualityScores", "printQScores");
143 /**************************************************************************************************/
145 void QualityScores::trimQScores(int start, int end){
150 hold = vector<int>(qScores.begin()+start, qScores.end());
154 hold = vector<int>(qScores.begin(), qScores.begin()+end); //not sure if indexing is correct
158 seqLength = qScores.size();
160 catch(exception& e) {
161 m->errorOut(e, "QualityScores", "trimQScores");
166 /**************************************************************************************************/
168 void QualityScores::flipQScores(){
171 vector<int> temp = qScores;
172 for(int i=0;i<seqLength;i++){
173 qScores[seqLength - i - 1] = temp[i];
177 catch(exception& e) {
178 m->errorOut(e, "QualityScores", "flipQScores");
183 /**************************************************************************************************/
185 bool QualityScores::stripQualThreshold(Sequence& sequence, double qThreshold){
187 string rawSequence = sequence.getUnaligned();
188 int seqLength = sequence.getNumBases();
190 if(seqName != sequence.getName()){
191 m->mothurOut("sequence name mismatch btwn fasta: " + sequence.getName() + " and qual file: " + seqName);
192 m->mothurOutEndLine();
196 for(int i=0;i<seqLength;i++){
198 if(qScores[i] < qThreshold){
203 sequence.setUnaligned(rawSequence.substr(0,end));
204 trimQScores(-1, end);
208 catch(exception& e) {
209 m->errorOut(e, "QualityScores", "flipQScores");
215 /**************************************************************************************************/
217 bool QualityScores::stripQualRollingAverage(Sequence& sequence, double qThreshold){
219 string rawSequence = sequence.getUnaligned();
220 int seqLength = sequence.getNumBases();
222 if(seqName != sequence.getName()){
223 m->mothurOut("sequence name mismatch btwn fasta: " + sequence.getName() + " and qual file: " + seqName);
224 m->mothurOutEndLine();
228 double rollingSum = 0.0000;
230 for(int i=0;i<seqLength;i++){
232 rollingSum += (double)qScores[i];
234 if(rollingSum / (double)(i+1) < qThreshold){
236 // cout << i+1 << '\t' << seqName << '\t' << rollingSum / (double)(i+1) << endl;
241 if(end == -1){ end = seqLength; }
243 sequence.setUnaligned(rawSequence.substr(0,end));
244 trimQScores(-1, end);
248 catch(exception& e) {
249 m->errorOut(e, "QualityScores", "flipQScores");
255 /**************************************************************************************************/
257 bool QualityScores::stripQualWindowAverage(Sequence& sequence, int stepSize, int windowSize, double qThreshold){
259 string rawSequence = sequence.getUnaligned();
260 int seqLength = sequence.getNumBases();
262 if(seqName != sequence.getName()){
263 m->mothurOut("sequence name mismatch between fasta: " + sequence.getName() + " and qual file: " + seqName);
264 m->mothurOutEndLine();
267 int end = windowSize;
271 while(start < seqLength){
272 double windowSum = 0.0000;
274 for(int i=start;i<end;i++){
275 windowSum += qScores[i];
277 double windowAverage = windowSum / (double)(end-start);
279 if(windowAverage < qThreshold){
280 end = end - stepSize;
284 end = start + windowSize;
285 if(end >= seqLength){ end = seqLength - 1; }
289 if(end == -1){ end = seqLength; }
291 sequence.setUnaligned(rawSequence.substr(0,end));
292 trimQScores(-1, end);
296 catch(exception& e) {
297 m->errorOut(e, "QualityScores", "flipQScores");
303 /**************************************************************************************************/
305 double QualityScores::calculateAverage(){
307 double aveQScore = 0.0000;
309 for(int i=0;i<seqLength;i++){
310 aveQScore += (double) qScores[i];
312 aveQScore /= (double) seqLength;
317 /**************************************************************************************************/
319 bool QualityScores::cullQualAverage(Sequence& sequence, double qAverage){
321 string rawSequence = sequence.getUnaligned();
322 bool success = 0; //guilty until proven innocent
324 if(seqName != sequence.getName()) {
325 m->mothurOut("sequence name mismatch btwn fasta: " + sequence.getName() + " and qual file: " + seqName);
326 m->mothurOutEndLine();
329 double aveQScore = calculateAverage();
331 if(aveQScore >= qAverage) { success = 1; }
332 else { success = 0; }
336 catch(exception& e) {
337 m->errorOut(e, "TrimSeqsCommand", "cullQualAverage");
342 /**************************************************************************************************/