matrixoutputcommand.cpp

   1 /*
   2  *  matrixoutputcommand.cpp
   3  *  Mothur
   4  *
   5  *  Created by Sarah Westcott on 5/20/09.
   6  *  Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
   7  *
   8  */
   9
  10 #include "matrixoutputcommand.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 //**********************************************************************************************************************
  24
  25 MatrixOutputCommand::MatrixOutputCommand(string option)  {
  26         try {
  27                 globaldata = GlobalData::getInstance();
  28                 abort = false;
  29                 allLines = 1;
  30                 labels.clear();
  31                 Groups.clear();
  32                 Estimators.clear();
  33
  34                 //allow user to run help
  35                 if(option == "help") { validCalculator = new ValidCalculators(); help(); abort = true; }
  36
  37                 else {
  38                         //valid paramters for this command
  39                         string Array[] =  {"label","calc","groups","outputdir","inputdir", "output"};
  40                         vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
  41
  42                         OptionParser parser(option);
  43                         map<string,string> parameters  = parser.getParameters();
  44
  45                         ValidParameters validParameter;
  46
  47                         //check to make sure all parameters are valid for command
  48                         for (map<string,string>::iterator it = parameters.begin(); it != parameters.end(); it++) {
  49                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
  50                         }
  51
  52                         //if the user changes the output directory command factory will send this info to us in the output parameter
  53                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){
  54                                 outputDir = "";
  55                                 outputDir += hasPath(globaldata->inputFileName); //if user entered a file with a path then preserve it
  56                         }
  57
  58                         //make sure the user has already run the read.otu command
  59                         if (globaldata->getSharedFile() == "") {
  60                                 if (globaldata->getListFile() == "") { m->mothurOut("You must read a list and a group, or a shared before you can use the dist.shared command."); m->mothurOutEndLine(); abort = true; }
  61                                 else if (globaldata->getGroupFile() == "") { m->mothurOut("You must read a list and a group, or a shared before you can use the dist.shared command."); m->mothurOutEndLine(); abort = true; }
  62                         }
  63
  64                         //check for optional parameter and set defaults
  65                         // ...at some point should added some additional type checking...
  66                         label = validParameter.validFile(parameters, "label", false);
  67                         if (label == "not found") { label = ""; }
  68                         else {
  69                                 if(label != "all") {  splitAtDash(label, labels);  allLines = 0;  }
  70                                 else { allLines = 1;  }
  71                         }
  72
  73                         output = validParameter.validFile(parameters, "output", false);         if(output == "not found"){      output = "lt"; }
  74                         if ((output != "lt") && (output != "square")) { m->mothurOut(output + " is not a valid output form. Options are lt and square. I will use lt."); m->mothurOutEndLine(); output = "lt"; }
  75
  76                         //if the user has not specified any labels use the ones from read.otu
  77                         if (label == "") {
  78                                 allLines = globaldata->allLines;
  79                                 labels = globaldata->labels;
  80                         }
  81
  82                         groups = validParameter.validFile(parameters, "groups", false);
  83                         if (groups == "not found") { groups = ""; }
  84                         else {
  85                                 splitAtDash(groups, Groups);
  86                                 globaldata->Groups = Groups;
  87                         }
  88
  89                         calc = validParameter.validFile(parameters, "calc", false);
  90                         if (calc == "not found") { calc = "jclass-thetayc";  }
  91                         else {
  92                                  if (calc == "default")  {  calc = "jclass-thetayc";  }
  93                         }
  94                         splitAtDash(calc, Estimators);
  95
  96                         if (abort == false) {
  97
  98                                 validCalculator = new ValidCalculators();
  99
 100                                 int i;
 101                                 for (i=0; i<Estimators.size(); i++) {
 102                                         if (validCalculator->isValidCalculator("matrix", Estimators[i]) == true) {
 103                                                 if (Estimators[i] == "jabund") {
 104                                                         matrixCalculators.push_back(new JAbund());
 105                                                 }else if (Estimators[i] == "sorabund") {
 106                                                         matrixCalculators.push_back(new SorAbund());
 107                                                 }else if (Estimators[i] == "jclass") {
 108                                                         matrixCalculators.push_back(new Jclass());
 109                                                 }else if (Estimators[i] == "sorclass") {
 110                                                         matrixCalculators.push_back(new SorClass());
 111                                                 }else if (Estimators[i] == "jest") {
 112                                                         matrixCalculators.push_back(new Jest());
 113                                                 }else if (Estimators[i] == "sorest") {
 114                                                         matrixCalculators.push_back(new SorEst());
 115                                                 }else if (Estimators[i] == "thetayc") {
 116                                                         matrixCalculators.push_back(new ThetaYC());
 117                                                 }else if (Estimators[i] == "thetan") {
 118                                                         matrixCalculators.push_back(new ThetaN());
 119                                                 }else if (Estimators[i] == "morisitahorn") {
 120                                                         matrixCalculators.push_back(new MorHorn());
 121                                                 }else if (Estimators[i] == "braycurtis") {
 122                                                         matrixCalculators.push_back(new BrayCurtis());
 123                                                 }
 124                                         }
 125                                 }
 126
 127                         }
 128                 }
 129
 130         }
 131         catch(exception& e) {
 132                 m->errorOut(e, "MatrixOutputCommand", "MatrixOutputCommand");
 133                 exit(1);
 134         }
 135 }
 136
 137 //**********************************************************************************************************************
 138
 139 void MatrixOutputCommand::help(){
 140         try {
 141                 m->mothurOut("The dist.shared command can only be executed after a successful read.otu command.\n");
 142                 m->mothurOut("The dist.shared command parameters are groups, calc, output and label.  No parameters are required.\n");
 143                 m->mothurOut("The groups parameter allows you to specify which of the groups in your groupfile you would like included used.\n");
 144                 m->mothurOut("The group names are separated by dashes. The label parameter allows you to select what distance levels you would like distance matrices created for, and is also separated by dashes.\n");
 145                 m->mothurOut("The dist.shared command should be in the following format: dist.shared(groups=yourGroups, calc=yourCalcs, label=yourLabels).\n");
 146                 m->mothurOut("The output parameter allows you to specify format of your distance matrix. Options are lt, and square. The default is lt.\n");
 147                 m->mothurOut("Example dist.shared(groups=A-B-C, calc=jabund-sorabund).\n");
 148                 m->mothurOut("The default value for groups is all the groups in your groupfile.\n");
 149                 m->mothurOut("The default value for calc is jclass and thetayc.\n");
 150                 validCalculator->printCalc("matrix", cout);
 151                 m->mothurOut("The dist.shared command outputs a .dist file for each calculator you specify at each distance you choose.\n");
 152                 m->mothurOut("Note: No spaces between parameter labels (i.e. groups), '=' and parameters (i.e.yourGroups).\n\n");
 153         }
 154         catch(exception& e) {
 155                 m->errorOut(e, "MatrixOutputCommand", "help");
 156                 exit(1);
 157         }
 158 }
 159
 160
 161 //**********************************************************************************************************************
 162
 163 MatrixOutputCommand::~MatrixOutputCommand(){
 164         if (abort == false) {
 165                 delete input; globaldata->ginput = NULL;
 166                 delete read;
 167                 delete validCalculator;
 168         }
 169 }
 170
 171 //**********************************************************************************************************************
 172
 173 int MatrixOutputCommand::execute(){
 174         try {
 175
 176                 if (abort == true) {    return 0;       }
 177
 178                 //if the users entered no valid calculators don't execute command
 179                 if (matrixCalculators.size() == 0) { m->mothurOut("No valid calculators."); m->mothurOutEndLine();  return 0; }
 180
 181                 //you have groups
 182                 read = new ReadOTUFile(globaldata->inputFileName);
 183                 read->read(&*globaldata);
 184
 185                 input = globaldata->ginput;
 186                 lookup = input->getSharedRAbundVectors();
 187                 string lastLabel = lookup[0]->getLabel();
 188
 189                 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
 190                 set<string> processedLabels;
 191                 set<string> userLabels = labels;
 192
 193                 if (lookup.size() < 2) { m->mothurOut("You have not provided enough valid groups.  I cannot run the command."); m->mothurOutEndLine(); return 0;}
 194
 195                 numGroups = lookup.size();
 196
 197                 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;  }
 198
 199                 //as long as you are not at the end of the file or done wih the lines you want
 200                 while((lookup[0] != NULL) && ((allLines == 1) || (userLabels.size() != 0))) {
 201
 202                         if (m->control_pressed) { delete read; delete input; globaldata->ginput = NULL; for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  } for (int i = 0; i < outputNames.size(); i++) {  remove(outputNames[i].c_str()); } globaldata->Groups.clear(); return 0;  }
 203
 204                         if(allLines == 1 || labels.count(lookup[0]->getLabel()) == 1){
 205                                 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 206                                 process(lookup);
 207
 208                                 processedLabels.insert(lookup[0]->getLabel());
 209                                 userLabels.erase(lookup[0]->getLabel());
 210                         }
 211
 212                         if ((anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
 213                                 string saveLabel = lookup[0]->getLabel();
 214
 215                                 for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 216                                 lookup = input->getSharedRAbundVectors(lastLabel);
 217
 218                                 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 219                                 process(lookup);
 220
 221                                 processedLabels.insert(lookup[0]->getLabel());
 222                                 userLabels.erase(lookup[0]->getLabel());
 223
 224                                 //restore real lastlabel to save below
 225                                 lookup[0]->setLabel(saveLabel);
 226                         }
 227
 228                         lastLabel = lookup[0]->getLabel();
 229
 230                         //get next line to process
 231                         for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 232                         lookup = input->getSharedRAbundVectors();
 233                 }
 234
 235                 if (m->control_pressed) { delete read; delete input; globaldata->ginput = NULL; for (int i = 0; i < outputNames.size(); i++) {  remove(outputNames[i].c_str()); } globaldata->Groups.clear(); return 0;  }
 236
 237                 //output error messages about any remaining user labels
 238                 set<string>::iterator it;
 239                 bool needToRun = false;
 240                 for (it = userLabels.begin(); it != userLabels.end(); it++) {
 241                         m->mothurOut("Your file does not include the label " + *it);
 242                         if (processedLabels.count(lastLabel) != 1) {
 243                                 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
 244                                 needToRun = true;
 245                         }else {
 246                                 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
 247                         }
 248                 }
 249
 250                 if (m->control_pressed) { delete read; delete input; globaldata->ginput = NULL;  for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } globaldata->Groups.clear(); return 0;  }
 251
 252                 //run last label if you need to
 253                 if (needToRun == true)  {
 254                         for (int i = 0; i < lookup.size(); i++) {  if (lookup[i] != NULL) {  delete lookup[i]; }  }
 255                         lookup = input->getSharedRAbundVectors(lastLabel);
 256
 257                         m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 258                         process(lookup);
 259                         for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 260                 }
 261
 262                 if (m->control_pressed) { delete read; delete input; globaldata->ginput = NULL;  for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } globaldata->Groups.clear(); return 0;  }
 263
 264                 //reset groups parameter
 265                 globaldata->Groups.clear();
 266
 267                 m->mothurOutEndLine();
 268                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 269                 for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
 270                 m->mothurOutEndLine();
 271
 272
 273                 return 0;
 274         }
 275         catch(exception& e) {
 276                 m->errorOut(e, "MatrixOutputCommand", "execute");
 277                 exit(1);
 278         }
 279 }
 280 /***********************************************************/
 281 void MatrixOutputCommand::printSims(ostream& out) {
 282         try {
 283
 284                 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
 285
 286                 //output num seqs
 287                 out << simMatrix.size() << endl;
 288
 289                 if (output == "lt") {
 290                         for (int m = 0; m < simMatrix.size(); m++)      {
 291                                 out << lookup[m]->getGroup() << '\t';
 292                                 for (int n = 0; n < m; n++)     {
 293                                         out << simMatrix[m][n] << '\t';
 294                                 }
 295                                 out << endl;
 296                         }
 297                 }else{
 298                         for (int m = 0; m < simMatrix.size(); m++)      {
 299                                 out << lookup[m]->getGroup() << '\t';
 300                                 for (int n = 0; n < simMatrix[m].size(); n++)   {
 301                                         out << simMatrix[m][n] << '\t';
 302                                 }
 303                                 out << endl;
 304                         }
 305                 }
 306         }
 307         catch(exception& e) {
 308                 m->errorOut(e, "MatrixOutputCommand", "printSims");
 309                 exit(1);
 310         }
 311 }
 312 /***********************************************************/
 313 int MatrixOutputCommand::process(vector<SharedRAbundVector*> thisLookup){
 314         try {
 315
 316                                 EstOutput data;
 317                                 vector<SharedRAbundVector*> subset;
 318
 319                                 //for each calculator
 320                                 for(int i = 0 ; i < matrixCalculators.size(); i++) {
 321
 322                                         //initialize simMatrix
 323                                         simMatrix.clear();
 324                                         simMatrix.resize(numGroups);
 325                                         for (int p = 0; p < simMatrix.size(); p++)      {
 326                                                 for (int j = 0; j < simMatrix.size(); j++)      {
 327                                                         simMatrix[p].push_back(0.0);
 328                                                 }
 329                                         }
 330
 331                                         for (int k = 0; k < thisLookup.size(); k++) {
 332                                                 for (int l = k; l < thisLookup.size(); l++) {
 333                                                         if (k != l) { //we dont need to similiarity of a groups to itself
 334                                                                 //get estimated similarity between 2 groups
 335
 336                                                                 if (m->control_pressed) { return 0; }
 337
 338                                                                 subset.clear(); //clear out old pair of sharedrabunds
 339                                                                 //add new pair of sharedrabunds
 340                                                                 subset.push_back(thisLookup[k]); subset.push_back(thisLookup[l]);
 341
 342                                                                 data = matrixCalculators[i]->getValues(subset); //saves the calculator outputs
 343                                                                 //save values in similarity matrix
 344                                                                 simMatrix[k][l] = 1.0 - data[0];  //convert similiarity to distance
 345                                                                 simMatrix[l][k] = 1.0 - data[0];  //convert similiarity to distance
 346                                                         }
 347                                                 }
 348                                         }
 349
 350                                         exportFileName = outputDir + getRootName(getSimpleName(globaldata->inputFileName)) + matrixCalculators[i]->getName() + "." + thisLookup[0]->getLabel() + "." + output + ".dist";
 351                                         openOutputFile(exportFileName, out);
 352                                         outputNames.push_back(exportFileName);
 353
 354                                         printSims(out);
 355                                         out.close();
 356
 357                                 }
 358
 359                                 return 0;
 360
 361         }
 362         catch(exception& e) {
 363                 m->errorOut(e, "MatrixOutputCommand", "process");
 364                 exit(1);
 365         }
 366 }
 367 /***********************************************************/
 368
 369
 370