]> git.donarmstrong.com Git - mothur.git/blobdiff - removerarecommand.cpp
added multiple processors option for Windows users to align.seqs, dist.seqs, summary...
[mothur.git] / removerarecommand.cpp
index 4ba493fb2c7daff5ccc400f48a4fc97047774385..351ace6db7ea98c033514086004091e6352f663e 100644 (file)
 #include "inputdata.h"
 
 //**********************************************************************************************************************
-vector<string> RemoveRareCommand::getValidParameters(){        
+vector<string> RemoveRareCommand::setParameters(){     
        try {
-               string Array[] =  {"rabund","sabund", "group", "list", "shared", "nseqs","groups","label","outputdir","inputdir"};
-               vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
+               CommandParameter plist("list", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(plist);
+               CommandParameter prabund("rabund", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(prabund);
+               CommandParameter psabund("sabund", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(psabund);
+               CommandParameter pshared("shared", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(pshared);
+               CommandParameter pgroup("group", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(pgroup);
+               CommandParameter pgroups("groups", "String", "", "", "", "", "",false,false); parameters.push_back(pgroups);
+               CommandParameter plabel("label", "String", "", "", "", "", "",false,false); parameters.push_back(plabel);
+               CommandParameter pnseqs("nseqs", "Number", "", "0", "", "", "",false,true); parameters.push_back(pnseqs);
+               CommandParameter pbygroup("bygroup", "Boolean", "", "f", "", "", "",false,true); parameters.push_back(pbygroup);
+               CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
+               CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
+               
+               vector<string> myArray;
+               for (int i = 0; i < parameters.size(); i++) {   myArray.push_back(parameters[i].name);          }
                return myArray;
        }
        catch(exception& e) {
-               m->errorOut(e, "RemoveRareCommand", "getValidParameters");
+               m->errorOut(e, "RemoveRareCommand", "setParameters");
+               exit(1);
+       }
+}
+//**********************************************************************************************************************
+string RemoveRareCommand::getHelpString(){     
+       try {
+               string helpString = "";
+               helpString += "The remove.rare command parameters are list, rabund, sabund, shared, group, label, groups, bygroup and nseqs.\n";
+               helpString += "The remove.rare command reads one of the following file types: list, rabund, sabund or shared file. It outputs a new file after removing the rare otus.\n";
+               helpString += "The groups parameter allows you to specify which of the groups you would like analyzed.  Default=all. You may separate group names with dashes.\n";
+               helpString += "The label parameter is used to analyze specific labels in your input. default=all. You may separate label names with dashes.\n";
+               helpString += "The bygroup parameter is only valid with the shared file. default=f, meaning remove any OTU that has nseqs or fewer sequences across all groups.\n";
+               helpString += "bygroups=T means remove any OTU that has nseqs or fewer sequences in each group (if groupA has 1 sequence and group B has 100 sequences in OTUZ and nseqs=1, then set the groupA count for OTUZ to 0 and keep groupB's count at 100.) \n";
+               helpString += "The nseqs parameter allows you to set the cutoff for an otu to be deemed rare. It is required.\n";
+               helpString += "The remove.rare command should be in the following format: remove.rare(shared=yourSharedFile, nseqs=yourRareCutoff).\n";
+               helpString += "Example remove.rare(shared=amazon.fn.shared, nseqs=2).\n";
+               helpString += "Note: No spaces between parameter labels (i.e. shared), '=' and parameters (i.e.yourSharedFile).\n";
+               return helpString;
+       }
+       catch(exception& e) {
+               m->errorOut(e, "RemoveRareCommand", "getHelpString");
                exit(1);
        }
 }
 //**********************************************************************************************************************
 RemoveRareCommand::RemoveRareCommand(){        
        try {
-               abort = true;
-               //initialize outputTypes
+               abort = true; calledHelp = true; 
+               setParameters();
                vector<string> tempOutNames;
                outputTypes["rabund"] = tempOutNames;
                outputTypes["sabund"] = tempOutNames;
@@ -43,41 +76,17 @@ RemoveRareCommand::RemoveRareCommand(){
        }
 }
 //**********************************************************************************************************************
-vector<string> RemoveRareCommand::getRequiredParameters(){     
-       try {
-               string Array[] =  {"nseqs"};
-               vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
-               return myArray;
-       }
-       catch(exception& e) {
-               m->errorOut(e, "RemoveRareCommand", "getRequiredParameters");
-               exit(1);
-       }
-}
-//**********************************************************************************************************************
-vector<string> RemoveRareCommand::getRequiredFiles(){  
-       try {
-               vector<string> myArray;
-               return myArray;
-       }
-       catch(exception& e) {
-               m->errorOut(e, "RemoveRareCommand", "getRequiredFiles");
-               exit(1);
-       }
-}
-//**********************************************************************************************************************
 RemoveRareCommand::RemoveRareCommand(string option)  {
        try {
-               abort = false;
+               abort = false; calledHelp = false;   
                allLines = 1;
                
                //allow user to run help
-               if(option == "help") { help(); abort = true; }
+               if(option == "help") { help(); abort = true; calledHelp = true; }
+               else if(option == "citation") { citation(); abort = true; calledHelp = true;}
                
                else {
-                       //valid paramters for this command
-                       string Array[] =  {"rabund","sabund", "group", "list", "shared", "nseqs","groups","label","outputdir","inputdir"};
-                       vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
+                       vector<string> myArray = setParameters();
                        
                        OptionParser parser(option);
                        map<string,string> parameters = parser.getParameters();
@@ -151,25 +160,52 @@ RemoveRareCommand::RemoveRareCommand(string option)  {
                        //check for file parameters
                        listfile = validParameter.validFile(parameters, "list", true);
                        if (listfile == "not open") { abort = true; }
-                       else if (listfile == "not found") {  listfile = "";  }  
+                       else if (listfile == "not found") {  listfile = "";  }
+                       else { m->setListFile(listfile); }
                        
                        sabundfile = validParameter.validFile(parameters, "sabund", true);
                        if (sabundfile == "not open") { abort = true; }
                        else if (sabundfile == "not found") {  sabundfile = "";  }      
+                       else { m->setSabundFile(sabundfile); }
                        
                        rabundfile = validParameter.validFile(parameters, "rabund", true);
                        if (rabundfile == "not open") { abort = true; }
                        else if (rabundfile == "not found") {  rabundfile = "";  }                              
+                       else { m->setRabundFile(rabundfile); }
                        
                        groupfile = validParameter.validFile(parameters, "group", true);
-                       if (groupfile == "not open") { abort = true; }
+                       if (groupfile == "not open") { groupfile = ""; abort = true; }
                        else if (groupfile == "not found") {  groupfile = "";  }        
+                       else { m->setGroupFile(groupfile); }
                        
                        sharedfile = validParameter.validFile(parameters, "shared", true);
-                       if (sharedfile == "not open") { abort = true; }
+                       if (sharedfile == "not open") { sharedfile = "";  abort = true; }
                        else if (sharedfile == "not found") {  sharedfile = "";  }
-                       
-                       if ((sabundfile == "") && (rabundfile == "")  && (sharedfile == "") && (listfile == ""))  { m->mothurOut("You must provide at least one of the following: rabund, sabund, shared or list."); m->mothurOutEndLine(); abort = true; }
+                       else { m->setSharedFile(sharedfile); }
+                       
+                       if ((sharedfile == "") && (listfile == "") && (rabundfile == "") && (sabundfile == "")) { 
+                               //is there are current file available for any of these?
+                               //give priority to shared, then list, then rabund, then sabund
+                               //if there is a current shared file, use it
+                               sharedfile = m->getSharedFile(); 
+                               if (sharedfile != "") {  m->mothurOut("Using " + sharedfile + " as input file for the shared parameter."); m->mothurOutEndLine(); }
+                               else { 
+                                       listfile = m->getListFile(); 
+                                       if (listfile != "") {  m->mothurOut("Using " + listfile + " as input file for the list parameter."); m->mothurOutEndLine(); }
+                                       else { 
+                                               rabundfile = m->getRabundFile(); 
+                                               if (rabundfile != "") {  m->mothurOut("Using " + rabundfile + " as input file for the rabund parameter."); m->mothurOutEndLine(); }
+                                               else { 
+                                                       sabundfile = m->getSabundFile(); 
+                                                       if (sabundfile != "") {  m->mothurOut("Using " + sabundfile + " as input file for the sabund parameter."); m->mothurOutEndLine(); }
+                                                       else { 
+                                                               m->mothurOut("No valid current files. You must provide a list, sabund, rabund or shared file."); m->mothurOutEndLine(); 
+                                                               abort = true;
+                                                       }
+                                               }
+                                       }
+                               }
+                       } 
                        
                        groups = validParameter.validFile(parameters, "groups", false);                 
                        if (groups == "not found") { groups = "all"; }
@@ -185,6 +221,13 @@ RemoveRareCommand::RemoveRareCommand(string option)  {
                        string temp = validParameter.validFile(parameters, "nseqs", false);      
                        if (temp == "not found") { m->mothurOut("nseqs is a required parameter."); m->mothurOutEndLine(); abort = true; }
                        else { convert(temp, nseqs); }
+                       
+                       temp = validParameter.validFile(parameters, "bygroup", false);   if (temp == "not found") { temp = "f"; }
+                       byGroup = m->isTrue(temp);
+                       
+                       if (byGroup && (sharedfile == "")) { m->mothurOut("The byGroup parameter is only valid with a shared file."); m->mothurOutEndLine(); }
+                       
+                       if ((groupfile != "") && (listfile == "")) { m->mothurOut("A groupfile is only valid with a list file."); m->mothurOutEndLine(); groupfile = ""; }
                }
                
        }
@@ -195,28 +238,10 @@ RemoveRareCommand::RemoveRareCommand(string option)  {
 }
 //**********************************************************************************************************************
 
-void RemoveRareCommand::help(){
-       try {
-               //m->mothurOut("The remove.seqs command reads an .accnos file and at least one of the following file types: fasta, name, group, list, taxonomy, quality or alignreport file.\n");
-               //m->mothurOut("It outputs a file containing the sequences NOT in the .accnos file.\n");
-               //m->mothurOut("The remove.seqs command parameters are accnos, fasta, name, group, list, taxonomy, qfile, alignreport and dups.  You must provide accnos and at least one of the file parameters.\n");
-               //m->mothurOut("The dups parameter allows you to remove the entire line from a name file if you remove any name from the line. default=true. \n");
-               //m->mothurOut("The remove.seqs command should be in the following format: remove.seqs(accnos=yourAccnos, fasta=yourFasta).\n");
-               //m->mothurOut("Example remove.seqs(accnos=amazon.accnos, fasta=amazon.fasta).\n");
-               //m->mothurOut("Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFasta).\n\n");
-       }
-       catch(exception& e) {
-               m->errorOut(e, "RemoveRareCommand", "help");
-               exit(1);
-       }
-}
-
-//**********************************************************************************************************************
-
 int RemoveRareCommand::execute(){
        try {
                
-               if (abort == true) { return 0; }
+               if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
                
                if (m->control_pressed) { return 0; }
                
@@ -224,15 +249,42 @@ int RemoveRareCommand::execute(){
                if (sabundfile != "")           {               processSabund();        }
                if (rabundfile != "")           {               processRabund();        }
                if (listfile != "")                     {               processList();          }
-               //if (sharedfile != "")         {               processShared();        }
+               if (sharedfile != "")           {               processShared();        }
                
-               if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) {        remove(outputNames[i].c_str()); } return 0; }
+               if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) {        m->mothurRemove(outputNames[i]); } return 0; }
                        
                if (outputNames.size() != 0) {
                        m->mothurOutEndLine();
                        m->mothurOut("Output File Names: "); m->mothurOutEndLine();
                        for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
                        m->mothurOutEndLine();
+                       
+                       //set rabund file as new current rabundfile
+                       string current = "";
+                       itTypes = outputTypes.find("rabund");
+                       if (itTypes != outputTypes.end()) {
+                               if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setRabundFile(current); }
+                       }
+                       
+                       itTypes = outputTypes.find("sabund");
+                       if (itTypes != outputTypes.end()) {
+                               if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setSabundFile(current); }
+                       }
+                       
+                       itTypes = outputTypes.find("group");
+                       if (itTypes != outputTypes.end()) {
+                               if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setGroupFile(current); }
+                       }
+                       
+                       itTypes = outputTypes.find("list");
+                       if (itTypes != outputTypes.end()) {
+                               if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setListFile(current); }
+                       }
+                       
+                       itTypes = outputTypes.find("shared");
+                       if (itTypes != outputTypes.end()) {
+                               if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setSharedFile(current); }
+                       }
                }
                
                return 0;               
@@ -301,7 +353,8 @@ int RemoveRareCommand::processList(){
                if (groupfile != "") { 
                        groupMap = new GroupMap(groupfile); groupMap->readMap(); 
                        SharedUtil util;
-                       util.setGroups(Groups, groupMap->namesOfGroups);
+                       vector<string> namesGroups = groupMap->getNamesOfGroups();
+                       util.setGroups(Groups, namesGroups);
                        m->openOutputFile(outputGroupFileName, outGroup);
                }
                
@@ -313,7 +366,7 @@ int RemoveRareCommand::processList(){
                        
                        //for each bin
                        for (int i = 0; i < list->getNumBins(); i++) {
-                               if (m->control_pressed) {  if (groupfile != "") { delete groupMap; outGroup.close(); remove(outputGroupFileName.c_str()); } out.close();  remove(outputFileName.c_str());  return 0; }
+                               if (m->control_pressed) {  if (groupfile != "") { delete groupMap; outGroup.close(); m->mothurRemove(outputGroupFileName); } out.close();  m->mothurRemove(outputFileName);  return 0; }
                                
                                //parse out names that are in accnos file
                                string binnames = list->get(i);
@@ -564,12 +617,170 @@ int RemoveRareCommand::processRabund(){
                
                return 0;
        }
+       catch(exception& e) {
+               m->errorOut(e, "RemoveRareCommand", "processRabund");
+               exit(1);
+       }
+}
+//**********************************************************************************************************************
+int RemoveRareCommand::processShared(){
+       try {
+               m->setGroups(Groups);
+               
+               string thisOutputDir = outputDir;
+               if (outputDir == "") {  thisOutputDir += m->hasPath(sharedfile);  }
+               string outputFileName = thisOutputDir + m->getRootName(m->getSimpleName(sharedfile)) + "pick" +  m->getExtension(sharedfile);
+               outputTypes["shared"].push_back(outputFileName); outputNames.push_back(outputFileName);
+               
+               ofstream out;
+               m->openOutputFile(outputFileName, out);
+               
+               //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
+               InputData input(sharedfile, "sharedfile");
+               vector<SharedRAbundVector*> lookup = input.getSharedRAbundVectors();
+               string lastLabel = lookup[0]->getLabel();
+               set<string> processedLabels;
+               set<string> userLabels = labels;
+               
+               while((lookup[0] != NULL) && ((allLines == 1) || (userLabels.size() != 0))) {
+                       
+                       if (m->control_pressed) { for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }  out.close(); return 0; }
+                       
+                       if(allLines == 1 || labels.count(lookup[0]->getLabel()) == 1){                  
+                               
+                               m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
+                               processedLabels.insert(lookup[0]->getLabel());
+                               userLabels.erase(lookup[0]->getLabel());
+                               
+                               if (!m->printedHeaders) { lookup[0]->printHeaders(out); }
+                               processLookup(lookup, out);
+                       }
+                       
+                       if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
+                               string saveLabel = lookup[0]->getLabel();
+                               
+                               for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  } 
+                               lookup = input.getSharedRAbundVectors(lastLabel);
+                               
+                               m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
+                               processedLabels.insert(lookup[0]->getLabel());
+                               userLabels.erase(lookup[0]->getLabel());
+                               
+                               if (!m->printedHeaders) { lookup[0]->printHeaders(out); }
+                               processLookup(lookup, out);                     
+                               
+                               //restore real lastlabel to save below
+                               lookup[0]->setLabel(saveLabel);
+                       }               
+                       
+                       lastLabel = lookup[0]->getLabel();                      
+                       
+                       for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  } 
+                       lookup = input.getSharedRAbundVectors();
+               }
+               
+               if (m->control_pressed) {  out.close(); return 0; }     
+               
+               //output error messages about any remaining user labels
+               set<string>::iterator it;
+               bool needToRun = false;
+               for (it = userLabels.begin(); it != userLabels.end(); it++) {  
+                       m->mothurOut("Your file does not include the label " + *it); 
+                       if (processedLabels.count(lastLabel) != 1) {
+                               m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
+                               needToRun = true;
+                       }else {
+                               m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
+                       }
+               }
+               
+               //run last label if you need to
+               if (needToRun == true)  {
+                       for (int i = 0; i < lookup.size(); i++) {  if (lookup[i] != NULL) {     delete lookup[i];       }  }
+                       lookup = input.getSharedRAbundVectors(lastLabel);
+                       
+                       m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
+                       
+                       if (!m->printedHeaders) { lookup[0]->printHeaders(out); }
+                       processLookup(lookup, out);     
+                       
+                       for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  } 
+               }
+               
+               return 0;
+       }
        catch(exception& e) {
                m->errorOut(e, "RemoveRareCommand", "processSabund");
                exit(1);
        }
 }
 //**********************************************************************************************************************
+int RemoveRareCommand::processLookup(vector<SharedRAbundVector*>& lookup, ofstream& out){
+       try {
+               
+               vector<SharedRAbundVector> newRabunds;  newRabunds.resize(lookup.size());
+               for (int i = 0; i < lookup.size(); i++) {  
+                       newRabunds[i].setGroup(lookup[i]->getGroup());
+                       newRabunds[i].setLabel(lookup[i]->getLabel());
+               }
+               
+               if (byGroup) {
+                       
+                       //for each otu
+                       for (int i = 0; i < lookup[0]->getNumBins(); i++) {
+                               bool allZero = true;
+                               
+                               if (m->control_pressed) { return 0; }
+                               
+                               //for each group
+                               for (int j = 0; j < lookup.size(); j++) {
+                                       
+                                       //are you rare?
+                                       if (lookup[j]->getAbundance(i) > nseqs) {
+                                               newRabunds[j].push_back(lookup[j]->getAbundance(i), newRabunds[j].getGroup());
+                                               allZero = false;
+                                       }else {
+                                               newRabunds[j].push_back(0, newRabunds[j].getGroup());
+                                       }
+                               }
+                               
+                               //eliminates zero otus
+                               if (allZero) { for (int j = 0; j < newRabunds.size(); j++) {  newRabunds[j].pop_back(); } }
+                       }
+               }else {
+                       //for each otu
+                       for (int i = 0; i < lookup[0]->getNumBins(); i++) {
+                               
+                               if (m->control_pressed) { return 0; }
+                               
+                               int totalAbund = 0;
+                               //get total otu abundance
+                               for (int j = 0; j < lookup.size(); j++) {
+                                       newRabunds[j].push_back(lookup[j]->getAbundance(i), newRabunds[j].getGroup());
+                                       totalAbund += lookup[j]->getAbundance(i);
+                               }
+                               
+                               //eliminates otus below rare cutoff
+                               if (totalAbund <= nseqs) { for (int j = 0; j < newRabunds.size(); j++) {  newRabunds[j].pop_back(); } }
+                       }
+               }
+               
+               //do we have any otus above the rare cutoff
+               if (newRabunds[0].getNumBins() != 0) { 
+                       for (int j = 0; j < newRabunds.size(); j++) { 
+                               out << newRabunds[j].getLabel() << '\t' << newRabunds[j].getGroup() << '\t';
+                               newRabunds[j].print(out); 
+                       }
+               }
+               
+               return 0;
+       }
+       catch(exception& e) {
+               m->errorOut(e, "RemoveRareCommand", "processLookup");
+               exit(1);
+       }
+}
+//**********************************************************************************************************************