--- /dev/null
+/*
+ * qualityscores.cpp
+ * Mothur
+ *
+ * Created by Pat Schloss on 7/12/10.
+ * Copyright 2010 Schloss Lab. All rights reserved.
+ *
+ */
+
+#include "qualityscores.h"
+
+/**************************************************************************************************/
+
+QualityScores::QualityScores(){
+ try {
+ m = MothurOut::getInstance();
+ seqName = "";
+ seqLength = -1;
+
+ }
+ catch(exception& e) {
+ m->errorOut(e, "QualityScores", "QualityScores");
+ exit(1);
+ }
+}
+
+/**************************************************************************************************/
+
+QualityScores::QualityScores(ifstream& qFile){
+ try {
+
+ m = MothurOut::getInstance();
+
+ seqName = "";
+ int score;
+
+ qFile >> seqName;
+ m->getline(qFile);
+ //cout << seqName << endl;
+ if (seqName == "") {
+ m->mothurOut("Error reading quality file, name blank at position, " + toString(qFile.tellg()));
+ m->mothurOutEndLine();
+ }
+ else{
+ seqName = seqName.substr(1);
+ }
+
+ string qScoreString = m->getline(qFile);
+ //cout << qScoreString << endl;
+ while(qFile.peek() != '>' && qFile.peek() != EOF){
+ if (m->control_pressed) { break; }
+ string temp = m->getline(qFile);
+ //cout << temp << endl;
+ qScoreString += ' ' + temp;
+ }
+ //cout << "done reading " << endl;
+ istringstream qScoreStringStream(qScoreString);
+ int count = 0;
+ while(!qScoreStringStream.eof()){
+ if (m->control_pressed) { break; }
+ string temp;
+ qScoreStringStream >> temp; m->gobble(qScoreStringStream);
+
+ //check temp to make sure its a number
+ if (!m->isContainingOnlyDigits(temp)) { m->mothurOut("[ERROR]: In sequence " + seqName + "'s quality scores, expected a number and got " + temp + ", setting score to 0."); m->mothurOutEndLine(); temp = "0"; }
+ convert(temp, score);
+
+ //cout << count << '\t' << score << endl;
+ qScores.push_back(score);
+ count++;
+ }
+ //qScores.pop_back();
+
+// string scores = "";
+//
+// while(!qFile.eof()){
+//
+// qFile >> seqName;
+//
+// //get name
+// if (seqName.length() != 0) {
+// seqName = seqName.substr(1);
+// while (!qFile.eof()) {
+// char c = qFile.get();
+// //gobble junk on line
+// if (c == 10 || c == 13){ break; }
+// }
+// m->gobble(qFile);
+// }
+//
+// //get scores
+// while(qFile){
+// char letter=qFile.get();
+// if((letter == '>')){ qFile.putback(letter); break; }
+// else if (isprint(letter)) { scores += letter; }
+// }
+// m->gobble(qFile);
+//
+// break;
+// }
+//
+// //convert scores string to qScores
+// istringstream qScoreStringStream(scores);
+//
+// int score;
+// while(!qScoreStringStream.eof()){
+//
+// if (m->control_pressed) { break; }
+//
+// qScoreStringStream >> score;
+// qScores.push_back(score);
+// }
+//
+// qScores.pop_back();
+
+ seqLength = qScores.size();
+ //cout << "seqlength = " << seqLength << '\t' << count << endl;
+
+ }
+ catch(exception& e) {
+ m->errorOut(e, "QualityScores", "QualityScores");
+ exit(1);
+ }
+
+}
+
+/**************************************************************************************************/
+
+string QualityScores::getName(){
+
+ try {
+ return seqName;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "QualityScores", "getName");
+ exit(1);
+ }
+}
+
+/**************************************************************************************************/
+
+void QualityScores::printQScores(ofstream& qFile){
+ try {
+
+ double aveQScore = calculateAverage();
+
+ qFile << '>' << seqName << '\t' << aveQScore << endl;
+
+ for(int i=0;i<seqLength;i++){
+ qFile << qScores[i] << ' ';
+ }
+ qFile << endl;
+
+ }
+ catch(exception& e) {
+ m->errorOut(e, "QualityScores", "printQScores");
+ exit(1);
+ }
+}
+
+/**************************************************************************************************/
+
+void QualityScores::trimQScores(int start, int end){
+ try {
+ vector<int> hold;
+
+
+ //cout << seqName << '\t' << start << '\t' << end << '\t' << qScores.size() << endl;
+ //for (int i = 0; i < qScores.size(); i++) { cout << qScores[i] << end; }
+ if(end == -1){
+ hold = vector<int>(qScores.begin()+start, qScores.end());
+ qScores = hold;
+ }
+ if(start == -1){
+ if(qScores.size() > end){
+ hold = vector<int>(qScores.begin(), qScores.begin()+end);
+ qScores = hold;
+ }
+ }
+
+ seqLength = qScores.size();
+ }
+ catch(exception& e) {
+ m->errorOut(e, "QualityScores", "trimQScores");
+ exit(1);
+ }
+}
+
+/**************************************************************************************************/
+
+void QualityScores::flipQScores(){
+ try {
+
+ vector<int> temp = qScores;
+ for(int i=0;i<seqLength;i++){
+ qScores[seqLength - i - 1] = temp[i];
+ }
+
+ }
+ catch(exception& e) {
+ m->errorOut(e, "QualityScores", "flipQScores");
+ exit(1);
+ }
+}
+
+/**************************************************************************************************/
+
+bool QualityScores::stripQualThreshold(Sequence& sequence, double qThreshold){
+ try {
+ string rawSequence = sequence.getUnaligned();
+ int seqLength = sequence.getNumBases();
+
+ if(seqName != sequence.getName()){
+ m->mothurOut("sequence name mismatch btwn fasta: " + sequence.getName() + " and qual file: " + seqName);
+ m->mothurOutEndLine();
+ }
+
+ int end;
+ for(int i=0;i<seqLength;i++){
+ end = i;
+ if(qScores[i] < qThreshold){
+ break;
+ }
+ }
+
+ //every score passed
+ if (end == (seqLength-1)) { end = seqLength; }
+
+ sequence.setUnaligned(rawSequence.substr(0,end));
+ trimQScores(-1, end);
+
+ return 1;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "QualityScores", "flipQScores");
+ exit(1);
+ }
+
+}
+
+/**************************************************************************************************/
+
+bool QualityScores::stripQualRollingAverage(Sequence& sequence, double qThreshold){
+ try {
+ string rawSequence = sequence.getUnaligned();
+ int seqLength = sequence.getNumBases();
+
+ if(seqName != sequence.getName()){
+ m->mothurOut("sequence name mismatch btwn fasta: " + sequence.getName() + " and qual file: " + seqName);
+ m->mothurOutEndLine();
+ }
+
+ int end = -1;
+ double rollingSum = 0.0000;
+
+ for(int i=0;i<seqLength;i++){
+
+ rollingSum += (double)qScores[i];
+
+ if(rollingSum / (double)(i+1) < qThreshold){
+ end = i;
+ break;
+ }
+ }
+
+ if(end == -1){ end = seqLength; }
+
+
+ sequence.setUnaligned(rawSequence.substr(0,end));
+ trimQScores(-1, end);
+
+
+ return 1;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "QualityScores", "flipQScores");
+ exit(1);
+ }
+
+}
+
+/**************************************************************************************************/
+
+bool QualityScores::stripQualWindowAverage(Sequence& sequence, int stepSize, int windowSize, double qThreshold){
+ try {
+ string rawSequence = sequence.getUnaligned();
+ int seqLength = sequence.getNumBases();
+
+ if(seqName != sequence.getName()){
+ m->mothurOut("sequence name mismatch between fasta: " + sequence.getName() + " and qual file: " + seqName);
+ m->mothurOutEndLine();
+ }
+
+ int end = windowSize;
+ int start = 0;
+
+ if(seqLength < windowSize) { return 0; }
+
+ while((start+windowSize) < seqLength){
+ double windowSum = 0.0000;
+
+ for(int i=start;i<end;i++){
+ windowSum += qScores[i];
+ }
+ double windowAverage = windowSum / (double)(end-start);
+
+ if(windowAverage < qThreshold){
+ end = end - stepSize;
+ break;
+ }
+
+ start += stepSize;
+ end = start + windowSize;
+
+ if(end >= seqLength){ end = seqLength; }
+
+ }
+
+ if(end == -1){ end = seqLength; }
+
+ //failed first window
+ if (end < windowSize) { return 0; }
+
+ sequence.setUnaligned(rawSequence.substr(0,end));
+ trimQScores(-1, end);
+
+ return 1;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "QualityScores", "stripQualWindowAverage");
+ exit(1);
+ }
+
+}
+
+/**************************************************************************************************/
+
+double QualityScores::calculateAverage(){
+
+ double aveQScore = 0.0000;
+
+ for(int i=0;i<seqLength;i++){
+ aveQScore += (double) qScores[i];
+ }
+ aveQScore /= (double) seqLength;
+
+ return aveQScore;
+}
+
+/**************************************************************************************************/
+
+bool QualityScores::cullQualAverage(Sequence& sequence, double qAverage){
+ try {
+ string rawSequence = sequence.getUnaligned();
+ bool success = 0; //guilty until proven innocent
+
+ if(seqName != sequence.getName()) {
+ m->mothurOut("sequence name mismatch btwn fasta: " + sequence.getName() + " and qual file: " + seqName);
+ m->mothurOutEndLine();
+ }
+
+ double aveQScore = calculateAverage();
+
+ if(aveQScore >= qAverage) { success = 1; }
+ else { success = 0; }
+
+ return success;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "QualityScores", "cullQualAverage");
+ exit(1);
+ }
+}
+
+/**************************************************************************************************/
+
+void QualityScores::updateQScoreErrorMap(map<char, vector<int> >& qualErrorMap, string errorSeq, int start, int stop, int weight){
+ try {
+
+ int seqLength = errorSeq.size();
+
+ int qIndex = start - 1;
+
+ for(int i=0;i<seqLength;i++){
+
+ if(errorSeq[i] == 'm') { qualErrorMap['m'][qScores[qIndex]] += weight; }
+ else if(errorSeq[i] == 's') { qualErrorMap['s'][qScores[qIndex]] += weight; }
+ else if(errorSeq[i] == 'i') { qualErrorMap['i'][qScores[qIndex]] += weight; }
+ else if(errorSeq[i] == 'a') { qualErrorMap['a'][qScores[qIndex]] += weight; /*if(qScores[qIndex] != 0){ cout << qIndex << '\t'; }*/ }
+ else if(errorSeq[i] == 'd') { /* there are no qScores for deletions */ }
+
+ if(errorSeq[i] != 'd') { qIndex++; }
+
+ if(qIndex > stop){ break; }
+ }
+ }
+ catch(exception& e) {
+ m->errorOut(e, "QualityScores", "updateQScoreErrorMap");
+ exit(1);
+ }
+}
+
+/**************************************************************************************************/
+
+void QualityScores::updateForwardMap(vector<vector<int> >& forwardMap, int start, int stop, int weight){
+ try {
+
+ int index = 0;
+ for(int i=start-1;i<stop;i++){
+ forwardMap[index++][qScores[i]] += weight;
+ }
+
+ }
+ catch(exception& e) {
+ m->errorOut(e, "QualityScores", "updateForwardMap");
+ exit(1);
+ }
+}
+
+/**************************************************************************************************/
+
+void QualityScores::updateReverseMap(vector<vector<int> >& reverseMap, int start, int stop, int weight){
+ try {
+
+ int index = 0;
+ for(int i=stop-1;i>=start;i--){
+ reverseMap[index++][qScores[i]] += weight;
+ }
+
+ }
+ catch(exception& e) {
+ m->errorOut(e, "QualityScores", "updateForwardMap");
+ exit(1);
+ }
+}
+
+/**************************************************************************************************/
projects / mothur.git / blobdiff