finished dereplicate changes to chimera commands. added taxonomy file to subsample...

diff --git a/abstractrandomforest.cpp b/abstractrandomforest.cpp
index ae60b773742c25e0af85f811167003641e4646a4..59e3bafe5fc53094ad41853a48e3a85347e1afc4 100644 (file)
--- a/abstractrandomforest.cpp
+++ b/abstractrandomforest.cpp
@@ -55,4 +55,5 @@ vector<int> AbstractRandomForest::getGlobalDiscardedFeatureIndices() {
        } 
 }
 
-/***********************************************************************/
\ No newline at end of file
+/***********************************************************************/
+

diff --git a/chimeraperseuscommand.cpp b/chimeraperseuscommand.cpp
index 1835862f6b970ffdf7469002866f33a6682257ad..842a65ee22fd54c730143cbed94f91688dba8006 100644 (file)
--- a/chimeraperseuscommand.cpp
+++ b/chimeraperseuscommand.cpp
@@ -649,6 +649,11 @@ int ChimeraPerseusCommand::execute(){
                if (itTypes != outputTypes.end()) {
                        if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setAccnosFile(current); }
                }
+        
+        itTypes = outputTypes.find("count");
+               if (itTypes != outputTypes.end()) {
+                       if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setCountTableFile(current); }
+               }
                
                m->mothurOutEndLine();
                m->mothurOut("Output File Names: "); m->mothurOutEndLine();

diff --git a/chimeraperseuscommand.h b/chimeraperseuscommand.h
index 664200fd85e486ddca190acea16228aeaf76ddc6..4fa4cfff55e9582ffd699d5a66409717df8fbd36 100644 (file)
--- a/chimeraperseuscommand.h
+++ b/chimeraperseuscommand.h
@@ -362,11 +362,11 @@ static DWORD WINAPI MyPerseusThreadFunction(LPVOID lpParam){
                     if (pDataArray->hasCount) {
                         while (!in.eof()) {
                             in >> name; pDataArray->m->gobble(in);
-                            outCountList << name << '\t' << groups[u] << endl;
+                            outCountList << name << '\t' << pDataArray->groups[u] << endl;
                         }
                         in.close();
                     }else {
-                        map<string, string> thisnamemap = parser->getNameMap(groups[u]);
+                        map<string, string> thisnamemap = parser->getNameMap(pDataArray->groups[u]);
                         map<string, string>::iterator itN;
                         ofstream out;
                         pDataArray->m->openOutputFile(accnosFileName+".temp", out);

diff --git a/chimeraslayercommand.cpp b/chimeraslayercommand.cpp
index d8e0952a8b9f89cd237946a27ea624f11fa42898..fceef2180163f60e3c33cfecb7e861839c4bd1c7 100644 (file)
--- a/chimeraslayercommand.cpp
+++ b/chimeraslayercommand.cpp
@@ -837,6 +837,11 @@ int ChimeraSlayerCommand::execute(){
                        }
                }
                
+        itTypes = outputTypes.find("count");
+               if (itTypes != outputTypes.end()) {
+                       if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setCountTableFile(current); }
+               }
+        
                m->mothurOutEndLine();
                m->mothurOut("Output File Names: "); m->mothurOutEndLine();
                for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }       

diff --git a/chimerauchimecommand.cpp b/chimerauchimecommand.cpp
index 1d7e252cb702f58138ae0e756460bb98273b8908..e82dc1bc8dce5ffd0887025d5816b5eec8d0eab9 100644 (file)
--- a/chimerauchimecommand.cpp
+++ b/chimerauchimecommand.cpp
@@ -820,6 +820,11 @@ int ChimeraUchimeCommand::execute(){
                if (itTypes != outputTypes.end()) {
                        if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setAccnosFile(current); }
                }
+        
+        itTypes = outputTypes.find("count");
+               if (itTypes != outputTypes.end()) {
+                       if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setCountTableFile(current); }
+               }
                
                m->mothurOutEndLine();
                m->mothurOut("Output File Names: "); m->mothurOutEndLine();

