lefsecommand.cpp

   1 //
   2 //  lefsecommand.cpp
   3 //  Mothur
   4 //
   5 //  Created by SarahsWork on 6/12/13.
   6 //  Copyright (c) 2013 Schloss Lab. All rights reserved.
   7 //
   8
   9 #include "lefsecommand.h"
  10 #include "linearalgebra.h"
  11
  12 //**********************************************************************************************************************
  13 vector<string> LefseCommand::setParameters(){
  14         try {
  15         CommandParameter pdesign("design", "InputTypes", "", "", "none", "none", "none","",false,true,true); parameters.push_back(pdesign);
  16         CommandParameter pshared("shared", "InputTypes", "", "", "none", "none", "none","summary",false,true,true); parameters.push_back(pshared);
  17         CommandParameter pclass("class", "String", "", "", "", "", "","",false,false); parameters.push_back(pclass);
  18         CommandParameter psubclass("subclass", "String", "", "", "", "", "","",false,false); parameters.push_back(psubclass);
  19                 CommandParameter plabel("label", "String", "", "", "", "", "","",false,false); parameters.push_back(plabel);
  20                 //CommandParameter pclasses("classes", "String", "", "", "", "", "","",false,false); parameters.push_back(pclasses);
  21         CommandParameter palpha("aalpha", "Number", "", "0.05", "", "", "","",false,false); parameters.push_back(palpha);
  22         CommandParameter pwalpha("walpha", "Number", "", "0.05", "", "", "","",false,false); parameters.push_back(pwalpha);
  23
  24         CommandParameter plda("lda", "Number", "", "2.0", "", "", "","",false,false); parameters.push_back(plda);
  25         CommandParameter pwilc("wilc", "Boolean", "", "T", "", "", "","",false,false); parameters.push_back(pwilc);
  26         CommandParameter pnormmillion("norm", "Boolean", "", "T", "", "", "","",false,false); parameters.push_back(pnormmillion);
  27         CommandParameter piters("iters", "Number", "", "30", "", "", "","",false,false); parameters.push_back(piters);
  28         //CommandParameter pwilcsamename("wilcsamename", "Boolean", "", "F", "", "", "","",false,false); parameters.push_back(pwilcsamename);
  29         CommandParameter pcurv("curv", "Boolean", "", "F", "", "", "","",false,false); parameters.push_back(pcurv);
  30         CommandParameter pfiters("fboots", "Number", "", "0.67", "", "", "","",false,false); parameters.push_back(pfiters);
  31         CommandParameter pstrict("strict", "Multiple", "0-1-2", "0", "", "", "","",false,false); parameters.push_back(pstrict);
  32         CommandParameter pminc("minc", "Number", "", "10", "", "", "","",false,false); parameters.push_back(pminc);
  33         CommandParameter pmulticlass_strat("multiclass", "Multiple", "onevone-onevall", "onevall", "", "", "","",false,false); parameters.push_back(pmulticlass_strat);
  34         //CommandParameter psubject("subject", "Boolean", "", "F", "", "", "","",false,false); parameters.push_back(psubject);
  35
  36
  37         //not used in their current code, but in parameters
  38         //CommandParameter pnlogs("nlogs", "Number", "", "3", "", "", "","",false,false); parameters.push_back(pnlogs);
  39         //CommandParameter pranktec("ranktec", "Multiple", "lda-svm", "lda", "", "", "","",false,false); parameters.push_back(pranktec); // svm not implemented in their source yet.
  40         //CommandParameter psvmnorm("svmnorm", "Boolean", "", "T", "", "", "","",false,false); parameters.push_back(psvmnorm); //not used because svm not implemented yet.
  41
  42
  43         //every command must have inputdir and outputdir.  This allows mothur users to redirect input and output files.
  44                 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(pinputdir);
  45                 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(poutputdir);
  46
  47                 vector<string> myArray;
  48                 for (int i = 0; i < parameters.size(); i++) {   myArray.push_back(parameters[i].name);          }
  49                 return myArray;
  50         }
  51         catch(exception& e) {
  52                 m->errorOut(e, "LefseCommand", "setParameters");
  53                 exit(1);
  54         }
  55 }
  56 //**********************************************************************************************************************
  57 string LefseCommand::getHelpString(){
  58         try {
  59                 string helpString = "";
  60                 helpString += "The lefse command allows you to ....\n";
  61                 helpString += "The lefse command parameters are: shared, design, class, subclass, label, walpha, aalpha, lda, wilc, iters, curv, fboots, strict, minc, multiclass and norm.\n";
  62                 helpString += "The class parameter is used to indicate the which category you would like used for the Kruskal Wallis analysis. If none is provided first category is used.\n";
  63         helpString += "The subclass parameter is used to indicate the .....If none is provided, second category is used, or if only one category subclass is ignored. \n";
  64         helpString += "The aalpha parameter is used to set the alpha value for the Krukal Wallis Anova test Default=0.05. \n";
  65         helpString += "The walpha parameter is used to set the alpha value for the Wilcoxon test. Default=0.05. \n";
  66         helpString += "The lda parameter is used to set the threshold on the absolute value of the logarithmic LDA score. Default=2.0. \n";
  67         helpString += "The wilc parameter is used to indicate whether to perform the Wilcoxon test. Default=T. \n";
  68         helpString += "The iters parameter is used to set the number of bootstrap iteration for LDA. Default=30. \n";
  69         //helpString += "The wilcsamename parameter is used to indicate whether perform the wilcoxon test only among the subclasses with the same name. Default=F. \n";
  70         helpString += "The curv parameter is used to set whether perform the wilcoxon testing the Curtis's approach [BETA VERSION] Default=F. \n";
  71         helpString += "The norm parameter is used to multiply relative abundances by 1000000. Recommended when very low values are present. Default=T. \n";
  72         helpString += "The fboots parameter is used to set the subsampling fraction value for each bootstrap iteration. Default=0.67. \n";
  73         helpString += "The strict parameter is used to set the multiple testing correction options. 0 no correction (more strict, default), 1 correction for independent comparisons, 2 correction for independent comparison. Options = 0,1,2. Default=0. \n";
  74         helpString += "The minc parameter is used to minimum number of samples per subclass for performing wilcoxon test. Default=10. \n";
  75         helpString += "The multiclass parameter is used to (for multiclass tasks) set whether the test is performed in a one-against-one ( onevone - more strict!) or in a one-against-all setting ( onevall - less strict). Default=onevall. \n";
  76         //helpString += "The classes parameter is used to indicate the classes you would like to use. Classes should be inputted in the following format: classes=label<value1|value2|value3>-label<value1|value2>. For example to include groups from treatment with value early or late and age= young or old.  class=treatment<Early|Late>-age<young|old>.\n";
  77         helpString += "The label parameter is used to indicate which distances in the shared file you would like to use. labels are separated by dashes.\n";
  78                 helpString += "The lefse command should be in the following format: lefse(shared=final.an.shared, design=final.design, class=treatment, subclass=age).\n";
  79         return helpString;
  80         }
  81         catch(exception& e) {
  82                 m->errorOut(e, "LefseCommand", "getHelpString");
  83                 exit(1);
  84         }
  85 }
  86 //**********************************************************************************************************************
  87 string LefseCommand::getOutputPattern(string type) {
  88     try {
  89         string pattern = "";
  90
  91         if (type == "summary") {  pattern = "[filename],[distance],lefse_summary"; }
  92         else { m->mothurOut("[ERROR]: No definition for type " + type + " output pattern.\n"); m->control_pressed = true;  }
  93
  94         return pattern;
  95     }
  96     catch(exception& e) {
  97         m->errorOut(e, "LefseCommand", "getOutputPattern");
  98         exit(1);
  99     }
 100 }
 101 //**********************************************************************************************************************
 102 LefseCommand::LefseCommand(){
 103         try {
 104                 abort = true; calledHelp = true;
 105                 setParameters();
 106         vector<string> tempOutNames;
 107                 outputTypes["summary"] = tempOutNames;
 108         }
 109         catch(exception& e) {
 110                 m->errorOut(e, "LefseCommand", "LefseCommand");
 111                 exit(1);
 112         }
 113 }
 114 //**********************************************************************************************************************
 115 LefseCommand::LefseCommand(string option)  {
 116         try {
 117                 abort = false; calledHelp = false;
 118         allLines = 1;
 119
 120                 //allow user to run help
 121                 if(option == "help") { help(); abort = true; calledHelp = true; }
 122                 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
 123
 124                 else {
 125                         //valid paramters for this command
 126                         vector<string> myArray = setParameters();
 127
 128                         OptionParser parser(option);
 129                         map<string,string> parameters = parser.getParameters();
 130
 131                         ValidParameters validParameter;
 132                         map<string,string>::iterator it;
 133                         //check to make sure all parameters are valid for command
 134                         for (it = parameters.begin(); it != parameters.end(); it++) {
 135                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
 136                         }
 137
 138                         vector<string> tempOutNames;
 139             outputTypes["summary"] = tempOutNames;
 140
 141                         //if the user changes the input directory command factory will send this info to us in the output parameter
 142                         string inputDir = validParameter.validFile(parameters, "inputdir", false);
 143                         if (inputDir == "not found"){   inputDir = "";          }
 144                         else {
 145
 146                 string path;
 147                                 it = parameters.find("design");
 148                                 //user has given a template file
 149                                 if(it != parameters.end()){
 150                                         path = m->hasPath(it->second);
 151                                         //if the user has not given a path then, add inputdir. else leave path alone.
 152                                         if (path == "") {       parameters["desing"] = inputDir + it->second;           }
 153                                 }
 154
 155                 it = parameters.find("shared");
 156                                 //user has given a template file
 157                                 if(it != parameters.end()){
 158                                         path = m->hasPath(it->second);
 159                                         //if the user has not given a path then, add inputdir. else leave path alone.
 160                                         if (path == "") {       parameters["shared"] = inputDir + it->second;           }
 161                                 }
 162             }
 163
 164             //get shared file, it is required
 165                         sharedfile = validParameter.validFile(parameters, "shared", true);
 166                         if (sharedfile == "not open") { sharedfile = ""; abort = true; }
 167                         else if (sharedfile == "not found") {
 168                                 //if there is a current shared file, use it
 169                                 sharedfile = m->getSharedFile();
 170                                 if (sharedfile != "") { m->mothurOut("Using " + sharedfile + " as input file for the shared parameter."); m->mothurOutEndLine(); }
 171                                 else {  m->mothurOut("You have no current sharedfile and the shared parameter is required."); m->mothurOutEndLine(); abort = true; }
 172                         }else { m->setSharedFile(sharedfile); }
 173
 174             //get shared file, it is required
 175                         designfile = validParameter.validFile(parameters, "design", true);
 176                         if (designfile == "not open") { designfile = ""; abort = true; }
 177                         else if (designfile == "not found") {
 178                                 //if there is a current shared file, use it
 179                                 designfile = m->getDesignFile();
 180                                 if (designfile != "") { m->mothurOut("Using " + designfile + " as input file for the design parameter."); m->mothurOutEndLine(); }
 181                                 else {  m->mothurOut("You have no current design file and the design parameter is required."); m->mothurOutEndLine(); abort = true; }
 182                         }else { m->setDesignFile(designfile); }
 183
 184             //if the user changes the output directory command factory will send this info to us in the output parameter
 185                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){
 186                                 outputDir = m->hasPath(sharedfile); //if user entered a file with a path then preserve it
 187                         }
 188
 189             string label = validParameter.validFile(parameters, "label", false);
 190                         if (label == "not found") { label = ""; }
 191                         else {
 192                                 if(label != "all") {  m->splitAtDash(label, labels);  allLines = 0;  }
 193                                 else { allLines = 1;  }
 194                         }
 195
 196             mclass = validParameter.validFile(parameters, "class", false);
 197                         if (mclass == "not found") { mclass = ""; }
 198
 199             subclass = validParameter.validFile(parameters, "subclass", false);
 200                         if (subclass == "not found") { subclass = mclass; }
 201
 202             string temp = validParameter.validFile(parameters, "aalpha", false);
 203                         if (temp == "not found") { temp = "0.05"; }
 204                         m->mothurConvert(temp, anovaAlpha);
 205
 206             temp = validParameter.validFile(parameters, "walpha", false);
 207                         if (temp == "not found") { temp = "0.05"; }
 208                         m->mothurConvert(temp, wilcoxonAlpha);
 209
 210             temp = validParameter.validFile(parameters, "wilc", false);
 211                         if (temp == "not found") { temp = "T"; }
 212                         wilc = m->isTrue(temp);
 213
 214             temp = validParameter.validFile(parameters, "norm", false);
 215                         if (temp == "not found") { temp = "T"; }
 216                         normMillion = m->isTrue(temp);
 217
 218             //temp = validParameter.validFile(parameters, "subject", false);
 219                         //if (temp == "not found") { temp = "F"; }
 220                         //subject = m->isTrue(temp);
 221
 222             temp = validParameter.validFile(parameters, "lda", false);
 223                         if (temp == "not found") { temp = "2.0"; }
 224                         m->mothurConvert(temp, ldaThreshold);
 225
 226             temp = validParameter.validFile(parameters, "iters", false);
 227                         if (temp == "not found") { temp = "30"; }
 228                         m->mothurConvert(temp, iters);
 229
 230             temp = validParameter.validFile(parameters, "fboots", false);
 231                         if (temp == "not found") { temp = "0.67"; }
 232                         m->mothurConvert(temp, fBoots);
 233
 234             //temp = validParameter.validFile(parameters, "wilcsamename", false);
 235                         //if (temp == "not found") { temp = "F"; }
 236                         //wilcsamename = m->isTrue(temp);
 237
 238             temp = validParameter.validFile(parameters, "curv", false);
 239                         if (temp == "not found") { temp = "F"; }
 240                         curv = m->isTrue(temp);
 241
 242             temp = validParameter.validFile(parameters, "strict", false);
 243             if (temp == "not found"){   temp = "0";             }
 244                         if ((temp != "0") && (temp != "1") && (temp != "2")) { m->mothurOut("Invalid strict option: choices are 0, 1 or 2."); m->mothurOutEndLine(); abort=true; }
 245             else {  m->mothurConvert(temp, strict); }
 246
 247             temp = validParameter.validFile(parameters, "minc", false);
 248                         if (temp == "not found") { temp = "10"; }
 249                         m->mothurConvert(temp, minC);
 250
 251             multiClassStrat = validParameter.validFile(parameters, "multiclass", false);
 252             if (multiClassStrat == "not found"){        multiClassStrat = "onevall";            }
 253                         if ((multiClassStrat != "onevall") && (multiClassStrat != "onevone")) { m->mothurOut("Invalid multiclass option: choices are onevone or onevall."); m->mothurOutEndLine(); abort=true; }
 254                 }
 255
 256         }
 257         catch(exception& e) {
 258                 m->errorOut(e, "LefseCommand", "LefseCommand");
 259                 exit(1);
 260         }
 261 }
 262 //**********************************************************************************************************************
 263
 264 int LefseCommand::execute(){
 265         try {
 266         srand(1982);
 267         //for reading lefse formatted file and running in mothur for testing - pass number of rows used for design file
 268         if (false) {  makeShared(1); exit(1); }
 269
 270                 if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
 271
 272         DesignMap designMap(designfile);
 273
 274         //if user did not select class use first column
 275         if (mclass == "") {  mclass = designMap.getDefaultClass(); m->mothurOut("\nYou did not provide a class, using " + mclass +".\n\n"); if (subclass == "") { subclass = mclass; } }
 276
 277         InputData input(sharedfile, "sharedfile");
 278         vector<SharedRAbundFloatVector*> lookup = input.getSharedRAbundFloatVectors();
 279         string lastLabel = lookup[0]->getLabel();
 280
 281         //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
 282         set<string> processedLabels;
 283         set<string> userLabels = labels;
 284
 285
 286         //as long as you are not at the end of the file or done wih the lines you want
 287         while((lookup[0] != NULL) && ((allLines == 1) || (userLabels.size() != 0))) {
 288
 289             if (m->control_pressed) { for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }  return 0; }
 290
 291             if(allLines == 1 || labels.count(lookup[0]->getLabel()) == 1){
 292
 293                 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 294
 295                 process(lookup, designMap);
 296
 297                 processedLabels.insert(lookup[0]->getLabel());
 298                 userLabels.erase(lookup[0]->getLabel());
 299             }
 300
 301             if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
 302                 string saveLabel = lookup[0]->getLabel();
 303
 304                 for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 305                 lookup = input.getSharedRAbundFloatVectors(lastLabel);
 306                 m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 307
 308                 process(lookup, designMap);
 309
 310                 processedLabels.insert(lookup[0]->getLabel());
 311                 userLabels.erase(lookup[0]->getLabel());
 312
 313                 //restore real lastlabel to save below
 314                 lookup[0]->setLabel(saveLabel);
 315             }
 316
 317             lastLabel = lookup[0]->getLabel();
 318             //prevent memory leak
 319             for (int i = 0; i < lookup.size(); i++) {  delete lookup[i]; lookup[i] = NULL; }
 320
 321             if (m->control_pressed) { return 0; }
 322
 323             //get next line to process
 324             lookup = input.getSharedRAbundFloatVectors();
 325         }
 326
 327         if (m->control_pressed) {  return 0; }
 328
 329         //output error messages about any remaining user labels
 330         set<string>::iterator it;
 331         bool needToRun = false;
 332         for (it = userLabels.begin(); it != userLabels.end(); it++) {
 333             m->mothurOut("Your file does not include the label " + *it);
 334             if (processedLabels.count(lastLabel) != 1) {
 335                 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
 336                 needToRun = true;
 337             }else {
 338                 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
 339             }
 340         }
 341
 342         //run last label if you need to
 343         if (needToRun == true)  {
 344             for (int i = 0; i < lookup.size(); i++) { if (lookup[i] != NULL) { delete lookup[i]; } }
 345             lookup = input.getSharedRAbundFloatVectors(lastLabel);
 346
 347             m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine();
 348             process(lookup, designMap);
 349
 350             for (int i = 0; i < lookup.size(); i++) {  delete lookup[i];  }
 351         }
 352
 353
 354         //output files created by command
 355                 m->mothurOutEndLine();
 356                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 357                 for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
 358                 m->mothurOutEndLine();
 359         srand(time(NULL));
 360         return 0;
 361
 362     }
 363         catch(exception& e) {
 364                 m->errorOut(e, "LefseCommand", "execute");
 365                 exit(1);
 366         }
 367 }
 368 //**********************************************************************************************************************
 369
 370 int LefseCommand::process(vector<SharedRAbundFloatVector*>& lookup, DesignMap& designMap) {
 371         try {
 372         vector<string> classes;
 373         vector<string> subclasses;
 374         map<string, string> subclass2Class;
 375         map<string, set<string> > class2SubClasses; //maps class name to vector of its subclasses
 376         map<string, vector<int> > subClass2GroupIndex; //maps subclass name to vector of indexes in lookup from that subclass. old -> 1,2,3 means groups in location 1,2,3 of lookup are from old.  Saves time below.
 377         map<string, vector<int> > class2GroupIndex; //maps subclass name to vector of indexes in lookup from that class. old -> 1,2,3 means groups in location 1,2,3 of lookup are from old.  Saves time below.
 378         if (normMillion) {  normalize(lookup);  }
 379         for (int j = 0; j < lookup.size(); j++) {
 380             string group = lookup[j]->getGroup();
 381             string treatment = designMap.get(group, mclass); //get value for this group in this category
 382             string thisSub = designMap.get(group, subclass);
 383             map<string, string>::iterator it = subclass2Class.find(thisSub);
 384             if (it == subclass2Class.end()) {
 385                 subclass2Class[thisSub] = treatment;
 386                 vector<int> temp; temp.push_back(j);
 387                 subClass2GroupIndex[thisSub] = temp;
 388             }
 389             else {
 390                 if (it->second != treatment) {
 391                     //m->mothurOut("[WARNING]: subclass " + thisSub + " has members in " + it->second + " and " + treatment + ". Subclass members must be from the same class for Wilcoxon. Changing " + thisSub + " to " + treatment + "_" + thisSub + ".\n");
 392                     thisSub = treatment + "_" + thisSub;
 393                     subclass2Class[thisSub] = treatment;
 394                     vector<int> temp; temp.push_back(j);
 395                     subClass2GroupIndex[thisSub] = temp;
 396                 }else { subClass2GroupIndex[thisSub].push_back(j); }
 397             }
 398
 399             map<string, set<string> >::iterator itClass = class2SubClasses.find(treatment);
 400             if (itClass == class2SubClasses.end()) {
 401                 set<string> temp; temp.insert(thisSub);
 402                 class2SubClasses[treatment] = temp;
 403                 vector<int> temp2; temp2.push_back(j);
 404                 class2GroupIndex[treatment] = temp2;
 405                 classes.push_back(treatment);
 406             }else{
 407                 class2SubClasses[treatment].insert(thisSub);
 408                 class2GroupIndex[treatment].push_back(j);
 409             }
 410         }
 411         //sort classes so order is right
 412         sort(classes.begin(), classes.end());
 413
 414         vector< vector<double> > means = getMeans(lookup, class2GroupIndex); //[numOTUs][classes] - classes in same order as class2GroupIndex
 415
 416         //run kruskal wallis on each otu
 417         map<int, double> significantOtuLabels = runKruskalWallis(lookup, designMap);
 418
 419         int numSigBeforeWilcox = significantOtuLabels.size();
 420
 421         if (m->debug) { m->mothurOut("[DEBUG]: completed Kruskal Wallis\n"); }
 422
 423         //check for subclass
 424         string wilcoxString = "";
 425         if ((subclass != "") && wilc) {  significantOtuLabels = runWilcoxon(lookup, designMap, significantOtuLabels, class2SubClasses, subClass2GroupIndex, subclass2Class);  wilcoxString += " ( " + toString(numSigBeforeWilcox) + " ) before internal wilcoxon"; }
 426
 427         int numSigAfterWilcox = significantOtuLabels.size();
 428
 429         if (m->debug) { m->mothurOut("[DEBUG]: completed Wilcoxon\n"); }
 430
 431         m->mothurOut("\nNumber of significantly discriminative features: " + toString(numSigAfterWilcox) + wilcoxString + ".\n");
 432
 433         map<int, double> sigOTUSLDA;
 434         if (numSigAfterWilcox > 0) {
 435             sigOTUSLDA = testLDA(lookup, significantOtuLabels, class2GroupIndex, subClass2GroupIndex);
 436             m->mothurOut("Number of discriminative features with abs LDA score > " + toString(ldaThreshold) + " : " + toString(significantOtuLabels.size()) + ".\n");
 437         }
 438         else { m->mothurOut("No features with significant differences between the classes.\n"); }
 439
 440         if (m->debug) { m->mothurOut("[DEBUG]: completed lda\n"); }
 441
 442         printResults(means, significantOtuLabels, sigOTUSLDA, lookup[0]->getLabel(), classes);
 443
 444         return 0;
 445     }
 446         catch(exception& e) {
 447                 m->errorOut(e, "LefseCommand", "process");
 448                 exit(1);
 449         }
 450 }
 451 //**********************************************************************************************************************
 452 int LefseCommand::normalize(vector<SharedRAbundFloatVector*>& lookup) {
 453         try {
 454         vector<double> mul;
 455         for (int i = 0; i < lookup.size(); i++) {
 456             double sum = 0.0;
 457             for (int j = 0; j < lookup[i]->getNumBins(); j++) { sum += lookup[i]->getAbundance(j); }
 458             mul.push_back(1000000.0/sum);
 459         }
 460
 461         for (int i = 0; i < lookup.size(); i++) {
 462             for (int j = 0; j < lookup[i]->getNumBins(); j++) {
 463                 lookup[i]->set(j, lookup[i]->getAbundance(j)*mul[i], lookup[i]->getGroup());
 464             }
 465         }
 466
 467         return 0;
 468     }
 469         catch(exception& e) {
 470                 m->errorOut(e, "LefseCommand", "normalize");
 471                 exit(1);
 472         }
 473 }
 474 //**********************************************************************************************************************
 475 map<int, double> LefseCommand::runKruskalWallis(vector<SharedRAbundFloatVector*>& lookup, DesignMap& designMap) {
 476         try {
 477         map<int, double> significantOtuLabels;
 478         int numBins = lookup[0]->getNumBins();
 479         //sanity check to make sure each treatment has a group in the shared file
 480         set<string> treatments;
 481         for (int j = 0; j < lookup.size(); j++) {
 482             string group = lookup[j]->getGroup();
 483             string treatment = designMap.get(group, mclass); //get value for this group in this category
 484             treatments.insert(treatment);
 485         }
 486         if (treatments.size() < 2) { m->mothurOut("[ERROR]: need at least 2 things to classes to compare, quitting.\n"); m->control_pressed = true; }
 487
 488         LinearAlgebra linear;
 489         for (int i = 0; i < numBins; i++) {
 490             if (m->control_pressed) { break; }
 491
 492             vector<spearmanRank> values;
 493             for (int j = 0; j < lookup.size(); j++) {
 494                 string group = lookup[j]->getGroup();
 495                 string treatment = designMap.get(group, mclass); //get value for this group in this category
 496                 spearmanRank temp(treatment, lookup[j]->getAbundance(i));
 497                 values.push_back(temp);
 498             }
 499
 500             double pValue = 0.0;
 501             double H = linear.calcKruskalWallis(values, pValue);
 502
 503             if (pValue < anovaAlpha) {  significantOtuLabels[i] = pValue;  }
 504         }
 505
 506         return significantOtuLabels;
 507     }
 508         catch(exception& e) {
 509                 m->errorOut(e, "LefseCommand", "runKruskalWallis");
 510                 exit(1);
 511         }
 512 }
 513 //**********************************************************************************************************************
 514 //assumes not neccessarily paired
 515 map<int, double> LefseCommand::runWilcoxon(vector<SharedRAbundFloatVector*>& lookup, DesignMap& designMap, map<int, double> bins, map<string, set<string> >& class2SubClasses, map<string, vector<int> >& subClass2GroupIndex, map<string, string> subclass2Class) {
 516     try {
 517         map<int, double> significantOtuLabels;
 518         map<int, double>::iterator it;
 519         //if it exists and meets the following requirements run Wilcoxon
 520         /*
 521          1. Subclass members all belong to same main class
 522          anything else
 523         */
 524
 525         int numBins = lookup[0]->getNumBins();
 526         for (int i = 0; i < numBins; i++) {
 527             if (m->control_pressed) { break; }
 528
 529             it = bins.find(i);
 530             if (it != bins.end()) { //flagged in Kruskal Wallis
 531
 532                 vector<float> abunds;  for (int j = 0; j < lookup.size(); j++) { abunds.push_back(lookup[j]->getAbundance(i)); }
 533
 534                 bool sig = testOTUWilcoxon(class2SubClasses, abunds, subClass2GroupIndex, subclass2Class);
 535                 if (sig) { significantOtuLabels[i] = it->second;  }
 536
 537             }//bins flagged from kw
 538         }//for bins
 539
 540         return significantOtuLabels;
 541     }
 542     catch(exception& e) {
 543         m->errorOut(e, "LefseCommand", "runWilcoxon");
 544         exit(1);
 545     }
 546 }
 547 //**********************************************************************************************************************
 548 //lefse.py - test_rep_wilcoxon_r function
 549 bool LefseCommand::testOTUWilcoxon(map<string, set<string> >& class2SubClasses, vector<float> abunds, map<string, vector<int> >& subClass2GroupIndex, map<string, string> subclass2Class) {
 550     try {
 551         int totalOk = 0;
 552         double alphaMtc = wilcoxonAlpha;
 553         vector< set<string> > allDiffs;
 554         LinearAlgebra linear;
 555
 556         //for each subclass comparision
 557         map<string, set<string> >::iterator itB;
 558         for(map<string, set<string> >::iterator it=class2SubClasses.begin();it!=class2SubClasses.end();it++){
 559             itB = it;itB++;
 560             for(itB;itB!=class2SubClasses.end();itB++){
 561                 if (m->control_pressed) { return false; }
 562                 bool first = true;
 563                 int dirCmp = 0; // not set?? dir_cmp = "not_set" # 0=notset or none, 1=true, 2=false.
 564                 int curv_sign = 0;
 565                 int ok = 0;
 566                 int count = 0;
 567                 for (set<string>::iterator itClass1 = (it->second).begin(); itClass1 != (it->second).end(); itClass1++) {
 568                     bool br = false;
 569                     for (set<string>::iterator itClass2 = (itB->second).begin(); itClass2 != (itB->second).end(); itClass2++) {
 570                         string subclass1 = *itClass1;
 571                         string subclass2 = *itClass2;
 572                         count++;
 573
 574                         if (m->debug) { m->mothurOut( "[DEBUG comparing " + it->first + "-" + *itClass1 + " to " + itB->first + "-" + *itClass2 + "\n"); }
 575
 576                         string treatment1 = subclass2Class[subclass1];
 577                         string treatment2 = subclass2Class[subclass2];
 578                         int numSubs1 = class2SubClasses[treatment1].size();
 579                         int numSubs2 = class2SubClasses[treatment2].size();
 580
 581                         //if mul_cor != 0: alpha_mtc = th*l_subcl1*l_subcl2 if mul_cor == 2 else 1.0-math.pow(1.0-th,l_subcl1*l_subcl2)
 582                         if (strict != 0) { alphaMtc = wilcoxonAlpha * numSubs1 * numSubs2 ; }
 583                         if (strict == 2) {}else{ alphaMtc = 1.0-pow((1.0-wilcoxonAlpha),(double)(numSubs1 * numSubs2)); }
 584
 585                         //fill x and y with this comparisons data
 586                         vector<double> x; vector<double> y;
 587
 588                         //fill x and y
 589                         vector<int> xIndexes = subClass2GroupIndex[subclass1]; //indexes in lookup for this subclass
 590                         vector<int> yIndexes = subClass2GroupIndex[subclass2]; //indexes in lookup for this subclass
 591                         for (int k = 0; k < yIndexes.size(); k++) { y.push_back(abunds[yIndexes[k]]);  }
 592                         for (int k = 0; k < xIndexes.size(); k++) { x.push_back(abunds[xIndexes[k]]);  }
 593
 594                         // med_comp = False
 595                         //if len(cl1) < min_c or len(cl2) < min_c:
 596                         //med_comp = True
 597                         bool medComp = false; // are there enough samples per subclass
 598                         if ((xIndexes.size() < minC) || (yIndexes.size() < minC)) { medComp = true; }
 599
 600                         double sx = m->median(x);
 601                         double sy = m->median(y);
 602
 603                         //if cl1[0] == cl2[0] and len(set(cl1)) == 1 and  len(set(cl2)) == 1:
 604                         //tres, first = False, False
 605                         double pValue = 0.0;
 606                         double H = 0.0;
 607                         bool tres = true; //don't think this is set in the python source.  Not sure how that is handled, but setting it here.
 608                         if ((x[0] == y[0]) && (x.size() == 1) && (y.size() == 1)) { tres = false; first = false; }
 609                         else if (!medComp) {
 610                             H = linear.calcWilcoxon(x, y, pValue);
 611                             if (pValue < (alphaMtc*2.0)) { tres = true; }
 612                             else { tres = false; }
 613                         }
 614                         /*if first:
 615                          first = False
 616                          if not curv and ( med_comp or tres ):
 617                          dir_cmp = sx < sy
 618                          if sx == sy: br = True
 619                          elif curv:
 620                          dir_cmp = None
 621                          if med_comp or tres:
 622                          curv_sign += 1
 623                          dir_cmp = sx < sy
 624                          else: br = True
 625                          elif not curv and med_comp:
 626                          if ((sx < sy) != dir_cmp or sx == sy): br = True
 627                          elif curv:
 628                          if tres and dir_cmp == None:
 629                          curv_sign += 1
 630                          dir_cmp = sx < sy
 631                          if tres and dir_cmp != (sx < sy):
 632                          br = True
 633                          curv_sign = -1
 634                          elif not tres or (sx < sy) != dir_cmp or sx == sy: br = True
 635                          */
 636                         int sxSy = 2; //false
 637                         if (sx<sy) {  sxSy = 1; } //true
 638
 639                         if (first) {
 640                             first = false;
 641                             if ((!curv) && (medComp || tres)) {
 642                                 dirCmp = 2; if (sx<sy) { dirCmp = 1; } //dir_cmp = sx < sy
 643                                 if (sx == sy) { br = true; }
 644                             }else if (curv) {
 645                                 dirCmp = 0;
 646                                 if (medComp || tres) {
 647                                     curv_sign++;
 648                                     dirCmp = 2; if (sx<sy) { dirCmp = 1; } //dir_cmp = sx < sy
 649                                 }
 650                             }else { br = true; }
 651                         }else if (!curv && medComp) {
 652                             if (sxSy != dirCmp || sx == sy) { br = true; }
 653                         }else if (curv) {
 654                             if (tres && dirCmp == 0) { curv_sign++; }
 655                             dirCmp = 2; if (sx<sy) { dirCmp = 1; } //dir_cmp = sx < sy
 656                             if (tres && dirCmp != sxSy) { //if tres and dir_cmp != (sx < sy):
 657                                 br = true;
 658                                 curv_sign = -1;
 659                             }
 660                         }else if (!tres || sxSy != dirCmp || sx == sy) { br = true; } //elif not tres or (sx < sy) != dir_cmp or sx == sy: br = True
 661                         if (br) { break; }
 662                         ok++;
 663                     }//for class2 subclasses
 664                     if (br) { break; }
 665                 }//for class1 subclasses
 666                 bool diff = false;
 667                 if (curv) { diff = false; if (curv_sign > 0) { diff = true; } } //if curv: diff = curv_sign > 0
 668                 else { //else: diff = (ok == len(cl_hie[pair[1]])*len(cl_hie[pair[0]]))
 669                     diff = false;
 670                     if (ok == count) { diff = true; }
 671                 }
 672                 if (diff) { totalOk++; }
 673                 if (!diff && (multiClassStrat == "onevone")) { return false; }
 674                 if (diff && (multiClassStrat == "onevall")) { //all_diff.append(pair)
 675                     set<string> pair; pair.insert(it->first); pair.insert(itB->first);
 676                     allDiffs.push_back(pair);
 677                 }
 678             }//classes
 679         }//classes
 680
 681         if (multiClassStrat == "onevall") {
 682             int tot_k = class2SubClasses.size();
 683             for(map<string, set<string> >::iterator it=class2SubClasses.begin();it!=class2SubClasses.end();it++){
 684                 if (m->control_pressed) { return false; }
 685                 int nk = 0;
 686                 //is this class okay in all comparisons
 687                 for (int h = 0; h < allDiffs.size(); h++) {
 688                     if (allDiffs[h].count(it->first) != 0) {  nk++; }
 689                 }
 690                 if (nk == (tot_k-1)) { return true;  }//if nk == tot_k-1: return True
 691             }
 692             return false;
 693         }
 694
 695         return true;
 696     }
 697     catch(exception& e) {
 698         m->errorOut(e, "LefseCommand", "testOTUWilcoxon");
 699         exit(1);
 700     }
 701 }
 702 //**********************************************************************************************************************
 703 //modelled after lefse.py test_lda_r function
 704 map<int, double> LefseCommand::testLDA(vector<SharedRAbundFloatVector*>& lookup, map<int, double> bins, map<string, vector<int> >& class2GroupIndex, map<string, vector<int> >& subClass2GroupIndex) {
 705     try {
 706         map<int, double> sigOTUS;
 707         map<int, double>::iterator it;
 708         LinearAlgebra linear;
 709
 710         int numBins = lookup[0]->getNumBins();
 711         vector< vector<double> > adjustedLookup;
 712
 713         for (int i = 0; i < numBins; i++) {
 714             if (m->control_pressed) { break; }
 715
 716             if (m->debug) { m->mothurOut("[DEBUG]: bin = " + toString(i) + "\n."); }
 717
 718             it = bins.find(i);
 719             if (it != bins.end()) { //flagged in Kruskal Wallis and Wilcoxon(if we ran it)
 720
 721                 if (m->debug) { m->mothurOut("[DEBUG]:flagged bin = " + toString(i) + "\n."); }
 722
 723                 //fill x with this OTUs abundances
 724                 vector<double> x;
 725                 for (int j = 0; j < lookup.size(); j++) {  x.push_back(lookup[j]->getAbundance(i));  }
 726
 727                 //go through classes
 728                 for (map<string, vector<int> >::iterator it = class2GroupIndex.begin(); it != class2GroupIndex.end(); it++) {
 729
 730                     if (m->debug) { m->mothurOut("[DEBUG]: class = " + it->first + "\n."); }
 731
 732                     //max(float(feats['class'].count(c))*0.5,4)
 733                     //max(numGroups in this class*0.5, 4.0)
 734                     double necessaryNum = ((double)((it->second).size())*0.5);
 735                     if (4.0 > necessaryNum) { necessaryNum = 4.0; }
 736
 737                     set<double> uniques;
 738                     for (int j = 0; j < (it->second).size(); j++) { uniques.insert(x[(it->second)[j]]); }
 739
 740                     //if len(set([float(v[1]) for v in ff if v[0] == c])) > max(float(feats['class'].count(c))*0.5,4): continue
 741                     if ((double)(uniques.size()) > necessaryNum) {  }
 742                     else {
 743                         //feats[k][i] = math.fabs(feats[k][i] + lrand.normalvariate(0.0,max(feats[k][i]*0.05,0.01)))
 744                         for (int j = 0; j < (it->second).size(); j++) { //(it->second) contains indexes of abundance for this class
 745                             double sigma = max((x[(it->second)[j]]*0.05), 0.01);
 746                             x[(it->second)[j]] = abs(x[(it->second)[j]] + linear.normalvariate(0.0, sigma));
 747                         }
 748                     }
 749                 }
 750                 adjustedLookup.push_back(x);
 751             }
 752         }
 753
 754         //go through classes
 755         int minCl = 1e6;
 756         map<int, string> indexToClass;
 757         vector<string> classes;
 758         for (map<string, vector<int> >::iterator it = class2GroupIndex.begin(); it != class2GroupIndex.end(); it++) {
 759             //class with minimum number of groups
 760             if ((it->second).size() < minCl) { minCl = (it->second).size(); }
 761             for (int i = 0; i < (it->second).size(); i++) { indexToClass[(it->second)[i]] = it->first; }
 762             classes.push_back(it->first);
 763         }
 764
 765         int numGroups = lookup.size(); //lfk
 766         int fractionNumGroups = numGroups * fBoots; //rfk
 767         minCl = (int)((float)(minCl*fBoots*fBoots*0.05));
 768         minCl = max(minCl, 1);
 769
 770         if (m->debug) { m->mothurOut("[DEBUG]: about to start iters. \n."); }
 771
 772         vector< vector< vector<double> > > results;//[iters][numComparison][numOTUs]
 773         for (int j = 0; j < iters; j++) {
 774             if (m->control_pressed) { return sigOTUS; }
 775
 776             if (m->debug) { m->mothurOut("[DEBUG]: iter = " + toString(j) + "\n."); }
 777
 778             //find "good" random vector
 779             vector<int> rand_s;
 780             for (int h = 0; h < 1000; h++) { //generate a vector of length fractionNumGroups with range 0 to numGroups-1
 781                 rand_s.clear();
 782                 for (int k = 0; k < fractionNumGroups; k++) {  rand_s.push_back(m->getRandomIndex(numGroups-1)); }
 783                 if (!contastWithinClassesOrFewPerClass(adjustedLookup, rand_s, minCl, class2GroupIndex, indexToClass)) { h+=1000; } //break out of loop
 784             }
 785             if (m->control_pressed) { return sigOTUS; }
 786
 787             //print data in R input format for testing
 788             if (false) {
 789                 vector<string> groups; for (int h = 0; h < rand_s.size(); h++) {  groups.push_back(lookup[rand_s[h]]->getGroup()); }
 790                 printToCoutForRTesting(adjustedLookup, rand_s, class2GroupIndex, bins, subClass2GroupIndex, groups);
 791             }
 792
 793             //for each pair of classes
 794             vector< vector<double> > temp = lda(adjustedLookup, rand_s, indexToClass, classes); //[numComparison][numOTUs]
 795             if (temp.size() != 0) { results.push_back(temp); }
 796         }
 797
 798         if (m->control_pressed) { return sigOTUS; }
 799
 800         //m = max([numpy.mean([means[k][kk][p] for kk in range(boots)]) for p in range(len(pairs))])
 801         int k = 0;
 802         for (it = bins.begin(); it != bins.end(); it++) { //[numOTUs] - need to go through bins so we can tie adjustedLookup back to the binNumber. adjustedLookup[0] ->bins entry[0].
 803             vector<double> averageForEachComparison; averageForEachComparison.resize(results[0].size(), 0.0);
 804             double maxM = 0.0; //max of averages for each comparison
 805             for (int j = 0; j < results[0].size(); j++) { //numComparisons
 806                 for (int i = 0; i < results.size(); i++) { //iters
 807                     averageForEachComparison[j]+= results[i][j][k];
 808                 }
 809                 averageForEachComparison[j] /= (double) results.size();
 810                 if (averageForEachComparison[j] > maxM) { maxM = averageForEachComparison[j]; }
 811             }
 812             //res[k] = math.copysign(1.0,m)*math.log(1.0+math.fabs(m),10)
 813             double multiple = 1.0; if (maxM < 0.0) { multiple = -1.0; }
 814             double resK = multiple * log10(1.0+abs(maxM));
 815             if (resK > ldaThreshold) { sigOTUS[it->first] = resK; }
 816             k++;
 817         }
 818
 819         return sigOTUS;
 820     }
 821     catch(exception& e) {
 822         m->errorOut(e, "LefseCommand", "testLDA");
 823         exit(1);
 824     }
 825 }
 826 //**********************************************************************************************************************
 827 vector< vector<double> > LefseCommand::getMeans(vector<SharedRAbundFloatVector*>& lookup, map<string, vector<int> >& class2GroupIndex) {
 828     try {
 829         int numBins = lookup[0]->getNumBins();
 830         int numClasses = class2GroupIndex.size();
 831         vector< vector<double> > means; //[numOTUS][classes]
 832         means.resize(numBins);
 833         for (int i = 0; i < means.size(); i++) {  means[i].resize(numClasses, 0.0); }
 834
 835         map<int, string> indexToClass;
 836         int count = 0;
 837         //shortcut for vectors below
 838         map<string, int> quickIndex;
 839         vector<int> classCounts;
 840         for (map<string, vector<int> >::iterator it = class2GroupIndex.begin(); it != class2GroupIndex.end(); it++) {
 841             for (int i = 0; i < (it->second).size(); i++) { indexToClass[(it->second)[i]] = it->first; }
 842             quickIndex[it->first] = count; count++;
 843             classCounts.push_back((it->second).size());
 844         }
 845
 846         for (int i = 0; i < numBins; i++) {
 847             for (int j = 0; j < lookup.size(); j++) {
 848                 if (m->control_pressed) { return means; }
 849                 means[i][quickIndex[indexToClass[j]]] += lookup[j]->getAbundance(i);
 850             }
 851         }
 852
 853         for (int i = 0; i < numBins; i++) {
 854             for (int j = 0; j < numClasses; j++) { means[i][j] /= (double) classCounts[j];  }
 855         }
 856
 857         return means;
 858     }
 859     catch(exception& e) {
 860         m->errorOut(e, "LefseCommand", "getMeans");
 861         exit(1);
 862     }
 863 }
 864 //**********************************************************************************************************************
 865 vector< vector<double> > LefseCommand::lda(vector< vector<double> >& adjustedLookup, vector<int> rand_s, map<int, string>& indexToClass, vector<string> classes) {
 866     try {
 867         //shortcut for vectors below
 868         map<string, int> quickIndex;
 869         for (int i = 0; i < classes.size(); i++) { quickIndex[classes[i]] = i; }
 870
 871         vector<string> randClass; //classes for rand sample
 872         vector<int> counts; counts.resize(classes.size(), 0);
 873         for (int i = 0; i < rand_s.size(); i++) {
 874             string thisClass = indexToClass[rand_s[i]];
 875             randClass.push_back(thisClass);
 876             counts[quickIndex[thisClass]]++;
 877         }
 878
 879         vector< vector<double> > a; //[numOTUs][numSampled]
 880         for (int i = 0; i < adjustedLookup.size(); i++) {
 881             vector<double> temp;
 882             for (int j = 0; j < rand_s.size(); j++) {
 883                 temp.push_back(adjustedLookup[i][rand_s[j]]);
 884             }
 885             a.push_back(temp);
 886         }
 887
 888         LinearAlgebra linear;
 889         vector< vector<double> > means; bool ignore;
 890         vector< vector<double> > scaling = linear.lda(a, randClass, means, ignore); //means are returned sorted, quickIndex sorts as well since it uses a map. means[class][otu] =
 891         if (ignore) { scaling.clear(); return scaling; }
 892         if (m->control_pressed) { return scaling; }
 893
 894         vector< vector<double> > w; w.resize(a.size()); //w.unit <- w/sqrt(sum(w^2))
 895         double denom = 0.0;
 896         for (int i = 0; i < scaling.size(); i++) { w[i].push_back(scaling[i][0]); denom += (w[i][0]*w[i][0]); }
 897         denom = sqrt(denom);
 898         for (int i = 0; i < w.size(); i++) {  w[i][0] /= denom;  } //[numOTUs][1] - w.unit
 899
 900         //robjects.r('LD <- xy.matrix%*%w.unit') [numSampled][numOtus] * [numOTUs][1]
 901         vector< vector<double> > LD = linear.matrix_mult(linear.transpose(a), w);
 902
 903         //find means for each groups LDs
 904         vector<double> LDMeans; LDMeans.resize(classes.size(), 0.0); //means[0] -> average for [group0].
 905         for (int i = 0; i < LD.size(); i++) {  LDMeans[quickIndex[randClass[i]]] += LD[i][0]; }
 906         for (int i = 0; i < LDMeans.size(); i++) { LDMeans[i] /= (double) counts[i]; }
 907
 908                 //calculate for each comparisons i.e. with groups A,B,C = AB, AC, BC = 3;
 909         vector< vector<double> > results;// [numComparison][numOTUs]
 910                 for (int i = 0; i < LDMeans.size(); i++) {
 911                         for (int l = 0; l < i; l++) {
 912                 if (m->control_pressed) { return scaling; }
 913                 //robjects.r('effect.size <- abs(mean(LD[sub_d[,"class"]=="'+p[0]+'"]) - mean(LD[sub_d[,"class"]=="'+p[1]+'"]))')
 914                 double effectSize = abs(LDMeans[i] - LDMeans[l]);
 915                 //scal = robjects.r('wfinal <- w.unit * effect.size')
 916                 vector<double> compResults;
 917                 for (int j = 0; j < w.size(); j++) { //[numOTUs][1]
 918                     //coeff = [abs(float(v)) if not math.isnan(float(v)) else 0.0 for v in scal]
 919                     double coeff = abs(w[j][0]*effectSize); if (isnan(coeff) || isinf(coeff)) { coeff = 0.0; }
 920                     //gm = abs(res[p[0]][j] - res[p[1]][j]) - res is the means for each group for each otu
 921                     double gm = abs(means[i][j] - means[l][j]);
 922                     //means[k][i].append((gm+coeff[j])*0.5)
 923                     compResults.push_back((gm+coeff)*0.5);
 924                 }
 925                 results.push_back(compResults);
 926             }
 927                 }
 928
 929         return results;
 930     }
 931     catch(exception& e) {
 932         m->errorOut(e, "LefseCommand", "lda");
 933         exit(1);
 934     }
 935 }
 936
 937 //**********************************************************************************************************************
 938 //modelled after lefse.py contast_within_classes_or_few_per_class function
 939 bool LefseCommand::contastWithinClassesOrFewPerClass(vector< vector<double> >& lookup, vector<int> rands, int minCl, map<string, vector<int> > class2GroupIndex, map<int, string> indexToClass) {
 940     try {
 941         set<string> cls;
 942         int countFound = 0;
 943
 944         for (int i = 0; i < rands.size(); i++) { //fill cls with the classes represented in the random selection
 945             for (map<string, vector<int> >::iterator it = class2GroupIndex.begin(); it != class2GroupIndex.end(); it++) {
 946                 if (m->inUsersGroups(rands[i], (it->second))) {
 947                     cls.insert(it->first);
 948                     countFound++;
 949                 }
 950             }
 951         }
 952
 953         //sanity check
 954         if (rands.size() != countFound) { m->mothurOut("oops, should never get here, missing something.\n"); }
 955
 956         if (cls.size() < class2GroupIndex.size()) { return true; } //some classes are not present in sampling
 957
 958         for (set<string>::iterator it = cls.begin(); it != cls.end(); it++) {
 959             if (cls.count(*it) < minCl) { return true; } //this sampling has class count below minimum
 960         }
 961
 962         //for this otu
 963         int numBins = lookup.size();
 964         for (int i = 0; i < numBins; i++) {
 965             if (m->control_pressed) { break; }
 966
 967             //break up random sampling by class
 968             map<string, set<double> > class2Values; //maps class name -> set of abunds present in random sampling. F003Early -> 0.001, 0.003...
 969             for (int j = 0; j < rands.size(); j++) {
 970                 class2Values[indexToClass[rands[j]]].insert(lookup[i][rands[j]]);
 971                 //rands[j] = index of randomly selected group in lookup, randIndex2Class[rands[j]] = class this group belongs to. lookup[rands[j]]->getAbundance(i) = abundance of this group for this OTU.
 972             }
 973             //are the unique values less than we want
 974             //if (len(set(col)) <= min_cl and min_cl > 1) or (min_cl == 1 and len(set(col)) <= 1):
 975             for (map<string, set<double> >::iterator it = class2Values.begin(); it != class2Values.end(); it++) {
 976                 if (((it->second).size() <= minCl && minCl > 1) || (minCl == 1 && (it->second).size() <= 1)) {  return true; }
 977             }
 978         }
 979
 980         return false;
 981     }
 982     catch(exception& e) {
 983         m->errorOut(e, "LefseCommand", "contastWithinClassesOrFewPerClass");
 984         exit(1);
 985     }
 986 }
 987 //**********************************************************************************************************************
 988 int LefseCommand::printResults(vector< vector<double> > means, map<int, double> sigKW, map<int, double> sigLDA, string label, vector<string> classes) {
 989     try {
 990         map<string, string> variables;
 991         variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(sharedfile));
 992         variables["[distance]"] = label;
 993         string outputFileName = getOutputFileName("summary",variables);
 994                 ofstream out;
 995                 m->openOutputFile(outputFileName, out);
 996                 outputNames.push_back(outputFileName); outputTypes["summary"].push_back(outputFileName);
 997
 998         //output headers
 999         out << "OTU\tLogMaxMean\tClass\tLDA\tpValue\n";
1000
1001         string temp = "";
1002         for (int i = 0; i < means.size(); i++) { //[numOTUs][classes]
1003             //find max mean of classes
1004             double maxMean = -1.0; string maxClass = "none";
1005             for (int j = 0; j < means[i].size(); j++) {   if (means[i][j] > maxMean) { maxMean = means[i][j]; maxClass = classes[j]; } }
1006
1007             //str(math.log(max(max(v),1.0),10.0))
1008             double logMaxMean = 1.0;
1009             if (maxMean > logMaxMean) { logMaxMean = maxMean; }
1010             logMaxMean = log10(logMaxMean);
1011
1012             out << m->currentSharedBinLabels[i] << '\t' << logMaxMean << '\t';
1013             if (m->debug) { temp = m->currentSharedBinLabels[i] + '\t' + toString(logMaxMean) + '\t'; }
1014
1015             map<int, double>::iterator it = sigLDA.find(i);
1016             if (it != sigLDA.end()) {
1017                 out << maxClass << '\t' << it->second << '\t' << sigKW[i] << endl; //sigLDA is a subset of sigKW so no need to look
1018                 if (m->debug) { temp += maxClass + '\t' + toString(it->second) + '\t' + toString(sigKW[i]) + '\n'; m->mothurOut(temp); temp = ""; }
1019             }else { out << '-' << endl; }
1020         }
1021
1022         out.close();
1023
1024         return 0;
1025     }
1026     catch(exception& e) {
1027         m->errorOut(e, "LefseCommand", "printResults");
1028         exit(1);
1029     }
1030 }
1031 //**********************************************************************************************************************
1032 //printToCoutForRTesting(adjustedLookup, rand_s, class2GroupIndex, numBins);
1033 bool LefseCommand::printToCoutForRTesting(vector< vector<double> >& adjustedLookup, vector<int> rand_s, map<string, vector<int> >& class2GroupIndex, map<int, double> bins, map<string, vector<int> >& subClass2GroupIndex, vector<string> groups) {
1034     try {
1035         cout << "rand_s = ";
1036         for (int h = 0; h < rand_s.size(); h++) { cout << rand_s[h] << '\t'; } cout << endl;
1037
1038         //print otu data
1039         int count = 0;
1040         for (map<int, double>::iterator it = bins.begin(); it != bins.end(); it++) {
1041             if (m->control_pressed) { break; }
1042
1043             cout << m->currentSharedBinLabels[it->first] << " <- c(";
1044             for (int h = 0; h < rand_s.size()-1; h++) {  cout << (adjustedLookup[count][rand_s[h]]) << ", "; }
1045             cout << (adjustedLookup[count][rand_s[rand_s.size()-1]]) << ")\n";
1046             count++;
1047         }
1048         /*
1049         string tempOutput = "";
1050         for (int h = 0; h < rand_s.size(); h++) {
1051             //find class this index is in
1052             for (map<string, vector<int> >::iterator it = class2GroupIndex.begin(); it!= class2GroupIndex.end(); it++) {
1053                 if (m->inUsersGroups(rand_s[h], (it->second)) ) {   cout << (h+1) << " <- c(\"" +it->first + "\")\n" ; }
1054             }
1055         }*/
1056
1057
1058         string tempOutput = "treatments <- c(";
1059         for (int h = 0; h < rand_s.size(); h++) {
1060             //find class this index is in
1061             for (map<string, vector<int> >::iterator it = class2GroupIndex.begin(); it!= class2GroupIndex.end(); it++) {
1062                 if (m->inUsersGroups(rand_s[h], (it->second)) ) {   tempOutput += "\"" +it->first + "\"" + ","; } //"\"" +it->first + "\""
1063             }
1064         }
1065         tempOutput = tempOutput.substr(0, tempOutput.length()-1);
1066         tempOutput += ")\n";
1067         cout << tempOutput;
1068
1069          /*
1070         if (subclass != "") {
1071             string tempOutput = "sub <- c(";
1072             for (int h = 0; h < rand_s.size(); h++) {
1073                 //find class this index is in
1074                 for (map<string, vector<int> >::iterator it = subClass2GroupIndex.begin(); it!= subClass2GroupIndex.end(); it++) {
1075                     if (m->inUsersGroups(rand_s[h], (it->second)) ) {   tempOutput += "\"" +it->first + "\"" + ','; }
1076                 }
1077             }
1078             tempOutput = tempOutput.substr(0, tempOutput.length()-1);
1079             tempOutput += ")\n";
1080             cout << tempOutput;
1081         }
1082
1083         if (subject) {
1084             string tempOutput = "group <- c(";
1085             for (int h = 0; h < groups.size(); h++) {
1086                 tempOutput += "\"" +groups[h] + "\"" + ',';
1087             }
1088             tempOutput = tempOutput.substr(0, tempOutput.length()-1);
1089             tempOutput += ")\n";
1090             cout << tempOutput;
1091         }*/
1092
1093
1094         //print data frame
1095         tempOutput = "dat <- data.frame(";
1096         for (map<int, double>::iterator it = bins.begin(); it != bins.end(); it++) {
1097             if (m->control_pressed) { break; }
1098
1099             tempOutput += "\"" + m->currentSharedBinLabels[it->first] + "\"=" + m->currentSharedBinLabels[it->first] + ",";
1100         }
1101         //tempOutput = tempOutput.substr(0, tempOutput.length()-1);
1102         tempOutput += " class=treatments";
1103         //if (subclass != "") { tempOutput += ", subclass=sub"; }
1104         //if (subject) { tempOutput += ", subject=group"; }
1105         tempOutput += ")\n";
1106         cout << tempOutput;
1107
1108         tempOutput = "z <- suppressWarnings(mylda(as.formula(class ~ ";
1109         for (map<int, double>::iterator it = bins.begin(); it != bins.end(); it++) {
1110             if (m->control_pressed) { break; }
1111
1112             tempOutput +=  m->currentSharedBinLabels[it->first] + "+";
1113         }
1114         tempOutput = tempOutput.substr(0, tempOutput.length()-1); //rip off extra plus sign
1115         tempOutput += "), data = dat, tol = 1e-10))";
1116         cout << tempOutput + "\nz\n";
1117         cout << "w <- z$scaling[,1]\n"; //robjects.r('w <- z$scaling[,1]')
1118         cout << "w.unit <- w/sqrt(sum(w^2))\n"; //robjects.r('w.unit <- w/sqrt(sum(w^2))')
1119         cout << "ss <- dat[,-match(\"class\",colnames(dat))]\n"; //robjects.r('ss <- sub_d[,-match("class",colnames(sub_d))]')
1120         //if (subclass != "") { cout << "ss <- ss[,-match(\"subclass\",colnames(ss))]\n";  }//robjects.r('ss <- ss[,-match("subclass",colnames(ss))]')
1121         //if (subject) { cout << "ss <- ss[,-match(\"subject\",colnames(ss))]\n";  }//robjects.r('ss <- ss[,-match("subject",colnames(ss))]')
1122         cout << "xy.matrix <- as.matrix(ss)\n"; //robjects.r('xy.matrix <- as.matrix(ss)')
1123         cout << "LD <- xy.matrix%*%w.unit\n"; //robjects.r('LD <- xy.matrix%*%w.unit')
1124         cout << "effect.size <- abs(mean(LD[dat[,\"class\"]==\"'+p[0]+'\"]) - mean(LD[dat[,\"class\"]==\"'+p[1]+'\"]))\n"; //robjects.r('effect.size <- abs(mean(LD[sub_d[,"class"]=="'+p[0]+'"]) - mean(LD[sub_d[,"class"]=="'+p[1]+'"]))')
1125         cout << "wfinal <- w.unit * effect.size\n"; //scal = robjects.r('wfinal <- w.unit * effect.size')
1126         cout << "mm <- z$means\n"; //rres = robjects.r('mm <- z$means')
1127
1128         return true;
1129     }
1130     catch(exception& e) {
1131         m->errorOut(e, "LefseCommand", "printToCoutForRTesting");
1132         exit(1);
1133     }
1134 }
1135 //**********************************************************************************************************************
1136 int LefseCommand::makeShared(int numDesignLines) {
1137     try {
1138         ifstream in;
1139         m->openInputFile(sharedfile, in);
1140         vector< vector<string> > lines;
1141         for(int i = 0; i < numDesignLines; i++) {
1142             if (m->control_pressed) { return 0; }
1143
1144             string line = m->getline(in);
1145             cout << line << endl;
1146             vector<string> pieces = m->splitWhiteSpace(line);
1147             lines.push_back(pieces);
1148         }
1149
1150         ofstream out;
1151         m->openOutputFile(sharedfile+".design", out); out << "group" << '\t';
1152         for (int j = 0; j < lines.size(); j++) { out << lines[j][0] << '\t'; } out << endl;
1153         for (int j = 1; j < lines[0].size(); j++) {
1154             out <<(j-1) << '\t';
1155             for (int i = 0; i < lines.size(); i++) {
1156                  out << lines[i][j] << '\t';
1157             }
1158             out << endl;
1159         }
1160         out.close();
1161         DesignMap design(sharedfile+".design");
1162
1163         vector<SharedRAbundFloatVector*> lookup;
1164         for (int k = 0; k < lines[0].size()-1; k++) {
1165             SharedRAbundFloatVector* temp = new SharedRAbundFloatVector();
1166             temp->setLabel("0.03");
1167             temp->setGroup(toString(k));
1168             lookup.push_back(temp);
1169         }
1170
1171         m->currentSharedBinLabels.clear();
1172         int count = 0;
1173         while (!in.eof()) {
1174             if (m->control_pressed) { return 0; }
1175
1176             string line = m->getline(in);
1177             vector<string> pieces = m->splitWhiteSpace(line);
1178
1179             float sum = 0.0;
1180             for (int i = 1; i < pieces.size(); i++) {
1181                 float value; m->mothurConvert(pieces[i], value);
1182                 sum += value;
1183             }
1184
1185             if (sum != 0.0) {
1186                 //cout << count << '\t';
1187                 for (int i = 1; i < pieces.size(); i++) {
1188                     float value; m->mothurConvert(pieces[i], value);
1189                     lookup[i-1]->push_back(value, toString(i-1));
1190                     //cout << pieces[i] << '\t';
1191                 }
1192                 m->currentSharedBinLabels.push_back(toString(count));
1193                 //m->currentBinLabels.push_back(pieces[0]);
1194                 //cout << line<< endl;
1195                 //cout << endl;
1196             }
1197             count++;
1198         }
1199         in.close();
1200
1201         for (int k = 0; k < lookup.size(); k++) {
1202             //cout << "0.03" << '\t' << toString(k) << endl; lookup[k]->print(cout);
1203         }
1204
1205         process(lookup, design);
1206
1207         return 0;
1208     }
1209     catch(exception& e) {
1210         m->errorOut(e, "LefseCommand", "printToCoutForRTesting");
1211         exit(1);
1212     }
1213 }
1214
1215 //**********************************************************************************************************************
1216
1217