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