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