diff --git a/cluster.cpp b/cluster.cpp
index 18133a29030cc47d10ab7a1c7f3572cec74e4c29..0a70fbfee99957de6f08315562f9f0d9ac14fee0 100644 (file)
--- a/cluster.cpp
+++ b/cluster.cpp
@@ -64,7 +64,7 @@ void Cluster::update(double& cutOFF){
         nRowCells = dMatrix->seqVec[smallRow].size();
         
                vector<int> foundCol(nColCells, 0);
-        //cout << dMatrix->getNNodes() << " small cell: " << smallRow << '\t' << smallCol << endl;  
+        //cout << dMatrix->getNNodes() << " small cell: " << smallRow << '\t' << smallCol << endl;
                int search;
                bool changed;
         
@@ -91,7 +91,8 @@ void Cluster::update(double& cutOFF){
                                //if not merged it you need it for warning 
                                if ((!merged) && (method == "average" || method == "weighted")) {  
                                        if (cutOFF > dMatrix->seqVec[smallRow][i].dist) {  
-                                               cutOFF = dMatrix->seqVec[smallRow][i].dist;  
+                                               cutOFF = dMatrix->seqVec[smallRow][i].dist;
+                        //cout << "changing cutoff to " << cutOFF << endl;
                                        }
                     
                                }

diff --git a/pairwiseseqscommand.cpp b/pairwiseseqscommand.cpp
index 767fdb52702bd1a3652d6bf2c1e2b244958e48ca..e937b1661b90b5f8e3340178222c21a9a9cd8d9e 100644 (file)
--- a/pairwiseseqscommand.cpp
+++ b/pairwiseseqscommand.cpp
@@ -643,7 +643,8 @@ int PairwiseSeqsCommand::driver(int startLine, int endLine, string dFileName, fl
                outFile << setprecision(4);
                
                if((output == "lt") && startLine == 0){ outFile << alignDB.getNumSeqs() << endl;        }
-               
+               int countSmall = 0;
+        int countAll = 0;
                for(int i=startLine;i<endLine;i++){
                        if(output == "lt")      {       
                                string name = alignDB.get(i).getName();
@@ -672,11 +673,18 @@ int PairwiseSeqsCommand::driver(int startLine, int endLine, string dFileName, fl
                                seqI.setAligned(alignment->getSeqAAln());
                                seqJ.setAligned(alignment->getSeqBAln());
 
-                               
+                               //cout << seqI.getName() << '\t' << seqJ.getName() << endl;
+                //cout << alignment->getSeqAAln() << endl << alignment->getSeqBAln() << endl;
+                
                                distCalculator->calcDist(seqI, seqJ);
                                double dist = distCalculator->getDist();
-                               
+                
+                //cout << "dist = " << dist << endl;
+                                               
                                if(dist <= cutoff){
+                    if (dist < 0.01) { countSmall++; }
+                    countAll++;
+
                                        if (output == "column") { outFile << alignDB.get(i).getName() << ' ' << alignDB.get(j).getName() << ' ' << dist << endl; }
                                }
                                if (output == "lt") {  outFile << dist << '\t'; }
@@ -690,7 +698,7 @@ int PairwiseSeqsCommand::driver(int startLine, int endLine, string dFileName, fl
                        
                }
                m->mothurOut(toString(endLine-1) + "\t" + toString(time(NULL) - startTime)); m->mothurOutEndLine();
-               
+               cout << "num less than 0.01 = " << countSmall << " of " << countAll << endl;
                outFile.close();
         delete alignment;
         delete distCalculator;

diff --git a/removeseqscommand.cpp b/removeseqscommand.cpp
index 5c0216366a98c83bc11c776ad31b0da1086120cd..b4fb467d6e395ce9b1655e043859cff0e2c77f67 100644 (file)
--- a/removeseqscommand.cpp
+++ b/removeseqscommand.cpp
@@ -813,7 +813,7 @@ int RemoveSeqsCommand::readTax(){
                while(!in.eof()){
                        if (m->control_pressed) { in.close();  out.close();  m->mothurRemove(outputFileName);  return 0; }
                        
-                       in >> name;                             //read from first column
+                       in >> name;     m->gobble(in);                  //read from first column
                        in >> tax;                      //read from second column
                        
             if (!dups) {//adjust name if needed

diff --git a/subsamplecommand.cpp b/subsamplecommand.cpp
index cfc7b1d9283779fdf9cc4a21a60c3f6d49d5177b..32745a25ee0ae66e6f6cbcb50cc4db4469cca21f 100644 (file)
--- a/subsamplecommand.cpp
+++ b/subsamplecommand.cpp
@@ -10,6 +10,7 @@
 #include "subsamplecommand.h"
 #include "sharedutilities.h"
 #include "deconvolutecommand.h"
+#include "getseqscommand.h"
 #include "subsample.h"
 
 //**********************************************************************************************************************
@@ -17,6 +18,7 @@ vector<string> SubSampleCommand::setParameters(){
        try {           
                CommandParameter pfasta("fasta", "InputTypes", "", "", "none", "FLSSR", "none","fasta",false,false,true); parameters.push_back(pfasta);
         CommandParameter pname("name", "InputTypes", "", "", "NameCount", "none", "none","name",false,false,true); parameters.push_back(pname);
+        CommandParameter ptaxonomy("taxonomy", "InputTypes", "", "", "none", "none", "none","taxonomy",false,false,true); parameters.push_back(ptaxonomy);
         CommandParameter pcount("count", "InputTypes", "", "", "NameCount-CountGroup", "none", "none","count",false,false,true); parameters.push_back(pcount);
                CommandParameter pgroup("group", "InputTypes", "", "", "CountGroup", "none", "none","group",false,false,true); parameters.push_back(pgroup);
                CommandParameter plist("list", "InputTypes", "", "", "none", "FLSSR", "none","list",false,false,true); parameters.push_back(plist);
@@ -44,7 +46,7 @@ string SubSampleCommand::getHelpString(){
        try {
                string helpString = "";
                helpString += "The sub.sample command is designed to be used as a way to normalize your data, or create a smaller set from your original set.\n";
-               helpString += "The sub.sample command parameters are fasta, name, list, group, count, rabund, sabund, shared, groups, size, persample and label.  You must provide a fasta, list, sabund, rabund or shared file as an input file.\n";
+               helpString += "The sub.sample command parameters are fasta, name, list, group, count, rabund, sabund, shared, taxonomy, groups, size, persample and label.  You must provide a fasta, list, sabund, rabund or shared file as an input file.\n";
                helpString += "The namefile is only used with the fasta file, not with the listfile, because the list file should contain all sequences.\n";
                helpString += "The groups parameter allows you to specify which of the groups in your groupfile you would like included. The group names are separated by dashes.\n";
                helpString += "The label parameter allows you to select what distance levels you would like, and are also separated by dashes.\n";
@@ -70,11 +72,12 @@ string SubSampleCommand::getOutputPattern(string type) {
         string pattern = "";
         
         if (type == "fasta")            {   pattern = "[filename],subsample,[extension]";    }
-        else if (type == "sabund")    {   pattern = "[filename],subsample,[extension]";    }
+        else if (type == "sabund")      {   pattern = "[filename],subsample,[extension]";    }
         else if (type == "name")        {   pattern = "[filename],subsample,[extension]";    }
         else if (type == "group")       {   pattern = "[filename],subsample,[extension]";    }
         else if (type == "count")       {   pattern = "[filename],subsample,[extension]";    }
         else if (type == "list")        {   pattern = "[filename],subsample,[extension]";    }
+        else if (type == "taxonomy")    {   pattern = "[filename],subsample,[extension]";    }
         else if (type == "shared")      {   pattern = "[filename],[distance],subsample,[extension]";    }
         else if (type == "rabund")      {   pattern = "[filename],subsample,[extension]";    }
         else { m->mothurOut("[ERROR]: No definition for type " + type + " output pattern.\n"); m->control_pressed = true;  }
@@ -100,6 +103,7 @@ SubSampleCommand::SubSampleCommand(){
                outputTypes["name"] = tempOutNames;
                outputTypes["group"] = tempOutNames;
         outputTypes["count"] = tempOutNames;
+        outputTypes["taxonomy"] = tempOutNames;
        }
        catch(exception& e) {
                m->errorOut(e, "SubSampleCommand", "GetRelAbundCommand");
@@ -140,6 +144,7 @@ SubSampleCommand::SubSampleCommand(string option) {
                        outputTypes["name"] = tempOutNames;
                        outputTypes["group"] = tempOutNames;
             outputTypes["count"] = tempOutNames;
+            outputTypes["taxonomy"] = tempOutNames;
                                        
                        //if the user changes the output directory command factory will send this info to us in the output parameter 
                        outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){  outputDir = ""; }
@@ -212,6 +217,14 @@ SubSampleCommand::SubSampleCommand(string option) {
                                        //if the user has not given a path then, add inputdir. else leave path alone.
                                        if (path == "") {       parameters["count"] = inputDir + it->second;            }
                                }
+                
+                it = parameters.find("taxonomy");
+                               //user has given a template file
+                               if(it != parameters.end()){
+                                       path = m->hasPath(it->second);
+                                       //if the user has not given a path then, add inputdir. else leave path alone.
+                                       if (path == "") {       parameters["taxonomy"] = inputDir + it->second;         }
+                               }
                        }
                        
                        //check for required parameters
@@ -249,6 +262,11 @@ SubSampleCommand::SubSampleCommand(string option) {
                        if (groupfile == "not open") { groupfile = ""; abort = true; }  
                        else if (groupfile == "not found") { groupfile = ""; }
                        else { m->setGroupFile(groupfile); }
+            
+            taxonomyfile = validParameter.validFile(parameters, "taxonomy", true);
+                       if (taxonomyfile == "not open") { taxonomyfile = ""; abort = true; }
+                       else if (taxonomyfile == "not found") { taxonomyfile = ""; }
+                       else { m->setTaxonomyFile(taxonomyfile); }
                        
             countfile = validParameter.validFile(parameters, "count", true);
                        if (countfile == "not open") { countfile = ""; abort = true; }
@@ -297,7 +315,10 @@ SubSampleCommand::SubSampleCommand(string option) {
                 }
             }
                        
-                       if ((namefile != "") && (fastafile == "")) { m->mothurOut("You may only use a namefile with a fastafile."); m->mothurOutEndLine(); abort = true; }
+                       if ((namefile != "") && ((fastafile == "") && (taxonomyfile == ""))) { m->mothurOut("You may only use a name file with a fasta file or taxonomy file."); m->mothurOutEndLine(); abort = true; }
+            
+            if ((taxonomyfile != "") && ((fastafile == "") && (listfile == ""))) { m->mothurOut("You may only use a taxonomyfile with a fastafile or listfile."); m->mothurOutEndLine(); abort = true; }
+            
                        
                        if ((fastafile == "") && (listfile == "") && (sabundfile == "") && (rabundfile == "") && (sharedfile == "")) {
                                m->mothurOut("You must provide a fasta, list, sabund, rabund or shared file as an input file."); m->mothurOutEndLine(); abort = true; }
@@ -389,6 +410,10 @@ int SubSampleCommand::execute(){
                        if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setCountTableFile(current); }
                }
                
+        itTypes = outputTypes.find("taxonomy");
+               if (itTypes != outputTypes.end()) {
+                       if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setTaxonomyFile(current); }
+               }
                
                m->mothurOutEndLine();
                m->mothurOut("Output File Names: "); m->mothurOutEndLine();
@@ -632,6 +657,7 @@ int SubSampleCommand::getSubSampleFasta() {
                }
                in.close();     
                out.close();
+        
                
                if (count != subset.size()) {
                        m->mothurOut("[ERROR]: The subset selected contained " + toString(subset.size()) + " sequences, but I only found " + toString(count) + " of those in the fastafile."); m->mothurOutEndLine();
@@ -665,7 +691,33 @@ int SubSampleCommand::getSubSampleFasta() {
                        m->mothurOut("/******************************************/"); m->mothurOutEndLine(); 
                        
                        m->mothurOut("Done."); m->mothurOutEndLine();
-               }
+            
+            if (taxonomyfile != "") {
+                set<string> tempSubset;
+                //get new unique names from fasta file
+                //read through fasta file outputting only the names on the subsample list after deconvolute
+                ifstream in2;
+                m->openInputFile(outputFileName, in2);
+                
+                while (!in2.eof()) {
+                    Sequence seq(in2); m->gobble(in2);
+                    if (seq.getName() != "") {
+                        tempSubset.insert(seq.getName());
+                    }
+                }
+                in2.close();
+                
+                //send that list to getTax
+                int tcount = getTax(tempSubset);
+                if (tcount != tempSubset.size()) { m->mothurOut("[ERROR]: subsampled fasta file contains " + toString(tempSubset.size()) + " sequences, but I only found " + toString(tcount) + " in your taxonomy file, please correct."); m->mothurOutEndLine(); }
+            }
+               }else {
+            if (taxonomyfile != "") {
+                int tcount = getTax(subset);
+                if (tcount != subset.size()) { m->mothurOut("[ERROR]: subsampled fasta file contains " + toString(subset.size()) + " sequences, but I only found " + toString(tcount) + " in your taxonomy file, please correct."); m->mothurOutEndLine(); }
+            
+            }  //should only contain uniques.
+        }
                
                outputTypes["fasta"].push_back(outputFileName);  outputNames.push_back(outputFileName);
                
@@ -938,7 +990,9 @@ int SubSampleCommand::processShared(vector<SharedRAbundVector*>& thislookup) {
 //**********************************************************************************************************************
 int SubSampleCommand::getSubSampleList() {
        try {
-               
+        
+               if (namefile != "") { m->readNames(namefile, nameMap); }
+        
                string thisOutputDir = outputDir;
                if (outputDir == "") {  thisOutputDir += m->hasPath(listfile);  }
                map<string, string> variables; 
@@ -1208,6 +1262,63 @@ int SubSampleCommand::getSubSampleList() {
                out.close();  
                if (list != NULL) { delete list; }
                delete input;
+        
+        if (taxonomyfile != "") {
+            if (namefile == "") {
+                //fake nameMap
+                for (set<string>::iterator it = subset.begin(); it != subset.end(); it++) {
+                    vector<string> temp; temp.push_back(*it);
+                    nameMap[*it] = temp;
+                }
+                int tcount = getTax(subset);
+                if (tcount != subset.size()) { m->mothurOut("[ERROR]: subsampled list file contains " + toString(subset.size()) + " sequences, but I only found " + toString(tcount) + " in your taxonomy file, did you forget a name file? Please correct."); m->mothurOutEndLine(); }
+            }else {
+                string tempAccnos = "temp.accnos";
+                ofstream outAccnos;
+                m->openOutputFile(tempAccnos, outAccnos);
+                for (set<string>::iterator it = subset.begin(); it != subset.end(); it++) { outAccnos << *it << endl; }
+                outAccnos.close();
+                
+                m->mothurOut("Sampling taxonomy and name file... "); m->mothurOutEndLine();
+                string thisNameOutputDir = outputDir;
+                if (outputDir == "") {  thisNameOutputDir += m->hasPath(namefile);  }
+                map<string, string> variables;
+                variables["[filename]"] = thisNameOutputDir + m->getRootName(m->getSimpleName(namefile));
+                variables["[extension]"] = m->getExtension(namefile);
+                string outputNameFileName = getOutputFileName("name", variables);
+                
+                string thisTaxOutputDir = outputDir;
+                if (outputDir == "") {  thisTaxOutputDir += m->hasPath(taxonomyfile);  }
+                variables["[filename]"] = thisTaxOutputDir + m->getRootName(m->getSimpleName(taxonomyfile));
+                variables["[extension]"] = m->getExtension(taxonomyfile);
+                string outputTaxFileName = getOutputFileName("taxonomy", variables);
+                
+                
+                //use unique.seqs to create new name and fastafile
+                string inputString = "dups=f, name=" + namefile + ", taxonomy=" + taxonomyfile + ", accnos=" + tempAccnos;
+                m->mothurOut("/******************************************/"); m->mothurOutEndLine();
+                m->mothurOut("Running command: get.seqs(" + inputString + ")"); m->mothurOutEndLine();
+                m->mothurCalling = true;
+                
+                Command* getCommand = new GetSeqsCommand(inputString);
+                getCommand->execute();
+                
+                map<string, vector<string> > filenames = getCommand->getOutputFiles();
+                
+                delete getCommand;
+                m->mothurCalling = false;
+                
+                m->renameFile(filenames["name"][0], outputNameFileName);
+                m->renameFile(filenames["taxonomy"][0], outputTaxFileName);
+                
+                outputTypes["name"].push_back(outputNameFileName);  outputNames.push_back(outputNameFileName);
+                outputNames.push_back(outputTaxFileName); outputTypes["taxonomy"].push_back(outputTaxFileName);
+                
+                m->mothurOut("/******************************************/"); m->mothurOutEndLine(); 
+                
+                m->mothurOut("Done."); m->mothurOutEndLine();
+            }
+        }
                                                
                return 0;
  
@@ -1579,7 +1690,63 @@ int SubSampleCommand::processSabund(SAbundVector*& sabund, ofstream& out) {
                m->errorOut(e, "SubSampleCommand", "processSabund");
                exit(1);
        }
-}                      
+}
+
+//**********************************************************************************************************************
+int SubSampleCommand::getTax(set<string>& subset) {
+       try {
+
+        string thisTaxOutputDir = outputDir;
+        if (outputDir == "") {  thisTaxOutputDir += m->hasPath(taxonomyfile);  }
+        map<string, string> variables;
+        variables["[filename]"] = thisTaxOutputDir + m->getRootName(m->getSimpleName(taxonomyfile));
+        variables["[extension]"] = m->getExtension(taxonomyfile);
+        string outputTaxFileName = getOutputFileName("taxonomy", variables);
+        ofstream outTax;
+        m->openOutputFile(outputTaxFileName, outTax);
+        outputNames.push_back(outputTaxFileName); outputTypes["taxonomy"].push_back(outputTaxFileName);
+        
+        //read through fasta file outputting only the names on the subsample list
+        ifstream inTax;
+        m->openInputFile(taxonomyfile, inTax);
+        
+        string tname, tax;
+        int tcount = 0;
+        map<string, vector<string> >::iterator itNameMap;
+        
+        while(!inTax.eof()){
+            
+            if (m->control_pressed) { inTax.close(); outTax.close();  return 0; }
+            
+            inTax >> tname;    m->gobble(inTax);                       //read from first column
+            inTax >> tax;      m->gobble(inTax);               //read from second column
+            
+            //does the subset contain a sequence that this sequence represents
+            itNameMap = nameMap.find(tname);
+            if (itNameMap != nameMap.end()) {
+                vector<string> nameRepresents = itNameMap->second;
+                
+            
+                for (int i = 0; i < nameRepresents.size(); i++){
+                    if (subset.count(nameRepresents[i]) != 0) {
+                        outTax << nameRepresents[i] << '\t' << tax << endl;
+                        tcount++;
+                       
+                    }
+                }
+            }else{ m->mothurOut("[ERROR]: " + tname + " is missing, please correct."); m->mothurOutEndLine(); }
+        }
+        inTax.close();
+        outTax.close();
+        
+        return tcount;
+    }
+    catch(exception& e) {
+        m->errorOut(e, "SubSampleCommand", "getTax");
+        exit(1);
+    }
+}
+
 //**********************************************************************************************************************
 
 

diff --git a/subsamplecommand.h b/subsamplecommand.h
index 7bafd9bfb573a03527b95d179130fc066e9234bb..7b251eee11c2fa7ad5be7baa6cea4427c90265c7 100644 (file)
--- a/subsamplecommand.h
+++ b/subsamplecommand.h
@@ -40,7 +40,7 @@ public:
        
 private:       
        bool abort, pickedGroups, allLines, persample;
-       string listfile, groupfile, countfile, sharedfile, rabundfile, sabundfile, fastafile, namefile;
+       string listfile, groupfile, countfile, sharedfile, rabundfile, sabundfile, fastafile, namefile, taxonomyfile;
        set<string> labels; //holds labels to be used
        string groups, label, outputDir;
        vector<string> Groups, outputNames;
@@ -60,6 +60,7 @@ private:
        int processList(ListVector*&, ofstream&, set<string>&);
        int getNames();
        int readNames();
+    int getTax(set<string>&);
        
 };