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 #include "readphylip.h"
  12 #include "readcolumn.h"
  13 #include "readmatrix.hpp"
  14 #include "clusterdoturcommand.h"
  15
  16 //**********************************************************************************************************************
  17 vector<string> ClusterCommand::setParameters(){
  18         try {
  19                 CommandParameter pphylip("phylip", "InputTypes", "", "", "PhylipColumn", "PhylipColumn", "none",false,false); parameters.push_back(pphylip);
  20                 CommandParameter pname("name", "InputTypes", "", "", "none", "none", "ColumnName",false,false); parameters.push_back(pname);
  21                 CommandParameter pcolumn("column", "InputTypes", "", "", "PhylipColumn", "PhylipColumn", "ColumnName",false,false); parameters.push_back(pcolumn);
  22                 CommandParameter pcutoff("cutoff", "Number", "", "10", "", "", "",false,false); parameters.push_back(pcutoff);
  23                 CommandParameter pprecision("precision", "Number", "", "100", "", "", "",false,false); parameters.push_back(pprecision);
  24                 CommandParameter pmethod("method", "Multiple", "furthest-nearest-average-weighted", "average", "", "", "",false,false); parameters.push_back(pmethod);
  25                 CommandParameter pshowabund("showabund", "Boolean", "", "T", "", "", "",false,false); parameters.push_back(pshowabund);
  26                 CommandParameter ptiming("timing", "Boolean", "", "F", "", "", "",false,false); parameters.push_back(ptiming);
  27                 CommandParameter psim("sim", "Boolean", "", "F", "", "", "",false,false); parameters.push_back(psim);
  28                 CommandParameter phard("hard", "Boolean", "", "T", "", "", "",false,false); parameters.push_back(phard);
  29                 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
  30                 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
  31
  32                 vector<string> myArray;
  33                 for (int i = 0; i < parameters.size(); i++) {   myArray.push_back(parameters[i].name);          }
  34                 return myArray;
  35         }
  36         catch(exception& e) {
  37                 m->errorOut(e, "ClusterCommand", "setParameters");
  38                 exit(1);
  39         }
  40 }
  41 //**********************************************************************************************************************
  42 string ClusterCommand::getHelpString(){
  43         try {
  44                 string helpString = "";
  45                 helpString += "The cluster command parameter options are phylip, column, name, method, cuttoff, hard, precision, sim, showabund and timing. Phylip or column and name are required, unless you have a valid current file.\n";
  46                 helpString += "The cluster command should be in the following format: \n";
  47                 helpString += "cluster(method=yourMethod, cutoff=yourCutoff, precision=yourPrecision) \n";
  48                 helpString += "The acceptable cluster methods are furthest, nearest, average and weighted.  If no method is provided then average is assumed.\n";
  49                 return helpString;
  50         }
  51         catch(exception& e) {
  52                 m->errorOut(e, "ClusterCommand", "getHelpString");
  53                 exit(1);
  54         }
  55 }
  56 //**********************************************************************************************************************
  57 ClusterCommand::ClusterCommand(){
  58         try {
  59                 abort = true; calledHelp = true;
  60                 setParameters();
  61                 vector<string> tempOutNames;
  62                 outputTypes["list"] = tempOutNames;
  63                 outputTypes["rabund"] = tempOutNames;
  64                 outputTypes["sabund"] = tempOutNames;
  65         }
  66         catch(exception& e) {
  67                 m->errorOut(e, "ClusterCommand", "ClusterCommand");
  68                 exit(1);
  69         }
  70 }
  71 //**********************************************************************************************************************
  72 //This function checks to make sure the cluster command has no errors and then clusters based on the method chosen.
  73 ClusterCommand::ClusterCommand(string option)  {
  74         try{
  75                 abort = false; calledHelp = false;
  76
  77                 //allow user to run help
  78                 if(option == "help") { help(); abort = true; calledHelp = true; }
  79                 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
  80
  81                 else {
  82                         vector<string> myArray = setParameters();
  83
  84                         OptionParser parser(option);
  85                         map<string,string> parameters = parser.getParameters();
  86                         map<string,string>::iterator it;
  87
  88                         ValidParameters validParameter;
  89
  90                         //check to make sure all parameters are valid for command
  91                         for (it = parameters.begin(); it != parameters.end(); it++) {
  92                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {
  93                                         abort = true;
  94                                 }
  95                         }
  96
  97                         //initialize outputTypes
  98                         vector<string> tempOutNames;
  99                         outputTypes["list"] = tempOutNames;
 100                         outputTypes["rabund"] = tempOutNames;
 101                         outputTypes["sabund"] = tempOutNames;
 102
 103                         //if the user changes the output directory command factory will send this info to us in the output parameter
 104                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){  outputDir = "";         }
 105
 106                         string inputDir = validParameter.validFile(parameters, "inputdir", false);
 107                         if (inputDir == "not found"){   inputDir = "";          }
 108                         else {
 109                                 string path;
 110                                 it = parameters.find("phylip");
 111                                 //user has given a template file
 112                                 if(it != parameters.end()){
 113                                         path = m->hasPath(it->second);
 114                                         //if the user has not given a path then, add inputdir. else leave path alone.
 115                                         if (path == "") {       parameters["phylip"] = inputDir + it->second;           }
 116                                 }
 117
 118                                 it = parameters.find("column");
 119                                 //user has given a template file
 120                                 if(it != parameters.end()){
 121                                         path = m->hasPath(it->second);
 122                                         //if the user has not given a path then, add inputdir. else leave path alone.
 123                                         if (path == "") {       parameters["column"] = inputDir + it->second;           }
 124                                 }
 125
 126                                 it = parameters.find("name");
 127                                 //user has given a template file
 128                                 if(it != parameters.end()){
 129                                         path = m->hasPath(it->second);
 130                                         //if the user has not given a path then, add inputdir. else leave path alone.
 131                                         if (path == "") {       parameters["name"] = inputDir + it->second;             }
 132                                 }
 133                         }
 134
 135                         //check for required parameters
 136                         phylipfile = validParameter.validFile(parameters, "phylip", true);
 137                         if (phylipfile == "not open") { phylipfile = ""; abort = true; }
 138                         else if (phylipfile == "not found") { phylipfile = ""; }
 139                         else {  distfile = phylipfile;  format = "phylip";      m->setPhylipFile(phylipfile); }
 140
 141                         columnfile = validParameter.validFile(parameters, "column", true);
 142                         if (columnfile == "not open") { columnfile = ""; abort = true; }
 143                         else if (columnfile == "not found") { columnfile = ""; }
 144                         else {  distfile = columnfile; format = "column"; m->setColumnFile(columnfile); }
 145
 146                         namefile = validParameter.validFile(parameters, "name", true);
 147                         if (namefile == "not open") { abort = true; }
 148                         else if (namefile == "not found") { namefile = ""; }
 149                         else { m->setNameFile(namefile); }
 150
 151                         if ((phylipfile == "") && (columnfile == "")) {
 152                                 //is there are current file available for either of these?
 153                                 //give priority to column, then phylip
 154                                 columnfile = m->getColumnFile();
 155                                 if (columnfile != "") {  distfile = columnfile; format = "column"; m->mothurOut("Using " + columnfile + " as input file for the column parameter."); m->mothurOutEndLine(); }
 156                                 else {
 157                                         phylipfile = m->getPhylipFile();
 158                                         if (phylipfile != "") { distfile = phylipfile;  format = "phylip"; m->mothurOut("Using " + phylipfile + " as input file for the phylip parameter."); m->mothurOutEndLine(); }
 159                                         else {
 160                                                 m->mothurOut("No valid current files. You must provide a phylip or column file before you can use the cluster command."); m->mothurOutEndLine();
 161                                                 abort = true;
 162                                         }
 163                                 }
 164                         }
 165                         else if ((phylipfile != "") && (columnfile != "")) { m->mothurOut("When executing a cluster command you must enter ONLY ONE of the following: phylip or column."); m->mothurOutEndLine(); abort = true; }
 166
 167                         if (columnfile != "") {
 168                                 if (namefile == "") {
 169                                         namefile = m->getNameFile();
 170                                         if (namefile != "") {  m->mothurOut("Using " + namefile + " as input file for the name parameter."); m->mothurOutEndLine(); }
 171                                         else {
 172                                                 m->mothurOut("You need to provide a namefile if you are going to use the column format."); m->mothurOutEndLine();
 173                                                 abort = true;
 174                                         }
 175                                 }
 176                         }
 177
 178                         //check for optional parameter and set defaults
 179                         // ...at some point should added some additional type checking...
 180                         //get user cutoff and precision or use defaults
 181                         string temp;
 182                         temp = validParameter.validFile(parameters, "precision", false);
 183                         if (temp == "not found") { temp = "100"; }
 184                         //saves precision legnth for formatting below
 185                         length = temp.length();
 186                         m->mothurConvert(temp, precision);
 187
 188                         temp = validParameter.validFile(parameters, "hard", false);                     if (temp == "not found") { temp = "T"; }
 189                         hard = m->isTrue(temp);
 190
 191                         temp = validParameter.validFile(parameters, "sim", false);                              if (temp == "not found") { temp = "F"; }
 192                         sim = m->isTrue(temp);
 193
 194                         temp = validParameter.validFile(parameters, "cutoff", false);
 195                         if (temp == "not found") { temp = "10"; }
 196                         m->mothurConvert(temp, cutoff);
 197                         cutoff += (5 / (precision * 10.0));
 198
 199                         method = validParameter.validFile(parameters, "method", false);
 200                         if (method == "not found") { method = "average"; }
 201
 202                         if ((method == "furthest") || (method == "nearest") || (method == "average") || (method == "weighted")) { }
 203                         else { m->mothurOut("Not a valid clustering method.  Valid clustering algorithms are furthest, nearest, average, and weighted."); m->mothurOutEndLine(); abort = true; }
 204
 205                         showabund = validParameter.validFile(parameters, "showabund", false);
 206                         if (showabund == "not found") { showabund = "T"; }
 207
 208                         timing = validParameter.validFile(parameters, "timing", false);
 209                         if (timing == "not found") { timing = "F"; }
 210
 211                 }
 212         }
 213         catch(exception& e) {
 214                 m->errorOut(e, "ClusterCommand", "ClusterCommand");
 215                 exit(1);
 216         }
 217 }
 218 //**********************************************************************************************************************
 219 ClusterCommand::~ClusterCommand(){}
 220 //**********************************************************************************************************************
 221
 222 int ClusterCommand::execute(){
 223         try {
 224
 225                 if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
 226
 227                 //phylip file given and cutoff not given - use cluster.classic because it uses less memory and is faster
 228                 if ((format == "phylip") && (cutoff > 10.0)) {
 229                         m->mothurOutEndLine(); m->mothurOut("You are using a phylip file and no cutoff.  I will run cluster.classic to save memory and time."); m->mothurOutEndLine();
 230
 231                         //run unique.seqs for deconvolute results
 232                         string inputString = "phylip=" + distfile;
 233                         if (namefile != "") { inputString += ", name=" + namefile; }
 234                         inputString += ", precision=" + toString(precision);
 235                         inputString += ", method=" + method;
 236                         if (hard)       { inputString += ", hard=T";    }
 237                         else            { inputString += ", hard=F";    }
 238                         if (sim)        { inputString += ", sim=T";             }
 239                         else            { inputString += ", sim=F";             }
 240
 241
 242                         m->mothurOutEndLine();
 243                         m->mothurOut("/------------------------------------------------------------/"); m->mothurOutEndLine();
 244                         m->mothurOut("Running command: cluster.classic(" + inputString + ")"); m->mothurOutEndLine();
 245
 246                         Command* clusterClassicCommand = new ClusterDoturCommand(inputString);
 247                         clusterClassicCommand->execute();
 248                         delete clusterClassicCommand;
 249
 250                         m->mothurOut("/------------------------------------------------------------/"); m->mothurOutEndLine();
 251
 252                         return 0;
 253                 }
 254
 255                 ReadMatrix* read;
 256                 if (format == "column") { read = new ReadColumnMatrix(columnfile, sim); }       //sim indicates whether its a similarity matrix
 257                 else if (format == "phylip") { read = new ReadPhylipMatrix(phylipfile, sim); }
 258
 259                 read->setCutoff(cutoff);
 260
 261                 NameAssignment* nameMap = NULL;
 262                 if(namefile != ""){
 263                         nameMap = new NameAssignment(namefile);
 264                         nameMap->readMap();
 265                 }
 266
 267                 read->read(nameMap);
 268                 list = read->getListVector();
 269                 matrix = read->getMatrix();
 270                 rabund = new RAbundVector(list->getRAbundVector());
 271                 delete read;
 272
 273                 if (m->control_pressed) { //clean up
 274                         delete list; delete matrix; delete rabund;
 275                         sabundFile.close();rabundFile.close();listFile.close();
 276                         for (int i = 0; i < outputNames.size(); i++) {  m->mothurRemove(outputNames[i]);        } outputTypes.clear();
 277                         return 0;
 278                 }
 279
 280                 //create cluster
 281                 if (method == "furthest")       {       cluster = new CompleteLinkage(rabund, list, matrix, cutoff, method); }
 282                 else if(method == "nearest"){   cluster = new SingleLinkage(rabund, list, matrix, cutoff, method); }
 283                 else if(method == "average"){   cluster = new AverageLinkage(rabund, list, matrix, cutoff, method);     }
 284                 else if(method == "weighted"){  cluster = new WeightedLinkage(rabund, list, matrix, cutoff, method);    }
 285                 tag = cluster->getTag();
 286
 287                 if (outputDir == "") { outputDir += m->hasPath(distfile); }
 288                 fileroot = outputDir + m->getRootName(m->getSimpleName(distfile));
 289
 290                 m->openOutputFile(fileroot+ tag + ".sabund",    sabundFile);
 291                 m->openOutputFile(fileroot+ tag + ".rabund",    rabundFile);
 292                 m->openOutputFile(fileroot+ tag + ".list",              listFile);
 293
 294                 outputNames.push_back(fileroot+ tag + ".sabund"); outputTypes["sabund"].push_back(fileroot+ tag + ".sabund");
 295                 outputNames.push_back(fileroot+ tag + ".rabund"); outputTypes["rabund"].push_back(fileroot+ tag + ".rabund");
 296                 outputNames.push_back(fileroot+ tag + ".list"); outputTypes["list"].push_back(fileroot+ tag + ".list");
 297
 298
 299                 time_t estart = time(NULL);
 300                 float previousDist = 0.00000;
 301                 float rndPreviousDist = 0.00000;
 302                 oldRAbund = *rabund;
 303                 oldList = *list;
 304
 305                 print_start = true;
 306                 start = time(NULL);
 307                 loops = 0;
 308                 double saveCutoff = cutoff;
 309
 310                 while (matrix->getSmallDist() < cutoff && matrix->getNNodes() > 0){
 311
 312                         if (m->control_pressed) { //clean up
 313                                 delete list; delete matrix; delete rabund; delete cluster;
 314                                 sabundFile.close();rabundFile.close();listFile.close();
 315                                 for (int i = 0; i < outputNames.size(); i++) {  m->mothurRemove(outputNames[i]);        } outputTypes.clear();
 316                                 return 0;
 317                         }
 318
 319                         if (print_start && m->isTrue(timing)) {
 320                                 m->mothurOut("Clustering (" + tag + ") dist " + toString(matrix->getSmallDist()) + "/"
 321                                         + toString(m->roundDist(matrix->getSmallDist(), precision))
 322                                         + "\t(precision: " + toString(precision) + ", Nodes: " + toString(matrix->getNNodes()) + ")");
 323                                 cout.flush();
 324                                 print_start = false;
 325                         }
 326
 327                         loops++;
 328
 329                         cluster->update(cutoff);
 330
 331                         float dist = matrix->getSmallDist();
 332                         float rndDist;
 333                         if (hard) {
 334                                 rndDist = m->ceilDist(dist, precision);
 335                         }else{
 336                                 rndDist = m->roundDist(dist, precision);
 337                         }
 338
 339                         if(previousDist <= 0.0000 && dist != previousDist){
 340                                 printData("unique");
 341                         }
 342                         else if(rndDist != rndPreviousDist){
 343                                 printData(toString(rndPreviousDist,  length-1));
 344                         }
 345
 346                         previousDist = dist;
 347                         rndPreviousDist = rndDist;
 348                         oldRAbund = *rabund;
 349                         oldList = *list;
 350                 }
 351
 352                 if (print_start && m->isTrue(timing)) {
 353                         m->mothurOut("Clustering (" + tag + ") for distance " + toString(previousDist) + "/" + toString(rndPreviousDist)
 354                                          + "\t(precision: " + toString(precision) + ", Nodes: " + toString(matrix->getNNodes()) + ")");
 355                         cout.flush();
 356                         print_start = false;
 357                 }
 358
 359                 if(previousDist <= 0.0000){
 360                         printData("unique");
 361                 }
 362                 else if(rndPreviousDist<cutoff){
 363                         printData(toString(rndPreviousDist, length-1));
 364                 }
 365
 366                 delete matrix;
 367                 delete list;
 368                 delete rabund;
 369                 delete cluster;
 370
 371                 sabundFile.close();
 372                 rabundFile.close();
 373                 listFile.close();
 374
 375                 if (saveCutoff != cutoff) {
 376                         if (hard)       {  saveCutoff = m->ceilDist(saveCutoff, precision);     }
 377                         else            {       saveCutoff = m->roundDist(saveCutoff, precision);  }
 378
 379                         m->mothurOut("changed cutoff to " + toString(cutoff)); m->mothurOutEndLine();
 380                 }
 381
 382                 //set list file as new current listfile
 383                 string current = "";
 384                 itTypes = outputTypes.find("list");
 385                 if (itTypes != outputTypes.end()) {
 386                         if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setListFile(current); }
 387                 }
 388
 389                 //set rabund file as new current rabundfile
 390                 itTypes = outputTypes.find("rabund");
 391                 if (itTypes != outputTypes.end()) {
 392                         if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setRabundFile(current); }
 393                 }
 394
 395                 //set sabund file as new current sabundfile
 396                 itTypes = outputTypes.find("sabund");
 397                 if (itTypes != outputTypes.end()) {
 398                         if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setSabundFile(current); }
 399                 }
 400
 401                 m->mothurOutEndLine();
 402                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 403                 for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
 404                 m->mothurOutEndLine();
 405
 406
 407                 //if (m->isTrue(timing)) {
 408                         m->mothurOut("It took " + toString(time(NULL) - estart) + " seconds to cluster"); m->mothurOutEndLine();
 409                 //}
 410
 411
 412                 return 0;
 413         }
 414         catch(exception& e) {
 415                 m->errorOut(e, "ClusterCommand", "execute");
 416                 exit(1);
 417         }
 418 }
 419
 420 //**********************************************************************************************************************
 421
 422 void ClusterCommand::printData(string label){
 423         try {
 424                 if (m->isTrue(timing)) {
 425                         m->mothurOut("\tTime: " + toString(time(NULL) - start) + "\tsecs for " + toString(oldRAbund.getNumBins())
 426                      + "\tclusters. Updates: " + toString(loops)); m->mothurOutEndLine();
 427                 }
 428                 print_start = true;
 429                 loops = 0;
 430                 start = time(NULL);
 431
 432                 oldRAbund.setLabel(label);
 433                 if (m->isTrue(showabund)) {
 434                         oldRAbund.getSAbundVector().print(cout);
 435                 }
 436                 oldRAbund.print(rabundFile);
 437                 oldRAbund.getSAbundVector().print(sabundFile);
 438
 439                 oldList.setLabel(label);
 440                 oldList.print(listFile);
 441         }
 442         catch(exception& e) {
 443                 m->errorOut(e, "ClusterCommand", "printData");
 444                 exit(1);
 445         }
 446
 447
 448 }
 449 //**********************************************************************************************************************