phylotree.cpp

   1 /*
   2  *  doTaxonomy.cpp
   3  *
   4  *
   5  *  Created by Pat Schloss on 6/17/09.
   6  *  Copyright 2009 Patrick D. Schloss. All rights reserved.
   7  *
   8  */
   9
  10 #include "phylotree.h"
  11
  12 /**************************************************************************************************/
  13
  14 PhyloTree::PhyloTree(){
  15         try {
  16                 m = MothurOut::getInstance();
  17                 numNodes = 1;
  18                 numSeqs = 0;
  19                 tree.push_back(TaxNode("Root"));
  20                 tree[0].heirarchyID = "0";
  21                 maxLevel = 0;
  22                 calcTotals = true;
  23         }
  24         catch(exception& e) {
  25                 m->errorOut(e, "PhyloTree", "PhyloTree");
  26                 exit(1);
  27         }
  28 }
  29 /**************************************************************************************************/
  30
  31 PhyloTree::PhyloTree(ifstream& in, string filename){
  32         try {
  33                 m = MothurOut::getInstance();
  34                 calcTotals = false;
  35
  36                 #ifdef USE_MPI
  37                         MPI_File inMPI;
  38                         MPI_Offset size;
  39                         MPI_Status status;
  40
  41                         char inFileName[1024];
  42                         strcpy(inFileName, filename.c_str());
  43
  44                         MPI_File_open(MPI_COMM_WORLD, inFileName, MPI_MODE_RDONLY, MPI_INFO_NULL, &inMPI);
  45                         MPI_File_get_size(inMPI, &size);
  46
  47                         char* buffer = new char[size];
  48                         MPI_File_read(inMPI, buffer, size, MPI_CHAR, &status);
  49
  50                         string tempBuf = buffer;
  51                         if (tempBuf.length() > size) { tempBuf = tempBuf.substr(0, size);  }
  52                         istringstream iss (tempBuf,istringstream::in);
  53                         delete buffer;
  54
  55                         iss >> numNodes; gobble(iss);
  56
  57                         tree.resize(numNodes);
  58
  59                         for (int i = 0; i < tree.size(); i++) {
  60                                 iss >> tree[i].name >> tree[i].level >> tree[i].parent; gobble(iss);
  61                         }
  62
  63                         //read genus nodes
  64                         int numGenus = 0;
  65                         iss >> numGenus; gobble(iss);
  66
  67                         int gnode, gsize;
  68                         totals.clear();
  69                         for (int i = 0; i < numGenus; i++) {
  70                                 iss >> gnode >> gsize; gobble(iss);
  71
  72                                 uniqueTaxonomies[gnode] = gnode;
  73                                 totals.push_back(gsize);
  74                         }
  75
  76                         MPI_File_close(&inMPI);
  77
  78                 #else
  79                         in >> numNodes; gobble(in);
  80
  81                         tree.resize(numNodes);
  82
  83                         for (int i = 0; i < tree.size(); i++) {
  84                                 in >> tree[i].name >> tree[i].level >> tree[i].parent; gobble(in);
  85                         }
  86
  87                         //read genus nodes
  88                         int numGenus = 0;
  89                         in >> numGenus; gobble(in);
  90
  91                         int gnode, gsize;
  92                         totals.clear();
  93                         for (int i = 0; i < numGenus; i++) {
  94                                 in >> gnode >> gsize; gobble(in);
  95
  96                                 uniqueTaxonomies[gnode] = gnode;
  97                                 totals.push_back(gsize);
  98                         }
  99
 100                         in.close();
 101
 102                 #endif
 103
 104         }
 105         catch(exception& e) {
 106                 m->errorOut(e, "PhyloTree", "PhyloTree");
 107                 exit(1);
 108         }
 109 }
 110 /**************************************************************************************************/
 111
 112 PhyloTree::PhyloTree(string tfile){
 113         try {
 114                 m = MothurOut::getInstance();
 115                 numNodes = 1;
 116                 numSeqs = 0;
 117                 tree.push_back(TaxNode("Root"));
 118                 tree[0].heirarchyID = "0";
 119                 maxLevel = 0;
 120                 calcTotals = true;
 121                 string name, tax;
 122
 123
 124                 #ifdef USE_MPI
 125                         int pid, num;
 126                         vector<long> positions;
 127
 128                         MPI_Status status;
 129                         MPI_File inMPI;
 130                         MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
 131
 132                         char inFileName[1024];
 133                         strcpy(inFileName, tfile.c_str());
 134
 135                         MPI_File_open(MPI_COMM_WORLD, inFileName, MPI_MODE_RDONLY, MPI_INFO_NULL, &inMPI);  //comm, filename, mode, info, filepointer
 136
 137                         if (pid == 0) {
 138                                 positions = setFilePosEachLine(tfile, num);
 139
 140                                 //send file positions to all processes
 141                                 MPI_Bcast(&num, 1, MPI_INT, 0, MPI_COMM_WORLD);  //send numSeqs
 142                                 MPI_Bcast(&positions[0], (num+1), MPI_LONG, 0, MPI_COMM_WORLD); //send file pos
 143                         }else{
 144                                 MPI_Bcast(&num, 1, MPI_INT, 0, MPI_COMM_WORLD); //get numSeqs
 145                                 positions.resize(num);
 146                                 MPI_Bcast(&positions[0], (num+1), MPI_LONG, 0, MPI_COMM_WORLD); //get file positions
 147                         }
 148
 149                         //read file
 150                         for(int i=0;i<num;i++){
 151                                 //read next sequence
 152                                 int length = positions[i+1] - positions[i];
 153                                 char* buf4 = new char[length];
 154
 155                                 MPI_File_read_at(inMPI, positions[i], buf4, length, MPI_CHAR, &status);
 156
 157                                 string tempBuf = buf4;
 158                                 if (tempBuf.length() > length) { tempBuf = tempBuf.substr(0, length); }
 159                                 delete buf4;
 160
 161                                 istringstream iss (tempBuf,istringstream::in);
 162                                 iss >> name >> tax;
 163                                 addSeqToTree(name, tax);
 164                         }
 165
 166                         MPI_File_close(&inMPI);
 167
 168                 #else
 169                         ifstream in;
 170                         openInputFile(tfile, in);
 171
 172                         //read in users taxonomy file and add sequences to tree
 173                         while(!in.eof()){
 174                                 in >> name >> tax; gobble(in);
 175
 176                                 addSeqToTree(name, tax);
 177                         }
 178                         in.close();
 179                 #endif
 180
 181                 assignHeirarchyIDs(0);
 182
 183                 //create file for summary if needed
 184                 setUp(tfile);
 185         }
 186         catch(exception& e) {
 187                 m->errorOut(e, "PhyloTree", "PhyloTree");
 188                 exit(1);
 189         }
 190 }
 191
 192 /**************************************************************************************************/
 193
 194 string PhyloTree::getNextTaxon(string& heirarchy){
 195         try {
 196                 string currentLevel = "";
 197                 if(heirarchy != ""){
 198                         int pos = heirarchy.find_first_of(';');
 199                         currentLevel=heirarchy.substr(0,pos);
 200                         if (pos != (heirarchy.length()-1)) {  heirarchy=heirarchy.substr(pos+1);  }
 201                         else { heirarchy = ""; }
 202                 }
 203                 return currentLevel;
 204         }
 205         catch(exception& e) {
 206                 m->errorOut(e, "PhyloTree", "getNextTaxon");
 207                 exit(1);
 208         }
 209 }
 210
 211 /**************************************************************************************************/
 212
 213 int PhyloTree::addSeqToTree(string seqName, string seqTaxonomy){
 214         try {
 215                 numSeqs++;
 216
 217                 map<string, int>::iterator childPointer;
 218
 219                 int currentNode = 0;
 220                 int level = 1;
 221
 222                 tree[0].accessions.push_back(seqName);
 223                 string taxon;// = getNextTaxon(seqTaxonomy);
 224
 225                 while(seqTaxonomy != ""){
 226
 227                         level++;
 228
 229                         if (m->control_pressed) { return 0; }
 230
 231                         //somehow the parent is getting one too many accnos
 232                         //use print to reassign the taxa id
 233                         taxon = getNextTaxon(seqTaxonomy);
 234
 235                         childPointer = tree[currentNode].children.find(taxon);
 236
 237                         if(childPointer != tree[currentNode].children.end()){   //if the node already exists, move on
 238                                 currentNode = childPointer->second;
 239                                 tree[currentNode].accessions.push_back(seqName);
 240                                 name2Taxonomy[seqName] = currentNode;
 241                         }
 242                         else{                                                                                   //otherwise, create it
 243                                 tree.push_back(TaxNode(taxon));
 244                                 numNodes++;
 245                                 tree[currentNode].children[taxon] = numNodes-1;
 246                                 tree[numNodes-1].parent = currentNode;
 247
 248                                 //                      int numChildren = tree[currentNode].children.size();
 249                                 //                      string heirarchyID = tree[currentNode].heirarchyID;
 250                                 //                      tree[currentNode].accessions.push_back(seqName);
 251
 252                                 currentNode = tree[currentNode].children[taxon];
 253                                 tree[currentNode].accessions.push_back(seqName);
 254                                 name2Taxonomy[seqName] = currentNode;
 255                                 //                      tree[currentNode].level = level;
 256                                 //                      tree[currentNode].childNumber = numChildren;
 257                                 //                      tree[currentNode].heirarchyID = heirarchyID + '.' + toString(tree[currentNode].childNumber);
 258                         }
 259
 260                         if (seqTaxonomy == "") {   uniqueTaxonomies[currentNode] = currentNode; }
 261                 }
 262
 263         }
 264         catch(exception& e) {
 265                 m->errorOut(e, "PhyloTree", "addSeqToTree");
 266                 exit(1);
 267         }
 268 }
 269 /**************************************************************************************************/
 270 vector<int> PhyloTree::getGenusNodes()  {
 271         try {
 272                 genusIndex.clear();
 273                 //generate genusIndexes
 274                 map<int, int>::iterator it2;
 275                 for (it2=uniqueTaxonomies.begin(); it2!=uniqueTaxonomies.end(); it2++) {  genusIndex.push_back(it2->first);     }
 276
 277                 return genusIndex;
 278         }
 279         catch(exception& e) {
 280                 m->errorOut(e, "PhyloTree", "getGenusNodes");
 281                 exit(1);
 282         }
 283 }
 284 /**************************************************************************************************/
 285 vector<int> PhyloTree::getGenusTotals() {
 286         try {
 287
 288                 if (calcTotals) {
 289                         totals.clear();
 290                         //reset counts because we are on a new word
 291                         for (int j = 0; j < genusIndex.size(); j++) {
 292                                 totals.push_back(tree[genusIndex[j]].accessions.size());
 293                         }
 294                         return totals;
 295                 }else{
 296                         return totals;
 297                 }
 298
 299         }
 300         catch(exception& e) {
 301                 m->errorOut(e, "PhyloTree", "getGenusNodes");
 302                 exit(1);
 303         }
 304 }
 305 /**************************************************************************************************/
 306
 307 void PhyloTree::assignHeirarchyIDs(int index){
 308         try {
 309                 map<string,int>::iterator it;
 310                 int counter = 1;
 311
 312                 for(it=tree[index].children.begin();it!=tree[index].children.end();it++){
 313                         tree[it->second].heirarchyID = tree[index].heirarchyID + '.' + toString(counter);
 314                         counter++;
 315                         tree[it->second].level = tree[index].level + 1;
 316
 317                         //save maxLevel for binning the unclassified seqs
 318                         if (tree[it->second].level > maxLevel) { maxLevel = tree[it->second].level; }
 319
 320                         assignHeirarchyIDs(it->second);
 321                 }
 322         }
 323         catch(exception& e) {
 324                 m->errorOut(e, "PhyloTree", "assignHeirarchyIDs");
 325                 exit(1);
 326         }
 327 }
 328 /**************************************************************************************************/
 329 void PhyloTree::setUp(string tfile){
 330         try{
 331                 string taxFileNameTest = tfile.substr(0,tfile.find_last_of(".")+1) + "tree.sum";
 332
 333                 #ifdef USE_MPI
 334                         int pid;
 335                         MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
 336
 337                         if (pid == 0) {  binUnclassified(taxFileNameTest);  }
 338
 339                 #else
 340                         //create file needed for summary if it doesn't exist
 341                         ifstream FileTest(taxFileNameTest.c_str());
 342
 343                         if (!FileTest) {
 344                                 binUnclassified(taxFileNameTest);
 345                         }
 346                 #endif
 347         }
 348         catch(exception& e) {
 349                 m->errorOut(e, "PhyloTree", "setUp");
 350                 exit(1);
 351         }
 352 }
 353 /**************************************************************************************************/
 354 void PhyloTree::binUnclassified(string file){
 355         try {
 356
 357                 ofstream out;
 358                 openOutputFile(file, out);
 359
 360                 map<string, int>::iterator itBin;
 361                 map<string, int>::iterator childPointer;
 362
 363                 vector<TaxNode> copy = tree;
 364                 int copyNodes = numNodes;
 365
 366                 //go through the seqs and if a sequence finest taxon is not the same level as the most finely defined taxon then classify it as unclassified where necessary
 367                 for (itBin = name2Taxonomy.begin(); itBin != name2Taxonomy.end(); itBin++) {
 368
 369                         if (m->control_pressed) {  out.close(); break;  }
 370
 371                         int level = copy[itBin->second].level;
 372                         int currentNode = itBin->second;
 373
 374                         //this sequence is unclassified at some levels
 375                         while(level != maxLevel){
 376
 377                                 level++;
 378
 379                                 string taxon = "unclassified";
 380
 381                                 //does the parent have a child names 'unclassified'?
 382                                 childPointer = copy[currentNode].children.find(taxon);
 383
 384                                 if(childPointer != copy[currentNode].children.end()){   //if the node already exists, move on
 385                                         currentNode = childPointer->second; //currentNode becomes 'unclassified'
 386                                         copy[currentNode].accessions.push_back(itBin->first);  //add this seq
 387                                 }
 388                                 else{                                                                                   //otherwise, create it
 389                                         copy.push_back(TaxNode(taxon));
 390                                         copyNodes++;
 391                                         copy[currentNode].children[taxon] = copyNodes-1;
 392                                         copy[copyNodes-1].parent = currentNode;
 393                                         copy[copyNodes-1].level = copy[currentNode].level + 1;
 394
 395                                         currentNode = copy[currentNode].children[taxon];
 396                                         copy[currentNode].accessions.push_back(itBin->first);
 397                                 }
 398                         }
 399                 }
 400
 401                 if (!m->control_pressed) {
 402                         //print copy tree
 403                         print(out, copy);
 404                 }
 405
 406         }
 407         catch(exception& e) {
 408                 m->errorOut(e, "PhyloTree", "binUnclassified");
 409                 exit(1);
 410         }
 411 }
 412 /**************************************************************************************************/
 413 string PhyloTree::getFullTaxonomy(string seqName) {
 414         try {
 415                 string tax = "";
 416
 417                 int currentNode = name2Taxonomy[seqName];
 418
 419                 while (tree[currentNode].parent != -1) {
 420                         tax = tree[currentNode].name + ";" + tax;
 421                         currentNode = tree[currentNode].parent;
 422                 }
 423
 424                 return tax;
 425         }
 426         catch(exception& e) {
 427                 m->errorOut(e, "PhyloTree", "getFullTaxonomy");
 428                 exit(1);
 429         }
 430 }
 431 /**************************************************************************************************/
 432
 433 void PhyloTree::print(ofstream& out, vector<TaxNode>& copy){
 434         try {
 435                 out << copy.size() << endl;
 436
 437                 for (int i = 0; i < copy.size(); i++) {
 438
 439                         out << copy[i].level << '\t'<< copy[i].name << '\t' << copy[i].children.size() << '\t';
 440
 441                         map<string,int>::iterator it;
 442                         for(it=copy[i].children.begin();it!=copy[i].children.end();it++){
 443                                 out << it->first << '\t' << it->second << '\t';
 444                         }
 445                         out << endl;
 446                 }
 447
 448                 out.close();
 449         }
 450         catch(exception& e) {
 451                 m->errorOut(e, "PhyloTree", "print");
 452                 exit(1);
 453         }
 454 }
 455 /**************************************************************************************************/
 456 void PhyloTree::printTreeNodes(string treefilename) {
 457         try {
 458
 459                 #ifdef USE_MPI
 460                         int pid;
 461                         MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
 462
 463                         if (pid == 0) {
 464
 465                 #endif
 466
 467                         ofstream outTree;
 468                         openOutputFile(treefilename, outTree);
 469
 470                         //print treenodes
 471                         outTree << tree.size() << endl;
 472                         for (int i = 0; i < tree.size(); i++) {
 473                                 outTree << tree[i].name << '\t' << tree[i].level << '\t' << tree[i].parent << endl;
 474                         }
 475
 476                         //print genus nodes
 477                         outTree << endl << uniqueTaxonomies.size() << endl;
 478                         map<int, int>::iterator it2;
 479                         for (it2=uniqueTaxonomies.begin(); it2!=uniqueTaxonomies.end(); it2++) {  outTree << it2->first << '\t' << tree[it2->first].accessions.size() << endl;  }
 480                         outTree << endl;
 481
 482
 483                         outTree.close();
 484
 485                 #ifdef USE_MPI
 486                         }
 487                 #endif
 488
 489
 490         }
 491         catch(exception& e) {
 492                 m->errorOut(e, "PhyloTree", "printTreeNodes");
 493                 exit(1);
 494         }
 495 }
 496 /**************************************************************************************************/
 497
 498
 499