small bug with trimming of quality scores over a window. lost last base if sequence...

diff --git a/indicatorcommand.cpp b/indicatorcommand.cpp
index 2d1c63e21dbcd82df19957cb533d1e420e434fdc..76816df40a7b962decb1ebcaaeb421546c77539a 100644 (file)
--- a/indicatorcommand.cpp
+++ b/indicatorcommand.cpp
@@ -41,7 +41,7 @@ string IndicatorCommand::getHelpString(){
                helpString += "The new tree contains labels at each internal node.  The label is the node number so you can relate the tree to the summary file.\n";
                helpString += "The summary file lists the indicator value for each OTU for each node.\n";
                helpString += "The indicator command parameters are tree, groups, shared, relabund, design and label. The tree parameter is required as well as either shared or relabund.\n";
-               helpString += "The design parameter allows you to provide a design file to relate the tree to the shared or relabund file.\n";          
+               helpString += "The design parameter allows you to relate the tree to the shared or relabund file, if your tree contains the grouping names, or if no tree is provided to group your groups into groupings.\n";                  
                helpString += "The groups parameter allows you to specify which of the groups in your shared or relabund you would like analyzed, or if you provide a design file the groups in your design file.  The groups may be entered separated by dashes.\n";
                helpString += "The label parameter indicates at what distance your tree relates to the shared or relabund.\n";
                helpString += "The iters parameter allows you to set number of randomization for the P value.  The default is 1000.";
@@ -379,6 +379,7 @@ int IndicatorCommand::GetIndicatorSpecies(){
                ofstream out;
                m->openOutputFile(outputFileName, out);
                out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
+               m->mothurOutEndLine(); m->mothurOut("Species\tIndicatorValue\tpValue\n");
                
                int numBins = 0;
                if (sharedfile != "") { numBins = lookup[0]->getNumBins(); }
@@ -472,8 +473,6 @@ int IndicatorCommand::GetIndicatorSpecies(){
                /******************************************************/
                //output indicator values to table form               //
                /*****************************************************/
-               int indicatorOTU;
-               double maxValue, pvalue; maxValue=0.0; pvalue=0.0;
                out << "OTU\tIndicator_Value\tpValue" << endl;
                for (int j = 0; j < indicatorValues.size(); j++) {
                                
@@ -484,20 +483,16 @@ int IndicatorCommand::GetIndicatorSpecies(){
                        if (pValues[j] > (1/(float)iters)) { out << pValues[j] << endl; } 
                        else { out << "<" << (1/(float)iters) << endl; }
                        
-                       if (maxValue < indicatorValues[j]) { 
-                               maxValue = indicatorValues[j];
-                               indicatorOTU = j;
+                       if (pValues[j] <= 0.05) {
+                               cout << "OTU" << j+1 << '\t' << indicatorValues[j]  << '\t';
+                               string pValueString = "<" + toString((1/(float)iters)); 
+                               if (pValues[j] > (1/(float)iters)) { pValueString = toString(pValues[j]); cout << pValues[j];} 
+                               else { cout << "<" << (1/(float)iters); }
+                               m->mothurOutJustToLog("OTU" + toString(j+1) + "\t" + toString(indicatorValues[j]) + "\t" + pValueString); 
+                               m->mothurOutEndLine(); 
                        }
                }
                
-               m->mothurOutEndLine(); m->mothurOut("Species\tIndicatorValue\tpValue\n");
-               cout << "OTU" << indicatorOTU+1 << '\t' << maxValue << '\t';
-               string pValueString = "<" + toString((1/(float)iters)); 
-               if (pValues[indicatorOTU] > (1/(float)iters)) { pValueString = toString(pValues[indicatorOTU]); cout << pValues[indicatorOTU];} 
-               else { cout << "<" << (1/(float)iters); }
-               m->mothurOutJustToLog("OTU" + toString(indicatorOTU+1) + "\t" + toString(maxValue) + "\t" + pValueString); 
-               m->mothurOutEndLine(); 
-               
                out.close();
                
                return 0;
@@ -554,7 +549,7 @@ int IndicatorCommand::GetIndicatorSpecies(Tree*& T){
                        
                //for each node
                for (int i = T->getNumLeaves(); i < T->getNumNodes(); i++) {
-                               
+                       //cout << endl << i+1 << endl;  
                        if (m->control_pressed) { out.close(); return 0; }
                        
                        /*****************************************************/
@@ -699,8 +694,6 @@ int IndicatorCommand::GetIndicatorSpecies(Tree*& T){
                        //output indicator values to table form + label tree  //
                        /*****************************************************/
                        out << (i+1) << '\t';
-                       int indicatorOTU;
-                       double maxValue, pvalue; maxValue=0.0; pvalue=0.0;
                        for (int j = 0; j < indicatorValues.size(); j++) {
                                
                                if (m->control_pressed) { out.close(); return 0; }
@@ -711,23 +704,18 @@ int IndicatorCommand::GetIndicatorSpecies(Tree*& T){
                                        out << indicatorValues[j] << '\t' << pValues[j] << '\t';
                                }
                                
-                               if (maxValue < indicatorValues[j]) { 
-                                       maxValue = indicatorValues[j];
-                                       indicatorOTU = j;
+                               if (pValues[j] <= 0.05) {
+                                       cout << i+1 << "\tOTU" << j+1 << '\t' << indicatorValues[j]  << '\t';
+                                       string pValueString = "<" + toString((1/(float)iters)); 
+                                       if (pValues[j] > (1/(float)iters)) { pValueString = toString(pValues[j]); cout << pValues[j];} 
+                                       else { cout << "<" << (1/(float)iters); }
+                                       m->mothurOutJustToLog(toString(i) + "\tOTU" + toString(j+1) + "\t" + toString(indicatorValues[j]) + "\t" + pValueString); 
+                                       m->mothurOutEndLine(); 
                                }
                        }
                        out << endl;
                        
                        T->tree[i].setLabel((i+1));
-                       
-                       cout << i+1 << "\tOTU" << indicatorOTU+1 << '\t' << maxValue << '\t';
-                       string pValueString = "<" + toString((1/(float)iters)); 
-                       if (pValues[indicatorOTU] > (1/(float)iters)) { pValueString = toString(pValues[indicatorOTU]); cout << pValues[indicatorOTU];} 
-                       else { cout << "<" << (1/(float)iters); }
-                       m->mothurOutJustToLog(toString(i) + "\tOTU" + toString(indicatorOTU+1) + "\t" + toString(maxValue) + "\t" + pValueString); 
-                       m->mothurOutEndLine(); 
-                       
-                       
                }
                out.close();
        

diff --git a/qualityscores.cpp b/qualityscores.cpp
index 520d0bb56eae9c66a936311ca71aa5330bf29aa1..ae85d86199c3ce9fbec09aac5525b17db3c4297a 100644 (file)
--- a/qualityscores.cpp
+++ b/qualityscores.cpp
@@ -286,23 +286,26 @@ bool QualityScores::stripQualWindowAverage(Sequence& sequence, int stepSize, int
 
                if(seqLength < windowSize) {    return 0;       }
                        
-               while(start < seqLength){
+               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 - 1;    }
+                               
+                       if(end >= seqLength){   end = seqLength;        }
+                               
                }
-               
+       
                if(end == -1){  end = seqLength;        }
                
                //failed first window

diff --git a/trimseqscommand.cpp b/trimseqscommand.cpp
index 082bc7950c590d9251100c3a6cb45dd18af6d2aa..a2c8021a57f6f3f076c8d1c9be7353eb8cffff93 100644 (file)
--- a/trimseqscommand.cpp
+++ b/trimseqscommand.cpp
@@ -604,7 +604,7 @@ int TrimSeqsCommand::driverCreateTrim(string filename, string qFileName, string
                                        else if(qRollAverage != 0)      {       success = currQual.stripQualRollingAverage(currSeq, qRollAverage);      }
                                        else if(qWindowAverage != 0){   success = currQual.stripQualWindowAverage(currSeq, qWindowStep, qWindowSize, qWindowAverage);   }
                                        else                                            {       success = 1;                            }
-                               
+                                       
                                        //you don't want to trim, if it fails above then scrap it
                                        if ((!qtrim) && (origLength != currSeq.getNumBases())) {  success = 0; }