phylosummary.cpp

   1 /*
   2  *  rawTrainingDataMaker.cpp
   3  *  Mothur
   4  *
   5  *  Created by westcott on 4/21/10.
   6  *  Copyright 2010 Schloss Lab. All rights reserved.
   7  *
   8  */
   9
  10 #include "phylosummary.h"
  11 #include "referencedb.h"
  12 /**************************************************************************************************/
  13
  14 PhyloSummary::PhyloSummary(string refTfile, CountTable* c, bool r){
  15         try {
  16                 m = MothurOut::getInstance();
  17                 maxLevel = 0;
  18                 ignore = false;
  19         numSeqs = 0;
  20         relabund = r;
  21
  22                 ct = c;
  23         groupmap = NULL;
  24
  25                 //check for necessary files
  26         if (refTfile == "saved") { ReferenceDB* rdb = ReferenceDB::getInstance(); refTfile = rdb->getSavedTaxonomy(); }
  27                 string taxFileNameTest = m->getFullPathName((refTfile.substr(0,refTfile.find_last_of(".")+1) + "tree.sum"));
  28                 ifstream FileTest(taxFileNameTest.c_str());
  29
  30                 if (!FileTest) {
  31                         m->mothurOut("Error: can't find " + taxFileNameTest + "."); m->mothurOutEndLine(); exit(1);
  32                 }else{
  33                         readTreeStruct(FileTest);
  34                 }
  35
  36                 tree[0].rank = "0";
  37                 assignRank(0);
  38
  39         }
  40         catch(exception& e) {
  41                 m->errorOut(e, "PhyloSummary", "PhyloSummary");
  42                 exit(1);
  43         }
  44 }
  45
  46 /**************************************************************************************************/
  47
  48 PhyloSummary::PhyloSummary(CountTable* c, bool r){
  49         try {
  50                 m = MothurOut::getInstance();
  51                 maxLevel = 0;
  52                 ignore = true;
  53         numSeqs = 0;
  54         relabund = r;
  55
  56                 ct = c;
  57         groupmap = NULL;
  58
  59                 tree.push_back(rawTaxNode("Root"));
  60                 tree[0].rank = "0";
  61         }
  62         catch(exception& e) {
  63                 m->errorOut(e, "PhyloSummary", "PhyloSummary");
  64                 exit(1);
  65         }
  66 }
  67 /**************************************************************************************************/
  68 PhyloSummary::PhyloSummary(string refTfile, GroupMap* g, bool r){
  69         try {
  70                 m = MothurOut::getInstance();
  71                 maxLevel = 0;
  72                 ignore = false;
  73         numSeqs = 0;
  74         relabund = r;
  75
  76                 groupmap = g;
  77         ct = NULL;
  78
  79                 //check for necessary files
  80         if (refTfile == "saved") { ReferenceDB* rdb = ReferenceDB::getInstance(); refTfile = rdb->getSavedTaxonomy(); }
  81                 string taxFileNameTest = m->getFullPathName((refTfile.substr(0,refTfile.find_last_of(".")+1) + "tree.sum"));
  82                 ifstream FileTest(taxFileNameTest.c_str());
  83
  84                 if (!FileTest) {
  85                         m->mothurOut("Error: can't find " + taxFileNameTest + "."); m->mothurOutEndLine(); exit(1);
  86                 }else{
  87                         readTreeStruct(FileTest);
  88                 }
  89
  90                 tree[0].rank = "0";
  91                 assignRank(0);
  92
  93         }
  94         catch(exception& e) {
  95                 m->errorOut(e, "PhyloSummary", "PhyloSummary");
  96                 exit(1);
  97         }
  98 }
  99
 100 /**************************************************************************************************/
 101
 102 PhyloSummary::PhyloSummary(GroupMap* g, bool r){
 103         try {
 104                 m = MothurOut::getInstance();
 105                 maxLevel = 0;
 106                 ignore = true;
 107         numSeqs = 0;
 108         relabund = r;
 109
 110                 groupmap = g;
 111         ct = NULL;
 112
 113                 tree.push_back(rawTaxNode("Root"));
 114                 tree[0].rank = "0";
 115         }
 116         catch(exception& e) {
 117                 m->errorOut(e, "PhyloSummary", "PhyloSummary");
 118                 exit(1);
 119         }
 120 }
 121 /**************************************************************************************************/
 122
 123 int PhyloSummary::summarize(string userTfile){
 124         try {
 125                 map<string, string> temp;
 126         m->readTax(userTfile, temp);
 127
 128         for (map<string, string>::iterator itTemp = temp.begin(); itTemp != temp.end();) {
 129             addSeqToTree(itTemp->first, itTemp->second);
 130             temp.erase(itTemp++);
 131         }
 132
 133         return numSeqs;
 134         }
 135         catch(exception& e) {
 136                 m->errorOut(e, "PhyloSummary", "summarize");
 137                 exit(1);
 138         }
 139 }
 140
 141 /**************************************************************************************************/
 142
 143 string PhyloSummary::getNextTaxon(string& heirarchy){
 144         try {
 145                 string currentLevel = "";
 146                 if(heirarchy != ""){
 147                         int pos = heirarchy.find_first_of(';');
 148                         currentLevel=heirarchy.substr(0,pos);
 149                         if (pos != (heirarchy.length()-1)) {  heirarchy=heirarchy.substr(pos+1);  }
 150                         else { heirarchy = ""; }
 151                 }
 152                 return currentLevel;
 153         }
 154         catch(exception& e) {
 155                 m->errorOut(e, "PhyloSummary", "getNextTaxon");
 156                 exit(1);
 157         }
 158 }
 159
 160 /**************************************************************************************************/
 161
 162 int PhyloSummary::addSeqToTree(string seqName, string seqTaxonomy){
 163         try {
 164
 165                 numSeqs++;
 166
 167                 map<string, int>::iterator childPointer;
 168
 169                 int currentNode = 0;
 170                 string taxon;
 171
 172                 int level = 0;
 173
 174                 //are there confidence scores, if so remove them
 175                 if (seqTaxonomy.find_first_of('(') != -1) {  m->removeConfidences(seqTaxonomy); }
 176
 177                 while (seqTaxonomy != "") {
 178
 179                         if (m->control_pressed) { return 0; }
 180
 181                         //somehow the parent is getting one too many accnos
 182                         //use print to reassign the taxa id
 183                         taxon = getNextTaxon(seqTaxonomy);
 184
 185                         childPointer = tree[currentNode].children.find(taxon);
 186
 187                         if(childPointer != tree[currentNode].children.end()){   //if the node already exists, update count and move on
 188                                 int thisCount = 1;
 189
 190                 if (groupmap != NULL) {
 191                                         //find out the sequences group
 192                                         string group = groupmap->getGroup(seqName);
 193
 194                                         if (group == "not found") {  m->mothurOut("[WARNING]: " + seqName + " is not in your groupfile, and will be included in the overall total, but not any group total."); m->mothurOutEndLine();  }
 195
 196                                         //do you have a count for this group?
 197                                         map<string, int>::iterator itGroup = tree[childPointer->second].groupCount.find(group);
 198
 199                                         //if yes, increment it - there should not be a case where we can't find it since we load group in read
 200                                         if (itGroup != tree[childPointer->second].groupCount.end()) {
 201                                                 tree[childPointer->second].groupCount[group]++;
 202                                         }
 203                                 }else if (ct != NULL) {
 204                     if (ct->hasGroupInfo()) {
 205                         vector<int> groupCounts = ct->getGroupCounts(seqName);
 206                         vector<string> groups = ct->getNamesOfGroups();
 207                         for (int i = 0; i < groups.size(); i++) {
 208
 209                             if (groupCounts[i] != 0) {
 210                                 //do you have a count for this group?
 211                                 map<string, int>::iterator itGroup = tree[childPointer->second].groupCount.find(groups[i]);
 212
 213                                 //if yes, increment it - there should not be a case where we can't find it since we load group in read
 214                                 if (itGroup != tree[childPointer->second].groupCount.end()) {
 215                                     tree[childPointer->second].groupCount[groups[i]] += groupCounts[i];
 216                                 }
 217                             }
 218                         }
 219                     }
 220                     thisCount = ct->getNumSeqs(seqName);
 221                 }
 222
 223                                 tree[childPointer->second].total += thisCount;
 224
 225                                 currentNode = childPointer->second;
 226                         }else{
 227                                 if (ignore) {
 228
 229                                         tree.push_back(rawTaxNode(taxon));
 230                                         int index = tree.size() - 1;
 231
 232                                         tree[index].parent = currentNode;
 233                                         tree[index].level = (level+1);
 234                                         tree[currentNode].children[taxon] = index;
 235                     int thisCount = 1;
 236
 237                                         //initialize groupcounts
 238                                         if (groupmap != NULL) {
 239                                                 vector<string> mGroups = groupmap->getNamesOfGroups();
 240                                                 for (int j = 0; j < mGroups.size(); j++) {
 241                                                         tree[index].groupCount[mGroups[j]] = 0;
 242                                                 }
 243
 244                                                 //find out the sequences group
 245                                                 string group = groupmap->getGroup(seqName);
 246
 247                                                 if (group == "not found") {  m->mothurOut("[WARNING]: " + seqName + " is not in your groupfile, and will be included in the overall total, but not any group total."); m->mothurOutEndLine();  }
 248
 249                                                 //do you have a count for this group?
 250                                                 map<string, int>::iterator itGroup = tree[index].groupCount.find(group);
 251
 252                                                 //if yes, increment it - there should not be a case where we can't find it since we load group in read
 253                                                 if (itGroup != tree[index].groupCount.end()) {
 254                                                         tree[index].groupCount[group]++;
 255                                                 }
 256                                         }else if (ct != NULL) {
 257                         if (ct->hasGroupInfo()) {
 258                             vector<string> mGroups = ct->getNamesOfGroups();
 259                             for (int j = 0; j < mGroups.size(); j++) {
 260                                 tree[index].groupCount[mGroups[j]] = 0;
 261                             }
 262                             vector<int> groupCounts = ct->getGroupCounts(seqName);
 263                             vector<string> groups = ct->getNamesOfGroups();
 264
 265                             for (int i = 0; i < groups.size(); i++) {
 266                                 if (groupCounts[i] != 0) {
 267
 268                                     //do you have a count for this group?
 269                                     map<string, int>::iterator itGroup = tree[index].groupCount.find(groups[i]);
 270
 271                                     //if yes, increment it - there should not be a case where we can't find it since we load group in read
 272                                     if (itGroup != tree[index].groupCount.end()) {
 273                                         tree[index].groupCount[groups[i]]+=groupCounts[i];
 274                                     }
 275                                 }
 276                             }
 277                         }
 278                         thisCount = ct->getNumSeqs(seqName);
 279                     }
 280
 281                     tree[index].total = thisCount;
 282                                         currentNode = index;
 283
 284                                 }else{ //otherwise, error
 285                                         m->mothurOut("Warning: cannot find taxon " + taxon + " in reference taxonomy tree at level " + toString(tree[currentNode].level) + " for " + seqName + ". This may cause totals of daughter levels not to add up in summary file."); m->mothurOutEndLine();
 286                                         break;
 287                                 }
 288                         }
 289
 290                         level++;
 291
 292                         if ((seqTaxonomy == "") && (level < maxLevel)) {  //if you think you are done and you are not.
 293                                 for (int k = level; k < maxLevel; k++) {  seqTaxonomy += "unclassified;";   }
 294                         }
 295                 }
 296                 return 0;
 297         }
 298         catch(exception& e) {
 299                 m->errorOut(e, "PhyloSummary", "addSeqToTree");
 300                 exit(1);
 301         }
 302 }
 303 /**************************************************************************************************/
 304
 305 int PhyloSummary::addSeqToTree(string seqTaxonomy, map<string, bool> containsGroup){
 306         try {
 307                 numSeqs++;
 308
 309                 map<string, int>::iterator childPointer;
 310
 311                 int currentNode = 0;
 312                 string taxon;
 313
 314                 int level = 0;
 315
 316                 //are there confidence scores, if so remove them
 317                 if (seqTaxonomy.find_first_of('(') != -1) {  m->removeConfidences(seqTaxonomy); }
 318
 319                 while (seqTaxonomy != "") {
 320
 321                         if (m->control_pressed) { return 0; }
 322
 323                         //somehow the parent is getting one too many accnos
 324                         //use print to reassign the taxa id
 325                         taxon = getNextTaxon(seqTaxonomy);
 326
 327                         childPointer = tree[currentNode].children.find(taxon);
 328
 329                         if(childPointer != tree[currentNode].children.end()){   //if the node already exists, update count and move on
 330                 for (map<string, bool>::iterator itGroup = containsGroup.begin(); itGroup != containsGroup.end(); itGroup++) {
 331                     if (itGroup->second == true) {
 332                         tree[childPointer->second].groupCount[itGroup->first]++;
 333                     }
 334                 }
 335
 336                                 tree[childPointer->second].total++;
 337
 338                                 currentNode = childPointer->second;
 339                         }else{
 340                                 if (ignore) {
 341
 342                                         tree.push_back(rawTaxNode(taxon));
 343                                         int index = tree.size() - 1;
 344
 345                                         tree[index].parent = currentNode;
 346                                         tree[index].level = (level+1);
 347                                         tree[index].total = 1;
 348                                         tree[currentNode].children[taxon] = index;
 349
 350                     for (map<string, bool>::iterator itGroup = containsGroup.begin(); itGroup != containsGroup.end(); itGroup++) {
 351                         if (itGroup->second == true) {
 352                             tree[index].groupCount[itGroup->first]++;
 353                         }
 354                     }
 355
 356                                         currentNode = index;
 357
 358                                 }else{ //otherwise, error
 359                                         m->mothurOut("Warning: cannot find taxon " + taxon + " in reference taxonomy tree at level " + toString(tree[currentNode].level) + ". This may cause totals of daughter levels not to add up in summary file."); m->mothurOutEndLine();
 360                                         break;
 361                                 }
 362                         }
 363
 364                         level++;
 365
 366                         if ((seqTaxonomy == "") && (level < maxLevel)) {  //if you think you are done and you are not.
 367                                 for (int k = level; k < maxLevel; k++) {  seqTaxonomy += "unclassified;";   }
 368                         }
 369                 }
 370                 return 0;
 371         }
 372         catch(exception& e) {
 373                 m->errorOut(e, "PhyloSummary", "addSeqToTree");
 374                 exit(1);
 375         }
 376 }
 377
 378 /**************************************************************************************************/
 379
 380 void PhyloSummary::assignRank(int index){
 381         try {
 382                 map<string,int>::iterator it;
 383                 int counter = 1;
 384
 385                 for(it=tree[index].children.begin();it!=tree[index].children.end();it++){
 386                         tree[it->second].rank = tree[index].rank + '.' + toString(counter);
 387                         counter++;
 388                         assignRank(it->second);
 389                 }
 390         }
 391         catch(exception& e) {
 392                 m->errorOut(e, "PhyloSummary", "assignRank");
 393                 exit(1);
 394         }
 395 }
 396 /**************************************************************************************************/
 397
 398 void PhyloSummary::print(ofstream& out){
 399         try {
 400
 401                 if (ignore)     {  assignRank(0); }
 402         vector<string> mGroups;
 403                 //print labels
 404                 out << "taxlevel\t rankID\t taxon\t daughterlevels\t total\t";
 405                 if (groupmap != NULL) {
 406                         //so the labels match the counts below, since the map sorts them automatically...
 407                         //sort(groupmap->namesOfGroups.begin(), groupmap->namesOfGroups.end());
 408             mGroups = groupmap->getNamesOfGroups();
 409                         for (int i = 0; i < mGroups.size(); i++) {
 410                                 out << mGroups[i] << '\t';
 411                         }
 412                 }else if (ct != NULL) {
 413             if (ct->hasGroupInfo()) {
 414                 mGroups = ct->getNamesOfGroups();
 415                 for (int i = 0; i < mGroups.size(); i++) {
 416                     out << mGroups[i] << '\t';
 417                 }
 418             }
 419         }
 420                 out << endl;
 421
 422                 int totalChildrenInTree = 0;
 423                 map<string, int>::iterator itGroup;
 424                 map<string,int>::iterator it;
 425                 for(it=tree[0].children.begin();it!=tree[0].children.end();it++){
 426                         if (tree[it->second].total != 0)  {
 427                                 totalChildrenInTree++;
 428                                 tree[0].total += tree[it->second].total;
 429
 430                                 if (groupmap != NULL) {
 431                                         for (int i = 0; i < mGroups.size(); i++) { tree[0].groupCount[mGroups[i]] += tree[it->second].groupCount[mGroups[i]]; }
 432                                 }else if ( ct != NULL) {
 433                     if (ct->hasGroupInfo()) { for (int i = 0; i < mGroups.size(); i++) { tree[0].groupCount[mGroups[i]] += tree[it->second].groupCount[mGroups[i]]; } }
 434                 }
 435                         }
 436                 }
 437
 438                 //print root
 439         if (relabund) {
 440             out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
 441             out << tree[0].level << "\t" << tree[0].rank << "\t" << tree[0].name << "\t" << totalChildrenInTree << "\t" << (tree[0].total/(double) tree[0].total) << '\t';
 442
 443             if (groupmap != NULL) {
 444                 for (int i = 0; i < mGroups.size(); i++) {
 445                     double thisNum = tree[0].groupCount[mGroups[i]];
 446                     thisNum /= (double) groupmap->getNumSeqs(mGroups[i]);
 447                     out << thisNum << '\t';
 448                 }
 449             }else if ( ct != NULL) {
 450                 if (ct->hasGroupInfo()) {
 451                     for (int i = 0; i < mGroups.size(); i++) {
 452                         double thisNum = tree[0].groupCount[mGroups[i]];
 453                         thisNum /= (double) ct->getGroupCount(mGroups[i]);
 454                         out << thisNum << '\t';
 455                     }
 456                 }
 457             }
 458             out << endl;
 459         }else {
 460             out << tree[0].level << "\t" << tree[0].rank << "\t" << tree[0].name << "\t" << totalChildrenInTree << "\t" << tree[0].total << "\t";
 461             if (groupmap != NULL) {
 462                 for (int i = 0; i < mGroups.size(); i++) {  out << tree[0].groupCount[mGroups[i]] << '\t'; }
 463             }else if ( ct != NULL) {
 464                 if (ct->hasGroupInfo()) { for (int i = 0; i < mGroups.size(); i++) {  out << tree[0].groupCount[mGroups[i]] << '\t'; } }
 465             }
 466             out << endl;
 467                 }
 468
 469                 //print rest
 470                 print(0, out);
 471
 472         }
 473         catch(exception& e) {
 474                 m->errorOut(e, "PhyloSummary", "print");
 475                 exit(1);
 476         }
 477 }
 478 /**************************************************************************************************/
 479
 480 void PhyloSummary::print(ofstream& out, bool relabund){
 481         try {
 482
 483                 if (ignore) { assignRank(0); }
 484
 485                 int totalChildrenInTree = 0;
 486                 map<string, int>::iterator itGroup;
 487
 488                 map<string,int>::iterator it;
 489                 for(it=tree[0].children.begin();it!=tree[0].children.end();it++){
 490                         if (tree[it->second].total != 0)  {
 491                                 totalChildrenInTree++;
 492                                 tree[0].total += tree[it->second].total;
 493
 494                                 if (groupmap != NULL) {
 495                     vector<string> mGroups = groupmap->getNamesOfGroups();
 496                                         for (int i = 0; i < mGroups.size(); i++) { tree[0].groupCount[mGroups[i]] += tree[it->second].groupCount[mGroups[i]]; }
 497                                 }else if ( ct != NULL) {
 498                     vector<string> mGroups = ct->getNamesOfGroups();
 499                     if (ct->hasGroupInfo()) { for (int i = 0; i < mGroups.size(); i++) { tree[0].groupCount[mGroups[i]] += tree[it->second].groupCount[mGroups[i]]; } }
 500                 }
 501                         }
 502                 }
 503
 504                 //print root
 505                 out << tree[0].name << "\t" << "1.0000" << "\t"; //root relative abundance is 1, everyone classifies to root
 506
 507                 /*
 508                 if (groupmap != NULL) {
 509                         for (int i = 0; i < mGroups.size(); i++) {  out << tree[0].groupCount[mGroups[i]] << '\t'; }
 510         }else if ( ct != NULL) {
 511             if (ct->hasGroupInfo()) { for (int i = 0; i < mGroups.size(); i++) {  out << tree[0].groupCount[mGroups[i]] << '\t'; } }
 512         }*/
 513
 514         if (groupmap != NULL) {
 515             vector<string> mGroups = groupmap->getNamesOfGroups();
 516                         for (int i = 0; i < mGroups.size(); i++) {  out << "1.0000" << '\t'; }
 517         }else if ( ct != NULL) {
 518             vector<string> mGroups = ct->getNamesOfGroups();
 519             if (ct->hasGroupInfo()) { for (int i = 0; i < mGroups.size(); i++) {  out << "1.0000" << '\t'; } }
 520         }
 521
 522                 out << endl;
 523
 524                 //print rest
 525                 print(0, out, relabund);
 526
 527         }
 528         catch(exception& e) {
 529                 m->errorOut(e, "PhyloSummary", "print");
 530                 exit(1);
 531         }
 532 }
 533 /**************************************************************************************************/
 534
 535 void PhyloSummary::print(int i, ofstream& out){
 536         try {
 537                 map<string,int>::iterator it;
 538                 for(it=tree[i].children.begin();it!=tree[i].children.end();it++){
 539
 540                         if (tree[it->second].total != 0)  {
 541
 542                                 int totalChildrenInTree = 0;
 543
 544                                 map<string,int>::iterator it2;
 545                                 for(it2=tree[it->second].children.begin();it2!=tree[it->second].children.end();it2++){
 546                                         if (tree[it2->second].total != 0)  {   totalChildrenInTree++; }
 547                                 }
 548
 549                 if (relabund) {
 550                     out << tree[it->second].level << "\t" << tree[it->second].rank << "\t" << tree[it->second].name << "\t" << totalChildrenInTree << "\t" << (tree[it->second].total/(double) tree[0].total) << "\t";
 551                 }else {
 552                     out << tree[it->second].level << "\t" << tree[it->second].rank << "\t" << tree[it->second].name << "\t" << totalChildrenInTree << "\t" << tree[it->second].total << "\t";
 553                                 }
 554
 555                 if (relabund) {
 556                     map<string, int>::iterator itGroup;
 557                     if (groupmap != NULL) {
 558                         vector<string> mGroups = groupmap->getNamesOfGroups();
 559                         for (int i = 0; i < mGroups.size(); i++) {  out << (tree[it->second].groupCount[mGroups[i]]/(double)groupmap->getNumSeqs(mGroups[i])) << '\t'; }
 560                     }else if (ct != NULL) {
 561                         if (ct->hasGroupInfo()) {
 562                             vector<string> mGroups = ct->getNamesOfGroups();
 563                             for (int i = 0; i < mGroups.size(); i++) {  out << (tree[it->second].groupCount[mGroups[i]]/(double)ct->getGroupCount(mGroups[i])) << '\t'; }
 564                         }
 565                     }
 566                 }else {
 567                     map<string, int>::iterator itGroup;
 568                     if (groupmap != NULL) {
 569                         vector<string> mGroups = groupmap->getNamesOfGroups();
 570                         for (int i = 0; i < mGroups.size(); i++) {  out << tree[it->second].groupCount[mGroups[i]] << '\t'; }
 571                     }else if (ct != NULL) {
 572                         if (ct->hasGroupInfo()) {
 573                             vector<string> mGroups = ct->getNamesOfGroups();
 574                             for (int i = 0; i < mGroups.size(); i++) {  out << tree[it->second].groupCount[mGroups[i]] << '\t'; }
 575                         }
 576                     }
 577
 578                 }
 579                 out << endl;
 580
 581                         }
 582
 583                         print(it->second, out);
 584                 }
 585         }
 586         catch(exception& e) {
 587                 m->errorOut(e, "PhyloSummary", "print");
 588                 exit(1);
 589         }
 590 }
 591
 592 /**************************************************************************************************/
 593
 594 void PhyloSummary::print(int i, ofstream& out, bool relabund){
 595         try {
 596                 map<string,int>::iterator it;
 597                 for(it=tree[i].children.begin();it!=tree[i].children.end();it++){
 598
 599                         if (tree[it->second].total != 0)  {
 600
 601                                 int totalChildrenInTree = 0;
 602
 603                                 map<string,int>::iterator it2;
 604                                 for(it2=tree[it->second].children.begin();it2!=tree[it->second].children.end();it2++){
 605                                         if (tree[it2->second].total != 0)  {   totalChildrenInTree++; }
 606                                 }
 607
 608                 string nodeName = "";
 609                 int thisNode = it->second;
 610                 while (tree[thisNode].rank != "0") { //while you are not at top
 611                     if (m->control_pressed) { break; }
 612                     nodeName = tree[thisNode].name + "|" + nodeName;
 613                     thisNode = tree[thisNode].parent;
 614                 }
 615                 if (nodeName != "") { nodeName = nodeName.substr(0, nodeName.length()-1); }
 616
 617                                 out << nodeName << "\t" << (tree[it->second].total / (float)tree[i].total) << "\t";
 618
 619                                 map<string, int>::iterator itGroup;
 620                                 if (groupmap != NULL) {
 621                                         vector<string> mGroups = groupmap->getNamesOfGroups();
 622                                         for (int j = 0; j < mGroups.size(); j++) {
 623                         if (tree[i].groupCount[mGroups[j]] == 0) {
 624                             out << 0 << '\t';
 625                         }else { out << (tree[it->second].groupCount[mGroups[j]] / (float)tree[i].groupCount[mGroups[j]]) << '\t'; }
 626                     }
 627                                 }else if (ct != NULL) {
 628                     if (ct->hasGroupInfo()) {
 629                         vector<string> mGroups = ct->getNamesOfGroups();
 630                         for (int j = 0; j < mGroups.size(); j++) {
 631                             if (tree[i].groupCount[mGroups[j]] == 0) {
 632                                 out << 0 << '\t';
 633                             }else { out << (tree[it->second].groupCount[mGroups[j]] / (float)tree[i].groupCount[mGroups[j]]) << '\t'; }
 634                         }
 635                     }
 636                 }
 637                                 out << endl;
 638
 639                         }
 640
 641                         print(it->second, out, relabund);
 642                 }
 643         }
 644         catch(exception& e) {
 645                 m->errorOut(e, "PhyloSummary", "print");
 646                 exit(1);
 647         }
 648 }
 649 /**************************************************************************************************/
 650 void PhyloSummary::readTreeStruct(ifstream& in){
 651         try {
 652
 653                 //read version
 654                 string line = m->getline(in); m->gobble(in);
 655
 656                 int num;
 657
 658                 in >> num; m->gobble(in);
 659
 660                 tree.resize(num);
 661
 662                 in >> maxLevel; m->gobble(in);
 663
 664                 //read the tree file
 665                 for (int i = 0; i < tree.size(); i++) {
 666
 667                         in >> tree[i].level >> tree[i].name >> num; //num contains the number of children tree[i] has
 668
 669                         //set children
 670                         string childName;
 671                         int childIndex;
 672                         for (int j = 0; j < num; j++) {
 673                                 in >> childName >> childIndex;
 674                                 tree[i].children[childName] = childIndex;
 675                         }
 676
 677                         //initialize groupcounts
 678                         if (groupmap != NULL) {
 679                                 for (int j = 0; j < (groupmap->getNamesOfGroups()).size(); j++) {
 680                                         tree[i].groupCount[(groupmap->getNamesOfGroups())[j]] = 0;
 681                                 }
 682                         }else if (ct != NULL) {
 683                 if (ct->hasGroupInfo()) {
 684                     for (int j = 0; j < (ct->getNamesOfGroups()).size(); j++) {
 685                         tree[i].groupCount[(ct->getNamesOfGroups())[j]] = 0;
 686                     }
 687                 }
 688             }
 689
 690                         tree[i].total = 0;
 691
 692                         m->gobble(in);
 693
 694                         //if (tree[i].level > maxLevel) {  maxLevel = tree[i].level;  }
 695                 }
 696
 697         }
 698         catch(exception& e) {
 699                 m->errorOut(e, "PhyloSummary", "readTreeStruct");
 700                 exit(1);
 701         }
 702 }
 703 /**************************************************************************************************/
 704
 705
 706