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";  return -1;
  36                         //exit(1);
  37                 }
  38                 if(d != c){
  39                         cerr << "Error: Expected " << name << " in input file.  Found " << d << ".\n";  return -1;
  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";  return -1;
  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";  return -1;
  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 int ReadNewickTree::read() {
 115         try {
 116                 int c, error;
 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                                 error = 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                                 error = readTreeString();
 162
 163                                 //save trees for later commands
 164                                 globaldata->gTree.push_back(T);
 165                         }
 166                 }
 167                 return readOk;
 168         }
 169         catch(exception& e) {
 170                 cout << "Standard Error: " << e.what() << " has occurred in the ReadNewickTree class Function read. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 171                 exit(1);
 172         }
 173         catch(...) {
 174                 cout << "An unknown error has occurred in the ReadNewickTree class function read. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 175                 exit(1);
 176         }
 177 }
 178 /**************************************************************************************************/
 179 //This function read the file through the translation of the sequences names and updates treemap.
 180 void ReadNewickTree::nexusTranslation() {
 181         try {
 182
 183                 holder = "";
 184                 int numSeqs = globaldata->gTreemap->getNumSeqs(); //must save this some when we clear old names we can still know how many sequences there were
 185                 int comment = 0;
 186
 187                 // get past comments
 188                 while(holder != "translate" && holder != "Translate"){
 189                         if(holder == "[" || holder == "[!"){
 190                                 comment = 1;
 191                         }
 192                         if(holder == "]"){
 193                                 comment = 0;
 194                         }
 195                         filehandle >> holder;
 196                         if(holder == "tree" && comment != 1){return;}
 197                 }
 198
 199                 //update treemap
 200                 globaldata->gTreemap->namesOfSeqs.clear();
 201                 for(int i=0;i<numSeqs;i++){
 202                         string number, name;
 203                         filehandle >> number;
 204                         filehandle >> name;
 205                         name.erase(name.end()-1);  //erase the comma
 206                         //insert new one with new name
 207                         globaldata->gTreemap->treemap[toString(number)].groupname = globaldata->gTreemap->treemap[name].groupname;
 208                         globaldata->gTreemap->treemap[toString(number)].vectorIndex = globaldata->gTreemap->treemap[name].vectorIndex;
 209                         //erase old one.  so treemap[sarah].groupnumber is now treemap[1].groupnumber. if number is 1 and name is sarah.
 210                         globaldata->gTreemap->treemap.erase(name);
 211                         globaldata->gTreemap->namesOfSeqs.push_back(number);
 212                 }
 213         }
 214         catch(exception& e) {
 215                 cout << "Standard Error: " << e.what() << " has occurred in the ReadNewickTree class Function nexus. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 216                 exit(1);
 217         }
 218         catch(...) {
 219                 cout << "An unknown error has occurred in the ReadNewickTree class function nexus. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 220                 exit(1);
 221         }
 222 }
 223
 224 /**************************************************************************************************/
 225 int ReadNewickTree::readTreeString() {
 226         try {
 227
 228                 int n = 0;
 229                 int lc, rc, error;
 230
 231                 int rooted = 0;
 232
 233                 int ch = filehandle.peek();
 234
 235                 if(ch == '('){
 236                         n = numLeaves;  //number of leaves / sequences, we want node 1 to start where the leaves left off
 237
 238                         lc = readNewickInt(filehandle, n, T);
 239                         if (lc == -1) { return -1; } //reports an error in reading
 240
 241                         if(filehandle.peek()==','){
 242                                 error = readSpecialChar(filehandle,',',"comma");
 243                                 if (error == -1) { readOk = -1; return -1; }
 244                         }
 245                         // ';' means end of tree.
 246                         else if((ch=filehandle.peek())==';' || ch=='['){
 247                                 rooted = 1;
 248                         }
 249                         if(rooted != 1){
 250                                 rc = readNewickInt(filehandle, n, T);
 251                                 if (rc == -1) { return -1; } //reports an error in reading
 252                                 if(filehandle.peek() == ')'){
 253                                         error = readSpecialChar(filehandle,')',"right parenthesis");
 254                                         if (error == -1) { readOk = -1; return -1; }
 255                                 }
 256                         }
 257                 }
 258                 //note: treeclimber had the code below added - not sure why?
 259                 else{
 260                         filehandle.putback(ch);
 261                         char name[MAX_LINE];
 262                         filehandle.get(name, MAX_LINE,'\n');
 263                         SKIPLINE(filehandle, ch);
 264
 265                         n = T->getIndex(name);
 266
 267                         if(n!=0){
 268                                 cerr << "Internal error: The only taxon is not taxon 0.\n";
 269                                 //exit(1);
 270                                 readOk = -1; return -1;
 271                         }
 272                         lc = rc = -1;
 273                 }
 274
 275                 while((ch=filehandle.get())!=';'){;}
 276                 if(rooted != 1){
 277                         T->tree[n].setChildren(lc,rc);
 278                         T->tree[n].setBranchLength(0);
 279                         T->tree[n].setParent(-1);
 280                         if(lc!=-1){             T->tree[lc].setParent(n);               }
 281                         if(rc!=-1){             T->tree[rc].setParent(n);               }
 282                 }
 283                 return 0;
 284
 285         }
 286         catch(exception& e) {
 287                 cout << "Standard Error: " << e.what() << " has occurred in the ReadNewickTree class Function readTreeString. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 288                 exit(1);
 289         }
 290         catch(...) {
 291                 cout << "An unknown error has occurred in the ReadNewickTree class function readTreeString. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 292                 exit(1);
 293         }
 294
 295 }
 296 /**************************************************************************************************/
 297
 298 int ReadNewickTree::readNewickInt(istream& f, int& n, Tree* T) {
 299         try {
 300                 int error;
 301
 302                 int c = readNodeChar(f);
 303                 if (c == -1) { readOk = -1; return -1; }
 304
 305                 if(c == '('){
 306                         int lc = readNewickInt(f, n, T);
 307                         if (lc == -1) { return -1; } //reports an error in reading
 308                         error = readSpecialChar(f,',',"comma");
 309                         if (error == -1) { readOk = -1; return -1; }
 310
 311                         int rc = readNewickInt(f, n, T);
 312                         if (rc == -1) { return -1; }  //reports an error in reading
 313                         if(f.peek()==')'){
 314                                 error = readSpecialChar(f,')',"right parenthesis");
 315                                 if (error == -1) { readOk = -1; return -1; }
 316                         }
 317
 318                         if(f.peek() == ':'){
 319                                 error = readSpecialChar(f,':',"colon");
 320                                 if (error == -1) { readOk = -1; return -1; }
 321                                 if(n >= numNodes){      cerr << "Error: Too many nodes in input tree\n";  readOk = -1; return -1; }
 322                                 error = readBranchLength(f);
 323                                 if (error == -1) { readOk = -1; return -1; }
 324                                 T->tree[n].setBranchLength(error);
 325                         }else{T->tree[n].setBranchLength(0.0); }
 326
 327                         T->tree[n].setChildren(lc,rc);
 328                         T->tree[lc].setParent(n);
 329                         T->tree[rc].setParent(n);
 330
 331                         return n++;
 332                 }else{
 333                         f.putback(c);
 334                         string name = "";
 335                         char d=f.get();
 336                         while(d != ':' && d != ',' && d!=')' && d!='\n'){
 337                                 name += d;
 338                                 d=f.get();
 339                         }
 340
 341                         int blen = 0;
 342                         if(d == ':')    {               blen = 1;                       }
 343
 344                         f.putback(d);
 345
 346                         //set group info
 347                         string group = globaldata->gTreemap->getGroup(name);
 348
 349                         //find index in tree of name
 350                         int n1 = T->getIndex(name);
 351
 352                         //adds sequence names that are not in group file to the "xxx" group
 353                         if(n1 == -1) {
 354                                 cerr << "Name: " << name << " not found in your groupfile. \n"; readOk = -1; return n1;
 355
 356                                 //globaldata->gTreemap->namesOfSeqs.push_back(name);
 357                                 //globaldata->gTreemap->treemap[name].groupname = "xxx";
 358                                 //globaldata->gTreemap->treemap[name].vectorIndex = (globaldata->gTreemap->namesOfSeqs.size() - 1);
 359
 360                                 //map<string, int>::iterator it;
 361                                 //it = globaldata->gTreemap->seqsPerGroup.find("xxx");
 362                                 //if (it == globaldata->gTreemap->seqsPerGroup.end()) { //its a new group
 363                                 //      globaldata->gTreemap->namesOfGroups.push_back("xxx");
 364                                 //      globaldata->gTreemap->seqsPerGroup["xxx"] = 1;
 365                                 //}else {
 366                                 //      globaldata->gTreemap->seqsPerGroup["xxx"]++;
 367                                 //}
 368
 369                                 //find index in tree of name
 370                                 //n1 = T->getIndex(name);
 371                                 //group = "xxx";
 372                                 //numLeaves++;
 373                                 //numNodes = 2*numLeaves - 1;
 374                         }
 375
 376                         T->tree[n1].setGroup(group);
 377                         T->tree[n1].setChildren(-1,-1);
 378
 379                         if(blen == 1){
 380                                 f.get();
 381                                 error = readBranchLength(f);
 382                                 if (error == -1) { readOk = -1; return -1; }
 383                                 T->tree[n1].setBranchLength(error);
 384                         }else{
 385                                 T->tree[n1].setBranchLength(0.0);
 386                         }
 387
 388                         while((c=f.get())!=0 && (c != ':' && c != ',' && c!=')') )              {;}
 389                         f.putback(c);
 390
 391                         return n1;
 392                 }
 393         }
 394         catch(exception& e) {
 395                 cout << "Standard Error: " << e.what() << " has occurred in the ReadNewickTree class Function readNewickInt. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 396                 exit(1);
 397         }
 398         catch(...) {
 399                 cout << "An unknown error has occurred in the ReadNewickTree class function readNewickInt. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 400                 exit(1);
 401         }
 402 }
 403 /**************************************************************************************************/
 404 /**************************************************************************************************/
 405