pcoacommand.cpp

   1
   2 /*
   3  *  pcacommand.cpp
   4  *  Mothur
   5  *
   6  *  Created by westcott on 1/4/10.
   7  *  Copyright 2010 Schloss Lab. All rights reserved.
   8  *
   9  */
  10
  11 #include "pcoacommand.h"
  12 #include "readphylipvector.h"
  13
  14 //**********************************************************************************************************************
  15 vector<string> PCOACommand::getValidParameters(){
  16         try {
  17                 string Array[] =  {"phylip", "metric","outputdir","inputdir"};
  18                 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
  19                 return myArray;
  20         }
  21         catch(exception& e) {
  22                 m->errorOut(e, "PCOACommand", "getValidParameters");
  23                 exit(1);
  24         }
  25 }
  26 //**********************************************************************************************************************
  27 PCOACommand::PCOACommand(){
  28         try {
  29                 abort = true;
  30                 //initialize outputTypes
  31                 vector<string> tempOutNames;
  32                 outputTypes["pcoa"] = tempOutNames;
  33                 outputTypes["loadings"] = tempOutNames;
  34         }
  35         catch(exception& e) {
  36                 m->errorOut(e, "PCOACommand", "PCOACommand");
  37                 exit(1);
  38         }
  39 }
  40 //**********************************************************************************************************************
  41 vector<string> PCOACommand::getRequiredParameters(){
  42         try {
  43                 string Array[] =  {"phylip"};
  44                 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
  45                 return myArray;
  46         }
  47         catch(exception& e) {
  48                 m->errorOut(e, "PCOACommand", "getRequiredParameters");
  49                 exit(1);
  50         }
  51 }
  52 //**********************************************************************************************************************
  53 vector<string> PCOACommand::getRequiredFiles(){
  54         try {
  55                 vector<string> myArray;
  56                 return myArray;
  57         }
  58         catch(exception& e) {
  59                 m->errorOut(e, "PCOACommand", "getRequiredFiles");
  60                 exit(1);
  61         }
  62 }
  63 //**********************************************************************************************************************
  64
  65 PCOACommand::PCOACommand(string option)  {
  66         try {
  67                 abort = false;
  68
  69                 //allow user to run help
  70                 if(option == "help") { help(); abort = true; }
  71
  72                 else {
  73                         //valid paramters for this command
  74                         string Array[] =  {"phylip","metric","outputdir", "inputdir"};
  75                         vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
  76
  77                         OptionParser parser(option);
  78                         map<string, string> parameters = parser. getParameters();
  79
  80                         ValidParameters validParameter;
  81                         map<string, string>::iterator it;
  82
  83                         //check to make sure all parameters are valid for command
  84                         for (it = parameters.begin(); it != parameters.end(); it++) {
  85                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
  86                         }
  87                         //if the user changes the input directory command factory will send this info to us in the output parameter
  88                         string inputDir = validParameter.validFile(parameters, "inputdir", false);
  89                         if (inputDir == "not found"){   inputDir = "";          }
  90                         else {
  91                                 string path;
  92                                 it = parameters.find("phylip");
  93                                 //user has given a template file
  94                                 if(it != parameters.end()){
  95                                         path = m->hasPath(it->second);
  96                                         //if the user has not given a path then, add inputdir. else leave path alone.
  97                                         if (path == "") {       parameters["phylip"] = inputDir + it->second;           }
  98                                 }
  99                         }
 100
 101                         //initialize outputTypes
 102                         vector<string> tempOutNames;
 103                         outputTypes["pcoa"] = tempOutNames;
 104                         outputTypes["loadings"] = tempOutNames;
 105
 106                         //required parameters
 107                         phylipfile = validParameter.validFile(parameters, "phylip", true);
 108                         if (phylipfile == "not open") { abort = true; }
 109                         else if (phylipfile == "not found") { phylipfile = ""; abort = true; }
 110                         else {  filename = phylipfile;  }
 111
 112                         //if the user changes the output directory command factory will send this info to us in the output parameter
 113                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){
 114                                 outputDir = "";
 115                                 outputDir += m->hasPath(phylipfile); //if user entered a file with a path then preserve it
 116                         }
 117
 118                         //error checking on files
 119                         if (phylipfile == "")   { m->mothurOut("You must provide a distance file before running the pcoa command."); m->mothurOutEndLine(); abort = true; }
 120
 121                         string temp = validParameter.validFile(parameters, "metric", false);    if (temp == "not found"){       temp = "T";                             }
 122                         metric = m->isTrue(temp);
 123                 }
 124
 125         }
 126         catch(exception& e) {
 127                 m->errorOut(e, "PCOACommand", "PCOACommand");
 128                 exit(1);
 129         }
 130 }
 131 //**********************************************************************************************************************
 132 void PCOACommand::help(){
 133         try {
 134
 135                 m->mothurOut("The pcoa command parameters are phylip and metric"); m->mothurOutEndLine();
 136                 m->mothurOut("The phylip parameter allows you to enter your distance file."); m->mothurOutEndLine();
 137                 m->mothurOut("The metric parameter allows indicate you if would like the pearson correlation coefficient calculated. Default=True"); m->mothurOutEndLine();
 138                 m->mothurOut("Example pcoa(phylip=yourDistanceFile).\n");
 139                 m->mothurOut("Note: No spaces between parameter labels (i.e. phylip), '=' and parameters (i.e.yourDistanceFile).\n\n");
 140         }
 141         catch(exception& e) {
 142                 m->errorOut(e, "PCOACommand", "help");
 143                 exit(1);
 144         }
 145 }
 146 //**********************************************************************************************************************
 147 PCOACommand::~PCOACommand(){}
 148 //**********************************************************************************************************************
 149 int PCOACommand::execute(){
 150         try {
 151
 152                 if (abort == true) { return 0; }
 153
 154                 cout.setf(ios::fixed, ios::floatfield);
 155                 cout.setf(ios::showpoint);
 156                 cerr.setf(ios::fixed, ios::floatfield);
 157                 cerr.setf(ios::showpoint);
 158
 159                 vector<string> names;
 160                 vector<vector<double> > D;
 161
 162                 fbase = outputDir + m->getRootName(m->getSimpleName(filename));
 163
 164                 ReadPhylipVector readFile(filename);
 165                 names = readFile.read(D);
 166
 167                 if (m->control_pressed) { return 0; }
 168
 169                 double offset = 0.0000;
 170                 vector<double> d;
 171                 vector<double> e;
 172                 vector<vector<double> > G = D;
 173                 vector<vector<double> > copy_G;
 174
 175                 m->mothurOut("\nProcessing...\n\n");
 176
 177                 for(int count=0;count<2;count++){
 178                         recenter(offset, D, G);                                 if (m->control_pressed) { return 0; }
 179                         linearCalc.tred2(G, d, e);                              if (m->control_pressed) { return 0; }
 180                         linearCalc.qtli(d, e, G);                               if (m->control_pressed) { return 0; }
 181                         offset = d[d.size()-1];
 182                         if(offset > 0.0) break;
 183                 }
 184
 185                 if (m->control_pressed) { return 0; }
 186
 187                 output(fbase, names, G, d);
 188
 189                 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) {        remove(outputNames[i].c_str());  } return 0; }
 190
 191                 if (metric) {
 192
 193                         for (int i = 1; i < 4; i++) {
 194
 195                                 vector< vector<double> > EuclidDists = linearCalc.calculateEuclidianDistance(G, i); //G is the pcoa file
 196
 197                                 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) {        remove(outputNames[i].c_str());  } return 0; }
 198
 199                                 double corr = linearCalc.calcPearson(EuclidDists, D); //G is the pcoa file, D is the users distance matrix
 200
 201                                 m->mothurOut("Pearson's coefficient using " + toString(i) + " axis: " + toString(corr)); m->mothurOutEndLine();
 202
 203                                 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) {        remove(outputNames[i].c_str());  } return 0; }
 204                         }
 205                 }
 206
 207                 m->mothurOutEndLine();
 208                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 209                 for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
 210                 m->mothurOutEndLine();
 211
 212                 return 0;
 213         }
 214         catch(exception& e) {
 215                 m->errorOut(e, "PCOACommand", "execute");
 216                 exit(1);
 217         }
 218 }
 219 /*********************************************************************************************************************************/
 220
 221 void PCOACommand::get_comment(istream& f, char begin, char end){
 222         try {
 223                 char d=f.get();
 224                 while(d != end){        d = f.get();    }
 225                 d = f.peek();
 226         }
 227         catch(exception& e) {
 228                 m->errorOut(e, "PCOACommand", "get_comment");
 229                 exit(1);
 230         }
 231 }
 232 /*********************************************************************************************************************************/
 233
 234 void PCOACommand::recenter(double offset, vector<vector<double> > D, vector<vector<double> >& G){
 235         try {
 236                 int rank = D.size();
 237
 238                 vector<vector<double> > A(rank);
 239                 vector<vector<double> > C(rank);
 240                 for(int i=0;i<rank;i++){
 241                         A[i].resize(rank);
 242                         C[i].resize(rank);
 243                 }
 244
 245                 double scale = -1.0000 / (double) rank;
 246
 247                 for(int i=0;i<rank;i++){
 248                         A[i][i] = 0.0000;
 249                         C[i][i] = 1.0000 + scale;
 250                         for(int j=i+1;j<rank;j++){
 251                                 A[i][j] = A[j][i] = -0.5 * D[i][j] * D[i][j] + offset;
 252                                 C[i][j] = C[j][i] = scale;
 253                         }
 254                 }
 255
 256                 A = linearCalc.matrix_mult(C,A);
 257                 G = linearCalc.matrix_mult(A,C);
 258         }
 259         catch(exception& e) {
 260                 m->errorOut(e, "PCOACommand", "recenter");
 261                 exit(1);
 262         }
 263
 264 }
 265
 266 /*********************************************************************************************************************************/
 267
 268 void PCOACommand::output(string fnameRoot, vector<string> name_list, vector<vector<double> >& G, vector<double> d) {
 269         try {
 270                 int rank = name_list.size();
 271                 double dsum = 0.0000;
 272                 for(int i=0;i<rank;i++){
 273                         dsum += d[i];
 274                         for(int j=0;j<rank;j++){
 275                                 if(d[j] >= 0)   {       G[i][j] *= pow(d[j],0.5);       }
 276                                 else                    {       G[i][j] = 0.00000;                      }
 277                         }
 278                 }
 279
 280                 ofstream pcaData((fnameRoot+"pcoa.axes").c_str(), ios::trunc);
 281                 pcaData.setf(ios::fixed, ios::floatfield);
 282                 pcaData.setf(ios::showpoint);
 283                 outputNames.push_back(fnameRoot+"pcoa.axes");
 284                 outputTypes["pcoa"].push_back(fnameRoot+"pcoa.axes");
 285
 286                 ofstream pcaLoadings((fnameRoot+"pcoa.loadings").c_str(), ios::trunc);
 287                 pcaLoadings.setf(ios::fixed, ios::floatfield);
 288                 pcaLoadings.setf(ios::showpoint);
 289                 outputNames.push_back(fnameRoot+"pcoa.loadings");
 290                 outputTypes["loadings"].push_back(fnameRoot+"pcoa.loadings");
 291
 292                 pcaLoadings << "axis\tloading\n";
 293                 for(int i=0;i<rank;i++){
 294                         pcaLoadings << i+1 << '\t' << d[i] * 100.0 / dsum << endl;
 295                 }
 296
 297                 pcaData << "group";
 298                 for(int i=0;i<rank;i++){
 299                         pcaData << '\t' << "axis" << i+1;
 300                 }
 301                 pcaData << endl;
 302
 303                 for(int i=0;i<rank;i++){
 304                         pcaData << name_list[i] << '\t';
 305                         for(int j=0;j<rank;j++){
 306                                 pcaData << G[i][j] << '\t';
 307                         }
 308                         pcaData << endl;
 309                 }
 310         }
 311         catch(exception& e) {
 312                 m->errorOut(e, "PCOACommand", "output");
 313                 exit(1);
 314         }
 315 }
 316
 317 /*********************************************************************************************************************************/
 318