]> git.donarmstrong.com Git - mothur.git/blobdiff - hclustercommand.cpp
changes while testing
[mothur.git] / hclustercommand.cpp
index 94082dfc74dfe0969e769572fc6df219ecaa69d2..b991ccd6289b9afd44628c12bf451eb86529fa44 100644 (file)
 #include "hclustercommand.h"
 
 //**********************************************************************************************************************
-vector<string> HClusterCommand::getValidParameters(){  
+vector<string> HClusterCommand::setParameters(){       
        try {
-               string Array[] =  {"cutoff","hard","precision","method","phylip","column","name","sorted","showabund","timing","outputdir","inputdir"};
-               vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
+               CommandParameter pphylip("phylip", "InputTypes", "", "", "PhylipColumn", "PhylipColumn", "none","list-rabund-sabund",false,false,true); parameters.push_back(pphylip);
+               CommandParameter pname("name", "InputTypes", "", "", "none", "none", "ColumnName","",false,false,true); parameters.push_back(pname);
+               CommandParameter pcolumn("column", "InputTypes", "", "", "PhylipColumn", "PhylipColumn", "ColumnName","list-rabund-sabund",false,false,true); parameters.push_back(pcolumn);
+               CommandParameter pcutoff("cutoff", "Number", "", "10", "", "", "","",false,false,true); parameters.push_back(pcutoff);
+               CommandParameter pprecision("precision", "Number", "", "100", "", "", "","",false,false); parameters.push_back(pprecision);
+               CommandParameter pmethod("method", "Multiple", "furthest-nearest-average-weighted", "average", "", "", "","",false,false); parameters.push_back(pmethod);
+               CommandParameter phard("hard", "Boolean", "", "T", "", "", "","",false,false); parameters.push_back(phard);
+               CommandParameter psorted("sorted", "Boolean", "", "F", "", "", "","",false,false); parameters.push_back(psorted);
+               CommandParameter pshowabund("showabund", "Boolean", "", "T", "", "", "","",false,false); parameters.push_back(pshowabund);
+               CommandParameter ptiming("timing", "Boolean", "", "F", "", "", "","",false,false); parameters.push_back(ptiming);               
+               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, "HClusterCommand", "getValidParameters");
+               m->errorOut(e, "HClusterCommand", "setParameters");
                exit(1);
        }
 }
 //**********************************************************************************************************************
