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