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 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;
 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                                 readSpecialChar(filehandle,',',"comma");
 243                         }
 244                         // ';' means end of tree.
 245                         else if((ch=filehandle.peek())==';' || ch=='['){
 246                                 rooted = 1;
 247                         }
 248                         if(rooted != 1){
 249                                 rc = readNewickInt(filehandle, n, T);
 250                                 if (rc == -1) { return -1; } //reports an error in reading
 251                                 if(filehandle.peek() == ')'){
 252                                         readSpecialChar(filehandle,')',"right parenthesis");
 253                                 }
 254                         }
 255                 }
 256                 //note: treeclimber had the code below added - not sure why?
 257                 else{
 258                         filehandle.putback(ch);
 259                         char name[MAX_LINE];
 260                         filehandle.get(name, MAX_LINE,'\n');
 261                         SKIPLINE(filehandle, ch);
 262
 263                         n = T->getIndex(name);
 264
 265                         if(n!=0){
 266                                 cerr << "Internal error: The only taxon is not taxon 0.\n";
 267                                 //exit(1);
 268                                 readOk = -1; return -1;
 269                         }
 270                         lc = rc = -1;
 271                 }
 272
 273                 while((ch=filehandle.get())!=';'){;}
 274                 if(rooted != 1){
 275                         T->tree[n].setChildren(lc,rc);
 276                         T->tree[n].setBranchLength(0);
 277                         T->tree[n].setParent(-1);
 278                         if(lc!=-1){             T->tree[lc].setParent(n);               }
 279                         if(rc!=-1){             T->tree[rc].setParent(n);               }
 280                 }
 281                 return 0;
 282
 283         }
 284         catch(exception& e) {
 285                 cout << "Standard Error: " << e.what() << " has occurred in the ReadNewickTree class Function readTreeString. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 286                 exit(1);
 287         }
 288         catch(...) {
 289                 cout << "An unknown error has occurred in the ReadNewickTree class function readTreeString. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 290                 exit(1);
 291         }
 292
 293 }
 294 /**************************************************************************************************/
 295
 296 int ReadNewickTree::readNewickInt(istream& f, int& n, Tree* T) {
 297         try {
 298                 int c = readNodeChar(f);
 299
 300                 if(c == '('){
 301                         int lc = readNewickInt(f, n, T);
 302                         if (lc == -1) { return -1; } //reports an error in reading
 303                         readSpecialChar(f,',',"comma");
 304
 305                         int rc = readNewickInt(f, n, T);
 306                         if (rc == -1) { return -1; }  //reports an error in reading
 307                         if(f.peek()==')'){
 308                                 readSpecialChar(f,')',"right parenthesis");
 309                         }
 310
 311                         if(f.peek() == ':'){
 312                                 readSpecialChar(f,':',"colon");
 313
 314                                 if(n >= numNodes){      cerr << "Error: Too many nodes in input tree\n";  readOk = -1; return -1; }
 315
 316                                 T->tree[n].setBranchLength(readBranchLength(f));
 317                         }else{T->tree[n].setBranchLength(0.0); }
 318
 319                         T->tree[n].setChildren(lc,rc);
 320                         T->tree[lc].setParent(n);
 321                         T->tree[rc].setParent(n);
 322
 323                         return n++;
 324                 }else{
 325                         f.putback(c);
 326                         string name = "";
 327                         char d=f.get();
 328                         while(d != ':' && d != ',' && d!=')' && d!='\n'){
 329                                 name += d;
 330                                 d=f.get();
 331                         }
 332
 333                         int blen = 0;
 334                         if(d == ':')    {               blen = 1;                       }
 335
 336                         f.putback(d);
 337
 338                         //set group info
 339                         string group = globaldata->gTreemap->getGroup(name);
 340
 341                         //find index in tree of name
 342                         int n1 = T->getIndex(name);
 343
 344                         //adds sequence names that are not in group file to the "xxx" group
 345                         if(n1 == -1) {
 346                                 cerr << "Name: " << name << " not found in your groupfile. \n"; readOk = -1; return n1;
 347
 348                                 //globaldata->gTreemap->namesOfSeqs.push_back(name);
 349                                 //globaldata->gTreemap->treemap[name].groupname = "xxx";
 350                                 //globaldata->gTreemap->treemap[name].vectorIndex = (globaldata->gTreemap->namesOfSeqs.size() - 1);
 351
 352                                 //map<string, int>::iterator it;
 353                                 //it = globaldata->gTreemap->seqsPerGroup.find("xxx");
 354                                 //if (it == globaldata->gTreemap->seqsPerGroup.end()) { //its a new group
 355                                 //      globaldata->gTreemap->namesOfGroups.push_back("xxx");
 356                                 //      globaldata->gTreemap->seqsPerGroup["xxx"] = 1;
 357                                 //}else {
 358                                 //      globaldata->gTreemap->seqsPerGroup["xxx"]++;
 359                                 //}
 360
 361                                 //find index in tree of name
 362                                 //n1 = T->getIndex(name);
 363                                 //group = "xxx";
 364                                 //numLeaves++;
 365                                 //numNodes = 2*numLeaves - 1;
 366                         }
 367
 368                         T->tree[n1].setGroup(group);
 369                         T->tree[n1].setChildren(-1,-1);
 370
 371                         if(blen == 1){
 372                                 f.get();
 373                                 T->tree[n1].setBranchLength(readBranchLength(f));
 374                         }else{
 375                                 T->tree[n1].setBranchLength(0.0);
 376                         }
 377
 378                         while((c=f.get())!=0 && (c != ':' && c != ',' && c!=')') )              {;}
 379                         f.putback(c);
 380
 381                         return n1;
 382                 }
 383         }
 384         catch(exception& e) {
 385                 cout << "Standard Error: " << e.what() << " has occurred in the ReadNewickTree class Function readNewickInt. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 386                 exit(1);
 387         }
 388         catch(...) {
 389                 cout << "An unknown error has occurred in the ReadNewickTree class function readNewickInt. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 390                 exit(1);
 391         }
 392 }
 393 /**************************************************************************************************/
 394 /**************************************************************************************************/
 395