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("Sokal, R. R., & Rohlf, F. J. (1995). Biometry, 3rd edn. New York: Freeman.\n");
 140                 m->mothurOut("The mantel command reads two distance matrices and calculates the mantel correlation coefficient.\n");
 141                 m->mothurOut("The mantel command parameters are phylip1, phylip2, iters and method.  The phylip1 and phylip2 parameters are required.  Matrices must be the same size and contain the same names.\n");
 142                 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");
 143                 m->mothurOut("The iters parameter allows you to set number of randomization for the P value.  The default is 1000. \n");
 144                 m->mothurOut("The mantel command should be in the following format: mantel(phylip1=veg.dist, phylip2=env.dist).\n");
 145                 m->mothurOut("The mantel command outputs a .mantel file.\n");
 146                 m->mothurOut("Note: No spaces between parameter labels (i.e. phylip1), '=' and parameters (i.e. veg.dist).\n\n");
 147         }
 148         catch(exception& e) {
 149                 m->errorOut(e, "MantelCommand", "help");
 150                 exit(1);
 151         }
 152 }
 153
 154 //**********************************************************************************************************************
 155
 156 MantelCommand::~MantelCommand(){}
 157
 158 //**********************************************************************************************************************
 159
 160 int MantelCommand::execute(){
 161         try {
 162
 163                 if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
 164
 165                 /***************************************************/
 166                 //    reading distance files                                       //
 167                 /***************************************************/
 168
 169                 //read phylip1
 170                 ReadPhylipVector readMatrix(phylipfile1);
 171                 vector< vector<double> > matrix1;
 172                 vector<string> names1 = readMatrix.read(matrix1);
 173
 174                 if (m->control_pressed) { return 0; }
 175
 176                 //read phylip2
 177                 ReadPhylipVector readMatrix2(phylipfile2);
 178                 vector< vector<double> > matrix2;
 179                 vector<string> names2 = readMatrix2.read(matrix2);
 180
 181                 if (m->control_pressed) { return 0; }
 182
 183                 //make sure matrix2 and matrix1 are in the same order
 184                 if (names1 == names2) { //then everything is in same order and same size
 185                 }else if (names1.size() != names2.size()) { //wrong size no need to order, abort
 186                         m->mothurOut("[ERROR]: distance matrices are not the same size, aborting."); m->mothurOutEndLine();
 187                         m->control_pressed = true;
 188                 }else { //sizes are the same, but either the names are different or they are in different order
 189                         m->mothurOut("[WARNING]: Names do not match between distance files. Comparing based on order in files."); m->mothurOutEndLine();
 190                 }
 191
 192                 if (m->control_pressed) { return 0; }
 193
 194                 /***************************************************/
 195                 //    calculating mantel and signifigance                  //
 196                 /***************************************************/
 197
 198                 //calc mantel coefficient
 199                 LinearAlgebra linear;
 200                 double mantel = 0.0;
 201                 if (method == "pearson")                {  mantel = linear.calcPearson(matrix1, matrix2);       }
 202                 else if (method == "spearman")  {  mantel = linear.calcSpearman(matrix1, matrix2);      }
 203                 else if (method == "kendall")   {  mantel = linear.calcKendall(matrix1, matrix2);       }
 204
 205
 206                 //calc signifigance
 207                 int count = 0;
 208                 for (int i = 0; i < iters; i++) {
 209
 210                         if (m->control_pressed) { return 0; }
 211
 212                         //randomize matrix2
 213                         vector< vector<double> > matrix2Copy = matrix2;
 214                         random_shuffle(matrix2Copy.begin(), matrix2Copy.end());
 215
 216                         //calc random mantel
 217                         double randomMantel = 0.0;
 218                         if (method == "pearson")                {  randomMantel = linear.calcPearson(matrix1, matrix2Copy);     }
 219                         else if (method == "spearman")  {  randomMantel = linear.calcSpearman(matrix1, matrix2Copy);    }
 220                         else if (method == "kendall")   {  randomMantel = linear.calcKendall(matrix1, matrix2Copy);     }
 221
 222                         if (randomMantel >= mantel) { count++; }
 223                 }
 224
 225                 double pValue = count / (float) iters;
 226
 227                 if (m->control_pressed) { return 0; }
 228
 229                 string outputFile = outputDir + m->getRootName(m->getSimpleName(phylipfile1)) + "mantel";
 230                 outputNames.push_back(outputFile); outputTypes["mantel"].push_back(outputFile);
 231                 ofstream out;
 232
 233                 m->openOutputFile(outputFile, out);
 234
 235                 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
 236                 cout.setf(ios::fixed, ios::floatfield); cout.setf(ios::showpoint);
 237
 238                 out << "Mantel\tpValue" << endl;
 239                 out << mantel << '\t' << pValue << endl;
 240
 241                 out.close();
 242
 243                 cout << "\nmantel = " << mantel << "\tpValue = " << pValue << endl;
 244                 m->mothurOutJustToLog("\nmantel = " + toString(mantel) + "\tpValue = " + toString(pValue) + "\n");
 245
 246                 m->mothurOutEndLine();
 247                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 248                 for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
 249                 m->mothurOutEndLine();
 250
 251                 return 0;
 252         }
 253         catch(exception& e) {
 254                 m->errorOut(e, "MantelCommand", "execute");
 255                 exit(1);
 256         }
 257 }
 258
 259 //**********************************************************************************************************************
 260
 261