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
  13 //**********************************************************************************************************************
  14 vector<string> PCOACommand::getValidParameters(){
  15         try {
  16                 string Array[] =  {"phylip", "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, "PCOACommand", "getValidParameters");
  22                 exit(1);
  23         }
  24 }
  25 //**********************************************************************************************************************
  26 PCOACommand::PCOACommand(){
  27         try {
  28                 abort = true;
  29                 //initialize outputTypes
  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;
  67
  68                 //allow user to run help
  69                 if(option == "help") { help(); abort = 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) { return 0; }
 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                 read(filename, names, D);
 164
 165                 if (m->control_pressed) { return 0; }
 166
 167                 double offset = 0.0000;
 168                 vector<double> d;
 169                 vector<double> e;
 170                 vector<vector<double> > G = D;
 171                 vector<vector<double> > copy_G;
 172
 173                 m->mothurOut("\nProcessing...\n\n");
 174
 175                 for(int count=0;count<2;count++){
 176                         recenter(offset, D, G);                                 if (m->control_pressed) { return 0; }
 177                         linearCalc.tred2(G, d, e);                              if (m->control_pressed) { return 0; }
 178                         linearCalc.qtli(d, e, G);                               if (m->control_pressed) { return 0; }
 179                         offset = d[d.size()-1];
 180                         if(offset > 0.0) break;
 181                 }
 182
 183                 if (m->control_pressed) { return 0; }
 184
 185                 output(fbase, names, G, d);
 186
 187                 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) {        remove(outputNames[i].c_str());  } return 0; }
 188
 189                 if (metric) {
 190
 191                         for (int i = 1; i < 4; i++) {
 192
 193                                 vector< vector<double> > EuclidDists = linearCalc.calculateEuclidianDistance(G, i); //G is the pcoa file
 194
 195                                 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) {        remove(outputNames[i].c_str());  } return 0; }
 196
 197                                 double corr = linearCalc.calcPearson(EuclidDists, D); //G is the pcoa file, D is the users distance matrix
 198
 199                                 m->mothurOut("Pearson's coefficient using " + toString(i) + " axis: " + toString(corr)); m->mothurOutEndLine();
 200
 201                                 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) {        remove(outputNames[i].c_str());  } return 0; }
 202                         }
 203                 }
 204
 205                 m->mothurOutEndLine();
 206                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 207                 for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
 208                 m->mothurOutEndLine();
 209
 210                 return 0;
 211         }
 212         catch(exception& e) {
 213                 m->errorOut(e, "PCOACommand", "execute");
 214                 exit(1);
 215         }
 216 }
 217 /*********************************************************************************************************************************/
 218
 219 void PCOACommand::get_comment(istream& f, char begin, char end){
 220         try {
 221                 char d=f.get();
 222                 while(d != end){        d = f.get();    }
 223                 d = f.peek();
 224         }
 225         catch(exception& e) {
 226                 m->errorOut(e, "PCOACommand", "get_comment");
 227                 exit(1);
 228         }
 229 }
 230
 231 /*********************************************************************************************************************************/
 232
 233 int PCOACommand::read_phylip(istream& f, int square_m, vector<string>& name_list, vector<vector<double> >& d){
 234         try {
 235                 //     int count1=0;
 236                 //     int count2=0;
 237
 238                 int rank;
 239                 f >> rank;
 240
 241                 name_list.resize(rank);
 242                 d.resize(rank);
 243                 if(square_m == 1){
 244                         for(int i=0;i<rank;i++)
 245                                 d[i].resize(rank);
 246                         for(int i=0;i<rank;i++) {
 247                                 f >> name_list[i];
 248                                 //                      cout << i << "\t" << name_list[i] << endl;
 249                                 for(int j=0;j<rank;j++) {
 250                                         if (m->control_pressed) { return 0; }
 251
 252                                         f >> d[i][j];
 253                                         if (d[i][j] == -0.0000)
 254                                                 d[i][j] = 0.0000;
 255                                 }
 256                         }
 257                 }
 258                 else if(square_m == 2){
 259                         for(int i=0;i<rank;i++){
 260                                 d[i].resize(rank);
 261                         }
 262                         d[0][0] = 0.0000;
 263                         f >> name_list[0];
 264                         for(int i=1;i<rank;i++){
 265                                 f >> name_list[i];
 266                                 d[i][i]=0.0000;
 267                                 for(int j=0;j<i;j++){
 268                                         if (m->control_pressed) { return 0; }
 269                                         f >> d[i][j];
 270                                         if (d[i][j] == -0.0000)
 271                                                 d[i][j] = 0.0000;
 272                                         d[j][i]=d[i][j];
 273                                 }
 274                         }
 275                 }
 276
 277                 return 0;
 278         }
 279         catch(exception& e) {
 280                 m->errorOut(e, "PCOACommand", "read_phylip");
 281                 exit(1);
 282         }
 283
 284 }
 285
 286 /*********************************************************************************************************************************/
 287
 288 void PCOACommand::read(string fname, vector<string>& names, vector<vector<double> >& D){
 289         try {
 290                 ifstream f;
 291                 m->openInputFile(fname, f);
 292
 293                 //check whether matrix is square
 294                 char d;
 295                 int q = 1;
 296                 int numSeqs;
 297                 string name;
 298
 299                 f >> numSeqs >> name;
 300
 301                 while((d=f.get()) != EOF){
 302
 303                         //is d a number meaning its square
 304                         if(isalnum(d)){
 305                                 q = 1;
 306                                 break;
 307                         }
 308
 309                         //is d a line return meaning its lower triangle
 310                         if(d == '\n'){
 311                                 q = 2;
 312                                 break;
 313                         }
 314                 }
 315                 f.close();
 316
 317                 //reopen to get back to beginning
 318                 m->openInputFile(fname, f);
 319                 read_phylip(f, q, names, D);
 320         }
 321                 catch(exception& e) {
 322                 m->errorOut(e, "PCOACommand", "read");
 323                 exit(1);
 324         }
 325 }
 326
 327 /*********************************************************************************************************************************/
 328
 329 void PCOACommand::recenter(double offset, vector<vector<double> > D, vector<vector<double> >& G){
 330         try {
 331                 int rank = D.size();
 332
 333                 vector<vector<double> > A(rank);
 334                 vector<vector<double> > C(rank);
 335                 for(int i=0;i<rank;i++){
 336                         A[i].resize(rank);
 337                         C[i].resize(rank);
 338                 }
 339
 340                 double scale = -1.0000 / (double) rank;
 341
 342                 for(int i=0;i<rank;i++){
 343                         A[i][i] = 0.0000;
 344                         C[i][i] = 1.0000 + scale;
 345                         for(int j=i+1;j<rank;j++){
 346                                 A[i][j] = A[j][i] = -0.5 * D[i][j] * D[i][j] + offset;
 347                                 C[i][j] = C[j][i] = scale;
 348                         }
 349                 }
 350
 351                 A = linearCalc.matrix_mult(C,A);
 352                 G = linearCalc.matrix_mult(A,C);
 353         }
 354         catch(exception& e) {
 355                 m->errorOut(e, "PCOACommand", "recenter");
 356                 exit(1);
 357         }
 358
 359 }
 360
 361 /*********************************************************************************************************************************/
 362
 363 void PCOACommand::output(string fnameRoot, vector<string> name_list, vector<vector<double> >& G, vector<double> d) {
 364         try {
 365                 int rank = name_list.size();
 366                 double dsum = 0.0000;
 367                 for(int i=0;i<rank;i++){
 368                         dsum += d[i];
 369                         for(int j=0;j<rank;j++){
 370                                 if(d[j] >= 0)   {       G[i][j] *= pow(d[j],0.5);       }
 371                                 else                    {       G[i][j] = 0.00000;                      }
 372                         }
 373                 }
 374
 375                 ofstream pcaData((fnameRoot+"pcoa").c_str(), ios::trunc);
 376                 pcaData.setf(ios::fixed, ios::floatfield);
 377                 pcaData.setf(ios::showpoint);
 378                 outputNames.push_back(fnameRoot+"pcoa");
 379                 outputTypes["pcoa"].push_back(fnameRoot+"pcoa");
 380
 381                 ofstream pcaLoadings((fnameRoot+"pcoa.loadings").c_str(), ios::trunc);
 382                 pcaLoadings.setf(ios::fixed, ios::floatfield);
 383                 pcaLoadings.setf(ios::showpoint);
 384                 outputNames.push_back(fnameRoot+"pcoa.loadings");
 385                 outputTypes["loadings"].push_back(fnameRoot+"pcoa.loadings");
 386
 387                 pcaLoadings << "axis\tloading\n";
 388                 for(int i=0;i<rank;i++){
 389                         pcaLoadings << i+1 << '\t' << d[i] * 100.0 / dsum << endl;
 390                 }
 391
 392                 pcaData << "group";
 393                 for(int i=0;i<rank;i++){
 394                         pcaData << '\t' << "axis" << i+1;
 395                 }
 396                 pcaData << endl;
 397
 398                 for(int i=0;i<rank;i++){
 399                         pcaData << name_list[i] << '\t';
 400                         for(int j=0;j<rank;j++){
 401                                 pcaData << G[i][j] << '\t';
 402                         }
 403                         pcaData << endl;
 404                 }
 405         }
 406         catch(exception& e) {
 407                 m->errorOut(e, "PCOACommand", "output");
 408                 exit(1);
 409         }
 410 }
 411
 412 /*********************************************************************************************************************************/
 413