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){
31 m = MothurOut::getInstance();
39 istringstream nameStream(line);
41 nameStream >> seqName;
42 seqName = seqName.substr(1);
45 istringstream qualStream(line);
49 qScores.push_back(score);
53 seqLength = qScores.size();
56 m->errorOut(e, "QualityScores", "QualityScores");
62 /**************************************************************************************************/
64 string QualityScores::getName(){
70 m->errorOut(e, "QualityScores", "getName");
75 /**************************************************************************************************/
77 void QualityScores::printQScores(ofstream& qFile){
80 double aveQScore = calculateAverage();
82 qFile << '>' << seqName << '\t' << aveQScore << endl;
84 for(int i=0;i<seqLength;i++){
85 qFile << qScores[i] << ' ';
91 m->errorOut(e, "QualityScores", "printQScores");
96 /**************************************************************************************************/
98 void QualityScores::trimQScores(int start, int end){
103 hold = vector<int>(qScores.begin()+start, qScores.end());
107 hold = vector<int>(qScores.begin(), qScores.begin()+end); //not sure if indexing is correct
111 seqLength = qScores.size();
113 catch(exception& e) {
114 m->errorOut(e, "QualityScores", "trimQScores");
119 /**************************************************************************************************/
121 void QualityScores::flipQScores(){
124 vector<int> temp = qScores;
125 for(int i=0;i<seqLength;i++){
126 qScores[seqLength - i - 1] = temp[i];
130 catch(exception& e) {
131 m->errorOut(e, "QualityScores", "flipQScores");
136 /**************************************************************************************************/
138 bool QualityScores::stripQualThreshold(Sequence& sequence, double qThreshold){
140 string rawSequence = sequence.getUnaligned();
141 int seqLength = sequence.getNumBases();
143 if(seqName != sequence.getName()){
144 m->mothurOut("sequence name mismatch btwn fasta: " + sequence.getName() + " and qual file: " + seqName);
145 m->mothurOutEndLine();
149 for(int i=0;i<seqLength;i++){
151 if(qScores[i] < qThreshold){
156 sequence.setUnaligned(rawSequence.substr(0,end));
157 trimQScores(-1, end);
161 catch(exception& e) {
162 m->errorOut(e, "QualityScores", "flipQScores");
168 /**************************************************************************************************/
170 bool QualityScores::stripQualRollingAverage(Sequence& sequence, double qThreshold){
172 string rawSequence = sequence.getUnaligned();
173 int seqLength = sequence.getNumBases();
175 if(seqName != sequence.getName()){
176 m->mothurOut("sequence name mismatch btwn fasta: " + sequence.getName() + " and qual file: " + seqName);
177 m->mothurOutEndLine();
181 double rollingSum = 0.0000;
183 for(int i=0;i<seqLength;i++){
185 rollingSum += (double)qScores[i];
187 if(rollingSum / (double)(i+1) < qThreshold){
189 // cout << i+1 << '\t' << seqName << '\t' << rollingSum / (double)(i+1) << endl;
194 if(end == -1){ end = seqLength; }
196 sequence.setUnaligned(rawSequence.substr(0,end));
197 trimQScores(-1, end);
201 catch(exception& e) {
202 m->errorOut(e, "QualityScores", "flipQScores");
208 /**************************************************************************************************/
210 bool QualityScores::stripQualWindowAverage(Sequence& sequence, int stepSize, int windowSize, double qThreshold){
212 string rawSequence = sequence.getUnaligned();
213 int seqLength = sequence.getNumBases();
215 if(seqName != sequence.getName()){
216 m->mothurOut("sequence name mismatch btwn fasta: " + sequence.getName() + " and qual file: " + seqName);
217 m->mothurOutEndLine();
220 int end = windowSize;
224 while(start < seqLength){
225 double windowSum = 0.0000;
227 for(int i=start;i<end;i++){
228 windowSum += qScores[i];
230 double windowAverage = windowSum / (double)(end-start);
232 if(windowAverage < qThreshold){
237 end = start + windowSize;
238 if(end >= seqLength){ end = seqLength - 1; }
242 if(end == -1){ end = seqLength; }
244 sequence.setUnaligned(rawSequence.substr(0,end));
245 trimQScores(-1, end);
249 catch(exception& e) {
250 m->errorOut(e, "QualityScores", "flipQScores");
256 /**************************************************************************************************/
258 double QualityScores::calculateAverage(){
260 double aveQScore = 0.0000;
262 for(int i=0;i<seqLength;i++){
263 aveQScore += (double) qScores[i];
265 aveQScore /= (double) seqLength;
270 /**************************************************************************************************/
272 bool QualityScores::cullQualAverage(Sequence& sequence, double qAverage){
274 string rawSequence = sequence.getUnaligned();
275 bool success = 0; //guilty until proven innocent
277 if(seqName != sequence.getName()) {
278 m->mothurOut("sequence name mismatch btwn fasta: " + sequence.getName() + " and qual file: " + seqName);
279 m->mothurOutEndLine();
282 double aveQScore = calculateAverage();
284 if(aveQScore >= qAverage) { success = 1; }
285 else { success = 0; }
289 catch(exception& e) {
290 m->errorOut(e, "TrimSeqsCommand", "cullQualAverage");
295 /**************************************************************************************************/