pcacommand.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 "pcacommand.h"
  12
  13 //**********************************************************************************************************************
  14 vector<string> PCACommand::getValidParameters(){
  15         try {
  16                 string Array[] =  {"phylip", "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                 //initialize outputTypes
  29                 vector<string> tempOutNames;
  30                 outputTypes["pcoa"] = 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                 string Array[] =  {"phylip"};
  42                 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
  43                 return myArray;
  44         }
  45         catch(exception& e) {
  46                 m->errorOut(e, "PCACommand", "getRequiredParameters");
  47                 exit(1);
  48         }
  49 }
  50 //**********************************************************************************************************************
  51 vector<string> PCACommand::getRequiredFiles(){
  52         try {
  53                 vector<string> myArray;
  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;
  66
  67                 //allow user to run help
  68                 if(option == "help") { help(); abort = true; }
  69
  70                 else {
  71                         //valid paramters for this command
  72                         string Array[] =  {"phylip","outputdir", "inputdir"};
  73                         vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
  74
  75                         OptionParser parser(option);
  76                         map<string, string> parameters = parser. getParameters();
  77
  78                         ValidParameters validParameter;
  79                         map<string, string>::iterator it;
  80
  81                         //check to make sure all parameters are valid for command
  82                         for (it = parameters.begin(); it != parameters.end(); it++) {
  83                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
  84                         }
  85                         //if the user changes the input directory command factory will send this info to us in the output parameter
  86                         string inputDir = validParameter.validFile(parameters, "inputdir", false);
  87                         if (inputDir == "not found"){   inputDir = "";          }
  88                         else {
  89                                 string path;
  90                                 it = parameters.find("phylip");
  91                                 //user has given a template file
  92                                 if(it != parameters.end()){
  93                                         path = m->hasPath(it->second);
  94                                         //if the user has not given a path then, add inputdir. else leave path alone.
  95                                         if (path == "") {       parameters["phylip"] = inputDir + it->second;           }
  96                                 }
  97                         }
  98
  99                         //initialize outputTypes
 100                         vector<string> tempOutNames;
 101                         outputTypes["pcoa"] = tempOutNames;
 102                         outputTypes["loadings"] = tempOutNames;
 103
 104                         //required parameters
 105                         phylipfile = validParameter.validFile(parameters, "phylip", true);
 106                         if (phylipfile == "not open") { abort = true; }
 107                         else if (phylipfile == "not found") { phylipfile = ""; abort = true; }
 108                         else {  filename = phylipfile;  }
 109
 110                         //if the user changes the output directory command factory will send this info to us in the output parameter
 111                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){
 112                                 outputDir = "";
 113                                 outputDir += m->hasPath(phylipfile); //if user entered a file with a path then preserve it
 114                         }
 115
 116                         //error checking on files
 117                         if (phylipfile == "")   { m->mothurOut("You must provide a distance file before running the pca command."); m->mothurOutEndLine(); abort = true; }
 118                 }
 119
 120         }
 121         catch(exception& e) {
 122                 m->errorOut(e, "PCACommand", "PCACommand");
 123                 exit(1);
 124         }
 125 }
 126 //**********************************************************************************************************************
 127 void PCACommand::help(){
 128         try {
 129
 130                 m->mothurOut("The pca command..."); m->mothurOutEndLine();
 131         }
 132         catch(exception& e) {
 133                 m->errorOut(e, "PCACommand", "help");
 134                 exit(1);
 135         }
 136 }
 137 //**********************************************************************************************************************
 138 PCACommand::~PCACommand(){}
 139 //**********************************************************************************************************************
 140 int PCACommand::execute(){
 141         try {
 142
 143                 if (abort == true) { return 0; }
 144
 145                 cout.setf(ios::fixed, ios::floatfield);
 146                 cout.setf(ios::showpoint);
 147                 cerr.setf(ios::fixed, ios::floatfield);
 148                 cerr.setf(ios::showpoint);
 149
 150                 vector<string> names;
 151                 vector<vector<double> > D;
 152
 153                 fbase = outputDir + m->getRootName(m->getSimpleName(filename));
 154
 155                 read(filename, names, D);
 156
 157                 if (m->control_pressed) { return 0; }
 158
 159                 double offset = 0.0000;
 160                 vector<double> d;
 161                 vector<double> e;
 162                 vector<vector<double> > G = D;
 163                 vector<vector<double> > copy_G;
 164                 //int rank = D.size();
 165
 166                 m->mothurOut("\nProcessing...\n");
 167
 168                 for(int count=0;count<2;count++){
 169                         recenter(offset, D, G);         if (m->control_pressed) { return 0; }
 170                         tred2(G, d, e);                         if (m->control_pressed) { return 0; }
 171                         qtli(d, e, G);                          if (m->control_pressed) { return 0; }
 172                         offset = d[d.size()-1];
 173                         if(offset > 0.0) break;
 174                 }
 175
 176                 if (m->control_pressed) { return 0; }
 177
 178                 output(fbase, names, G, d);
 179
 180                 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) {        remove(outputNames[i].c_str());  } return 0; }
 181
 182                 m->mothurOutEndLine();
 183                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 184                 for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
 185                 m->mothurOutEndLine();
 186
 187                 return 0;
 188         }
 189         catch(exception& e) {
 190                 m->errorOut(e, "PCACommand", "execute");
 191                 exit(1);
 192         }
 193 }
 194 /*********************************************************************************************************************************/
 195
 196 inline double SIGN(const double a, const double b)
 197 {
 198     return b>=0 ? (a>=0 ? a:-a) : (a>=0 ? -a:a);
 199 }
 200
 201 /*********************************************************************************************************************************/
 202
 203 void PCACommand::get_comment(istream& f, char begin, char end){
 204         try {
 205                 char d=f.get();
 206                 while(d != end){        d = f.get();    }
 207                 d = f.peek();
 208         }
 209         catch(exception& e) {
 210                 m->errorOut(e, "PCACommand", "get_comment");
 211                 exit(1);
 212         }
 213 }
 214
 215 /*********************************************************************************************************************************/
 216
 217 int PCACommand::read_phylip(istream& f, int square_m, vector<string>& name_list, vector<vector<double> >& d){
 218         try {
 219                 //     int count1=0;
 220                 //     int count2=0;
 221
 222                 int rank;
 223                 f >> rank;
 224
 225                 name_list.resize(rank);
 226                 d.resize(rank);
 227                 if(square_m == 1){
 228                         for(int i=0;i<rank;i++)
 229                                 d[i].resize(rank);
 230                         for(int i=0;i<rank;i++) {
 231                                 f >> name_list[i];
 232                                 //                      cout << i << "\t" << name_list[i] << endl;
 233                                 for(int j=0;j<rank;j++) {
 234                                         if (m->control_pressed) { return 0; }
 235
 236                                         f >> d[i][j];
 237                                         if (d[i][j] == -0.0000)
 238                                                 d[i][j] = 0.0000;
 239                                 }
 240                         }
 241                 }
 242                 else if(square_m == 2){
 243                         for(int i=0;i<rank;i++){
 244                                 d[i].resize(rank);
 245                         }
 246                         d[0][0] = 0.0000;
 247                         f >> name_list[0];
 248                         for(int i=1;i<rank;i++){
 249                                 f >> name_list[i];
 250                                 d[i][i]=0.0000;
 251                                 for(int j=0;j<i;j++){
 252                                         if (m->control_pressed) { return 0; }
 253                                         f >> d[i][j];
 254                                         if (d[i][j] == -0.0000)
 255                                                 d[i][j] = 0.0000;
 256                                         d[j][i]=d[i][j];
 257                                 }
 258                         }
 259                 }
 260
 261                 return 0;
 262         }
 263         catch(exception& e) {
 264                 m->errorOut(e, "PCACommand", "read_phylip");
 265                 exit(1);
 266         }
 267
 268 }
 269
 270 /*********************************************************************************************************************************/
 271
 272 void PCACommand::read(string fname, vector<string>& names, vector<vector<double> >& D){
 273         try {
 274                 ifstream f;
 275                 m->openInputFile(fname, f);
 276
 277                 //check whether matrix is square
 278                 char d;
 279                 int q = 1;
 280                 int numSeqs;
 281                 string name;
 282
 283                 f >> numSeqs >> name;
 284
 285                 while((d=f.get()) != EOF){
 286
 287                         //is d a number meaning its square
 288                         if(isalnum(d)){
 289                                 q = 1;
 290                                 break;
 291                         }
 292
 293                         //is d a line return meaning its lower triangle
 294                         if(d == '\n'){
 295                                 q = 2;
 296                                 break;
 297                         }
 298                 }
 299                 f.close();
 300
 301                 //reopen to get back to beginning
 302                 m->openInputFile(fname, f);
 303                 read_phylip(f, q, names, D);
 304         }
 305                 catch(exception& e) {
 306                 m->errorOut(e, "PCACommand", "read");
 307                 exit(1);
 308         }
 309 }
 310
 311 /*********************************************************************************************************************************/
 312
 313 double PCACommand::pythag(double a, double b)   {       return(pow(a*a+b*b,0.5));       }
 314
 315 /*********************************************************************************************************************************/
 316
 317 void PCACommand::matrix_mult(vector<vector<double> > first, vector<vector<double> > second, vector<vector<double> >& product){
 318         try {
 319                 int first_rows = first.size();
 320                 int first_cols = first[0].size();
 321                 int second_cols = second[0].size();
 322
 323                 product.resize(first_rows);
 324                 for(int i=0;i<first_rows;i++){
 325                         product[i].resize(second_cols);
 326                 }
 327
 328                 for(int i=0;i<first_rows;i++){
 329                         for(int j=0;j<second_cols;j++){
 330                                 product[i][j] = 0.0;
 331                                 for(int k=0;k<first_cols;k++){
 332                                         product[i][j] += first[i][k] * second[k][j];
 333                                 }
 334                         }
 335                 }
 336         }
 337         catch(exception& e) {
 338                 m->errorOut(e, "PCACommand", "matrix_mult");
 339                 exit(1);
 340         }
 341
 342 }
 343
 344 /*********************************************************************************************************************************/
 345
 346 void PCACommand::recenter(double offset, vector<vector<double> > D, vector<vector<double> >& G){
 347         try {
 348                 int rank = D.size();
 349
 350                 vector<vector<double> > A(rank);
 351                 vector<vector<double> > C(rank);
 352                 for(int i=0;i<rank;i++){
 353                         A[i].resize(rank);
 354                         C[i].resize(rank);
 355                 }
 356
 357                 double scale = -1.0000 / (double) rank;
 358
 359                 for(int i=0;i<rank;i++){
 360                         A[i][i] = 0.0000;
 361                         C[i][i] = 1.0000 + scale;
 362                         for(int j=i+1;j<rank;j++){
 363                                 A[i][j] = A[j][i] = -0.5 * D[i][j] * D[i][j] + offset;
 364                                 C[i][j] = C[j][i] = scale;
 365                         }
 366                 }
 367
 368                 matrix_mult(C,A,A);
 369                 matrix_mult(A,C,G);
 370         }
 371         catch(exception& e) {
 372                 m->errorOut(e, "PCACommand", "recenter");
 373                 exit(1);
 374         }
 375
 376 }
 377
 378 /*********************************************************************************************************************************/
 379
 380 //  This function is taken from Numerical Recipes in C++ by Press et al., 2nd edition, pg. 479
 381
 382 void PCACommand::tred2(vector<vector<double> >& a, vector<double>& d, vector<double>& e){
 383         try {
 384                 double scale, hh, h, g, f;
 385
 386                 int n = a.size();
 387
 388                 d.resize(n);
 389                 e.resize(n);
 390
 391                 for(int i=n-1;i>0;i--){
 392                         int l=i-1;
 393                         h = scale = 0.0000;
 394                         if(l>0){
 395                                 for(int k=0;k<l+1;k++){
 396                                         scale += fabs(a[i][k]);
 397                                 }
 398                                 if(scale == 0.0){
 399                                         e[i] = a[i][l];
 400                                 }
 401                                 else{
 402                                         for(int k=0;k<l+1;k++){
 403                                                 a[i][k] /= scale;
 404                                                 h += a[i][k] * a[i][k];
 405                                         }
 406                                         f = a[i][l];
 407                                         g = (f >= 0.0 ? -sqrt(h) : sqrt(h));
 408                                         e[i] = scale * g;
 409                                         h -= f * g;
 410                                         a[i][l] = f - g;
 411                                         f = 0.0;
 412                                         for(int j=0;j<l+1;j++){
 413                                                 a[j][i] = a[i][j] / h;
 414                                                 g = 0.0;
 415                                                 for(int k=0;k<j+1;k++){
 416                                                         g += a[j][k] * a[i][k];
 417                                                 }
 418                                                 for(int k=j+1;k<l+1;k++){
 419                                                         g += a[k][j] * a[i][k];
 420                                                 }
 421                                                 e[j] = g / h;
 422                                                 f += e[j] * a[i][j];
 423                                         }
 424                                         hh = f / (h + h);
 425                                         for(int j=0;j<l+1;j++){
 426                                                 f = a[i][j];
 427                                                 e[j] = g = e[j] - hh * f;
 428                                                 for(int k=0;k<j+1;k++){
 429                                                         a[j][k] -= (f * e[k] + g * a[i][k]);
 430                                                 }
 431                                         }
 432                                 }
 433                         }
 434                         else{
 435                                 e[i] = a[i][l];
 436                         }
 437
 438                         d[i] = h;
 439                 }
 440
 441                 d[0] = 0.0000;
 442                 e[0] = 0.0000;
 443
 444                 for(int i=0;i<n;i++){
 445                         int l = i;
 446                         if(d[i] != 0.0){
 447                                 for(int j=0;j<l;j++){
 448                                         g = 0.0000;
 449                                         for(int k=0;k<l;k++){
 450                                                 g += a[i][k] * a[k][j];
 451                                         }
 452                                         for(int k=0;k<l;k++){
 453                                                 a[k][j] -= g * a[k][i];
 454                                         }
 455                                 }
 456                         }
 457                         d[i] = a[i][i];
 458                         a[i][i] = 1.0000;
 459                         for(int j=0;j<l;j++){
 460                                 a[j][i] = a[i][j] = 0.0;
 461                         }
 462                 }
 463         }
 464         catch(exception& e) {
 465                 m->errorOut(e, "PCACommand", "tred2");
 466                 exit(1);
 467         }
 468
 469 }
 470
 471 /*********************************************************************************************************************************/
 472
 473 //  This function is taken from Numerical Recipes in C++ by Press et al., 2nd edition, pg. 479
 474
 475 void PCACommand::qtli(vector<double>& d, vector<double>& e, vector<vector<double> >& z) {
 476         try {
 477                 int m, i, iter;
 478                 double s, r, p, g, f, dd, c, b;
 479
 480                 int n = d.size();
 481                 for(int i=1;i<=n;i++){
 482                         e[i-1] = e[i];
 483                 }
 484                 e[n-1] = 0.0000;
 485
 486                 for(int l=0;l<n;l++){
 487                         iter = 0;
 488                         do {
 489                                 for(m=l;m<n-1;m++){
 490                                         dd = fabs(d[m]) + fabs(d[m+1]);
 491                                         if(fabs(e[m])+dd == dd) break;
 492                                 }
 493                                 if(m != l){
 494                                         if(iter++ == 30) cerr << "Too many iterations in tqli\n";
 495                                         g = (d[l+1]-d[l]) / (2.0 * e[l]);
 496                                         r = pythag(g, 1.0);
 497                                         g = d[m] - d[l] + e[l] / (g + SIGN(r,g));
 498                                         s = c = 1.0;
 499                                         p = 0.0000;
 500                                         for(i=m-1;i>=l;i--){
 501                                                 f = s * e[i];
 502                                                 b = c * e[i];
 503                                                 e[i+1] = (r=pythag(f,g));
 504                                                 if(r==0.0){
 505                                                         d[i+1] -= p;
 506                                                         e[m] = 0.0000;
 507                                                         break;
 508                                                 }
 509                                                 s = f / r;
 510                                                 c = g / r;
 511                                                 g = d[i+1] - p;
 512                                                 r = (d[i] - g) * s + 2.0 * c * b;
 513                                                 d[i+1] = g + ( p = s * r);
 514                                                 g = c * r - b;
 515                                                 for(int k=0;k<n;k++){
 516                                                         f = z[k][i+1];
 517                                                         z[k][i+1] = s * z[k][i] + c * f;
 518                                                         z[k][i] = c * z[k][i] - s * f;
 519                                                 }
 520                                         }
 521                                         if(r == 0.00 && i >= l) continue;
 522                                         d[l] -= p;
 523                                         e[l] = g;
 524                                         e[m] = 0.0;
 525                                 }
 526                         } while (m != l);
 527                 }
 528
 529                 int k;
 530                 for(int i=0;i<n;i++){
 531                         p=d[k=i];
 532                         for(int j=i;j<n;j++){
 533                                 if(d[j] >= p){
 534                                         p=d[k=j];
 535                                 }
 536                         }
 537                         if(k!=i){
 538                                 d[k]=d[i];
 539                                 d[i]=p;
 540                                 for(int j=0;j<n;j++){
 541                                         p=z[j][i];
 542                                         z[j][i] = z[j][k];
 543                                         z[j][k] = p;
 544                                 }
 545                         }
 546                 }
 547         }
 548         catch(exception& e) {
 549                 m->errorOut(e, "PCACommand", "qtli");
 550                 exit(1);
 551         }
 552 }
 553
 554 /*********************************************************************************************************************************/
 555
 556 void PCACommand::output(string fnameRoot, vector<string> name_list, vector<vector<double> > G, vector<double> d) {
 557         try {
 558                 int rank = name_list.size();
 559                 double dsum = 0.0000;
 560                 for(int i=0;i<rank;i++){
 561                         dsum += d[i];
 562                         for(int j=0;j<rank;j++){
 563                                 if(d[j] >= 0)   {       G[i][j] *= pow(d[j],0.5);       }
 564                                 else                    {       G[i][j] = 0.00000;                      }
 565                         }
 566                 }
 567
 568                 ofstream pcaData((fnameRoot+"pcoa").c_str(), ios::trunc);
 569                 pcaData.setf(ios::fixed, ios::floatfield);
 570                 pcaData.setf(ios::showpoint);
 571                 outputNames.push_back(fnameRoot+"pcoa");
 572                 outputTypes["pcoa"].push_back(fnameRoot+"pcoa");
 573
 574                 ofstream pcaLoadings((fnameRoot+"pcoa.loadings").c_str(), ios::trunc);
 575                 pcaLoadings.setf(ios::fixed, ios::floatfield);
 576                 pcaLoadings.setf(ios::showpoint);
 577                 outputNames.push_back(fnameRoot+"pcoa.loadings");
 578                 outputTypes["loadings"].push_back(fnameRoot+"pcoa.loadings");
 579
 580                 pcaLoadings << "axis\tloading\n";
 581                 for(int i=0;i<rank;i++){
 582                         pcaLoadings << i+1 << '\t' << d[i] * 100.0 / dsum << endl;
 583                 }
 584
 585                 pcaData << "group";
 586                 for(int i=0;i<rank;i++){
 587                         pcaData << '\t' << "axis" << i+1;
 588                 }
 589                 pcaData << endl;
 590
 591                 for(int i=0;i<rank;i++){
 592                         pcaData << name_list[i] << '\t';
 593                         for(int j=0;j<rank;j++){
 594                                 pcaData << G[i][j] << '\t';
 595                         }
 596                         pcaData << endl;
 597                 }
 598         }
 599         catch(exception& e) {
 600                 m->errorOut(e, "PCACommand", "output");
 601                 exit(1);
 602         }
 603 }
 604
 605 /*********************************************************************************************************************************/
 606