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