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