heatmapsimcommand.cpp

   1 /*
   2  *  heatmapsimcommand.cpp
   3  *  Mothur
   4  *
   5  *  Created by Sarah Westcott on 6/8/09.
   6  *  Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
   7  *
   8  */
   9
  10 #include "heatmapsimcommand.h"
  11 #include "sharedjabund.h"
  12 #include "sharedsorabund.h"
  13 #include "sharedjclass.h"
  14 #include "sharedsorclass.h"
  15 #include "sharedjest.h"
  16 #include "sharedsorest.h"
  17 #include "sharedthetayc.h"
  18 #include "sharedthetan.h"
  19 #include "sharedmorisitahorn.h"
  20 #include "sharedbraycurtis.h"
  21
  22 //**********************************************************************************************************************
  23 vector<string> HeatMapSimCommand::getValidParameters(){
  24         try {
  25                 string Array[] =  {"groups","label", "calc","phylip","column","name","outputdir","inputdir"};
  26                 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
  27                 return myArray;
  28         }
  29         catch(exception& e) {
  30                 m->errorOut(e, "HeatMapSimCommand", "getValidParameters");
  31                 exit(1);
  32         }
  33 }
  34 //**********************************************************************************************************************
  35 HeatMapSimCommand::HeatMapSimCommand(){
  36         try {
  37                 abort = true;
  38                 //initialize outputTypes
  39                 vector<string> tempOutNames;
  40                 outputTypes["svg"] = tempOutNames;
  41         }
  42         catch(exception& e) {
  43                 m->errorOut(e, "HeatMapSimCommand", "HeatMapSimCommand");
  44                 exit(1);
  45         }
  46 }
  47 //**********************************************************************************************************************
  48 vector<string> HeatMapSimCommand::getRequiredParameters(){
  49         try {
  50                 vector<string> myArray;
  51                 return myArray;
  52         }
  53         catch(exception& e) {
  54                 m->errorOut(e, "HeatMapSimCommand", "getRequiredParameters");
  55                 exit(1);
  56         }
  57 }
  58 //**********************************************************************************************************************
  59 vector<string> HeatMapSimCommand::getRequiredFiles(){
  60         try {
  61                 vector<string> myArray;
  62                 return myArray;
  63         }
  64         catch(exception& e) {
  65                 m->errorOut(e, "HeatMapSimCommand", "getRequiredFiles");
  66                 exit(1);
  67         }
  68 }
  69 //**********************************************************************************************************************
  70
  71 HeatMapSimCommand::HeatMapSimCommand(string option)  {
  72         try {
  73                 globaldata = GlobalData::getInstance();
  74                 abort = false;
  75                 allLines = 1;
  76                 labels.clear();
  77                 Groups.clear();
  78                 Estimators.clear();
  79
  80                 //allow user to run help
  81                 if(option == "help") { validCalculator = new ValidCalculators(); help(); abort = true; }
  82
  83                 else {
  84                         //valid paramters for this command
  85                         string AlignArray[] =  {"groups","label", "calc","phylip","column","name","outputdir","inputdir"};
  86                         vector<string> myArray (AlignArray, AlignArray+(sizeof(AlignArray)/sizeof(string)));
  87
  88                         OptionParser parser(option);
  89                         map<string,string> parameters = parser.getParameters();
  90
  91                         ValidParameters validParameter;
  92                         map<string,string>::iterator it;
  93
  94                         //check to make sure all parameters are valid for command
  95                         for (it = parameters.begin(); it != parameters.end(); it++) {
  96                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
  97                         }
  98
  99                         //initialize outputTypes
 100                         vector<string> tempOutNames;
 101                         outputTypes["svg"] = tempOutNames;
 102
 103                         format = "";
 104                         //if the user changes the output directory command factory will send this info to us in the output parameter
 105                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){  outputDir = "";         }
 106
 107                         //if the user changes the input directory command factory will send this info to us in the output parameter
 108                         string inputDir = validParameter.validFile(parameters, "inputdir", false);
 109                         if (inputDir == "not found"){   inputDir = "";          }
 110                         else {
 111                                 string path;
 112                                 it = parameters.find("phylip");
 113                                 //user has given a template file
 114                                 if(it != parameters.end()){
 115                                         path = m->hasPath(it->second);
 116                                         //if the user has not given a path then, add inputdir. else leave path alone.
 117                                         if (path == "") {       parameters["phylip"] = inputDir + it->second;           }
 118                                 }
 119
 120                                 it = parameters.find("column");
 121                                 //user has given a template file
 122                                 if(it != parameters.end()){
 123                                         path = m->hasPath(it->second);
 124                                         //if the user has not given a path then, add inputdir. else leave path alone.
 125                                         if (path == "") {       parameters["column"] = inputDir + it->second;           }
 126                                 }
 127
 128                                 it = parameters.find("name");
 129                                 //user has given a template file
 130                                 if(it != parameters.end()){
 131                                         path = m->hasPath(it->second);
 132                                         //if the user has not given a path then, add inputdir. else leave path alone.
 133                                         if (path == "") {       parameters["name"] = inputDir + it->second;             }
 134                                 }
 135                         }
 136
 137                         //required parameters
 138                         phylipfile = validParameter.validFile(parameters, "phylip", true);
 139                         if (phylipfile == "not open") { abort = true; }
 140                         else if (phylipfile == "not found") { phylipfile = ""; }
 141                         else {  format = "phylip";      if (outputDir == "") { outputDir += m->hasPath(phylipfile); }  }
 142
 143                         columnfile = validParameter.validFile(parameters, "column", true);
 144                         if (columnfile == "not open") { abort = true; }
 145                         else if (columnfile == "not found") { columnfile = ""; }
 146                         else {  format = "column";      if (outputDir == "") { outputDir += m->hasPath(columnfile); } }
 147
 148                         namefile = validParameter.validFile(parameters, "name", true);
 149                         if (namefile == "not open") { abort = true; }
 150                         else if (namefile == "not found") { namefile = ""; }
 151
 152
 153                         //error checking on files
 154                         if ((globaldata->getSharedFile() == "") && ((phylipfile == "") && (columnfile == "")))  { m->mothurOut("You must run the read.otu command or provide a distance file before running the heatmap.sim command."); m->mothurOutEndLine(); abort = true; }
 155                         else if ((phylipfile != "") && (columnfile != "")) { m->mothurOut("When running the heatmap.sim command with a distance file you may not use both the column and the phylip parameters."); m->mothurOutEndLine(); abort = true; }
 156
 157                         if (columnfile != "") {
 158                                 if (namefile == "") {  m->mothurOut("You need to provide a namefile if you are going to use the column format."); m->mothurOutEndLine(); abort = true; }
 159                         }
 160
 161                         if (format == "") { format = "shared"; }
 162
 163                         //check for optional parameter and set defaults
 164                         // ...at some point should added some additional type checking...
 165                         if (format == "shared") {
 166                                 if (outputDir == "") { outputDir += m->hasPath(globaldata->getSharedFile()); }
 167
 168                                 label = validParameter.validFile(parameters, "label", false);
 169                                 if (label == "not found") { label = ""; }
 170                                 else {
 171                                         if(label != "all") {  m->splitAtDash(label, labels);  allLines = 0;  }
 172                                         else { allLines = 1;  }
 173                                 }
 174
 175                                 //if the user has not specified any labels use the ones from read.otu
 176                                 if (label == "") {
 177                                         allLines = globaldata->allLines;
 178                                         labels = globaldata->labels;
 179                                 }
 180
 181                                 calc = validParameter.validFile(parameters, "calc", false);
 182                                 if (calc == "not found") { calc = "jest-thetayc";  }
 183                                 else {
 184                                         if (calc == "default")  {  calc = "jest-thetayc";  }
 185                                 }
 186                                 m->splitAtDash(calc, Estimators);
 187
 188                                 groups = validParameter.validFile(parameters, "groups", false);
 189                                 if (groups == "not found") { groups = ""; }
 190                                 else {
 191                                         m->splitAtDash(groups, Groups);
 192                                         globaldata->Groups = Groups;
 193                                 }
 194                         }
 195
 196                         if (abort == false) {
 197                                 validCalculator = new ValidCalculators();
 198
 199                                 int i;
 200                                 for (i=0; i<Estimators.size(); i++) {
 201                                         if (validCalculator->isValidCalculator("heat", Estimators[i]) == true) {
 202                                                 if (Estimators[i] == "jabund") {
 203                                                         heatCalculators.push_back(new JAbund());
 204                                                 }else if (Estimators[i] == "sorabund") {
 205                                                         heatCalculators.push_back(new SorAbund());
 206                                                 }else if (Estimators[i] == "jclass") {
 207                                                         heatCalculators.push_back(new Jclass());
 208                                                 }else if (Estimators[i] == "sorclass") {
 209                                                         heatCalculators.push_back(new SorClass());
 210                                                 }else if (Estimators[i] == "jest") {
 211                                                         heatCalculators.push_back(new Jest());
 212                                                 }else if (Estimators[i] == "sorest") {
 213                                                         heatCalculators.push_back(new SorEst());
 214                                                 }else if (Estimators[i] == "thetayc") {
 215                                                         heatCalculators.push_back(new ThetaYC());
 216                                                 }else if (Estimators[i] == "thetan") {
 217                                                         heatCalculators.push_back(new ThetaN());
 218                                                 }else if (Estimators[i] == "morisitahorn") {
 219                                                         heatCalculators.push_back(new MorHorn());
 220                                                 }else if (Estimators[i] == "braycurtis") {
 221                                                         heatCalculators.push_back(new BrayCurtis());
 222                                                 }
 223                                         }
 224                                 }
 225
 226                         }
 227                 }
 228
 229
 230
 231         }
 232         catch(exception& e) {
 233                 m->errorOut(e, "HeatMapSimCommand", "HeatMapSimCommand");
 234                 exit(1);
 235         }
 236 }
 237
 238 //**********************************************************************************************************************
 239
 240 void HeatMapSimCommand::help(){
 241         try {
 242                 m->mothurOut("The heatmap.sim command can only be executed after a successful read.otu command, or by providing a distance file.\n");
 243                 m->mothurOut("The heatmap.sim command parameters are phylip, column, name, groups, calc and label.  No parameters are required.\n");
 244                 m->mothurOut("There are two ways to use the heatmap.sim command. The first is with the read.otu command. \n");
 245                 m->mothurOut("With the read.otu command you may use the groups, label and calc parameters. \n");
 246                 m->mothurOut("The groups parameter allows you to specify which of the groups in your groupfile you would like included in your heatmap.\n");
 247                 m->mothurOut("The group names are separated by dashes. The label parameter allows you to select what distance levels you would like a heatmap created for, and is also separated by dashes.\n");
 248                 m->mothurOut("The heatmap.sim command should be in the following format: heatmap.sim(groups=yourGroups, calc=yourCalc, label=yourLabels).\n");
 249                 m->mothurOut("Example heatmap.sim(groups=A-B-C, calc=jabund).\n");
 250                 m->mothurOut("The default value for groups is all the groups in your groupfile, and all labels in your inputfile will be used.\n");
 251                 validCalculator->printCalc("heat", cout);
 252                 m->mothurOut("The default value for calc is jclass-thetayc.\n");
 253                 m->mothurOut("The heatmap.sim command outputs a .svg file for each calculator you choose at each label you specify.\n");
 254                 m->mothurOut("The second way to use the heatmap.sim command is with a distance file representing the distance bewteen your groups. \n");
 255                 m->mothurOut("Using the command this way, the phylip or column parameter are required, and only one may be used.  If you use a column file the name filename is required. \n");
 256                 m->mothurOut("The heatmap.sim command should be in the following format: heatmap.sim(phylip=yourDistanceFile).\n");
 257                 m->mothurOut("Example heatmap.sim(phylip=amazonGroups.dist).\n");
 258                 m->mothurOut("Note: No spaces between parameter labels (i.e. groups), '=' and parameters (i.e.yourGroups).\n\n");
 259
 260         }
 261         catch(exception& e) {
 262                 m->errorOut(e, "HeatMapSimCommand", "help");
 263                 exit(1);
 264         }
 265 }
 266
 267 //**********************************************************************************************************************
 268
 269 HeatMapSimCommand::~HeatMapSimCommand(){}
 270
 271 //**********************************************************************************************************************
 272
 273 int HeatMapSimCommand::execute(){
 274         try {
 275
 276                 if (abort == true)  { return 0; }
 277
 278                 heatmap = new HeatMapSim(outputDir);
 279
 280                 if (format == "shared") {
 281                         runCommandShared();
 282                 }else if (format == "phylip") {
 283                         globaldata->inputFileName = phylipfile;
 284                         runCommandDist();
 285                 }else if (format == "column") {
 286                         globaldata->inputFileName = columnfile;
 287                         runCommandDist();
 288                 }
 289
 290                 delete heatmap;
 291                 delete validCalculator;
 292
 293                 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) {        remove(outputNames[i].c_str());  } outputTypes.clear(); return 0; }
 294
 295                 m->mothurOutEndLine();
 296                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 297                 for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
 298                 m->mothurOutEndLine();
 299
 300                 return 0;
 301         }
 302         catch(exception& e) {
 303                 m->errorOut(e, "HeatMapSimCommand", "execute");
 304                 exit(1);
 305         }
 306 }
 307
 308 //**********************************************************************************************************************
 309 int HeatMapSimCommand::runCommandShared() {
 310         try {
 311                 //if the users entered no valid calculators don't execute command
 312                 if (heatCalculators.size() == 0) { m->mothurOut("No valid calculators."); m->mothurOutEndLine(); return 0; }
 313
 314                 //you have groups
 315                 read = new ReadOTUFile(globaldata->inputFileName);
 316                 read->read(&*globaldata);
 317
 318                 input = globaldata->ginput;
 319                 lookup = input->getSharedRAbundVectors();
 320                 string lastLabel = lookup[0]->getLabel();
 321
 322                 if (lookup.size() < 2) { m->mothurOut("You have not provided enough valid groups.  I cannot run the command."); m->mothurOutEndLine(); return 0;}
 323
 324                 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
 325                 set<string> processedLabels;
 326                 set<string> userLabels = labels;
 327
 328                 if (m->control_pressed) { delete read; delete input; globaldata->ginput = NULL; for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }  globaldata->Groups.clear(); return 0; }
 329
 330                 //as long as you are not at the end of the file or done wih the lines you want
 331                 while((lookup[0] != NULL) && ((allLines == 1) || (userLabels.size() != 0))) {
 332
 333                         if (m->control_pressed) { delete read; delete input; globaldata->ginput = NULL; for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  } globaldata->Groups.clear(); return 0; }
 334
 335                         if(allLines == 1 || labels.count(lookup[0]->getLabel()) == 1){
 336
 337                                 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 338                                 vector<string> outfilenames = heatmap->getPic(lookup, heatCalculators);
 339                                 for(int i = 0; i < outfilenames.size(); i++) { outputNames.push_back(outfilenames[i]);  outputTypes["svg"].push_back(outfilenames[i]); }
 340
 341                                 processedLabels.insert(lookup[0]->getLabel());
 342                                 userLabels.erase(lookup[0]->getLabel());
 343                         }
 344
 345                         if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
 346                                 string saveLabel = lookup[0]->getLabel();
 347
 348                                 for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 349                                 lookup = input->getSharedRAbundVectors(lastLabel);
 350
 351                                 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 352                                 vector<string> outfilenames = heatmap->getPic(lookup, heatCalculators);
 353                                 for(int i = 0; i < outfilenames.size(); i++) { outputNames.push_back(outfilenames[i]); outputTypes["svg"].push_back(outfilenames[i]);  }
 354
 355                                 processedLabels.insert(lookup[0]->getLabel());
 356                                 userLabels.erase(lookup[0]->getLabel());
 357
 358                                 //restore real lastlabel to save below
 359                                 lookup[0]->setLabel(saveLabel);
 360                         }
 361
 362                         //prevent memory leak
 363
 364                         lastLabel = lookup[0]->getLabel();
 365
 366                         //get next line to process
 367                         for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 368                         lookup = input->getSharedRAbundVectors();
 369
 370                 }
 371
 372
 373                 if (m->control_pressed) { delete read; delete input; globaldata->ginput = NULL; globaldata->Groups.clear();  return 0; }
 374
 375                 //output error messages about any remaining user labels
 376                 set<string>::iterator it;
 377                 bool needToRun = false;
 378                 for (it = userLabels.begin(); it != userLabels.end(); it++) {
 379                         m->mothurOut("Your file does not include the label " + *it);
 380                         if (processedLabels.count(lastLabel) != 1) {
 381                                 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
 382                                 needToRun = true;
 383                         }else {
 384                                 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
 385                         }
 386                 }
 387
 388                 if (m->control_pressed) { delete read; delete input; globaldata->ginput = NULL;  globaldata->Groups.clear(); return 0; }
 389
 390                 //run last label if you need to
 391                 if (needToRun == true)  {
 392                         for (int i = 0; i < lookup.size(); i++) {  if (lookup[i] != NULL) { delete lookup[i]; } }
 393                         lookup = input->getSharedRAbundVectors(lastLabel);
 394
 395                         m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 396                         vector<string> outfilenames = heatmap->getPic(lookup, heatCalculators);
 397                         for(int i = 0; i < outfilenames.size(); i++) { outputNames.push_back(outfilenames[i]); outputTypes["svg"].push_back(outfilenames[i]);  }
 398
 399                         for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 400                 }
 401
 402                 if (m->control_pressed) { delete read; delete input; globaldata->ginput = NULL;  globaldata->Groups.clear(); return 0; }
 403
 404                 //reset groups parameter
 405                 globaldata->Groups.clear();
 406
 407                 delete input;  globaldata->ginput = NULL;
 408                 delete read;
 409
 410                 return 0;
 411         }
 412         catch(exception& e) {
 413                 m->errorOut(e, "HeatMapSimCommand", "runCommandShared");
 414                 exit(1);
 415         }
 416 }
 417 //**********************************************************************************************************************
 418 int HeatMapSimCommand::runCommandDist() {
 419         try {
 420
 421                 vector< vector<double> > matrix;
 422                 vector<string> names;
 423                 ifstream in;
 424
 425                 //read distance file and create distance vector and names vector
 426                 if (format == "phylip") {
 427                         //read phylip file
 428                         m->openInputFile(phylipfile, in);
 429
 430                         string name;
 431                         int numSeqs;
 432                         in >> numSeqs >> name;
 433
 434                         //save name
 435                         names.push_back(name);
 436
 437                         //resize the matrix and fill with zeros
 438                         matrix.resize(numSeqs);
 439                         for(int i = 0; i < numSeqs; i++) {
 440                                 matrix[i].resize(numSeqs, 0.0);
 441                         }
 442
 443                         //determine if matrix is square or lower triangle
 444                         //if it is square read the distances for the first sequence
 445                         char d;
 446                         bool square;
 447                         while((d=in.get()) != EOF){
 448
 449                                 //is d a number meaning its square
 450                                 if(isalnum(d)){
 451                                         square = true;
 452                                         in.putback(d);
 453
 454                                         for(int i=0;i<numSeqs;i++){
 455                                                 in >> matrix[0][i];
 456                                         }
 457                                         break;
 458                                 }
 459
 460                                 //is d a line return meaning its lower triangle
 461                                 if(d == '\n'){
 462                                         square = false;
 463                                         break;
 464                                 }
 465                         }
 466
 467                         //read rest of matrix
 468                         if (square == true) {
 469                                 for(int i=1;i<numSeqs;i++){
 470                                         in >> name;
 471                                         names.push_back(name);
 472
 473                                         if (m->control_pressed) { return 0; }
 474
 475                                         for(int j=0;j<numSeqs;j++) { in >> matrix[i][j];  }
 476                                         m->gobble(in);
 477                                 }
 478                         }else {
 479                                 double dist;
 480                                 for(int i=1;i<numSeqs;i++){
 481                                         in >> name;
 482                                         names.push_back(name);
 483
 484                                         if (m->control_pressed) { return 0; }
 485
 486                                         for(int j=0;j<i;j++){
 487                                                 in >> dist;
 488                                                 matrix[i][j] = dist;  matrix[j][i] = dist;
 489                                         }
 490                                         m->gobble(in);
 491                                 }
 492                         }
 493                         in.close();
 494                 }else {
 495                         //read names file
 496                         NameAssignment* nameMap = new NameAssignment(namefile);
 497                         nameMap->readMap();
 498
 499                         //put names in order in vector
 500                         for (int i = 0; i < nameMap->size(); i++) {
 501                                 names.push_back(nameMap->get(i));
 502                         }
 503
 504                         //resize matrix
 505                         matrix.resize(nameMap->size());
 506                         for (int i = 0; i < nameMap->size(); i++) {
 507                                 matrix[i].resize(nameMap->size(), 0.0);
 508                         }
 509
 510                         //read column file
 511                         string first, second;
 512                         double dist;
 513                         m->openInputFile(columnfile, in);
 514
 515                         while (!in.eof()) {
 516                                 in >> first >> second >> dist; m->gobble(in);
 517
 518                                 if (m->control_pressed) { return 0; }
 519
 520                                 map<string, int>::iterator itA = nameMap->find(first);
 521                                 map<string, int>::iterator itB = nameMap->find(second);
 522
 523                                 if(itA == nameMap->end()){  cerr << "AAError: Sequence '" << first << "' was not found in the names file, please correct\n"; exit(1);  }
 524                                 if(itB == nameMap->end()){  cerr << "ABError: Sequence '" << second << "' was not found in the names file, please correct\n"; exit(1);  }
 525
 526                                 //save distance
 527                                 matrix[itA->second][itB->second] = dist;
 528                                 matrix[itB->second][itA->second] = dist;
 529                         }
 530                         in.close();
 531
 532                         delete nameMap;
 533                 }
 534
 535
 536                 string outputFileName = heatmap->getPic(matrix, names);
 537                 outputNames.push_back(outputFileName); //vector<vector<double>>, vector<string>
 538                 outputTypes["svg"].push_back(outputFileName);
 539
 540                 return 0;
 541         }
 542         catch(exception& e) {
 543                 m->errorOut(e, "HeatMapSimCommand", "runCommandDist");
 544                 exit(1);
 545         }
 546 }
 547 //**********************************************************************************************************************
 548
 549
 550
 551
 552
 553