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