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