readtree.cpp

   1 /*
   2  *  readtree.cpp
   3  *  Mothur
   4  *
   5  *  Created by Sarah Westcott on 1/22/09.
   6  *  Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
   7  *
   8  */
   9
  10 #include "readtree.h"
  11
  12 /***********************************************************************/
  13 ReadTree::ReadTree() {
  14         try {
  15                 globaldata = GlobalData::getInstance();
  16                 globaldata->gTree.clear();
  17         }
  18         catch(exception& e) {
  19                 cout << "Standard Error: " << e.what() << " has occurred in the ReadTree class Function ReadTree. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  20                 exit(1);
  21         }
  22         catch(...) {
  23                 cout << "An unknown error has occurred in the ReadTree class function ReadTree. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  24                 exit(1);
  25         }
  26 }
  27 /***********************************************************************/
  28 int ReadTree::readSpecialChar(istream& f, char c, string name) {
  29     try {
  30
  31                 gobble(f);
  32                 char d = f.get();
  33
  34                 if(d == EOF){
  35                         cerr << "Error: Input file ends prematurely, expecting a " << name << "\n";
  36                         exit(1);
  37                 }
  38                 if(d != c){
  39                         cerr << "Error: Expected " << name << " in input file.  Found " << d << ".\n";
  40                         exit(1);
  41                 }
  42                 if(d == ')' && f.peek() == '\n'){
  43                         gobble(f);
  44                 }
  45                 return d;
  46         }
  47         catch(exception& e) {
  48                 cout << "Standard Error: " << e.what() << " has occurred in the ReadTree class Function readSpecialChar. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  49                 exit(1);
  50         }
  51         catch(...) {
  52                 cout << "An unknown error has occurred in the ReadTree class function readSpecialChar. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  53                 exit(1);
  54         }
  55 }
  56 /**************************************************************************************************/
  57
  58 int ReadTree::readNodeChar(istream& f) {
  59         try {
  60 //              while(isspace(d=f.get()))               {;}
  61                 gobble(f);
  62                 char d = f.get();
  63
  64                 if(d == EOF){
  65                         cerr << "Error: Input file ends prematurely, expecting a left parenthesis\n";
  66                         exit(1);
  67                 }
  68                 return d;
  69         }
  70         catch(exception& e) {
  71                 cout << "Standard Error: " << e.what() << " has occurred in the ReadTree class Function readNodeChar. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  72                 exit(1);
  73         }
  74         catch(...) {
  75                 cout << "An unknown error has occurred in the ReadTree class function readNodeChar. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  76                 exit(1);
  77         }
  78 }
  79
  80 /**************************************************************************************************/
  81
  82 float ReadTree::readBranchLength(istream& f) {
  83     try {
  84                 float b;
  85
  86                 if(!(f >> b)){
  87                         cerr << "Error: Missing branch length in input tree.\n";
  88                         exit(1);
  89                 }
  90                 gobble(f);
  91                 return b;
  92         }
  93         catch(exception& e) {
  94                 cout << "Standard Error: " << e.what() << " has occurred in the ReadTree class Function readBranchLength. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  95                 exit(1);
  96         }
  97         catch(...) {
  98                 cout << "An unknown error has occurred in the ReadTree class function readBranchLength. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  99                 exit(1);
 100         }
 101 }
 102
 103
 104 /***********************************************************************/
 105 /***********************************************************************/
 106
 107
 108 //Child Classes Below
 109
 110 /***********************************************************************/
 111 /***********************************************************************/
 112 //This class reads a file in Newick form and stores it in a tree.
 113
 114 void ReadNewickTree::read() {
 115         try {
 116                 int c;
 117                 int comment = 0;
 118
 119                 //if you are not a nexus file
 120                 if ((c = filehandle.peek()) != '#') {
 121                         while((c = filehandle.peek()) != EOF) {
 122                                 //make new tree
 123                                 T = new Tree();
 124                                 numNodes = T->getNumNodes();
 125                                 numLeaves = T->getNumLeaves();
 126
 127                                 readTreeString();
 128
 129                                 //save trees for later commands
 130                                 globaldata->gTree.push_back(T);
 131                                 gobble(filehandle);
 132                         }
 133                 //if you are a nexus file
 134                 }else if ((c = filehandle.peek()) == '#') {
 135                         nexusTranslation();  //reads file through the translation and updates treemap
 136                         while((c = filehandle.peek()) != EOF) {
 137                                 // get past comments
 138                                 while ((c = filehandle.peek()) != EOF) {
 139                                         if(holder == "[" || holder == "[!"){
 140                                                 comment = 1;
 141                                         }
 142                                         if(holder == "]"){
 143                                                 comment = 0;
 144                                         }
 145                                         if((holder == "tree" || holder == "end;") && comment != 1){ holder = ""; comment = 0; break;}
 146                                         filehandle >> holder;
 147                                 }
 148
 149                                 //pass over the "tree rep.6878900 = "
 150                                 while (((c = filehandle.get()) != '(') && ((c = filehandle.peek()) != EOF) ) {;}
 151
 152                                 if (c == EOF ) { break; }
 153                                 filehandle.putback(c);  //put back first ( of tree.
 154
 155                                 //make new tree
 156                                 T = new Tree();
 157                                 numNodes = T->getNumNodes();
 158                                 numLeaves = T->getNumLeaves();
 159
 160                                 //read tree info
 161                                 readTreeString();
 162
 163                                 //save trees for later commands
 164                                 globaldata->gTree.push_back(T);
 165                         }
 166                 }
 167         }
 168         catch(exception& e) {
 169                 cout << "Standard Error: " << e.what() << " has occurred in the ReadNewickTree class Function read. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 170                 exit(1);
 171         }
 172         catch(...) {
 173                 cout << "An unknown error has occurred in the ReadNewickTree class function read. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 174                 exit(1);
 175         }
 176 }
 177 /**************************************************************************************************/
 178 //This function read the file through the translation of the sequences names and updates treemap.
 179 void ReadNewickTree::nexusTranslation() {
 180         try {
 181
 182                 holder = "";
 183                 int numSeqs = globaldata->gTreemap->getNumSeqs(); //must save this some when we clear old names we can still know how many sequences there were
 184                 int comment = 0;
 185
 186                 // get past comments
 187                 while(holder != "translate" && holder != "Translate"){
 188                         if(holder == "[" || holder == "[!"){
 189                                 comment = 1;
 190                         }
 191                         if(holder == "]"){
 192                                 comment = 0;
 193                         }
 194                         filehandle >> holder;
 195                         if(holder == "tree" && comment != 1){return;}
 196                 }
 197
 198                 //update treemap
 199                 globaldata->gTreemap->namesOfSeqs.clear();
 200                 for(int i=0;i<numSeqs;i++){
 201                         string number, name;
 202                         filehandle >> number;
 203                         filehandle >> name;
 204                         name.erase(name.end()-1);  //erase the comma
 205                         //insert new one with new name
 206                         globaldata->gTreemap->treemap[toString(number)].groupname = globaldata->gTreemap->treemap[name].groupname;
 207                         globaldata->gTreemap->treemap[toString(number)].vectorIndex = globaldata->gTreemap->treemap[name].vectorIndex;
 208                         //erase old one.  so treemap[sarah].groupnumber is now treemap[1].groupnumber. if number is 1 and name is sarah.
 209                         globaldata->gTreemap->treemap.erase(name);
 210                         globaldata->gTreemap->namesOfSeqs.push_back(number);
 211                 }
 212         }
 213         catch(exception& e) {
 214                 cout << "Standard Error: " << e.what() << " has occurred in the ReadNewickTree class Function nexus. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 215                 exit(1);
 216         }
 217         catch(...) {
 218                 cout << "An unknown error has occurred in the ReadNewickTree class function nexus. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 219                 exit(1);
 220         }
 221 }
 222
 223 /**************************************************************************************************/
 224 void ReadNewickTree::readTreeString() {
 225         try {
 226
 227                 int n = 0;
 228                 int lc, rc;
 229
 230                 int rooted = 0;
 231
 232                 int ch = filehandle.peek();
 233
 234                 if(ch == '('){
 235                         n = numLeaves;  //number of leaves / sequences, we want node 1 to start where the leaves left off
 236
 237                         lc = readNewickInt(filehandle, n, T);
 238
 239                         if(filehandle.peek()==','){
 240                                 readSpecialChar(filehandle,',',"comma");
 241                         }
 242                         // ';' means end of tree.
 243                         else if((ch=filehandle.peek())==';' || ch=='['){
 244                                 rooted = 1;
 245                         }
 246                         if(rooted != 1){
 247                                 rc = readNewickInt(filehandle, n, T);
 248                                 if(filehandle.peek() == ')'){
 249                                         readSpecialChar(filehandle,')',"right parenthesis");
 250                                 }
 251                         }
 252                 }
 253                 //note: treeclimber had the code below added - not sure why?
 254                 else{
 255                         filehandle.putback(ch);
 256                         char name[MAX_LINE];
 257                         filehandle.get(name, MAX_LINE,'\n');
 258                         SKIPLINE(filehandle, ch);
 259
 260                         n = T->getIndex(name);
 261
 262                         if(n!=0){
 263                                 cerr << "Internal error: The only taxon is not taxon 0.\n";
 264                                 exit(1);
 265                         }
 266                         lc = rc = -1;
 267                 }
 268
 269                 while((ch=filehandle.get())!=';'){;}
 270                 if(rooted != 1){
 271                         T->tree[n].setChildren(lc,rc);
 272                         T->tree[n].setBranchLength(0);
 273                         T->tree[n].setParent(-1);
 274                         if(lc!=-1){             T->tree[lc].setParent(n);               }
 275                         if(rc!=-1){             T->tree[rc].setParent(n);               }
 276                         cout << "new loop "<< endl;
 277                 }
 278
 279         }
 280         catch(exception& e) {
 281                 cout << "Standard Error: " << e.what() << " has occurred in the ReadNewickTree class Function readTreeString. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 282                 exit(1);
 283         }
 284         catch(...) {
 285                 cout << "An unknown error has occurred in the ReadNewickTree class function readTreeString. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 286                 exit(1);
 287         }
 288
 289 }
 290 /**************************************************************************************************/
 291
 292 int ReadNewickTree::readNewickInt(istream& f, int& n, Tree* T) {
 293         try {
 294                 int c = readNodeChar(f);
 295
 296                 if(c == '('){
 297                         int lc = readNewickInt(f, n, T);
 298                         readSpecialChar(f,',',"comma");
 299
 300                         int rc = readNewickInt(f, n, T);
 301                         if(f.peek()==')'){
 302                                 readSpecialChar(f,')',"right parenthesis");
 303                         }
 304
 305                         if(f.peek() == ':'){
 306                                 readSpecialChar(f,':',"colon");
 307                                 if(n >= numNodes){      cerr << "Error: Too many nodes in input tree\n";  exit(1); }
 308                                 T->tree[n].setBranchLength(readBranchLength(f));
 309                         }else{T->tree[n].setBranchLength(0.0); }
 310
 311                         T->tree[n].setChildren(lc,rc);
 312                         T->tree[lc].setParent(n);
 313                         T->tree[rc].setParent(n);
 314
 315                         return n++;
 316                 }else{
 317                         f.putback(c);
 318                         string name = "";
 319                         char d=f.get();
 320                         while(d != ':' && d != ',' && d!=')' && d!='\n'){
 321                                 name += d;
 322                                 d=f.get();
 323                         }
 324
 325                         int blen = 0;
 326                         if(d == ':')    {               blen = 1;                       }
 327
 328                         f.putback(d);
 329
 330                         //set group info
 331                         string group = globaldata->gTreemap->getGroup(name);
 332
 333                         //find index in tree of name
 334                         int n1 = T->getIndex(name);
 335
 336                         //adds sequence names that are not in group file to the "xxx" group
 337                         if(n1 == -1) {
 338                                 cerr << "Name: " << name << " not found in your groupfile and it will be ignored. \n";
 339
 340                                 globaldata->gTreemap->namesOfSeqs.push_back(name);
 341                                 globaldata->gTreemap->treemap[name].groupname = "xxx";
 342                                 globaldata->gTreemap->treemap[name].vectorIndex = (globaldata->gTreemap->namesOfSeqs.size() - 1);
 343
 344                                 map<string, int>::iterator it;
 345                                 it = globaldata->gTreemap->seqsPerGroup.find("xxx");
 346                                 if (it == globaldata->gTreemap->seqsPerGroup.end()) { //its a new group
 347                                         globaldata->gTreemap->namesOfGroups.push_back("xxx");
 348                                         globaldata->gTreemap->seqsPerGroup["xxx"] = 1;
 349                                 }else {
 350                                         globaldata->gTreemap->seqsPerGroup["xxx"]++;
 351                                 }
 352
 353                                 //find index in tree of name
 354                                 n1 = T->getIndex(name);
 355                                 group = "xxx";
 356                                 numLeaves++;
 357                                 numNodes = 2*numLeaves - 1;
 358                                 T->resetTree();
 359                         }
 360
 361                         T->tree[n1].setGroup(group);
 362                         T->printTree();
 363                         T->tree[n1].setChildren(-1,-1);
 364
 365                         if(blen == 1){
 366                                 f.get();
 367                                 T->tree[n1].setBranchLength(readBranchLength(f));
 368                         }else{
 369                                 T->tree[n1].setBranchLength(0.0);
 370                         }
 371
 372                         while((c=f.get())!=0 && (c != ':' && c != ',' && c!=')') )              {;}
 373                         f.putback(c);
 374
 375                         return n1;
 376                 }
 377         }
 378         catch(exception& e) {
 379                 cout << "Standard Error: " << e.what() << " has occurred in the ReadNewickTree class Function readNewickInt. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 380                 exit(1);
 381         }
 382         catch(...) {
 383                 cout << "An unknown error has occurred in the ReadNewickTree class function readNewickInt. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 384                 exit(1);
 385         }
 386 }
 387 /**************************************************************************************************/
 388 /**************************************************************************************************/
 389