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();
60 for(int i=0;i<seqLength;i++){
62 qScores.push_back(score);
68 m->errorOut(e, "QualityScores", "QualityScores");
74 /**************************************************************************************************/
76 string QualityScores::getName(){
82 m->errorOut(e, "QualityScores", "getName");
87 /**************************************************************************************************/
89 void QualityScores::printQScores(ofstream& qFile){
92 double aveQScore = calculateAverage();
94 qFile << '>' << seqName << '\t' << aveQScore << endl;
96 for(int i=0;i<seqLength;i++){
97 qFile << qScores[i] << ' ';
102 catch(exception& e) {
103 m->errorOut(e, "QualityScores", "printQScores");
108 /**************************************************************************************************/
110 void QualityScores::trimQScores(int start, int end){
115 hold = vector<int>(qScores.begin()+start, qScores.end());
119 hold = vector<int>(qScores.begin(), qScores.begin()+end); //not sure if indexing is correct
123 seqLength = qScores.size();
125 catch(exception& e) {
126 m->errorOut(e, "QualityScores", "trimQScores");
131 /**************************************************************************************************/
133 void QualityScores::flipQScores(){
136 vector<int> temp = qScores;
137 for(int i=0;i<seqLength;i++){
138 qScores[seqLength - i - 1] = temp[i];
142 catch(exception& e) {
143 m->errorOut(e, "QualityScores", "flipQScores");
148 /**************************************************************************************************/
150 bool QualityScores::stripQualThreshold(Sequence& sequence, double qThreshold){
152 string rawSequence = sequence.getUnaligned();
153 int seqLength = sequence.getNumBases();
155 if(seqName != sequence.getName()){
156 m->mothurOut("sequence name mismatch btwn fasta: " + sequence.getName() + " and qual file: " + seqName);
157 m->mothurOutEndLine();
161 for(int i=0;i<seqLength;i++){
163 if(qScores[i] < qThreshold){
168 sequence.setUnaligned(rawSequence.substr(0,end));
169 trimQScores(-1, end);
173 catch(exception& e) {
174 m->errorOut(e, "QualityScores", "flipQScores");
180 /**************************************************************************************************/
182 bool QualityScores::stripQualRollingAverage(Sequence& sequence, double qThreshold){
184 string rawSequence = sequence.getUnaligned();
185 int seqLength = sequence.getNumBases();
187 if(seqName != sequence.getName()){
188 m->mothurOut("sequence name mismatch btwn fasta: " + sequence.getName() + " and qual file: " + seqName);
189 m->mothurOutEndLine();
193 double rollingSum = 0.0000;
195 for(int i=0;i<seqLength;i++){
197 rollingSum += (double)qScores[i];
199 if(rollingSum / (double)(i+1) < qThreshold){
201 // cout << i+1 << '\t' << seqName << '\t' << rollingSum / (double)(i+1) << endl;
206 if(end == -1){ end = seqLength; }
208 sequence.setUnaligned(rawSequence.substr(0,end));
209 trimQScores(-1, end);
213 catch(exception& e) {
214 m->errorOut(e, "QualityScores", "flipQScores");
220 /**************************************************************************************************/
222 bool QualityScores::stripQualWindowAverage(Sequence& sequence, int stepSize, int windowSize, double qThreshold){
224 string rawSequence = sequence.getUnaligned();
225 int seqLength = sequence.getNumBases();
227 if(seqName != sequence.getName()){
228 m->mothurOut("sequence name mismatch between fasta: " + sequence.getName() + " and qual file: " + seqName);
229 m->mothurOutEndLine();
232 int end = windowSize;
236 while(start < seqLength){
237 double windowSum = 0.0000;
239 for(int i=start;i<end;i++){
240 windowSum += qScores[i];
242 double windowAverage = windowSum / (double)(end-start);
244 if(windowAverage < qThreshold){
245 end = end - stepSize;
249 end = start + windowSize;
250 if(end >= seqLength){ end = seqLength - 1; }
254 if(end == -1){ end = seqLength; }
256 sequence.setUnaligned(rawSequence.substr(0,end));
257 trimQScores(-1, end);
261 catch(exception& e) {
262 m->errorOut(e, "QualityScores", "flipQScores");
268 /**************************************************************************************************/
270 double QualityScores::calculateAverage(){
272 double aveQScore = 0.0000;
274 for(int i=0;i<seqLength;i++){
275 aveQScore += (double) qScores[i];
277 aveQScore /= (double) seqLength;
282 /**************************************************************************************************/
284 bool QualityScores::cullQualAverage(Sequence& sequence, double qAverage){
286 string rawSequence = sequence.getUnaligned();
287 bool success = 0; //guilty until proven innocent
289 if(seqName != sequence.getName()) {
290 m->mothurOut("sequence name mismatch btwn fasta: " + sequence.getName() + " and qual file: " + seqName);
291 m->mothurOutEndLine();
294 double aveQScore = calculateAverage();
296 if(aveQScore >= qAverage) { success = 1; }
297 else { success = 0; }
301 catch(exception& e) {
302 m->errorOut(e, "TrimSeqsCommand", "cullQualAverage");
307 /**************************************************************************************************/