fixed trim.seqs bug with qtrim parameter and added num=1 special case to database...

[mothur.git] / trimseqscommand.cpp

diff --git a/trimseqscommand.cpp b/trimseqscommand.cpp
index 9fc3abb8422341e755a1e54da38453002cf01dce..2dc07794be5abb2f239e4b228c3aacbf51c51b6e 100644 (file)
--- a/trimseqscommand.cpp
+++ b/trimseqscommand.cpp
@@ -11,6 +11,7 @@
 #include "needlemanoverlap.hpp"
 
 //**********************************************************************************************************************
+
 vector<string> TrimSeqsCommand::getValidParameters(){  
        try {
                string Array[] =  {"fasta", "flip", "oligos", "maxambig", "maxhomop", "group","minlength", "maxlength", "qfile", 
@@ -25,7 +26,9 @@ vector<string> TrimSeqsCommand::getValidParameters(){
                exit(1);
        }
 }
+
 //**********************************************************************************************************************
+
 TrimSeqsCommand::TrimSeqsCommand(){    
        try {
                abort = true;
@@ -40,7 +43,9 @@ TrimSeqsCommand::TrimSeqsCommand(){
                exit(1);
        }
 }
+
 //**********************************************************************************************************************
+
 vector<string> TrimSeqsCommand::getRequiredParameters(){       
        try {
                string Array[] =  {"fasta"};
@@ -52,7 +57,9 @@ vector<string> TrimSeqsCommand::getRequiredParameters(){
                exit(1);
        }
 }
+
 //**********************************************************************************************************************
+
 vector<string> TrimSeqsCommand::getRequiredFiles(){    
        try {
                vector<string> myArray;
@@ -63,6 +70,7 @@ vector<string> TrimSeqsCommand::getRequiredFiles(){
                exit(1);
        }
 }
+
 //***************************************************************************************************************
 
 TrimSeqsCommand::TrimSeqsCommand(string option)  {
@@ -254,6 +262,7 @@ TrimSeqsCommand::TrimSeqsCommand(string option)  {
                exit(1);
        }
 }
+
 //**********************************************************************************************************************
 
 void TrimSeqsCommand::help(){
@@ -446,7 +455,7 @@ int TrimSeqsCommand::driverCreateTrim(string filename, string qFileName, string
        try {
                
                ofstream outFASTA;
-               int able = m->openOutputFile(trimFile, outFASTA);
+               m->openOutputFile(trimFile, outFASTA);
                
                ofstream scrapFASTA;
                m->openOutputFile(scrapFile, scrapFASTA);
@@ -547,16 +556,16 @@ int TrimSeqsCommand::driverCreateTrim(string filename, string qFileName, string
 
                                
                                if(qFileName != ""){
+                                       int origLength = currSeq.getNumBases();
                                        
                                        if(qThreshold != 0)                     {       success = currQual.stripQualThreshold(currSeq, qThreshold);                     }
                                        else if(qAverage != 0)          {       success = currQual.cullQualAverage(currSeq, qAverage);                          }
                                        else if(qRollAverage != 0)      {       success = currQual.stripQualRollingAverage(currSeq, qRollAverage);      }
                                        else if(qWindowAverage != 0){   success = currQual.stripQualWindowAverage(currSeq, qWindowStep, qWindowSize, qWindowAverage);   }
                                        else                                            {       success = 1;                            }
-
-//                                     if (qtrim == 1 && (origSeq.length() != currSeq.getUnaligned().length())) { 
-//                                             success = 0; //if you don't want to trim and the sequence does not meet quality requirements, move to scrap
-//                                     }
+                                       
+                                       //you don't want to trim, if it fails above then scrap it
+                                       if ((!qtrim) && (origLength != currSeq.getNumBases())) { success = 0; }
                                        
                                        if(!success)                            {       trashCode += 'q';       }
                                }                               
@@ -874,6 +883,7 @@ int TrimSeqsCommand::setLines(string filename, string qfilename, vector<unsigned
                exit(1);
        }
 }
+
 //***************************************************************************************************************
 
 void TrimSeqsCommand::getOligos(vector<string>& outFASTAVec, vector<string>& outQualVec){
@@ -1085,7 +1095,6 @@ int TrimSeqsCommand::stripBarcode(Sequence& seq, QualityScores& qual, int& group
                                oligo = oligo.substr(0,alnLength);
                                temp = temp.substr(0,alnLength);
                                
-                               int newStart=0;
                                int numDiff = countDiffs(oligo, temp);
                                
 //                             cout << oligo << '\t' << temp << '\t' << numDiff << endl;                               
@@ -1210,7 +1219,6 @@ int TrimSeqsCommand::stripForward(Sequence& seq, QualityScores& qual, int& group
                                oligo = oligo.substr(0,alnLength);
                                temp = temp.substr(0,alnLength);
                                
-                               int newStart=0;
                                int numDiff = countDiffs(oligo, temp);
                                
 //                             cout << oligo << '\t' << temp << '\t' << numDiff << endl;                               
@@ -1433,6 +1441,7 @@ bool TrimSeqsCommand::compareDNASeq(string oligo, string seq){
        }
 
 }
+
 //***************************************************************************************************************
 
 int TrimSeqsCommand::countDiffs(string oligo, string seq){