m = MothurOut::getInstance();
int score;
- seqName = getSequenceName(qFile);
+ seqName = getSequenceName(qFile); m->gobble(qFile);
+ if (m->debug) { m->mothurOut("[DEBUG]: name = '" + seqName + "'\n."); }
+
if (!m->control_pressed) {
- string qScoreString = m->getline(qFile);
- //cout << qScoreString << endl;
+ string qScoreString = m->getline(qFile); m->gobble(qFile);
+
+ if (m->debug) { m->mothurOut("[DEBUG]: scores = '" + qScoreString + "'\n."); }
+
while(qFile.peek() != '>' && qFile.peek() != EOF){
if (m->control_pressed) { break; }
- string temp = m->getline(qFile);
- //cout << temp << endl;
+ string temp = m->getline(qFile); m->gobble(qFile);
+ //if (m->debug) { m->mothurOut("[DEBUG]: scores = '" + temp + "'\n."); }
qScoreString += ' ' + temp;
}
- //cout << "done reading " << endl;
+ //cout << "done reading " << endl;
istringstream qScoreStringStream(qScoreString);
int count = 0;
while(!qScoreStringStream.eof()){
string temp;
qScoreStringStream >> temp; m->gobble(qScoreStringStream);
+ //if (m->debug) { m->mothurOut("[DEBUG]: score " + toString(qScores.size()) + " = '" + temp + "'\n."); }
+
//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);
}
seqLength = qScores.size();
- //cout << "seqlength = " << seqLength << '\t' << count << endl;
+ //cout << "seqlength = " << seqLength << endl;
}
catch(exception& e) {
name = name.substr(1);
- for (int i = 0; i < name.length(); i++) {
- if (name[i] == ':') { name[i] = '_'; m->changedSeqNames = true; }
- }
+ m->checkName(name);
}else{ m->mothurOut("Error in reading your qfile, at position " + toString(qFile.tellg()) + ". Blank name."); m->mothurOutEndLine(); m->control_pressed = true; }
void QualityScores::setName(string name) {
try {
- for (int i = 0; i < name.length(); i++) {
- if (name[i] == ':') { name[i] = '_'; m->changedSeqNames = true; }
- }
-
+ m->checkName(name);
seqName = name;
}
catch(exception& e) {
void QualityScores::printQScores(ofstream& qFile){
try {
- double aveQScore = calculateAverage();
+ double aveQScore = calculateAverage(false);
qFile << '>' << seqName << '\t' << aveQScore << endl;
if(seqName != sequence.getName()){
m->mothurOut("sequence name mismatch btwn fasta: " + sequence.getName() + " and qual file: " + seqName);
- m->mothurOutEndLine();
+ m->mothurOutEndLine(); m->control_pressed = true;
}
int end;
/**************************************************************************************************/
-bool QualityScores::stripQualRollingAverage(Sequence& sequence, double qThreshold){
+bool QualityScores::stripQualRollingAverage(Sequence& sequence, double qThreshold, bool logTransform){
try {
string rawSequence = sequence.getUnaligned();
int seqLength = sequence.getNumBases();
int end = -1;
double rollingSum = 0.0000;
+ double value = 0.0;
for(int i=0;i<seqLength;i++){
-
- rollingSum += (double)qScores[i];
+
+ if (logTransform) {
+ rollingSum += (double)pow(10.0, qScores[i]);
+ value = log10(rollingSum / (double)(i+1));
+
+ } //Sum 10^Q
+ else {
+ rollingSum += (double)qScores[i];
+ value = rollingSum / (double)(i+1);
+ }
+
- if(rollingSum / (double)(i+1) < qThreshold){
+ if(value < qThreshold){
end = i;
break;
}
/**************************************************************************************************/
-bool QualityScores::stripQualWindowAverage(Sequence& sequence, int stepSize, int windowSize, double qThreshold){
+bool QualityScores::stripQualWindowAverage(Sequence& sequence, int stepSize, int windowSize, double qThreshold, bool logTransform){
try {
string rawSequence = sequence.getUnaligned();
int seqLength = sequence.getNumBases();
double windowSum = 0.0000;
for(int i=start;i<end;i++){
- windowSum += qScores[i];
+ if (logTransform) { windowSum += pow(10.0, qScores[i]); }
+ else { windowSum += qScores[i]; }
}
- double windowAverage = windowSum / (double)(end-start);
+ double windowAverage = 0.0;
+ if (logTransform) { windowAverage = log10(windowSum / (double)(end-start)); }
+ else { windowAverage = windowSum / (double)(end-start); }
if(windowAverage < qThreshold){
end = end - stepSize;
/**************************************************************************************************/
-double QualityScores::calculateAverage(){
+double QualityScores::calculateAverage(bool logTransform){
double aveQScore = 0.0000;
for(int i=0;i<seqLength;i++){
- aveQScore += (double) qScores[i];
+ if (logTransform) { aveQScore += pow(10.0, qScores[i]); }
+ else { aveQScore += qScores[i]; }
}
aveQScore /= (double) seqLength;
+
+ if (logTransform) { aveQScore = log10(aveQScore /(double) seqLength); }
+ else { aveQScore /= (double) seqLength; }
return aveQScore;
}
/**************************************************************************************************/
-bool QualityScores::cullQualAverage(Sequence& sequence, double qAverage){
+bool QualityScores::cullQualAverage(Sequence& sequence, double qAverage, bool logTransform){
try {
string rawSequence = sequence.getUnaligned();
bool success = 0; //guilty until proven innocent
m->mothurOutEndLine();
}
- double aveQScore = calculateAverage();
+ double aveQScore = calculateAverage(logTransform);
if(aveQScore >= qAverage) { success = 1; }
else { success = 0; }