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