]> git.donarmstrong.com Git - mothur.git/blobdiff - clustersplitcommand.cpp
added pipeline commands which involved change to command factory and command class...
[mothur.git] / clustersplitcommand.cpp
index 38e839496536b6ae525c94efd0c5b024c7eee384..3c957a6ce84b5e5a489faaccf18a1dc40f99a02b 100644 (file)
 #include "readmatrix.hpp"
 #include "inputdata.h"
 
+
+//**********************************************************************************************************************
+vector<string> ClusterSplitCommand::getValidParameters(){      
+       try {
+               string AlignArray[] =  {"fasta","phylip","column","name","cutoff","precision","method","splitmethod","taxonomy","taxlevel","large","showabund","timing","hard","processors","outputdir","inputdir"};
+               vector<string> myArray (AlignArray, AlignArray+(sizeof(AlignArray)/sizeof(string)));
+               return myArray;
+       }
+       catch(exception& e) {
+               m->errorOut(e, "ClusterSplitCommand", "getValidParameters");
+               exit(1);
+       }
+}
+//**********************************************************************************************************************
+ClusterSplitCommand::ClusterSplitCommand(){    
+       try {
+               //initialize outputTypes
+               vector<string> tempOutNames;
+               outputTypes["list"] = tempOutNames;
+               outputTypes["rabund"] = tempOutNames;
+               outputTypes["sabund"] = tempOutNames;
+       }
+       catch(exception& e) {
+               m->errorOut(e, "ClusterSplitCommand", "ClusterSplitCommand");
+               exit(1);
+       }
+}
+//**********************************************************************************************************************
+vector<string> ClusterSplitCommand::getRequiredParameters(){   
+       try {
+               string Array[] =  {"fasta","phylip","column","or"};
+               vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
+               return myArray;
+       }
+       catch(exception& e) {
+               m->errorOut(e, "ClusterSplitCommand", "getRequiredParameters");
+               exit(1);
+       }
+}
+//**********************************************************************************************************************
+vector<string> ClusterSplitCommand::getRequiredFiles(){        
+       try {
+               vector<string> myArray;
+               return myArray;
+       }
+       catch(exception& e) {
+               m->errorOut(e, "ClusterSplitCommand", "getRequiredFiles");
+               exit(1);
+       }
+}
 //**********************************************************************************************************************
 //This function checks to make sure the cluster command has no errors and then clusters based on the method chosen.
 ClusterSplitCommand::ClusterSplitCommand(string option)  {
@@ -34,7 +84,7 @@ ClusterSplitCommand::ClusterSplitCommand(string option)  {
                        OptionParser parser(option);
                        map<string,string> parameters = parser.getParameters();
                        
-                       ValidParameters validParameter;
+                       ValidParameters validParameter("cluster.split");
                
                        //check to make sure all parameters are valid for command
                        map<string,string>::iterator it;
@@ -44,6 +94,12 @@ ClusterSplitCommand::ClusterSplitCommand(string option)  {
                                }
                        }
                        
+                       //initialize outputTypes
+                       vector<string> tempOutNames;
+                       outputTypes["list"] = tempOutNames;
+                       outputTypes["rabund"] = tempOutNames;
+                       outputTypes["sabund"] = tempOutNames;
+                       
                        globaldata->newRead();
                        
                        //if the user changes the output directory command factory will send this info to us in the output parameter 
@@ -57,7 +113,7 @@ ClusterSplitCommand::ClusterSplitCommand(string option)  {
                                it = parameters.find("phylip");
                                //user has given a template file
                                if(it != parameters.end()){ 
-                                       path = hasPath(it->second);
+                                       path = m->hasPath(it->second);
                                        //if the user has not given a path then, add inputdir. else leave path alone.
                                        if (path == "") {       parameters["phylip"] = inputDir + it->second;           }
                                }
@@ -65,7 +121,7 @@ ClusterSplitCommand::ClusterSplitCommand(string option)  {
                                it = parameters.find("column");
                                //user has given a template file
                                if(it != parameters.end()){ 
-                                       path = hasPath(it->second);
+                                       path = m->hasPath(it->second);
                                        //if the user has not given a path then, add inputdir. else leave path alone.
                                        if (path == "") {       parameters["column"] = inputDir + it->second;           }
                                }
@@ -73,7 +129,7 @@ ClusterSplitCommand::ClusterSplitCommand(string option)  {
                                it = parameters.find("name");
                                //user has given a template file
                                if(it != parameters.end()){ 
-                                       path = hasPath(it->second);
+                                       path = m->hasPath(it->second);
                                        //if the user has not given a path then, add inputdir. else leave path alone.
                                        if (path == "") {       parameters["name"] = inputDir + it->second;             }
                                }
@@ -81,7 +137,7 @@ ClusterSplitCommand::ClusterSplitCommand(string option)  {
                                it = parameters.find("taxonomy");
                                //user has given a template file
                                if(it != parameters.end()){ 
-                                       path = hasPath(it->second);
+                                       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;         }
                                }
@@ -89,7 +145,7 @@ ClusterSplitCommand::ClusterSplitCommand(string option)  {
                                it = parameters.find("fasta");
                                //user has given a template file
                                if(it != parameters.end()){ 
-                                       path = hasPath(it->second);
+                                       path = m->hasPath(it->second);
                                        //if the user has not given a path then, add inputdir. else leave path alone.
                                        if (path == "") {       parameters["fasta"] = inputDir + it->second;            }
                                }
@@ -142,10 +198,10 @@ ClusterSplitCommand::ClusterSplitCommand(string option)  {
                        convert(temp, precision); 
                        
                        temp = validParameter.validFile(parameters, "hard", false);                     if (temp == "not found") { temp = "F"; }
-                       hard = isTrue(temp);
+                       hard = m->isTrue(temp);
                        
                        temp = validParameter.validFile(parameters, "large", false);                    if (temp == "not found") { temp = "F"; }
-                       large = isTrue(temp);
+                       large = m->isTrue(temp);
                        
                        temp = validParameter.validFile(parameters, "processors", false);       if (temp == "not found"){       temp = "1";                             }
                        convert(temp, processors); 
@@ -195,7 +251,7 @@ void ClusterSplitCommand::help(){
                m->mothurOut("The cluster.split command can split your files in 3 ways. Splitting by distance file, by classification, or by classification also using a fasta file. \n");
                m->mothurOut("For the distance file method, you need only provide your distance file and mothur will split the file into distinct groups. \n");
                m->mothurOut("For the classification method, you need to provide your distance file and taxonomy file, and set the splitmethod to classify.  \n");
-               m->mothurOut("You will also need to set the taxlevel you want to split by. mothur will split the sequence into distinct taxonomy groups, and split the distance file based on those groups. \n");
+               m->mothurOut("You will also need to set the taxlevel you want to split by. mothur will split the sequences into distinct taxonomy groups, and split the distance file based on those groups. \n");
                m->mothurOut("For the classification method using a fasta file, you need to provide your fasta file, names file and taxonomy file.  \n");
                m->mothurOut("You will also need to set the taxlevel you want to split by. mothur will split the sequence into distinct taxonomy groups, and create distance files for each grouping. \n");
                m->mothurOut("The phylip and column parameter allow you to enter your distance file. \n");
@@ -206,7 +262,7 @@ void ClusterSplitCommand::help(){
                m->mothurOut("The method allows you to specify what clustering algorythm you want to use, default=furthest, option furthest, nearest, or average. \n");
                m->mothurOut("The splitmethod parameter allows you to specify how you want to split your distance file before you cluster, default=distance, options distance, classify or fasta. \n");
                m->mothurOut("The taxonomy parameter allows you to enter the taxonomy file for your sequences, this is only valid if you are using splitmethod=classify. Be sure your taxonomy file does not include the probability scores. \n");
-               m->mothurOut("The taxlevel parameter allows you to specify the taxonomy level you want to use to split the distance file, default=1. \n");
+               m->mothurOut("The taxlevel parameter allows you to specify the taxonomy level you want to use to split the distance file, default=1, meaning use the first taxon in each list. \n");
                m->mothurOut("The large parameter allows you to indicate that your distance matrix is too large to fit in RAM.  The default value is false.\n");
                #ifdef USE_MPI
                m->mothurOut("When using MPI, the processors parameter is set to the number of MPI processes running. \n");
@@ -237,6 +293,7 @@ int ClusterSplitCommand::execute(){
                vector<string> listFileNames;
                set<string> labels;
                string singletonName = "";
+               double saveCutoff = cutoff;
 
                //****************** file prep work ******************************//
                #ifdef USE_MPI
@@ -270,7 +327,7 @@ int ClusterSplitCommand::execute(){
                        if (namefile == "") {  //you need to make a namefile for split matrix
                                ofstream out;
                                namefile = phylipfile + ".names";
-                               openOutputFile(namefile, out);
+                               m->openOutputFile(namefile, out);
                                for (int i = 0; i < listToMakeNameFile->getNumBins(); i++) {
                                        string bin = listToMakeNameFile->get(i);
                                        out << bin << '\t' << bin << endl;
@@ -291,7 +348,7 @@ int ClusterSplitCommand::execute(){
                SplitMatrix* split;
                if (splitmethod == "distance")                  {       split = new SplitMatrix(distfile, namefile, taxFile, cutoff, splitmethod, large);                                                       }
                else if (splitmethod == "classify")             {       split = new SplitMatrix(distfile, namefile, taxFile, taxLevelCutoff, splitmethod, large);                                       }
-               else if (splitmethod == "fasta")                {       split = new SplitMatrix(fastafile, namefile, taxFile, taxLevelCutoff, splitmethod, processors, outputDir);      }
+               else if (splitmethod == "fasta")                {       split = new SplitMatrix(fastafile, namefile, taxFile, taxLevelCutoff, cutoff, splitmethod, processors, outputDir);      }
                else { m->mothurOut("Not a valid splitting method.  Valid splitting algorithms are distance, classify or fasta."); m->mothurOutEndLine(); return 0;             }
                
                split->split();
@@ -361,6 +418,10 @@ int ClusterSplitCommand::execute(){
                        for(int i = 1; i < processors; i++) { 
                                int num = dividedNames[i].size();
                                
+                               double tempCutoff;
+                               MPI_Recv(&tempCutoff, 1, MPI_DOUBLE, i, tag, MPI_COMM_WORLD, &status);
+                               if (tempCutoff < cutoff) { cutoff = tempCutoff; }
+                               
                                //send list filenames to root process
                                for (int j = 0; j < num; j++) {  
                                        int lengthList = 0;
@@ -429,6 +490,9 @@ int ClusterSplitCommand::execute(){
                        //process them
                        listFileNames = cluster(myNames, labels);
                        
+                       //send cutoff
+                       MPI_Send(&cutoff, 1, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD);
+                       
                        //send list filenames to root process
                        for (int j = 0; j < num; j++) {  
                                char tempListFileName[1024];
@@ -489,13 +553,13 @@ int ClusterSplitCommand::execute(){
                                        for(int i=0;i<processors;i++){
                                                string filename = toString(processIDS[i]) + ".temp";
                                                ifstream in;
-                                               openInputFile(filename, in);
+                                               m->openInputFile(filename, in);
                                                
-                                               in >> tag; gobble(in);
+                                               in >> tag; m->gobble(in);
                                                
                                                while(!in.eof()) {
                                                        string tempName;
-                                                       in >> tempName; gobble(in);
+                                                       in >> tempName; m->gobble(in);
                                                        listFileNames.push_back(tempName);
                                                }
                                                in.close();
@@ -504,11 +568,15 @@ int ClusterSplitCommand::execute(){
                                                //get labels
                                                filename = toString(processIDS[i]) + ".temp.labels";
                                                ifstream in2;
-                                               openInputFile(filename, in2);
+                                               m->openInputFile(filename, in2);
+                                               
+                                               float tempCutoff;
+                                               in2 >> tempCutoff; m->gobble(in2);
+                                               if (tempCutoff < cutoff) { cutoff = tempCutoff; }
                                                
                                                while(!in2.eof()) {
                                                        string tempName;
-                                                       in2 >> tempName; gobble(in2);
+                                                       in2 >> tempName; m->gobble(in2);
                                                        if (labels.count(tempName) == 0) { labels.insert(tempName); }
                                                }
                                                in2.close();
@@ -521,6 +589,8 @@ int ClusterSplitCommand::execute(){
        #endif  
                if (m->control_pressed) { for (int i = 0; i < listFileNames.size(); i++) { remove(listFileNames[i].c_str()); } return 0; }
                
+               if (saveCutoff != cutoff) { m->mothurOut("Cutoff was " + toString(saveCutoff) + " changed cutoff to " + toString(cutoff)); m->mothurOutEndLine();  }
+               
                m->mothurOut("It took " + toString(time(NULL) - estart) + " seconds to cluster"); m->mothurOutEndLine();
                
                //****************** merge list file and create rabund and sabund files ******************************//
@@ -572,12 +642,12 @@ map<float, int> ClusterSplitCommand::completeListFile(vector<string> listNames,
                //read in singletons
                if (singleton != "none") {
                        ifstream in;
-                       openInputFile(singleton, in);
+                       m->openInputFile(singleton, in);
                                
                        string firstCol, secondCol;
                        listSingle = new ListVector();
                        while (!in.eof()) {
-                               in >> firstCol >> secondCol; gobble(in);
+                               in >> firstCol >> secondCol; m->gobble(in);
                                listSingle->push_back(secondCol);
                        }
                        in.close();
@@ -592,9 +662,11 @@ map<float, int> ClusterSplitCommand::completeListFile(vector<string> listNames,
 
                        if ((*it != "unique") && (convertTestFloat(*it, temp) == true)) {       convert(*it, temp);     }
                        else if (*it == "unique")                                                                               {       temp = -1.0;            }
-                                               
-                       orderFloat.push_back(temp);
-                       labelBin[temp] = numSingleBins; //initialize numbins 
+                       
+                       if (temp <= cutoff) {
+                               orderFloat.push_back(temp);
+                               labelBin[temp] = numSingleBins; //initialize numbins 
+                       }
                }
        
                //sort order
@@ -616,7 +688,7 @@ map<float, int> ClusterSplitCommand::completeListFile(vector<string> listNames,
                        
                        string filledInList = listNames[k] + "filledInTemp";
                        ofstream outFilled;
-                       openOutputFile(filledInList, outFilled);
+                       m->openOutputFile(filledInList, outFilled);
        
                        //for each label needed
                        for(int l = 0; l < orderFloat.size(); l++){
@@ -673,16 +745,16 @@ map<float, int> ClusterSplitCommand::completeListFile(vector<string> listNames,
 //**********************************************************************************************************************
 int ClusterSplitCommand::mergeLists(vector<string> listNames, map<float, int> userLabels, ListVector* listSingle){
        try {
-               if (outputDir == "") { outputDir += hasPath(distfile); }
-               fileroot = outputDir + getRootName(getSimpleName(distfile));
+               if (outputDir == "") { outputDir += m->hasPath(distfile); }
+               fileroot = outputDir + m->getRootName(m->getSimpleName(distfile));
                
-               openOutputFile(fileroot+ tag + ".sabund",       outSabund);
-               openOutputFile(fileroot+ tag + ".rabund",       outRabund);
-               openOutputFile(fileroot+ tag + ".list",         outList);
+               m->openOutputFile(fileroot+ tag + ".sabund",    outSabund);
+               m->openOutputFile(fileroot+ tag + ".rabund",    outRabund);
+               m->openOutputFile(fileroot+ tag + ".list",              outList);
                                
-               outputNames.push_back(fileroot+ tag + ".sabund");
-               outputNames.push_back(fileroot+ tag + ".rabund");
-               outputNames.push_back(fileroot+ tag + ".list");
+               outputNames.push_back(fileroot+ tag + ".sabund");  outputTypes["list"].push_back(fileroot+ tag + ".list");
+               outputNames.push_back(fileroot+ tag + ".rabund");  outputTypes["rabund"].push_back(fileroot+ tag + ".rabund");
+               outputNames.push_back(fileroot+ tag + ".list");    outputTypes["sabund"].push_back(fileroot+ tag + ".sabund");
                
                map<float, int>::iterator itLabel;
 
@@ -702,7 +774,7 @@ int ClusterSplitCommand::mergeLists(vector<string> listNames, map<float, int> us
                        if (listSingle != NULL) {
                                for (int j = 0; j < listSingle->getNumBins(); j++) {
                                        outList << listSingle->get(j) << '\t';
-                                       rabund->push_back(getNumNames(listSingle->get(j)));
+                                       rabund->push_back(m->getNumNames(listSingle->get(j)));
                                }
                        }
                        
@@ -719,7 +791,7 @@ int ClusterSplitCommand::mergeLists(vector<string> listNames, map<float, int> us
                                else {          
                                        for (int j = 0; j < list->getNumBins(); j++) {
                                                outList << list->get(j) << '\t';
-                                               rabund->push_back(getNumNames(list->get(j)));
+                                               rabund->push_back(m->getNumNames(list->get(j)));
                                        }
                                        delete list;
                                }
@@ -759,7 +831,7 @@ void ClusterSplitCommand::printData(ListVector* oldList){
                RAbundVector oldRAbund = oldList->getRAbundVector();
                
                oldRAbund.setLabel(label);
-               if (isTrue(showabund)) {
+               if (m->isTrue(showabund)) {
                        oldRAbund.getSAbundVector().print(cout);
                }
                oldRAbund.print(outRabund);
@@ -795,7 +867,7 @@ int ClusterSplitCommand::createProcesses(vector < vector < map<string, string> >
                                //write out names to file
                                string filename = toString(getpid()) + ".temp";
                                ofstream out;
-                               openOutputFile(filename, out);
+                               m->openOutputFile(filename, out);
                                out << tag << endl;
                                for (int j = 0; j < listFileNames.size(); j++) { out << listFileNames[j] << endl;  }
                                out.close();
@@ -803,8 +875,9 @@ int ClusterSplitCommand::createProcesses(vector < vector < map<string, string> >
                                //print out labels
                                ofstream outLabels;
                                filename = toString(getpid()) + ".temp.labels";
-                               openOutputFile(filename, outLabels);
-               
+                               m->openOutputFile(filename, outLabels);
+                               
+                               outLabels << cutoff << endl;
                                for (set<string>::iterator it = labels.begin(); it != labels.end(); it++) {
                                        outLabels << (*it) << endl;
                                }
@@ -841,8 +914,11 @@ vector<string> ClusterSplitCommand::cluster(vector< map<string, string> > distNa
                
                vector<string> listFileNames;
                
+               double smallestCutoff = cutoff;
+               
                //cluster each distance file
                for (int i = 0; i < distNames.size(); i++) {
+                       if (m->control_pressed) { return listFileNames; }
                        
                        string thisNamefile = distNames[i].begin()->second;
                        string thisDistFile = distNames[i].begin()->first;
@@ -897,11 +973,11 @@ vector<string> ClusterSplitCommand::cluster(vector< map<string, string> > distNa
                        else if(method == "average"){   cluster = new AverageLinkage(rabund, list, matrix, cutoff, method);     }
                        tag = cluster->getTag();
                
-                       if (outputDir == "") { outputDir += hasPath(thisDistFile); }
-                       fileroot = outputDir + getRootName(getSimpleName(thisDistFile));
+                       if (outputDir == "") { outputDir += m->hasPath(thisDistFile); }
+                       fileroot = outputDir + m->getRootName(m->getSimpleName(thisDistFile));
                        
                        ofstream listFile;
-                       openOutputFile(fileroot+ tag + ".list", listFile);
+                       m->openOutputFile(fileroot+ tag + ".list",      listFile);
                
                        listFileNames.push_back(fileroot+ tag + ".list");
                
@@ -925,14 +1001,14 @@ vector<string> ClusterSplitCommand::cluster(vector< map<string, string> > distNa
                                        listFileNames.clear(); return listFileNames;
                                }
                
-                               cluster->update(cutoff);
+                               cluster->update(saveCutoff);
        
                                float dist = matrix->getSmallDist();
                                float rndDist;
                                if (hard) {
-                                       rndDist = ceilDist(dist, precision); 
+                                       rndDist = m->ceilDist(dist, precision); 
                                }else{
-                                       rndDist = roundDist(dist, precision); 
+                                       rndDist = m->roundDist(dist, precision); 
                                }
 
                                if(previousDist <= 0.0000 && dist != previousDist){
@@ -973,8 +1049,18 @@ vector<string> ClusterSplitCommand::cluster(vector< map<string, string> > distNa
                        
                        remove(thisDistFile.c_str());
                        remove(thisNamefile.c_str());
+                       
+                       if (saveCutoff != cutoff) { 
+                               if (hard)       {  saveCutoff = m->ceilDist(saveCutoff, precision);     }
+                               else            {       saveCutoff = m->roundDist(saveCutoff, precision);  }
+                       
+                               m->mothurOut("Cutoff was " + toString(cutoff) + " changed cutoff to " + toString(saveCutoff)); m->mothurOutEndLine();  
+                       }
+                       
+                       if (saveCutoff < smallestCutoff) { smallestCutoff = saveCutoff;  }
                }
                
+               cutoff = smallestCutoff;
                                        
                return listFileNames;