clustercommand.cpp

   1 /*
   2  *  clustercommand.cpp
   3  *  Dotur
   4  *
   5  *  Created by Sarah Westcott on 1/2/09.
   6  *  Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
   7  *
   8  */
   9
  10 #include "clustercommand.h"
  11
  12 //**********************************************************************************************************************
  13 vector<string> ClusterCommand::getValidParameters(){
  14         try {
  15                 string AlignArray[] =  {"cutoff","precision","method","showabund","timing","hard","outputdir","inputdir"};
  16                 vector<string> myArray (AlignArray, AlignArray+(sizeof(AlignArray)/sizeof(string)));
  17                 return myArray;
  18         }
  19         catch(exception& e) {
  20                 m->errorOut(e, "ClusterCommand", "getValidParameters");
  21                 exit(1);
  22         }
  23 }
  24 //**********************************************************************************************************************
  25 ClusterCommand::ClusterCommand(){
  26         try {
  27                 abort = true; calledHelp = true;
  28                 vector<string> tempOutNames;
  29                 outputTypes["list"] = tempOutNames;
  30                 outputTypes["rabund"] = tempOutNames;
  31                 outputTypes["sabund"] = tempOutNames;
  32         }
  33         catch(exception& e) {
  34                 m->errorOut(e, "ClusterCommand", "ClusterCommand");
  35                 exit(1);
  36         }
  37 }
  38 //**********************************************************************************************************************
  39 vector<string> ClusterCommand::getRequiredParameters(){
  40         try {
  41                 vector<string> myArray;
  42                 return myArray;
  43         }
  44         catch(exception& e) {
  45                 m->errorOut(e, "ClusterCommand", "getRequiredParameters");
  46                 exit(1);
  47         }
  48 }
  49 //**********************************************************************************************************************
  50 vector<string> ClusterCommand::getRequiredFiles(){
  51         try {
  52                 string Array[] =  {"phylip","column","or"};
  53                 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
  54                 return myArray;
  55         }
  56         catch(exception& e) {
  57                 m->errorOut(e, "ClusterCommand", "getRequiredFiles");
  58                 exit(1);
  59         }
  60 }
  61 //**********************************************************************************************************************
  62 //This function checks to make sure the cluster command has no errors and then clusters based on the method chosen.
  63 ClusterCommand::ClusterCommand(string option)  {
  64         try{
  65                 globaldata = GlobalData::getInstance();
  66
  67                 abort = false; calledHelp = false;
  68
  69                 //allow user to run help
  70                 if(option == "help") { help(); abort = true; calledHelp = true; }
  71
  72                 else {
  73                         //valid paramters for this command
  74                         string Array[] =  {"cutoff","precision","method","showabund","timing","hard","outputdir","inputdir"};
  75                         vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
  76
  77                         OptionParser parser(option);
  78                         map<string,string> parameters = parser.getParameters();
  79
  80                         ValidParameters validParameter;
  81
  82                         //check to make sure all parameters are valid for command
  83                         for (map<string,string>::iterator it = parameters.begin(); it != parameters.end(); it++) {
  84                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {
  85                                         abort = true;
  86                                 }
  87                         }
  88
  89                         //initialize outputTypes
  90                         vector<string> tempOutNames;
  91                         outputTypes["list"] = tempOutNames;
  92                         outputTypes["rabund"] = tempOutNames;
  93                         outputTypes["sabund"] = tempOutNames;
  94
  95                         //if the user changes the output directory command factory will send this info to us in the output parameter
  96                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){  outputDir = "";         }
  97
  98                         //error checking to make sure they read a distance file
  99                         if ((globaldata->gSparseMatrix == NULL) || (globaldata->gListVector == NULL)) {
 100                                 m->mothurOut("Before you use the cluster command, you first need to read in a distance matrix."); m->mothurOutEndLine();
 101                                 abort = true;
 102                         }
 103
 104                         //check for optional parameter and set defaults
 105                         // ...at some point should added some additional type checking...
 106                         //get user cutoff and precision or use defaults
 107                         string temp;
 108                         temp = validParameter.validFile(parameters, "precision", false);
 109                         if (temp == "not found") { temp = "100"; }
 110                         //saves precision legnth for formatting below
 111                         length = temp.length();
 112                         convert(temp, precision);
 113
 114                         temp = validParameter.validFile(parameters, "hard", false);                     if (temp == "not found") { temp = "F"; }
 115                         hard = m->isTrue(temp);
 116
 117                         temp = validParameter.validFile(parameters, "cutoff", false);
 118                         if (temp == "not found") { temp = "10"; }
 119                         convert(temp, cutoff);
 120                         cutoff += (5 / (precision * 10.0));
 121
 122                         method = validParameter.validFile(parameters, "method", false);
 123                         if (method == "not found") { method = "furthest"; }
 124
 125                         if ((method == "furthest") || (method == "nearest") || (method == "average") || (method == "weighted")) { }
 126                         else { m->mothurOut("Not a valid clustering method.  Valid clustering algorithms are furthest, nearest, average, and weighted."); m->mothurOutEndLine(); abort = true; }
 127
 128                         showabund = validParameter.validFile(parameters, "showabund", false);
 129                         if (showabund == "not found") { showabund = "T"; }
 130
 131                         timing = validParameter.validFile(parameters, "timing", false);
 132                         if (timing == "not found") { timing = "F"; }
 133
 134                         if (abort == false) {
 135
 136
 137                                                         //get matrix, list and rabund for execute
 138                                 if(globaldata->gSparseMatrix != NULL)   {       matrix = globaldata->gSparseMatrix;             }
 139
 140                                 if(globaldata->gListVector != NULL){
 141                                         list = globaldata->gListVector;
 142                                         rabund = new RAbundVector(list->getRAbundVector());
 143                                 }
 144
 145                                 //create cluster
 146                                 if (method == "furthest")       {       cluster = new CompleteLinkage(rabund, list, matrix, cutoff, method); }
 147                                 else if(method == "nearest"){   cluster = new SingleLinkage(rabund, list, matrix, cutoff, method); }
 148                                 else if(method == "average"){   cluster = new AverageLinkage(rabund, list, matrix, cutoff, method);     }
 149                                 else if(method == "weighted"){  cluster = new WeightedLinkage(rabund, list, matrix, cutoff, method);    }
 150                                 tag = cluster->getTag();
 151
 152                                 if (outputDir == "") { outputDir += m->hasPath(globaldata->inputFileName); }
 153                                 fileroot = outputDir + m->getRootName(m->getSimpleName(globaldata->inputFileName));
 154
 155                                 m->openOutputFile(fileroot+ tag + ".sabund",    sabundFile);
 156                                 m->openOutputFile(fileroot+ tag + ".rabund",    rabundFile);
 157                                 m->openOutputFile(fileroot+ tag + ".list",              listFile);
 158
 159                                 outputNames.push_back(fileroot+ tag + ".sabund"); outputTypes["sabund"].push_back(fileroot+ tag + ".sabund");
 160                                 outputNames.push_back(fileroot+ tag + ".rabund"); outputTypes["rabund"].push_back(fileroot+ tag + ".rabund");
 161                                 outputNames.push_back(fileroot+ tag + ".list"); outputTypes["list"].push_back(fileroot+ tag + ".list");
 162                         }
 163                 }
 164         }
 165         catch(exception& e) {
 166                 m->errorOut(e, "ClusterCommand", "ClusterCommand");
 167                 exit(1);
 168         }
 169 }
 170
 171 //**********************************************************************************************************************
 172
 173 void ClusterCommand::help(){
 174         try {
 175                 m->mothurOut("The cluster command can only be executed after a successful read.dist command.\n");
 176                 m->mothurOut("The cluster command parameter options are method, cuttoff, hard, precision, showabund and timing. No parameters are required.\n");
 177                 m->mothurOut("The cluster command should be in the following format: \n");
 178                 m->mothurOut("cluster(method=yourMethod, cutoff=yourCutoff, precision=yourPrecision) \n");
 179                 m->mothurOut("The acceptable cluster methods are furthest, nearest and average.  If no method is provided then furthest is assumed.\n\n");
 180         }
 181         catch(exception& e) {
 182                 m->errorOut(e, "ClusterCommand", "help");
 183                 exit(1);
 184         }
 185 }
 186
 187 //**********************************************************************************************************************
 188
 189 ClusterCommand::~ClusterCommand(){
 190         if (abort == false) {
 191                 delete cluster;
 192                 delete rabund;
 193         }
 194 }
 195
 196 //**********************************************************************************************************************
 197
 198 int ClusterCommand::execute(){
 199         try {
 200
 201                 if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
 202
 203                 time_t estart = time(NULL);
 204                 //int ndist = matrix->getNNodes();
 205                 float previousDist = 0.00000;
 206                 float rndPreviousDist = 0.00000;
 207                 oldRAbund = *rabund;
 208                 oldList = *list;
 209
 210                 print_start = true;
 211                 start = time(NULL);
 212                 loops = 0;
 213                 double saveCutoff = cutoff;
 214
 215                 while (matrix->getSmallDist() < cutoff && matrix->getNNodes() > 0){
 216
 217                         if (m->control_pressed) { //clean up
 218                                 delete globaldata->gSparseMatrix;  globaldata->gSparseMatrix = NULL;
 219                                 delete globaldata->gListVector;  globaldata->gListVector = NULL;
 220                                 if (globaldata->getFormat() == "phylip") { globaldata->setPhylipFile(""); }
 221                                 else if (globaldata->getFormat() == "column") { globaldata->setColumnFile(""); }
 222                                 sabundFile.close();rabundFile.close();listFile.close();
 223                                 for (int i = 0; i < outputNames.size(); i++) {  remove(outputNames[i].c_str());         } outputTypes.clear();
 224                                 return 0;
 225                         }
 226
 227                         if (print_start && m->isTrue(timing)) {
 228                                 m->mothurOut("Clustering (" + tag + ") dist " + toString(matrix->getSmallDist()) + "/"
 229                                         + toString(m->roundDist(matrix->getSmallDist(), precision))
 230                                         + "\t(precision: " + toString(precision) + ", Nodes: " + toString(matrix->getNNodes()) + ")");
 231                                 cout.flush();
 232                                 print_start = false;
 233                         }
 234
 235                         loops++;
 236
 237                         cluster->update(cutoff);
 238
 239                         float dist = matrix->getSmallDist();
 240                         float rndDist;
 241                         if (hard) {
 242                                 rndDist = m->ceilDist(dist, precision);
 243                         }else{
 244                                 rndDist = m->roundDist(dist, precision);
 245                         }
 246
 247                         if(previousDist <= 0.0000 && dist != previousDist){
 248                                 printData("unique");
 249                         }
 250                         else if(rndDist != rndPreviousDist){
 251                                 printData(toString(rndPreviousDist,  length-1));
 252                         }
 253
 254                         previousDist = dist;
 255                         rndPreviousDist = rndDist;
 256                         oldRAbund = *rabund;
 257                         oldList = *list;
 258                 }
 259
 260                 if (print_start && m->isTrue(timing)) {
 261                         m->mothurOut("Clustering (" + tag + ") for distance " + toString(previousDist) + "/" + toString(rndPreviousDist)
 262                                          + "\t(precision: " + toString(precision) + ", Nodes: " + toString(matrix->getNNodes()) + ")");
 263                         cout.flush();
 264                         print_start = false;
 265                 }
 266
 267                 if(previousDist <= 0.0000){
 268                         printData("unique");
 269                 }
 270                 else if(rndPreviousDist<cutoff){
 271                         printData(toString(rndPreviousDist, length-1));
 272                 }
 273
 274                 //delete globaldata's copy of the sparsematrix and listvector to free up memory
 275                 delete globaldata->gSparseMatrix;  globaldata->gSparseMatrix = NULL;
 276                 delete globaldata->gListVector;  globaldata->gListVector = NULL;
 277
 278                 //saves .list file so you can do the collect, rarefaction and summary commands without doing a read.list
 279                 if (globaldata->getFormat() == "phylip") { globaldata->setPhylipFile(""); }
 280                 else if (globaldata->getFormat() == "column") { globaldata->setColumnFile(""); }
 281
 282                 globaldata->setListFile(fileroot+ tag + ".list");
 283                 globaldata->setNameFile("");
 284                 globaldata->setFormat("list");
 285
 286                 sabundFile.close();
 287                 rabundFile.close();
 288                 listFile.close();
 289
 290                 if (saveCutoff != cutoff) {
 291                         if (hard)       {  saveCutoff = m->ceilDist(saveCutoff, precision);     }
 292                         else            {       saveCutoff = m->roundDist(saveCutoff, precision);  }
 293
 294                         m->mothurOut("changed cutoff to " + toString(cutoff)); m->mothurOutEndLine();
 295                 }
 296
 297                 //set list file as new current listfile
 298                 string current = "";
 299                 itTypes = outputTypes.find("list");
 300                 if (itTypes != outputTypes.end()) {
 301                         if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setListFile(current); }
 302                 }
 303
 304                 //set rabund file as new current rabundfile
 305                 itTypes = outputTypes.find("rabund");
 306                 if (itTypes != outputTypes.end()) {
 307                         if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setRabundFile(current); }
 308                 }
 309
 310                 //set sabund file as new current sabundfile
 311                 itTypes = outputTypes.find("sabund");
 312                 if (itTypes != outputTypes.end()) {
 313                         if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setSabundFile(current); }
 314                 }
 315
 316                 m->mothurOutEndLine();
 317                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 318                 for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
 319                 m->mothurOutEndLine();
 320
 321
 322                 //if (m->isTrue(timing)) {
 323                         m->mothurOut("It took " + toString(time(NULL) - estart) + " seconds to cluster"); m->mothurOutEndLine();
 324                 //}
 325
 326
 327                 return 0;
 328         }
 329         catch(exception& e) {
 330                 m->errorOut(e, "ClusterCommand", "execute");
 331                 exit(1);
 332         }
 333 }
 334
 335 //**********************************************************************************************************************
 336
 337 void ClusterCommand::printData(string label){
 338         try {
 339                 if (m->isTrue(timing)) {
 340                         m->mothurOut("\tTime: " + toString(time(NULL) - start) + "\tsecs for " + toString(oldRAbund.getNumBins())
 341                      + "\tclusters. Updates: " + toString(loops)); m->mothurOutEndLine();
 342                 }
 343                 print_start = true;
 344                 loops = 0;
 345                 start = time(NULL);
 346
 347                 oldRAbund.setLabel(label);
 348                 if (m->isTrue(showabund)) {
 349                         oldRAbund.getSAbundVector().print(cout);
 350                 }
 351                 oldRAbund.print(rabundFile);
 352                 oldRAbund.getSAbundVector().print(sabundFile);
 353
 354                 oldList.setLabel(label);
 355                 oldList.print(listFile);
 356         }
 357         catch(exception& e) {
 358                 m->errorOut(e, "ClusterCommand", "printData");
 359                 exit(1);
 360         }
 361
 362
 363 }
 364 //**********************************************************************************************************************