pcacommand.cpp

   1 /*
   2  *  pcacommand.cpp
   3  *  mothur
   4  *
   5  *  Created by westcott on 1/7/11.
   6  *  Copyright 2011 Schloss Lab. All rights reserved.
   7  *
   8  */
   9
  10 #include "pcacommand.h"
  11 #include "inputdata.h"
  12
  13 //**********************************************************************************************************************
  14 vector<string> PCACommand::getValidParameters(){
  15         try {
  16                 string Array[] =  {"label", "groups","metric","outputdir","inputdir"};
  17                 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
  18                 return myArray;
  19         }
  20         catch(exception& e) {
  21                 m->errorOut(e, "PCACommand", "getValidParameters");
  22                 exit(1);
  23         }
  24 }
  25 //**********************************************************************************************************************
  26 PCACommand::PCACommand(){
  27         try {
  28                 abort = true; calledHelp = true;
  29                 vector<string> tempOutNames;
  30                 outputTypes["pca"] = tempOutNames;
  31                 outputTypes["loadings"] = tempOutNames;
  32         }
  33         catch(exception& e) {
  34                 m->errorOut(e, "PCACommand", "PCACommand");
  35                 exit(1);
  36         }
  37 }
  38 //**********************************************************************************************************************
  39 vector<string> PCACommand::getRequiredParameters(){
  40         try {
  41                 vector<string> myArray;
  42                 return myArray;
  43         }
  44         catch(exception& e) {
  45                 m->errorOut(e, "PCACommand", "getRequiredParameters");
  46                 exit(1);
  47         }
  48 }
  49 //**********************************************************************************************************************
  50 vector<string> PCACommand::getRequiredFiles(){
  51         try {
  52                 string Array[] =  {"shared","relabund","or"};
  53                 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
  54                 return myArray;
  55         }
  56         catch(exception& e) {
  57                 m->errorOut(e, "PCACommand", "getRequiredFiles");
  58                 exit(1);
  59         }
  60 }
  61 //**********************************************************************************************************************
  62
  63 PCACommand::PCACommand(string option)  {
  64         try {
  65                 abort = false; calledHelp = false;
  66
  67                 globaldata = GlobalData::getInstance();
  68
  69                 //allow user to run help
  70                 if(option == "help") { help(); abort = true; calledHelp = true; }
  71
  72                 else {
  73                         //valid paramters for this command
  74                         string Array[] =  {"label","groups","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);              if (inputDir == "not found"){   inputDir = "";          }
  89
  90                         //initialize outputTypes
  91                         vector<string> tempOutNames;
  92                         outputTypes["pca"] = tempOutNames;
  93                         outputTypes["loadings"] = tempOutNames;
  94
  95                         //make sure the user has already run the read.otu command
  96                         if ((globaldata->getSharedFile() == "") && (globaldata->getRelAbundFile() == "")) {
  97                                 m->mothurOut("You must read a list and a group, shared or relabund file before you can use the pca command."); m->mothurOutEndLine(); abort = true;
  98                         }
  99
 100                         if (globaldata->getSharedFile() != "")          { mode = "shared"; inputFile = globaldata->getSharedFile();             }
 101                         if (globaldata->getRelAbundFile() != "")        { mode = "relabund"; inputFile = globaldata->getRelAbundFile(); }
 102
 103                         //if the user changes the output directory command factory will send this info to us in the output parameter
 104                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){
 105                                 outputDir = "";
 106                                 outputDir += m->hasPath(inputFile); //if user entered a file with a path then preserve it
 107                         }
 108
 109                         string temp = validParameter.validFile(parameters, "metric", false);    if (temp == "not found"){       temp = "T";                             }
 110                         metric = m->isTrue(temp);
 111
 112                         label = validParameter.validFile(parameters, "label", false);
 113                         if (label == "not found") { label = ""; labels = globaldata->labels; if(labels.size() == 0) {  m->mothurOut("You did not provide a label, I will use the first label in your inputfile."); m->mothurOutEndLine(); } }
 114                         else { m->splitAtDash(label, labels); }
 115
 116                         groups = validParameter.validFile(parameters, "groups", false);
 117                         if (groups == "not found") { groups = "";  }
 118                         else { m->splitAtDash(groups, Groups);  }
 119                         globaldata->Groups = Groups;
 120
 121                 }
 122
 123         }
 124         catch(exception& e) {
 125                 m->errorOut(e, "PCACommand", "PCACommand");
 126                 exit(1);
 127         }
 128 }
 129 //**********************************************************************************************************************
 130 void PCACommand::help(){
 131         try {
 132                 m->mothurOut("The pca command can only be run after a successful read.otu command of a shared or relabund file."); m->mothurOutEndLine();
 133                 m->mothurOut("The pca command parameters are label, groups and metric. No parameters are required."); m->mothurOutEndLine();
 134                 m->mothurOut("The label parameter is used to analyze specific labels in your input. Default is the first label in your shared or relabund file. Multiple labels may be separated by dashes.\n");
 135                 m->mothurOut("The groups parameter allows you to specify which groups you would like analyzed. Groupnames are separated by dashes.\n");
 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 pca(groups=yourGroups).\n");
 138                 m->mothurOut("Example pca(groups=A-B-C).\n");
 139                 m->mothurOut("Note: No spaces between parameter labels (i.e. groups), '=' and parameters (i.e.yourGroups).\n\n");
 140         }
 141         catch(exception& e) {
 142                 m->errorOut(e, "PCACommand", "help");
 143                 exit(1);
 144         }
 145 }
 146 //**********************************************************************************************************************
 147 PCACommand::~PCACommand(){}
 148 //**********************************************************************************************************************
 149 int PCACommand::execute(){
 150         try {
 151
 152                 if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
 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                 //get first line of shared file
 160                 vector< vector<double> > matrix;
 161                 InputData* input;
 162                 if (mode == "shared")                   {
 163                         input = new InputData(inputFile, "sharedfile");
 164                 }else if (mode == "relabund")   {
 165                         input = new InputData(inputFile, "relabund");
 166                 }else {  m->mothurOut("[ERROR]: filetype not recognized."); m->mothurOutEndLine();  return 0; }
 167
 168                 vector<SharedRAbundFloatVector*> lookupFloat = input->getSharedRAbundFloatVectors();
 169                 string lastLabel = lookupFloat[0]->getLabel();
 170
 171                 set<string> processedLabels;
 172                 set<string> userLabels = labels;
 173
 174                 //if the user gave no labels, then use the first one read
 175                 if (labels.size() == 0) {
 176                         label = lastLabel;
 177
 178                         process(lookupFloat);
 179                 }
 180
 181                 //as long as you are not at the end of the file or done wih the lines you want
 182                 while((lookupFloat[0] != NULL) && (userLabels.size() != 0)) {
 183
 184                         if (m->control_pressed) {  for (int i = 0; i < outputNames.size(); i++) {       remove(outputNames[i].c_str());  } delete input; for (int i = 0; i < lookupFloat.size(); i++) {  delete lookupFloat[i];  }  lookupFloat.clear(); return 0;  }
 185
 186                         if(labels.count(lookupFloat[0]->getLabel()) == 1){
 187                                 processedLabels.insert(lookupFloat[0]->getLabel());
 188                                 userLabels.erase(lookupFloat[0]->getLabel());
 189
 190                                 process(lookupFloat);
 191                         }
 192
 193                         if ((m->anyLabelsToProcess(lookupFloat[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
 194                                 string saveLabel = lookupFloat[0]->getLabel();
 195
 196                                 for (int i = 0; i < lookupFloat.size(); i++) {  delete lookupFloat[i];  }  lookupFloat.clear();
 197                                 lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
 198
 199                                 process(lookupFloat);
 200
 201                                 processedLabels.insert(lookupFloat[0]->getLabel());
 202                                 userLabels.erase(lookupFloat[0]->getLabel());
 203
 204                                 //restore real lastlabel to save below
 205                                 lookupFloat[0]->setLabel(saveLabel);
 206                         }
 207
 208                         lastLabel = lookupFloat[0]->getLabel();
 209
 210                         //get next line to process
 211                         //prevent memory leak
 212                         for (int i = 0; i < lookupFloat.size(); i++) {  delete lookupFloat[i];  } lookupFloat.clear();
 213                         lookupFloat = input->getSharedRAbundFloatVectors();
 214                 }
 215
 216
 217                 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) {        remove(outputNames[i].c_str());  } delete input; for (int i = 0; i < lookupFloat.size(); i++) {  delete lookupFloat[i];  } lookupFloat.clear(); return 0;  }
 218
 219                 //output error messages about any remaining user labels
 220                 set<string>::iterator it;
 221                 bool needToRun = false;
 222                 for (it = userLabels.begin(); it != userLabels.end(); it++) {
 223                         m->mothurOut("Your file does not include the label " + *it);
 224                         if (processedLabels.count(lastLabel) != 1) {
 225                                 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
 226                                 needToRun = true;
 227                         }else {
 228                                 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
 229                         }
 230                 }
 231
 232                 //run last label if you need to
 233                 if (needToRun == true)  {
 234                         for (int i = 0; i < lookupFloat.size(); i++) {  if (lookupFloat[i] != NULL) {   delete lookupFloat[i];  } }  lookupFloat.clear();
 235                         lookupFloat = input->getSharedRAbundFloatVectors(lastLabel);
 236
 237                         process(lookupFloat);
 238
 239                         for (int i = 0; i < lookupFloat.size(); i++) {  if (lookupFloat[i] != NULL) {   delete lookupFloat[i];  } } lookupFloat.clear();
 240                 }
 241
 242                 for (int i = 0; i < lookupFloat.size(); i++) {  if (lookupFloat[i] != NULL) {   delete lookupFloat[i];  } } lookupFloat.clear();
 243                 delete input;
 244
 245                 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) {        remove(outputNames[i].c_str());  } return 0; }
 246
 247                 m->mothurOutEndLine();
 248                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 249                 for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
 250                 m->mothurOutEndLine();
 251
 252                 return 0;
 253         }
 254         catch(exception& e) {
 255                 m->errorOut(e, "PCACommand", "execute");
 256                 exit(1);
 257         }
 258 }
 259 /**********************************************************************************************************************
 260 vector< vector<double> > PCACommand::createMatrix(vector<SharedRAbundFloatVector*> lookupFloat){
 261         try {
 262                 vector< vector<double> > matrix; matrix.resize(lookupFloat.size());
 263
 264                 //fill matrix with shared files relative abundances
 265                 for (int i = 0; i < lookupFloat.size(); i++) {
 266                         for (int j = 0; j < lookupFloat[i]->getNumBins(); j++) {
 267                                 matrix[i].push_back(lookupFloat[i]->getAbundance(j));
 268                         }
 269                 }
 270
 271                 vector< vector<double> > transposeMatrix; transposeMatrix.resize(matrix[0].size());
 272                 for (int i = 0; i < transposeMatrix.size(); i++) {
 273                         for (int j = 0; j < matrix.size(); j++) {
 274                                 transposeMatrix[i].push_back(matrix[j][i]);
 275                         }
 276                 }
 277
 278                 matrix = linearCalc.matrix_mult(matrix, transposeMatrix);
 279
 280                 return matrix;
 281         }
 282         catch(exception& e) {
 283                 m->errorOut(e, "PCACommand", "createMatrix");
 284                 exit(1);
 285         }
 286 }*/
 287 //**********************************************************************************************************************
 288 int PCACommand::process(vector<SharedRAbundFloatVector*>& lookupFloat){
 289         try {
 290                 m->mothurOut("\nProcessing " + lookupFloat[0]->getLabel()); m->mothurOutEndLine();
 291
 292                 vector< vector<double> > matrix; matrix.resize(lookupFloat.size());
 293
 294                 ofstream out;
 295                 string temp = outputDir + "matrix.transpose.out";
 296                 m->openOutputFile(temp, out);
 297                 out << "matrix" << endl;
 298
 299                 //fill matrix with shared files relative abundances
 300                 for (int i = 0; i < lookupFloat.size(); i++) {
 301                         for (int j = 0; j < lookupFloat[i]->getNumBins(); j++) {
 302                                 matrix[i].push_back(lookupFloat[i]->getAbundance(j));
 303                                 out << lookupFloat[i]->getAbundance(j) << '\t';
 304                         }
 305                         out << endl;
 306                 }
 307                 out << endl << endl << "transpose" << endl;
 308                 vector< vector<double> > transposeMatrix; transposeMatrix.resize(matrix[0].size());
 309                 for (int i = 0; i < transposeMatrix.size(); i++) {
 310                         for (int j = 0; j < matrix.size(); j++) {
 311                                 transposeMatrix[i].push_back(matrix[j][i]);
 312                                 out << matrix[j][i] << '\t';
 313                         }
 314                         out << endl;
 315                 }
 316
 317                 matrix = linearCalc.matrix_mult(matrix, transposeMatrix);
 318
 319                 out << endl << endl << "matrix mult" << endl;
 320                 for (int i = 0; i < matrix.size(); i++) {
 321                         for (int j = 0; j < matrix[i].size(); j++) {
 322                                 out << matrix[i][j] << '\t';
 323                         }
 324                     out << endl;
 325                 }
 326                 out.close();
 327
 328
 329                 double offset = 0.0000;
 330                 vector<double> d;
 331                 vector<double> e;
 332                 vector<vector<double> > G = matrix;
 333                 //vector<vector<double> > copy_G;
 334
 335                 for(int count=0;count<2;count++){
 336                         linearCalc.recenter(offset, matrix, G);         if (m->control_pressed) { return 0; }
 337                         linearCalc.tred2(G, d, e);                                      if (m->control_pressed) { return 0; }
 338                         linearCalc.qtli(d, e, G);                                       if (m->control_pressed) { return 0; }
 339                         offset = d[d.size()-1];
 340                         if(offset > 0.0) break;
 341                 }
 342
 343                 if (m->control_pressed) { return 0; }
 344
 345                 string fbase = outputDir + m->getRootName(m->getSimpleName(inputFile));
 346                 string outputFileName = fbase + lookupFloat[0]->getLabel();
 347                 output(outputFileName, globaldata->Groups, G, d);
 348
 349                 if (metric) {
 350
 351                         for (int i = 1; i < 4; i++) {
 352
 353                                 vector< vector<double> > EuclidDists = linearCalc.calculateEuclidianDistance(G, i); //G is the pca file
 354
 355                                 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) {        remove(outputNames[i].c_str());  } return 0; }
 356
 357                                 double corr = linearCalc.calcPearson(EuclidDists, matrix); //G is the pca file, D is the users distance matrix
 358
 359                                 m->mothurOut("Pearson's coefficient using " + toString(i) + " axis: " + toString(corr)); m->mothurOutEndLine();
 360
 361                                 m->mothurOut("Rsq " + toString(i) + " axis: " + toString(corr * corr)); m->mothurOutEndLine();
 362
 363                                 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) {        remove(outputNames[i].c_str());  } return 0; }
 364                         }
 365                 }
 366
 367                 return 0;
 368         }
 369         catch(exception& e) {
 370                 m->errorOut(e, "PCACommand", "process");
 371                 exit(1);
 372         }
 373 }
 374 /*********************************************************************************************************************************/
 375
 376 void PCACommand::output(string fnameRoot, vector<string> name_list, vector<vector<double> >& G, vector<double> d) {
 377         try {
 378                 int rank = name_list.size();
 379                 double dsum = 0.0000;
 380                 for(int i=0;i<rank;i++){
 381                         dsum += d[i];
 382                         for(int j=0;j<rank;j++){
 383                                 if(d[j] >= 0)   {       G[i][j] *= pow(d[j],0.5);       }
 384                                 else                    {       G[i][j] = 0.00000;                      }
 385                         }
 386                 }
 387
 388                 ofstream pcaData((fnameRoot+".pca.axes").c_str(), ios::trunc);
 389                 pcaData.setf(ios::fixed, ios::floatfield);
 390                 pcaData.setf(ios::showpoint);
 391                 outputNames.push_back(fnameRoot+".pca.axes");
 392                 outputTypes["pca"].push_back(fnameRoot+".pca.axes");
 393
 394                 ofstream pcaLoadings((fnameRoot+".pca.loadings").c_str(), ios::trunc);
 395                 pcaLoadings.setf(ios::fixed, ios::floatfield);
 396                 pcaLoadings.setf(ios::showpoint);
 397                 outputNames.push_back(fnameRoot+".pca.loadings");
 398                 outputTypes["loadings"].push_back(fnameRoot+".pca.loadings");
 399
 400                 pcaLoadings << "axis\tloading\n";
 401                 for(int i=0;i<rank;i++){
 402                         pcaLoadings << i+1 << '\t' << d[i] * 100.0 / dsum << endl;
 403                 }
 404
 405                 pcaData << "group";
 406                 for(int i=0;i<rank;i++){
 407                         pcaData << '\t' << "axis" << i+1;
 408                 }
 409                 pcaData << endl;
 410
 411                 for(int i=0;i<rank;i++){
 412                         pcaData << name_list[i] << '\t';
 413                         for(int j=0;j<rank;j++){
 414                                 pcaData << G[i][j] << '\t';
 415                         }
 416                         pcaData << endl;
 417                 }
 418         }
 419         catch(exception& e) {
 420                 m->errorOut(e, "PCACommand", "output");
 421                 exit(1);
 422         }
 423 }
 424 /*********************************************************************************************************************************/
 425
 426