-HClusterCommand::HClusterCommand(){    
+string HClusterCommand::getHelpString(){       
        try {
-               abort = true; calledHelp = true; 
-               vector<string> tempOutNames;
-               outputTypes["list"] = tempOutNames;
-               outputTypes["rabund"] = tempOutNames;
-               outputTypes["sabund"] = tempOutNames;
+               string helpString = "";
+               helpString += "The hcluster command parameter options are cutoff, precision, method, phylip, column, name, showabund, timing and sorted. Phylip or column and name are required, unless you have valid current files.\n";
+               helpString += "The phylip and column parameter allow you to enter your distance file, and sorted indicates whether your column distance file is already sorted. \n";
+               helpString += "The name parameter allows you to enter your name file and is required if your distance file is in column format. \n";
+               helpString += "The hcluster command should be in the following format: \n";
+               helpString += "hcluster(column=youDistanceFile, name=yourNameFile, method=yourMethod, cutoff=yourCutoff, precision=yourPrecision) \n";
+               helpString += "The acceptable hcluster methods are furthest, nearest, weighted and average.\n"; 
+               return helpString;
        }
        catch(exception& e) {
-               m->errorOut(e, "HClusterCommand", "HClusterCommand");
+               m->errorOut(e, "HClusterCommand", "getHelpString");
                exit(1);
        }
 }
 //**********************************************************************************************************************
-vector<string> HClusterCommand::getRequiredParameters(){       
-       try {
-               string Array[] =  {"phylip","column","or"};
-               vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
-               return myArray;
-       }
-       catch(exception& e) {
-               m->errorOut(e, "HClusterCommand", "getRequiredParameters");
-               exit(1);
-       }
+string HClusterCommand::getOutputPattern(string type) {
+    try {
+        string pattern = "";
+        
+        if (type == "list") {  pattern = "[filename],[clustertag],list"; } 
+        else if (type == "rabund") {  pattern = "[filename],[clustertag],rabund"; } 
+        else if (type == "sabund") {  pattern = "[filename],[clustertag],sabund"; }
+        else { m->mothurOut("[ERROR]: No definition for type " + type + " output pattern.\n"); m->control_pressed = true;  }
+        
+        return pattern;
+    }
+    catch(exception& e) {
+        m->errorOut(e, "HClusterCommand", "getOutputPattern");
+        exit(1);
+    }
 }
+
 //**********************************************************************************************************************
-vector<string> HClusterCommand::getRequiredFiles(){    
+HClusterCommand::HClusterCommand(){    
        try {
-               vector<string> myArray;
-               return myArray;
+               abort = true; calledHelp = true; 
+               setParameters();
+               vector<string> tempOutNames;
+               outputTypes["list"] = tempOutNames;
+               outputTypes["rabund"] = tempOutNames;
+               outputTypes["sabund"] = tempOutNames;
        }
        catch(exception& e) {
-               m->errorOut(e, "HClusterCommand", "getRequiredFiles");
+               m->errorOut(e, "HClusterCommand", "HClusterCommand");
                exit(1);
        }
 }
@@ -62,16 +88,14 @@ vector<string> HClusterCommand::getRequiredFiles(){
 //This function checks to make sure the cluster command has no errors and then clusters based on the method chosen.
 HClusterCommand::HClusterCommand(string option)  {
        try{
-               globaldata = GlobalData::getInstance();
                abort = false; calledHelp = false;   
                
                //allow user to run help
                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[] =  {"cutoff","hard","precision","method","phylip","column","name","sorted","showabund","timing","outputdir","inputdir"};
-                       vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
+                       vector<string> myArray = setParameters();
                        
                        OptionParser parser(option);
                        map<string,string> parameters = parser.getParameters();
@@ -92,8 +116,6 @@ HClusterCommand::HClusterCommand(string option)  {
                        outputTypes["rabund"] = tempOutNames;
                        outputTypes["sabund"] = tempOutNames;
                
-                       globaldata->newRead();
-                       
                        //if the user changes the input directory command factory will send this info to us in the output parameter 
                        string inputDir = validParameter.validFile(parameters, "inputdir", false);              
                        if (inputDir == "not found"){   inputDir = "";          }
@@ -131,22 +153,43 @@ HClusterCommand::HClusterCommand(string option)  {
                        phylipfile = validParameter.validFile(parameters, "phylip", true);
                        if (phylipfile == "not open") { abort = true; }
                        else if (phylipfile == "not found") { phylipfile = ""; }        
-                       else {  distfile = phylipfile;  format = "phylip";      }
+                       else {  distfile = phylipfile;  format = "phylip";      m->setPhylipFile(phylipfile); }
                        
                        columnfile = validParameter.validFile(parameters, "column", true);
                        if (columnfile == "not open") { abort = true; } 
                        else if (columnfile == "not found") { columnfile = ""; }
-                       else {  distfile = columnfile; format = "column";       }
+                       else {  distfile = columnfile; format = "column";       m->setColumnFile(columnfile); }
                        
                        namefile = validParameter.validFile(parameters, "name", true);
                        if (namefile == "not open") { abort = true; }   
                        else if (namefile == "not found") { namefile = ""; }
+                       else { m->setNameFile(namefile); }
                        
-                       if ((phylipfile == "") && (columnfile == "")) { m->mothurOut("When executing a hcluster command you must enter a phylip or a column."); m->mothurOutEndLine(); abort = true; }
+                       if ((phylipfile == "") && (columnfile == "")) { 
+                               //is there are current file available for either of these?
+                               //give priority to column, then phylip
+                               columnfile = m->getColumnFile(); 
+                               if (columnfile != "") {  m->mothurOut("Using " + columnfile + " as input file for the column parameter."); m->mothurOutEndLine(); }
+                               else { 
+                                       phylipfile = m->getPhylipFile(); 
+                                       if (phylipfile != "") {  m->mothurOut("Using " + phylipfile + " as input file for the phylip parameter."); m->mothurOutEndLine(); }
+                                       else { 
+                                               m->mothurOut("No valid current files. You must provide a phylip or column file before you can use the hcluster command."); m->mothurOutEndLine(); 
+                                               abort = true;
+                                       }
+                               }
+                       }
                        else if ((phylipfile != "") && (columnfile != "")) { m->mothurOut("When executing a hcluster command you must enter ONLY ONE of the following: phylip or column."); m->mothurOutEndLine(); abort = true; }
                
                        if (columnfile != "") {
-                               if (namefile == "") {  cout << "You need to provide a namefile if you are going to use the column format." << endl; abort = true; }
+                               if (namefile == "") { 
+                                       namefile = m->getNameFile(); 
+                                       if (namefile != "") {  m->mothurOut("Using " + namefile + " as input file for the name parameter."); m->mothurOutEndLine(); }
+                                       else { 
+                                               m->mothurOut("You need to provide a namefile if you are going to use the column format."); m->mothurOutEndLine(); 
+                                               abort = true; 
+                                       }       
+                               }
                        }
                        
                        //check for optional parameter and set defaults
@@ -157,18 +200,18 @@ HClusterCommand::HClusterCommand(string option)  {
                        if (temp == "not found") { temp = "100"; }
                        //saves precision legnth for formatting below
                        length = temp.length();
-                       convert(temp, precision); 
+                       m->mothurConvert(temp, precision); 
                        
-                       temp = validParameter.validFile(parameters, "hard", false);                     if (temp == "not found") { temp = "F"; }
+                       temp = validParameter.validFile(parameters, "hard", false);                     if (temp == "not found") { temp = "T"; }
                        hard = m->isTrue(temp);
                        
                        temp = validParameter.validFile(parameters, "cutoff", false);
                        if (temp == "not found") { temp = "10"; }
-                       convert(temp, cutoff); 
+                       m->mothurConvert(temp, cutoff); 
                        cutoff += (5 / (precision * 10.0)); 
                        
                        method = validParameter.validFile(parameters, "method", false);
-                       if (method == "not found") { method = "furthest"; }
+                       if (method == "not found") { method = "average"; }
                        
                        if ((method == "furthest") || (method == "nearest") || (method == "average") || (method == "weighted")) { }
                        else { m->mothurOut("Not a valid clustering method.  Valid clustering algorithms are furthest, nearest, average or weighted."); m->mothurOutEndLine(); abort = true; }
@@ -193,14 +236,22 @@ HClusterCommand::HClusterCommand(string option)  {
                                else if (method == "nearest")   { tag = "nn";  }
                                else if (method == "weighted")  { tag = "wn";  }
                                else                                                    { tag = "an";  }
+                
+                map<string, string> variables;
+                variables["[filename]"] = fileroot;
+                variables["[clustertag]"] = tag;
                        
-                               m->openOutputFile(fileroot+ tag + ".sabund",    sabundFile);
-                               m->openOutputFile(fileroot+ tag + ".rabund",    rabundFile);
-                               m->openOutputFile(fileroot+ tag + ".list",              listFile);
-                               
-                               outputNames.push_back(fileroot+ tag + ".sabund"); outputTypes["sabund"].push_back(fileroot+ tag + ".sabund");
-                               outputNames.push_back(fileroot+ tag + ".rabund"); outputTypes["rabund"].push_back(fileroot+ tag + ".rabund");
-                               outputNames.push_back(fileroot+ tag + ".list"); outputTypes["list"].push_back(fileroot+ tag + ".list");
+                               string sabundFileName = getOutputFileName("sabund",variables);
+                string rabundFileName = getOutputFileName("rabund",variables);
+                string listFileName = getOutputFileName("list", variables);
+                
+                m->openOutputFile(sabundFileName,      sabundFile);
+                m->openOutputFile(rabundFileName,      rabundFile);
+                m->openOutputFile(listFileName,        listFile);
+                
+                outputNames.push_back(sabundFileName); outputTypes["sabund"].push_back(sabundFileName);
+                outputNames.push_back(rabundFileName); outputTypes["rabund"].push_back(rabundFileName);
+                outputNames.push_back(listFileName); outputTypes["list"].push_back(listFileName);
                        }
                }
        }
@@ -212,52 +263,30 @@ HClusterCommand::HClusterCommand(string option)  {
 
 //**********************************************************************************************************************
 
-void HClusterCommand::help(){
-       try {
-               m->mothurOut("The hcluster command parameter options are cutoff, precision, method, phylip, column, name, showabund, timing and sorted. Phylip or column and name are required.\n");
-               m->mothurOut("The phylip and column parameter allow you to enter your distance file, and sorted indicates whether your column distance file is already sorted. \n");
-               m->mothurOut("The name parameter allows you to enter your name file and is required if your distance file is in column format. \n");
-               m->mothurOut("The hcluster command should be in the following format: \n");
-               m->mothurOut("hcluster(column=youDistanceFile, name=yourNameFile, method=yourMethod, cutoff=yourCutoff, precision=yourPrecision) \n");
-               m->mothurOut("The acceptable hcluster methods are furthest, nearest, weighted and average.\n\n");       
-       }
-       catch(exception& e) {
-               m->errorOut(e, "HClusterCommand", "help");
-               exit(1);
-       }
-}
-
-//**********************************************************************************************************************
-
-HClusterCommand::~HClusterCommand(){}
-
-//**********************************************************************************************************************
-
 int HClusterCommand::execute(){
        try {
        
                if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
                
+               NameAssignment* nameMap = NULL;
                if(namefile != ""){     
-                       globaldata->nameMap = new NameAssignment(namefile);
-                       globaldata->nameMap->readMap();
-               }else{
-                       globaldata->nameMap = NULL;
-               }
+                       nameMap = new NameAssignment(namefile);
+                       nameMap->readMap();
+               }               
                
                time_t estart = time(NULL);
                
                if (!sorted) {
                        read = new ReadCluster(distfile, cutoff, outputDir, true);      
                        read->setFormat(format);
-                       read->read(globaldata->nameMap);
+                       read->read(nameMap);
                        
                        if (m->control_pressed) {  
                                delete read; 
                                sabundFile.close();
                                rabundFile.close();
                                listFile.close();
-                               for (int i = 0; i < outputNames.size(); i++) {  remove(outputNames[i].c_str());  } outputTypes.clear();
+                               for (int i = 0; i < outputNames.size(); i++) {  m->mothurRemove(outputNames[i]);  } outputTypes.clear();
                                return 0;  
                        }
                        
@@ -266,14 +295,14 @@ int HClusterCommand::execute(){
                        list = read->getListVector();
                        delete read;
                }else {
-                       list = new ListVector(globaldata->nameMap->getListVector());
+                       list = new ListVector(nameMap->getListVector());
                }
                
                if (m->control_pressed) {  
                        sabundFile.close();
                        rabundFile.close();
                        listFile.close();
-                       for (int i = 0; i < outputNames.size(); i++) {  remove(outputNames[i].c_str());  } outputTypes.clear();
+                       for (int i = 0; i < outputNames.size(); i++) {  m->mothurRemove(outputNames[i]);  } outputTypes.clear();
                        return 0;  
                }
 
@@ -295,7 +324,7 @@ int HClusterCommand::execute(){
                print_start = true;
                start = time(NULL);
                                
-               cluster = new HCluster(rabund, list, method, distfile, globaldata->nameMap, cutoff);
+               cluster = new HCluster(rabund, list, method, distfile, nameMap, cutoff);
                vector<seqDist> seqs; seqs.resize(1); // to start loop
                
                if (m->control_pressed) {  
@@ -303,35 +332,41 @@ int HClusterCommand::execute(){
                                sabundFile.close();
                                rabundFile.close();
                                listFile.close();
-                               for (int i = 0; i < outputNames.size(); i++) {  remove(outputNames[i].c_str());  } outputTypes.clear();
+                               for (int i = 0; i < outputNames.size(); i++) {  m->mothurRemove(outputNames[i]);  } outputTypes.clear();
                                return 0;  
                }
 
+               float saveCutoff = cutoff;
                
                while (seqs.size() != 0){
                
                        seqs = cluster->getSeqs();
                        
+                       //to account for cutoff change in average neighbor
+                       if (seqs.size() != 0) {
+                               if (seqs[0].dist > cutoff) { break; }
+                       }
+                       
                        if (m->control_pressed) {  
                                delete cluster;
                                sabundFile.close();
                                rabundFile.close();
                                listFile.close();
-                               for (int i = 0; i < outputNames.size(); i++) {  remove(outputNames[i].c_str());  } outputTypes.clear();
+                               for (int i = 0; i < outputNames.size(); i++) {  m->mothurRemove(outputNames[i]);  } outputTypes.clear();
                                return 0;  
                        }
 
                        for (int i = 0; i < seqs.size(); i++) {  //-1 means skip me
                                
                                if (seqs[i].seq1 != seqs[i].seq2) {
-                                       cluster->update(seqs[i].seq1, seqs[i].seq2, seqs[i].dist);
+                                       cutoff = cluster->update(seqs[i].seq1, seqs[i].seq2, seqs[i].dist);
                                        
                                        if (m->control_pressed) {  
                                                delete cluster;
                                                sabundFile.close();
                                                rabundFile.close();
                                                listFile.close();
-                                               for (int i = 0; i < outputNames.size(); i++) {  remove(outputNames[i].c_str());  } outputTypes.clear();
+                                               for (int i = 0; i < outputNames.size(); i++) {  m->mothurRemove(outputNames[i]);  } outputTypes.clear();
                                                return 0;  
                                        }
 
@@ -364,7 +399,7 @@ int HClusterCommand::execute(){
                        sabundFile.close();
                        rabundFile.close();
                        listFile.close();
-                       for (int i = 0; i < outputNames.size(); i++) {  remove(outputNames[i].c_str());  } outputTypes.clear();
+                       for (int i = 0; i < outputNames.size(); i++) {  m->mothurRemove(outputNames[i]);  } outputTypes.clear();
                        return 0;  
                }
                                        
@@ -374,28 +409,44 @@ int HClusterCommand::execute(){
                else if(rndPreviousDist<cutoff){
                        printData(toString(rndPreviousDist, length-1));
                }
-               
-               //delete globaldata's copy of the sparsematrix and listvector to free up memory
-               delete globaldata->gListVector;  globaldata->gListVector = NULL;
-               
-               //saves .list file so you can do the collect, rarefaction and summary commands without doing a read.list
-               if (globaldata->getFormat() == "phylip") { globaldata->setPhylipFile(""); }
-               else if (globaldata->getFormat() == "column") { globaldata->setColumnFile(""); }
-               
-               globaldata->setListFile(fileroot+ tag + ".list");
-               globaldata->setNameFile("");
-               globaldata->setFormat("list");
-               
+                               
                sabundFile.close();
                rabundFile.close();
                listFile.close();
                delete cluster;
                
                if (m->control_pressed) {  
-                       for (int i = 0; i < outputNames.size(); i++) {  remove(outputNames[i].c_str());  } outputTypes.clear();
+                       for (int i = 0; i < outputNames.size(); i++) {  m->mothurRemove(outputNames[i]);  } outputTypes.clear();
                        return 0;  
                }
-
+               
+               
+               if (saveCutoff != cutoff) { 
+                       if (hard)       {  saveCutoff = m->ceilDist(saveCutoff, precision);     }
+                       else            {       saveCutoff = m->roundDist(saveCutoff, precision);  }
+                       
+                       m->mothurOut("changed cutoff to " + toString(cutoff)); m->mothurOutEndLine(); 
+               }
+               
+               //set list file as new current listfile
+               string current = "";
+               itTypes = outputTypes.find("list");
+               if (itTypes != outputTypes.end()) {
+                       if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setListFile(current); }
+               }
+               
+               //set rabund file as new current rabundfile
+               itTypes = outputTypes.find("rabund");
+               if (itTypes != outputTypes.end()) {
+                       if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setRabundFile(current); }
+               }
+               
+               //set sabund file as new current sabundfile
+               itTypes = outputTypes.find("sabund");
+               if (itTypes != outputTypes.end()) {
+                       if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setSabundFile(current); }
+               }
+               
                
                m->mothurOutEndLine();
                m->mothurOut("Output File Names: "); m->mothurOutEndLine();