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