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