mantelcommand.cpp

   1 /*
   2  *  mantelcommand.cpp
   3  *  mothur
   4  *
   5  *  Created by westcott on 2/9/11.
   6  *  Copyright 2011 Schloss Lab. All rights reserved.
   7  *
   8  */
   9
  10 #include "mantelcommand.h"
  11 #include "readphylipvector.h"
  12
  13 //**********************************************************************************************************************
  14 vector<string> MantelCommand::getValidParameters(){
  15         try {
  16                 string Array[] =  {"phylip1","phylip2","method","iters","outputdir","inputdir"};
  17                 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
  18                 return myArray;
  19         }
  20         catch(exception& e) {
  21                 m->errorOut(e, "MantelCommand", "getValidParameters");
  22                 exit(1);
  23         }
  24 }
  25 //**********************************************************************************************************************
  26 vector<string> MantelCommand::getRequiredParameters(){
  27         try {
  28                 string Array[] =  {"phylip1", "phylip2"};
  29                 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
  30                 return myArray;
  31         }
  32         catch(exception& e) {
  33                 m->errorOut(e, "MantelCommand", "getRequiredParameters");
  34                 exit(1);
  35         }
  36 }
  37 //**********************************************************************************************************************
  38 MantelCommand::MantelCommand(){
  39         try {
  40                 abort = true; calledHelp = true;
  41                 vector<string> tempOutNames;
  42                 outputTypes["mantel"] = tempOutNames;
  43         }
  44         catch(exception& e) {
  45                 m->errorOut(e, "MantelCommand", "MantelCommand");
  46                 exit(1);
  47         }
  48 }
  49
  50 //**********************************************************************************************************************
  51 vector<string> MantelCommand::getRequiredFiles(){
  52         try {
  53                 vector<string> myArray;
  54                 return myArray;
  55         }
  56         catch(exception& e) {
  57                 m->errorOut(e, "MantelCommand", "getRequiredFiles");
  58                 exit(1);
  59         }
  60 }
  61 //**********************************************************************************************************************
  62 MantelCommand::MantelCommand(string option)  {
  63         try {
  64                 abort = false; calledHelp = false;
  65
  66                 //allow user to run help
  67                 if(option == "help") { help(); abort = true; calledHelp = true; }
  68
  69                 else {
  70                         //valid paramters for this command
  71                         string Array[] =  {"phylip1","phylip2","method","iters","outputdir","inputdir"};
  72                         vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
  73
  74                         OptionParser parser(option);
  75                         map<string, string> parameters = parser.getParameters();
  76
  77                         ValidParameters validParameter;
  78                         map<string, string>::iterator it;
  79
  80                         //check to make sure all parameters are valid for command
  81                         for (it = parameters.begin(); it != parameters.end(); it++) {
  82                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
  83                         }
  84
  85                         vector<string> tempOutNames;
  86                         outputTypes["mantel"] = tempOutNames;
  87
  88                         //if the user changes the input directory command factory will send this info to us in the output parameter
  89                         string inputDir = validParameter.validFile(parameters, "inputdir", false);
  90                         if (inputDir == "not found"){   inputDir = "";          }
  91                         else {
  92                                 string path;
  93                                 it = parameters.find("phylip1");
  94                                 //user has given a template file
  95                                 if(it != parameters.end()){
  96                                         path = m->hasPath(it->second);
  97                                         //if the user has not given a path then, add inputdir. else leave path alone.
  98                                         if (path == "") {       parameters["phylip1"] = inputDir + it->second;          }
  99                                 }
 100
 101                                 it = parameters.find("phylip2");
 102                                 //user has given a template file
 103                                 if(it != parameters.end()){
 104                                         path = m->hasPath(it->second);
 105                                         //if the user has not given a path then, add inputdir. else leave path alone.
 106                                         if (path == "") {       parameters["phylip2"] = inputDir + it->second;          }
 107                                 }
 108                         }
 109
 110
 111                         //check for required parameters
 112                         phylipfile1 = validParameter.validFile(parameters, "phylip1", true);
 113                         if (phylipfile1 == "not open") { phylipfile1 = ""; abort = true; }
 114                         else if (phylipfile1 == "not found") { phylipfile1 = ""; m->mothurOut("phylip1 is a required parameter for the mantel command."); m->mothurOutEndLine(); abort = true;  }
 115
 116                         phylipfile2 = validParameter.validFile(parameters, "phylip2", true);
 117                         if (phylipfile2 == "not open") { phylipfile2 = ""; abort = true; }
 118                         else if (phylipfile2 == "not found") { phylipfile2 = ""; m->mothurOut("phylip2 is a required parameter for the mantel command."); m->mothurOutEndLine(); abort = true;  }
 119
 120                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){  outputDir = m->hasPath(phylipfile1);    }
 121
 122                         method = validParameter.validFile(parameters, "method", false);         if (method == "not found"){     method = "pearson";             }
 123
 124                         string temp = validParameter.validFile(parameters, "iters", false);                     if (temp == "not found") { temp = "1000"; }
 125                         convert(temp, iters);
 126
 127                         if ((method != "pearson") && (method != "spearman") && (method != "kendall")) { m->mothurOut(method + " is not a valid method. Valid methods are pearson, spearman, and kendall."); m->mothurOutEndLine(); abort = true; }
 128                 }
 129         }
 130         catch(exception& e) {
 131                 m->errorOut(e, "MantelCommand", "MantelCommand");
 132                 exit(1);
 133         }
 134 }
 135 //**********************************************************************************************************************
 136
 137 void MantelCommand::help(){
 138         try {
 139                 m->mothurOut("The mantel command reads two distance matrices and calculates the mantel correlation coefficient.\n");
 140                 m->mothurOut("The mantel command parameters are phylip1, phylip2 and method.  The phylip1 and phylip2 parameters are required.  Matrices must be the same size and contain the same names.\n");
 141                 m->mothurOut("The method parameter allows you to select what method you would like to use. Options are pearson, spearman and kendall. Default=pearson.\n");
 142                 m->mothurOut("The mantel command should be in the following format: mantel(phylip1=veg.dist, phylip2=env.dist).\n");
 143                 m->mothurOut("The mantel command outputs a .mantel file.\n");
 144                 m->mothurOut("Note: No spaces between parameter labels (i.e. phylip1), '=' and parameters (i.e. veg.dist).\n\n");
 145         }
 146         catch(exception& e) {
 147                 m->errorOut(e, "MantelCommand", "help");
 148                 exit(1);
 149         }
 150 }
 151
 152 //**********************************************************************************************************************
 153
 154 MantelCommand::~MantelCommand(){}
 155
 156 //**********************************************************************************************************************
 157
 158 int MantelCommand::execute(){
 159         try {
 160
 161                 if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
 162
 163                 /***************************************************/
 164                 //    reading distance files                                       //
 165                 /***************************************************/
 166
 167                 //read phylip1
 168                 ReadPhylipVector readMatrix(phylipfile1);
 169                 vector< vector<double> > matrix1;
 170                 vector<string> names1 = readMatrix.read(matrix1);
 171
 172                 if (m->control_pressed) { return 0; }
 173
 174                 //read phylip2
 175                 ReadPhylipVector readMatrix2(phylipfile2);
 176                 vector<seqDist> temp; //seqDist - int, int, float
 177                 vector<string> names2 = readMatrix2.read(temp);
 178
 179                 if (m->control_pressed) { return 0; }
 180
 181                 //fill matrix2 making sure to make sure the distances are in the same order as matrix1
 182                 vector< vector<double> > matrix2;
 183                 if (names1 == names2) { //then everything is in same order and same size
 184
 185                         //initialize space
 186                         matrix2.resize(names2.size());
 187                         for (int i = 0; i < matrix2.size(); i++) { matrix2[i].resize(names2.size(), 0.0); }
 188
 189                         //fill matrix2
 190                         for (int i = 0; i < temp.size(); i++) {
 191                                 matrix2[temp[i].seq1][temp[i].seq2] = temp[i].dist;
 192                                 matrix2[temp[i].seq2][temp[i].seq1] = temp[i].dist;
 193                         }
 194
 195                 }else if (names1.size() != names2.size()) { //wrong size no need to order, abort
 196                         m->mothurOut("[ERROR]: distance matrices are not the same size, aborting."); m->mothurOutEndLine();
 197                         m->control_pressed = true;
 198                 }else { //sizes are the same, but either the names are different or they are in different order
 199
 200                         //map location of name in names1 to location of name in names2
 201                         map<string, int> names1Map;
 202                         map<string, int>::iterator it;
 203                         for (int i = 0; i < names1.size(); i++) {  names1Map[names1[i]] = i; }
 204
 205                         map<string, int> names2Map;
 206                         bool nameError = false;
 207                         for (int i = 0; i < names2.size(); i++) {
 208
 209                                 //if you find one name error stop looking
 210                                 if (!nameError) {
 211                                         it = names1Map.find(names2[i]);
 212                                         if (it == names1Map.end()) { nameError = true; }
 213                                 }
 214
 215                                 //are names different
 216                                 names2Map[names2[i]] = i;
 217                         }
 218
 219                         //initialize space
 220                         matrix2.resize(names2.size());
 221                         for (int i = 0; i < matrix2.size(); i++) { matrix2[i].resize(names2.size(), 0.0); }
 222
 223                         //fill matrix2
 224                         //are we comparing apples to apples?
 225                         if (nameError) {
 226                                 m->mothurOut("[WARNING]: Names do not match between distance files. Comparing based on order in files."); m->mothurOutEndLine();
 227
 228                                 for (int i = 0; i < temp.size(); i++) {
 229                                         matrix2[temp[i].seq1][temp[i].seq2] = temp[i].dist;
 230                                         matrix2[temp[i].seq2][temp[i].seq1] = temp[i].dist;
 231                                 }
 232
 233                         }else { //no name error just different orders, so reorder
 234
 235                                 for (int i = 0; i < temp.size(); i++) {
 236
 237                                         //what's the location of this distance comparison in matrix1
 238                                         string matrix2NameI = names2[temp[i].seq1];
 239                                         string matrix2NameJ = names2[temp[i].seq2];
 240                                         int locationI = names1Map[matrix2NameI];
 241                                         int locationJ = names1Map[matrix2NameJ];
 242
 243                                         matrix2[locationI][locationJ] = temp[i].dist;
 244                                         matrix2[locationJ][locationI] = temp[i].dist;
 245                                 }
 246                         }
 247
 248                 }
 249
 250                 //frees up space
 251                 temp.clear();
 252
 253                 if (m->control_pressed) { return 0; }
 254
 255                 /***************************************************/
 256                 //    calculating mantel and signifigance                  //
 257                 /***************************************************/
 258
 259                 //calc mantel coefficient
 260                 LinearAlgebra linear;
 261                 double mantel = 0.0;
 262                 if (method == "pearson")                {  mantel = linear.calcPearson(matrix1, matrix2);       }
 263                 else if (method == "spearman")  {  mantel = linear.calcSpearman(matrix1, matrix2);      }
 264                 else if (method == "kendall")   {  mantel = linear.calcKendall(matrix1, matrix2);       }
 265
 266
 267                 //calc signifigance
 268                 int count = 0;
 269                 for (int i = 0; i < iters; i++) {
 270
 271                         if (m->control_pressed) { return 0; }
 272
 273                         //randomize matrix2
 274                         vector< vector<double> > matrix2Copy = matrix2;
 275                         random_shuffle(matrix2Copy.begin(), matrix2Copy.end());
 276
 277                         //calc random mantel
 278                         double randomMantel = 0.0;
 279                         if (method == "pearson")                {  randomMantel = linear.calcPearson(matrix1, matrix2Copy);     }
 280                         else if (method == "spearman")  {  randomMantel = linear.calcSpearman(matrix1, matrix2Copy);    }
 281                         else if (method == "kendall")   {  randomMantel = linear.calcKendall(matrix1, matrix2Copy);     }
 282
 283                         if (randomMantel >= mantel) { count++; }
 284                 }
 285
 286                 double pValue = count / (float) iters;
 287
 288                 if (m->control_pressed) { return 0; }
 289
 290                 string outputFile = outputDir + m->getRootName(m->getSimpleName(phylipfile1)) + "mantel";
 291                 outputNames.push_back(outputFile); outputTypes["mantel"].push_back(outputFile);
 292                 ofstream out;
 293
 294                 m->openOutputFile(outputFile, out);
 295
 296                 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
 297                 cout.setf(ios::fixed, ios::floatfield); cout.setf(ios::showpoint);
 298
 299                 out << "Mantel\tpValue" << endl;
 300                 out << mantel << '\t' << pValue << endl;
 301
 302                 out.close();
 303
 304                 cout << "\nmantel = " << mantel << "\tpValue = " << pValue << endl;
 305                 m->mothurOutJustToLog("\nmantel = " + toString(mantel) + "\tpValue = " + toString(pValue) + "\n");
 306
 307                 m->mothurOutEndLine();
 308                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 309                 for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
 310                 m->mothurOutEndLine();
 311
 312                 return 0;
 313         }
 314         catch(exception& e) {
 315                 m->errorOut(e, "MantelCommand", "execute");
 316                 exit(1);
 317         }
 318 }
 319
 320 //**********************************************************************************************************************
 321
 322