fixed seq. error align=f issue. added std and ave calcs to mothurout. working on...

diff --git a/mothurout.cpp b/mothurout.cpp
index 468c063cb5c7e74a8d5901a1a38203f47409e2d5..240e7ec2e5ac24ec3e637b78ae28dd7608a0389d 100644 (file)
--- a/mothurout.cpp
+++ b/mothurout.cpp
@@ -2925,6 +2925,79 @@ bool MothurOut::checkReleaseVersion(ifstream& file, string version) {
        }
 }
 /**************************************************************************************************/
+vector<double> MothurOut::getAverages(vector< vector<double> >& dists) {
+       try{
+        vector<double> averages; //averages.resize(numComp, 0.0);
+        for (int i = 0; i < dists[0].size(); i++) { averages.push_back(0.0); }
+      
+        for (int thisIter = 0; thisIter < dists.size(); thisIter++) {
+            for (int i = 0; i < dists[thisIter].size(); i++) {  
+                averages[i] += dists[thisIter][i];
+            }
+        }
+        
+        //finds average.
+        for (int i = 0; i < averages.size(); i++) {  averages[i] /= (double) dists.size(); }
+        
+        return averages;
+    }
+       catch(exception& e) {
+               errorOut(e, "MothurOut", "getAverages");                
+               exit(1);
+       }
+}
+/**************************************************************************************************/
+vector<double> MothurOut::getStandardDeviation(vector< vector<double> >& dists) {
+       try{
+        
+        vector<double> averages = getAverages(dists);
+        
+        //find standard deviation
+        vector<double> stdDev; //stdDev.resize(numComp, 0.0);
+        for (int i = 0; i < dists[0].size(); i++) { stdDev.push_back(0.0); }
+        
+        for (int thisIter = 0; thisIter < dists.size(); thisIter++) { //compute the difference of each dist from the mean, and square the result of each
+            for (int j = 0; j < dists[thisIter].size(); j++) {
+                stdDev[j] += ((dists[thisIter][j] - averages[j]) * (dists[thisIter][j] - averages[j]));
+            }
+        }
+        for (int i = 0; i < stdDev.size(); i++) {  
+            stdDev[i] /= (double) dists.size(); 
+            stdDev[i] = sqrt(stdDev[i]);
+        }
+        
+        return stdDev;
+    }
+       catch(exception& e) {
+               errorOut(e, "MothurOut", "getAverages");                
+               exit(1);
+       }
+}
+/**************************************************************************************************/
+vector<double> MothurOut::getStandardDeviation(vector< vector<double> >& dists, vector<double>& averages) {
+       try{
+        //find standard deviation
+        vector<double> stdDev; //stdDev.resize(numComp, 0.0);
+        for (int i = 0; i < dists[0].size(); i++) { stdDev.push_back(0.0); }
+        
+        for (int thisIter = 0; thisIter < dists.size(); thisIter++) { //compute the difference of each dist from the mean, and square the result of each
+            for (int j = 0; j < dists[thisIter].size(); j++) {
+                stdDev[j] += ((dists[thisIter][j] - averages[j]) * (dists[thisIter][j] - averages[j]));
+            }
+        }
+        for (int i = 0; i < stdDev.size(); i++) {  
+            stdDev[i] /= (double) dists.size(); 
+            stdDev[i] = sqrt(stdDev[i]);
+        }
+        
+        return stdDev;
+    }
+       catch(exception& e) {
+               errorOut(e, "MothurOut", "getAverages");                
+               exit(1);
+       }
+}
+/**************************************************************************************************/
 bool MothurOut::isContainingOnlyDigits(string input) {
        try{
                

diff --git a/mothurout.h b/mothurout.h
index 0d2e86f24cb258b24e78367e1427fc7ee2bc7d31..84f9be31f82ebc2c27aa19563fb5366519e1cba8 100644 (file)
--- a/mothurout.h
+++ b/mothurout.h
@@ -155,6 +155,9 @@ class MothurOut {
                unsigned int fromBase36(string);
                int getRandomIndex(int); //highest
         double getStandardDeviation(vector<int>&);
+        vector<double> getStandardDeviation(vector< vector<double> >&);
+        vector<double> getStandardDeviation(vector< vector<double> >&, vector<double>&);
+        vector<double> getAverages(vector< vector<double> >&);
 
                int control_pressed;
                bool executing, runParse, jumble, gui, mothurCalling, debug;

diff --git a/refchimeratest.cpp b/refchimeratest.cpp
index 701fb9eb6a711fd82108d6694790b02992bc6816..c084bffd40c94ef9c14a43d05de39e768dc96bea 100644 (file)
--- a/refchimeratest.cpp
+++ b/refchimeratest.cpp
@@ -24,7 +24,8 @@ RefChimeraTest::RefChimeraTest(vector<Sequence>& refs, bool aligned) : aligned(a
        referenceSeqs.resize(numRefSeqs);
        referenceNames.resize(numRefSeqs);
        for(int i=0;i<numRefSeqs;i++){
-               referenceSeqs[i] = refs[i].getAligned();
+               if (aligned) { referenceSeqs[i] = refs[i].getAligned(); }
+        else { referenceSeqs[i] = refs[i].getUnaligned(); }
                referenceNames[i] = refs[i].getName();
        }
        

diff --git a/seqerrorcommand.cpp b/seqerrorcommand.cpp
index 1fe60e8954acd14dc1a1e02d37c0d8bd987cc1c0..67e43aa0d4b7f6a6ff251bc45a8254ca49b8b162 100644 (file)
--- a/seqerrorcommand.cpp
+++ b/seqerrorcommand.cpp
@@ -298,7 +298,7 @@ SeqErrorCommand::SeqErrorCommand(string option)  {
                        }
             else{
                 if(reportFileName != ""){
-                    m->mothurOut("we are ignoring the report file if your sequences are not aligned.  we will check that the sequences in your fasta and and qual fileare the same length.");
+                    m->mothurOut("we are ignoring the report file if your sequences are not aligned.  we will check that the sequences in your fasta and and qual file are the same length.");
                     m->mothurOutEndLine();
                 }
             }
@@ -759,7 +759,10 @@ int SeqErrorCommand::driver(string filename, string qFileName, string rFileName,
             int numParentSeqs = -1;
             int closestRefIndex = -1;
                         
-            numParentSeqs = chimeraTest.analyzeQuery(query.getName(), query.getAligned(), outChimeraReport);
+            string querySeq = query.getAligned();
+            if (!aligned) {  querySeq = query.getUnaligned();  }
+            
+            numParentSeqs = chimeraTest.analyzeQuery(query.getName(), querySeq, outChimeraReport);
             
             closestRefIndex = chimeraTest.getClosestRefIndex();
             

diff --git a/unifracunweightedcommand.cpp b/unifracunweightedcommand.cpp
index 5f3cffc2b1ab807b279ba3185444dc96a0a3b863..19a2ece3bb4253e6da6199515d9c203d2d656400 100644 (file)
--- a/unifracunweightedcommand.cpp
+++ b/unifracunweightedcommand.cpp
@@ -481,20 +481,29 @@ int UnifracUnweightedCommand::execute() {
 int UnifracUnweightedCommand::getAverageSTDMatrices(vector< vector<double> >& dists, int treeNum) {
        try {
         //we need to find the average distance and standard deviation for each groups distance
-        
         //finds sum
-        vector<double> averages; averages.resize(numComp, 0); 
+        vector<double> averages; //averages.resize(numComp, 0.0);
+        for (int i = 0; i < numComp; i++) { averages.push_back(0.0); }
+        
+        if (m->debug) { m->mothurOut("[DEBUG]: numcomparisons = " + toString(numComp) + ", subsampleIters = " + toString(subsampleIters) + "\n"); }
+        
         for (int thisIter = 0; thisIter < subsampleIters; thisIter++) {
             for (int i = 0; i < dists[thisIter].size(); i++) {  
                 averages[i] += dists[thisIter][i];
             }
         }
         
+        if (m->debug) { m->mothurOut("[DEBUG]: numcomparisons = " + toString(numComp) + ", subsampleIters = " + toString(subsampleIters) + "\n"); }
+        
         //finds average.
-        for (int i = 0; i < averages.size(); i++) {  averages[i] /= (float) subsampleIters; }
+        for (int i = 0; i < averages.size(); i++) {  
+            averages[i] /= (float) subsampleIters; 
+            if (m->debug) { m->mothurOut("[DEBUG]: i = " + toString(i) + ", averages[i] = " + toString(averages[i]) + "\n"); }
+        }
         
         //find standard deviation
-        vector<double> stdDev; stdDev.resize(numComp, 0);
+        vector<double> stdDev; //stdDev.resize(numComp, 0.0);
+        for (int i = 0; i < numComp; i++) { stdDev.push_back(0.0); }
         
         for (int thisIter = 0; thisIter < subsampleIters; thisIter++) { //compute the difference of each dist from the mean, and square the result of each
             for (int j = 0; j < dists[thisIter].size(); j++) {
@@ -504,20 +513,28 @@ int UnifracUnweightedCommand::getAverageSTDMatrices(vector< vector<double> >& di
         for (int i = 0; i < stdDev.size(); i++) {  
             stdDev[i] /= (float) subsampleIters; 
             stdDev[i] = sqrt(stdDev[i]);
+            if (m->debug) { m->mothurOut("[DEBUG]: i = " + toString(i) + ", stdDev[i] = " + toString(stdDev[i]) + "\n"); }
         }
         
         //make matrix with scores in it
-        vector< vector<double> > avedists;     avedists.resize(m->getNumGroups());
+        vector< vector<double> > avedists;     //avedists.resize(m->getNumGroups());
         for (int i = 0; i < m->getNumGroups(); i++) {
-            avedists[i].resize(m->getNumGroups(), 0.0);
+            vector<double> temp;
+            for (int j = 0; j < m->getNumGroups(); j++) { temp.push_back(0.0); }
+            avedists.push_back(temp);
         }
         
         //make matrix with scores in it
-        vector< vector<double> > stddists;     stddists.resize(m->getNumGroups());
+        vector< vector<double> > stddists;     //stddists.resize(m->getNumGroups());
         for (int i = 0; i < m->getNumGroups(); i++) {
-            stddists[i].resize(m->getNumGroups(), 0.0);
+            vector<double> temp;
+            for (int j = 0; j < m->getNumGroups(); j++) { temp.push_back(0.0); }
+            //stddists[i].resize(m->getNumGroups(), 0.0);
+            stddists.push_back(temp);
         }
         
+        if (m->debug) { m->mothurOut("[DEBUG]: about to fill matrix.\n"); }
+        
         //flip it so you can print it
         int count = 0;
         for (int r=0; r<m->getNumGroups(); r++) { 
@@ -530,6 +547,8 @@ int UnifracUnweightedCommand::getAverageSTDMatrices(vector< vector<double> >& di
             }
         }
         
+        if (m->debug) { m->mothurOut("[DEBUG]: done filling matrix.\n"); }
+        
         map<string, string> variables; 
                variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(treefile));
         variables["[tag]"] = toString(treeNum+1);