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