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