5 * Created by Sarah Westcott on 1/22/09.
6 * Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
12 /*****************************************************************/
13 Tree::Tree(int num, TreeMap* t) : tmap(t) {
15 m = MothurOut::getInstance();
18 numNodes = 2*numLeaves - 1;
20 tree.resize(numNodes);
23 m->errorOut(e, "Tree", "Tree - numNodes");
27 /*****************************************************************/
28 Tree::Tree(string g) { //do not use tree generated by this its just to extract the treenames, its a chicken before the egg thing that needs to be revisited.
30 m = MothurOut::getInstance();
31 parseTreeFile(); m->runParse = false;
34 m->errorOut(e, "Tree", "Tree - just parse");
38 /*****************************************************************/
39 Tree::Tree(TreeMap* t) : tmap(t) {
41 m = MothurOut::getInstance();
43 if (m->runParse == true) { parseTreeFile(); m->runParse = false; }
44 //for(int i = 0; i < globaldata->Treenames.size(); i++) { cout << i << '\t' << globaldata->Treenames[i] << endl; }
45 numLeaves = m->Treenames.size();
46 numNodes = 2*numLeaves - 1;
48 tree.resize(numNodes);
50 //initialize groupNodeInfo
51 for (int i = 0; i < (tmap->getNamesOfGroups()).size(); i++) {
52 groupNodeInfo[(tmap->getNamesOfGroups())[i]].resize(0);
55 //initialize tree with correct number of nodes, name and group info.
56 for (int i = 0; i < numNodes; i++) {
57 //initialize leaf nodes
58 if (i <= (numLeaves-1)) {
59 tree[i].setName(m->Treenames[i]);
62 string group = tmap->getGroup(m->Treenames[i]);
64 vector<string> tempGroups; tempGroups.push_back(group);
65 tree[i].setGroup(tempGroups);
66 groupNodeInfo[group].push_back(i);
68 //set pcount and pGroup for groupname to 1.
69 tree[i].pcount[group] = 1;
70 tree[i].pGroups[group] = 1;
72 //Treemap knows name, group and index to speed up search
73 tmap->setIndex(m->Treenames[i], i);
75 //intialize non leaf nodes
76 }else if (i > (numLeaves-1)) {
78 vector<string> tempGroups;
79 tree[i].setGroup(tempGroups);
85 m->errorOut(e, "Tree", "Tree");
89 /*****************************************************************/
90 Tree::Tree(TreeMap* t, vector< vector<double> >& sims) : tmap(t) {
92 m = MothurOut::getInstance();
94 if (m->runParse == true) { parseTreeFile(); m->runParse = false; }
95 //for(int i = 0; i < globaldata->Treenames.size(); i++) { cout << i << '\t' << globaldata->Treenames[i] << endl; }
96 numLeaves = m->Treenames.size();
97 numNodes = 2*numLeaves - 1;
99 tree.resize(numNodes);
101 //initialize groupNodeInfo
102 for (int i = 0; i < (tmap->getNamesOfGroups()).size(); i++) {
103 groupNodeInfo[(tmap->getNamesOfGroups())[i]].resize(0);
106 //initialize tree with correct number of nodes, name and group info.
107 for (int i = 0; i < numNodes; i++) {
108 //initialize leaf nodes
109 if (i <= (numLeaves-1)) {
110 tree[i].setName(m->Treenames[i]);
113 string group = tmap->getGroup(m->Treenames[i]);
115 vector<string> tempGroups; tempGroups.push_back(group);
116 tree[i].setGroup(tempGroups);
117 groupNodeInfo[group].push_back(i);
119 //set pcount and pGroup for groupname to 1.
120 tree[i].pcount[group] = 1;
121 tree[i].pGroups[group] = 1;
123 //Treemap knows name, group and index to speed up search
124 tmap->setIndex(m->Treenames[i], i);
126 //intialize non leaf nodes
127 }else if (i > (numLeaves-1)) {
129 vector<string> tempGroups;
130 tree[i].setGroup(tempGroups);
134 //build tree from matrix
136 map<int, int> indexes; //maps row in simMatrix to vector index in the tree
137 int numGroups = (tmap->getNamesOfGroups()).size();
138 for (int g = 0; g < numGroups; g++) { indexes[g] = g; }
140 //do merges and create tree structure by setting parents and children
141 //there are numGroups - 1 merges to do
142 for (int i = 0; i < (numGroups - 1); i++) {
143 float largest = -1000.0;
145 if (m->control_pressed) { break; }
148 //find largest value in sims matrix by searching lower triangle
149 for (int j = 1; j < sims.size(); j++) {
150 for (int k = 0; k < j; k++) {
151 if (sims[j][k] > largest) { largest = sims[j][k]; row = j; column = k; }
155 //set non-leaf node info and update leaves to know their parents
157 tree[numGroups + i].setChildren(indexes[row], indexes[column]);
160 tree[indexes[row]].setParent(numGroups + i);
161 tree[indexes[column]].setParent(numGroups + i);
163 //blength = distance / 2;
164 float blength = ((1.0 - largest) / 2);
167 tree[indexes[row]].setBranchLength(blength - tree[indexes[row]].getLengthToLeaves());
168 tree[indexes[column]].setBranchLength(blength - tree[indexes[column]].getLengthToLeaves());
170 //set your length to leaves to your childs length plus branchlength
171 tree[numGroups + i].setLengthToLeaves(tree[indexes[row]].getLengthToLeaves() + tree[indexes[row]].getBranchLength());
175 indexes[row] = numGroups+i;
176 indexes[column] = numGroups+i;
178 //remove highest value that caused the merge.
179 sims[row][column] = -1000.0;
180 sims[column][row] = -1000.0;
182 //merge values in simsMatrix
183 for (int n = 0; n < sims.size(); n++) {
184 //row becomes merge of 2 groups
185 sims[row][n] = (sims[row][n] + sims[column][n]) / 2;
186 sims[n][row] = sims[row][n];
188 sims[column][n] = -1000.0;
189 sims[n][column] = -1000.0;
193 //adjust tree to make sure root to tip length is .5
194 int root = findRoot();
195 tree[root].setBranchLength((0.5 - tree[root].getLengthToLeaves()));
198 catch(exception& e) {
199 m->errorOut(e, "Tree", "Tree");
203 /*****************************************************************/
205 /*****************************************************************/
206 void Tree::addNamesToCounts(map<string, string> nameMap) {
208 //ex. seq1 seq2,seq3,se4
214 //before this function seq1.pcount = pasture -> 1
215 //after seq1.pcount = pasture -> 2, forest -> 1, ocean -> 1
217 //before this function seq1.pgroups = pasture -> 1
218 //after seq1.pgroups = pasture -> 1 since that is the dominant group
221 //go through each leaf and update its pcounts and pgroups
225 for (int i = 0; i < numLeaves; i++) {
227 string name = tree[i].getName();
229 map<string, string>::iterator itNames = nameMap.find(name);
231 if (itNames == nameMap.end()) { m->mothurOut(name + " is not in your name file, please correct."); m->mothurOutEndLine(); exit(1); }
233 vector<string> dupNames;
234 m->splitAtComma(nameMap[name], dupNames);
236 map<string, int>::iterator itCounts;
238 set<string> groupsAddedForThisNode;
239 for (int j = 0; j < dupNames.size(); j++) {
241 string group = tmap->getGroup(dupNames[j]);
243 if (dupNames[j] != name) {//you already added yourself in the constructor
245 if (groupsAddedForThisNode.count(group) == 0) { groupNodeInfo[group].push_back(i); groupsAddedForThisNode.insert(group); } //if you have not already added this node for this group, then add it
248 itCounts = tree[i].pcount.find(group);
249 if (itCounts == tree[i].pcount.end()) { //new group, add it
250 tree[i].pcount[group] = 1;
252 tree[i].pcount[group]++;
256 itCounts = tree[i].pGroups.find(group);
257 if (itCounts == tree[i].pGroups.end()) { //new group, add it
258 tree[i].pGroups[group] = 1;
260 tree[i].pGroups[group]++;
264 if(tree[i].pGroups[group] > maxPars){
265 maxPars = tree[i].pGroups[group];
267 }else { groupsAddedForThisNode.insert(group); } //add it so you don't add it to groupNodeInfo again
270 if (maxPars > 1) { //then we have some more dominant groups
271 //erase all the groups that are less than maxPars because you found a more dominant group.
272 for(it=tree[i].pGroups.begin();it!=tree[i].pGroups.end();){
273 if(it->second < maxPars){
274 tree[i].pGroups.erase(it++);
277 //set one remaining groups to 1
278 for(it=tree[i].pGroups.begin();it!=tree[i].pGroups.end();it++){
279 tree[i].pGroups[it->first] = 1;
283 //update groups to reflect all the groups this node represents
284 vector<string> nodeGroups;
285 map<string, int>::iterator itGroups;
286 for (itGroups = tree[i].pcount.begin(); itGroups != tree[i].pcount.end(); itGroups++) {
287 nodeGroups.push_back(itGroups->first);
289 tree[i].setGroup(nodeGroups);
295 //cout << "addNamesToCounts\t" << (B - A) / CLOCKS_PER_SEC << endl;
298 catch(exception& e) {
299 m->errorOut(e, "Tree", "addNamesToCounts");
303 /*****************************************************************/
304 int Tree::getIndex(string searchName) {
306 //Treemap knows name, group and index to speed up search
307 // getIndex function will return the vector index or -1 if seq is not found.
308 int index = tmap->getIndex(searchName);
312 catch(exception& e) {
313 m->errorOut(e, "Tree", "getIndex");
317 /*****************************************************************/
319 void Tree::setIndex(string searchName, int index) {
321 //set index in treemap
322 tmap->setIndex(searchName, index);
324 catch(exception& e) {
325 m->errorOut(e, "Tree", "setIndex");
329 /*****************************************************************/
330 int Tree::assembleTree(map<string, string> nameMap) {
335 //if user has given a names file we want to include that info in the pgroups and pcount info.
336 if(nameMap.size() != 0) { addNamesToCounts(nameMap); }
338 //build the pGroups in non leaf nodes to be used in the parsimony calcs.
339 for (int i = numLeaves; i < numNodes; i++) {
340 if (m->control_pressed) { return 1; }
342 tree[i].pGroups = (mergeGroups(i));
343 tree[i].pcount = (mergeGcounts(i));
348 catch(exception& e) {
349 m->errorOut(e, "Tree", "assembleTree");
353 /*****************************************************************
354 int Tree::assembleTree(string n) {
357 //build the pGroups in non leaf nodes to be used in the parsimony calcs.
358 for (int i = numLeaves; i < numNodes; i++) {
359 if (m->control_pressed) { return 1; }
361 tree[i].pGroups = (mergeGroups(i));
362 tree[i].pcount = (mergeGcounts(i));
365 //cout << "assembleTree\t" << (B-A) / CLOCKS_PER_SEC << endl;
368 catch(exception& e) {
369 m->errorOut(e, "Tree", "assembleTree");
373 /*****************************************************************/
374 //assumes leaf node names are in groups and no names file - used by indicator command
375 void Tree::getSubTree(Tree* Ctree, vector<string> Groups) {
378 //copy Tree since we are going to destroy it
379 Tree* copy = new Tree(tmap);
380 copy->getCopy(Ctree);
381 map<string, string> empty;
382 copy->assembleTree(empty);
384 //we want to select some of the leaf nodes to create the output tree
385 //go through the input Tree starting at parents of leaves
386 for (int i = 0; i < numNodes; i++) {
388 //initialize leaf nodes
389 if (i <= (numLeaves-1)) {
390 tree[i].setName(Groups[i]);
393 string group = tmap->getGroup(Groups[i]);
394 vector<string> tempGroups; tempGroups.push_back(group);
395 tree[i].setGroup(tempGroups);
396 groupNodeInfo[group].push_back(i);
398 //set pcount and pGroup for groupname to 1.
399 tree[i].pcount[group] = 1;
400 tree[i].pGroups[group] = 1;
402 //Treemap knows name, group and index to speed up search
403 tmap->setIndex(Groups[i], i);
405 //intialize non leaf nodes
406 }else if (i > (numLeaves-1)) {
408 vector<string> tempGroups;
409 tree[i].setGroup(tempGroups);
413 set<int> removedLeaves;
414 for (int i = 0; i < copy->getNumLeaves(); i++) {
416 if (removedLeaves.count(i) == 0) {
419 int parent = copy->tree[i].getParent();
423 if (m->inUsersGroups(copy->tree[i].getName(), Groups)) {
424 //find my siblings name
425 int parentRC = copy->tree[parent].getRChild();
426 int parentLC = copy->tree[parent].getLChild();
428 //if I am the right child, then my sib is the left child
429 int sibIndex = parentRC;
430 if (parentRC == i) { sibIndex = parentLC; }
432 string sibsName = copy->tree[sibIndex].getName();
434 //if yes, is my sibling
435 if ((m->inUsersGroups(sibsName, Groups)) || (sibsName == "")) {
436 //we both are okay no trimming required
438 //i am, my sib is not, so remove sib by setting my parent to my grandparent
439 int grandparent = copy->tree[parent].getParent();
440 int grandparentLC = copy->tree[grandparent].getLChild();
441 int grandparentRC = copy->tree[grandparent].getRChild();
443 //whichever of my granparents children was my parent now equals me
444 if (grandparentLC == parent) { grandparentLC = i; }
445 else { grandparentRC = i; }
447 copy->tree[i].setParent(grandparent);
448 copy->tree[i].setBranchLength((copy->tree[i].getBranchLength()+copy->tree[parent].getBranchLength()));
449 if (grandparent != -1) {
450 copy->tree[grandparent].setChildren(grandparentLC, grandparentRC);
452 removedLeaves.insert(sibIndex);
455 //find my siblings name
456 int parentRC = copy->tree[parent].getRChild();
457 int parentLC = copy->tree[parent].getLChild();
459 //if I am the right child, then my sib is the left child
460 int sibIndex = parentRC;
461 if (parentRC == i) { sibIndex = parentLC; }
463 string sibsName = copy->tree[sibIndex].getName();
465 //if no is my sibling
466 if ((m->inUsersGroups(sibsName, Groups)) || (sibsName == "")) {
467 //i am not, but my sib is
468 int grandparent = copy->tree[parent].getParent();
469 int grandparentLC = copy->tree[grandparent].getLChild();
470 int grandparentRC = copy->tree[grandparent].getRChild();
472 //whichever of my granparents children was my parent now equals my sib
473 if (grandparentLC == parent) { grandparentLC = sibIndex; }
474 else { grandparentRC = sibIndex; }
476 copy->tree[sibIndex].setParent(grandparent);
477 copy->tree[sibIndex].setBranchLength((copy->tree[sibIndex].getBranchLength()+copy->tree[parent].getBranchLength()));
478 if (grandparent != -1) {
479 copy->tree[grandparent].setChildren(grandparentLC, grandparentRC);
481 removedLeaves.insert(i);
483 //neither of us are, so we want to eliminate ourselves and our parent
484 //so set our parents sib to our great-grandparent
485 int parent = copy->tree[i].getParent();
486 int grandparent = copy->tree[parent].getParent();
488 if (grandparent != -1) {
489 int greatgrandparent = copy->tree[grandparent].getParent();
490 int greatgrandparentLC, greatgrandparentRC;
491 if (greatgrandparent != -1) {
492 greatgrandparentLC = copy->tree[greatgrandparent].getLChild();
493 greatgrandparentRC = copy->tree[greatgrandparent].getRChild();
496 int grandparentLC = copy->tree[grandparent].getLChild();
497 int grandparentRC = copy->tree[grandparent].getRChild();
499 parentsSibIndex = grandparentLC;
500 if (grandparentLC == parent) { parentsSibIndex = grandparentRC; }
502 //whichever of my greatgrandparents children was my grandparent
503 if (greatgrandparentLC == grandparent) { greatgrandparentLC = parentsSibIndex; }
504 else { greatgrandparentRC = parentsSibIndex; }
506 copy->tree[parentsSibIndex].setParent(greatgrandparent);
507 copy->tree[parentsSibIndex].setBranchLength((copy->tree[parentsSibIndex].getBranchLength()+copy->tree[grandparent].getBranchLength()));
508 if (greatgrandparent != -1) {
509 copy->tree[greatgrandparent].setChildren(greatgrandparentLC, greatgrandparentRC);
512 copy->tree[parent].setParent(-1);
513 //cout << "issues with making subtree" << endl;
515 removedLeaves.insert(sibIndex);
516 removedLeaves.insert(i);
524 for (int i = 0; i < copy->getNumNodes(); i++) {
526 if (copy->tree[i].getParent() == -1) { root = i; break; }
529 int nextSpot = numLeaves;
530 populateNewTree(copy->tree, root, nextSpot);
534 catch(exception& e) {
535 m->errorOut(e, "Tree", "getSubTree");
539 /*****************************************************************/
540 //assumes nameMap contains unique names as key or is empty.
541 //assumes numLeaves defined in tree constructor equals size of seqsToInclude and seqsToInclude only contains unique seqs.
542 int Tree::getSubTree(Tree* copy, vector<string> seqsToInclude, map<string, string> nameMap) {
545 if (numLeaves != seqsToInclude.size()) { m->mothurOut("[ERROR]: numLeaves does not equal numUniques, cannot create subtree.\n"); m->control_pressed = true; return 0; }
547 getSubTree(copy, seqsToInclude);
548 if (nameMap.size() != 0) { addNamesToCounts(nameMap); }
550 //build the pGroups in non leaf nodes to be used in the parsimony calcs.
551 for (int i = numLeaves; i < numNodes; i++) {
552 if (m->control_pressed) { return 1; }
554 tree[i].pGroups = (mergeGroups(i));
555 tree[i].pcount = (mergeGcounts(i));
560 catch(exception& e) {
561 m->errorOut(e, "Tree", "getSubTree");
565 /*****************************************************************/
566 int Tree::populateNewTree(vector<Node>& oldtree, int node, int& index) {
569 if (oldtree[node].getLChild() != -1) {
570 int rc = populateNewTree(oldtree, oldtree[node].getLChild(), index);
571 int lc = populateNewTree(oldtree, oldtree[node].getRChild(), index);
573 tree[index].setChildren(lc, rc);
574 tree[rc].setParent(index);
575 tree[lc].setParent(index);
577 tree[index].setBranchLength(oldtree[node].getBranchLength());
578 tree[rc].setBranchLength(oldtree[oldtree[node].getLChild()].getBranchLength());
579 tree[lc].setBranchLength(oldtree[oldtree[node].getRChild()].getBranchLength());
582 }else { //you are a leaf
583 int indexInNewTree = tmap->getIndex(oldtree[node].getName());
584 return indexInNewTree;
587 catch(exception& e) {
588 m->errorOut(e, "Tree", "populateNewTree");
592 /*****************************************************************/
593 void Tree::getCopy(Tree* copy) {
596 //for each node in the tree copy its info
597 for (int i = 0; i < numNodes; i++) {
599 tree[i].setName(copy->tree[i].getName());
602 tree[i].setGroup(copy->tree[i].getGroup());
605 tree[i].setBranchLength(copy->tree[i].getBranchLength());
608 tree[i].setParent(copy->tree[i].getParent());
611 tree[i].setChildren(copy->tree[i].getLChild(), copy->tree[i].getRChild());
613 //copy index in node and tmap
614 tree[i].setIndex(copy->tree[i].getIndex());
615 setIndex(copy->tree[i].getName(), getIndex(copy->tree[i].getName()));
618 tree[i].pGroups = copy->tree[i].pGroups;
621 tree[i].pcount = copy->tree[i].pcount;
624 groupNodeInfo = copy->groupNodeInfo;
627 catch(exception& e) {
628 m->errorOut(e, "Tree", "getCopy");
632 /*****************************************************************/
633 //returns a map with a groupname and the number of times that group was seen in the children
634 //for instance if your children are white and black then it would return a map with 2 entries
635 // p[white] = 1 and p[black] = 1. Now go up a level and merge that with a node who has p[white] = 1
636 //and you get p[white] = 2, p[black] = 1, but you erase the p[black] because you have a p value higher than 1.
638 map<string, int> Tree::mergeGroups(int i) {
640 int lc = tree[i].getLChild();
641 int rc = tree[i].getRChild();
643 //set parsimony groups to left child
644 map<string,int> parsimony = tree[lc].pGroups;
648 //look at right child groups and update maxPars if right child has something higher for that group.
649 for(it=tree[rc].pGroups.begin();it!=tree[rc].pGroups.end();it++){
650 it2 = parsimony.find(it->first);
651 if (it2 != parsimony.end()) {
652 parsimony[it->first]++;
654 parsimony[it->first] = 1;
657 if(parsimony[it->first] > maxPars){
658 maxPars = parsimony[it->first];
662 // this is true if right child had a greater parsimony for a certain group
664 //erase all the groups that are only 1 because you found something with 2.
665 for(it=parsimony.begin();it!=parsimony.end();){
667 parsimony.erase(it++);
670 //set one remaining groups to 1
671 //so with our above example p[white] = 2 would be left and it would become p[white] = 1
672 for(it=parsimony.begin();it!=parsimony.end();it++){
673 parsimony[it->first] = 1;
680 catch(exception& e) {
681 m->errorOut(e, "Tree", "mergeGroups");
685 /*****************************************************************/
686 //returns a map with a groupname and the number of times that group was seen in the children
687 //for instance if your children are white and black then it would return a map with 2 entries
688 // p[white] = 1 and p[black] = 1. Now go up a level and merge that with a node who has p[white] = 1
689 //and you get p[white] = 2, p[black] = 1, but you erase the p[black] because you have a p value higher than 1.
691 map<string, int> Tree::mergeUserGroups(int i, vector<string> g) {
694 int lc = tree[i].getLChild();
695 int rc = tree[i].getRChild();
697 //loop through nodes groups removing the ones the user doesn't want
698 for(it=tree[lc].pGroups.begin();it!=tree[lc].pGroups.end();){
699 if (m->inUsersGroups(it->first, g) != true) {
700 tree[lc].pGroups.erase(it++);
704 //loop through nodes groups removing the ones the user doesn't want
705 for(it=tree[rc].pGroups.begin();it!=tree[rc].pGroups.end();){
706 if (m->inUsersGroups(it->first, g) != true) {
707 tree[rc].pGroups.erase(it++);
711 //set parsimony groups to left child
712 map<string,int> parsimony = tree[lc].pGroups;
716 //look at right child groups and update maxPars if right child has something higher for that group.
717 for(it=tree[rc].pGroups.begin();it!=tree[rc].pGroups.end();it++){
718 it2 = parsimony.find(it->first);
719 if (it2 != parsimony.end()) {
720 parsimony[it->first]++;
722 parsimony[it->first] = 1;
725 if(parsimony[it->first] > maxPars){
726 maxPars = parsimony[it->first];
730 // this is true if right child had a greater parsimony for a certain group
732 //erase all the groups that are only 1 because you found something with 2.
733 for(it=parsimony.begin();it!=parsimony.end();){
735 parsimony.erase(it++);
739 for(it=parsimony.begin();it!=parsimony.end();it++){
740 parsimony[it->first] = 1;
746 catch(exception& e) {
747 m->errorOut(e, "Tree", "mergeUserGroups");
753 /**************************************************************************************************/
755 map<string,int> Tree::mergeGcounts(int position) {
757 map<string,int>::iterator pos;
759 int lc = tree[position].getLChild();
760 int rc = tree[position].getRChild();
762 map<string,int> sum = tree[lc].pcount;
764 for(it=tree[rc].pcount.begin();it!=tree[rc].pcount.end();it++){
765 sum[it->first] += it->second;
769 catch(exception& e) {
770 m->errorOut(e, "Tree", "mergeGcounts");
774 /**************************************************************************************************/
775 void Tree::randomLabels(vector<string> g) {
778 //initialize groupNodeInfo
779 for (int i = 0; i < (tmap->getNamesOfGroups()).size(); i++) {
780 groupNodeInfo[(tmap->getNamesOfGroups())[i]].resize(0);
783 for(int i = 0; i < numLeaves; i++){
785 //get random index to switch with
786 z = int((float)(i+1) * (float)(rand()) / ((float)RAND_MAX+1.0));
788 //you only want to randomize the nodes that are from a group the user wants analyzed, so
789 //if either of the leaf nodes you are about to switch are not in the users groups then you don't want to switch them.
792 treez = m->inUsersGroups(tree[z].getGroup(), g);
793 treei = m->inUsersGroups(tree[i].getGroup(), g);
795 if ((treez == true) && (treei == true)) {
796 //switches node i and node z's info.
797 map<string,int> lib_hold = tree[z].pGroups;
798 tree[z].pGroups = (tree[i].pGroups);
799 tree[i].pGroups = (lib_hold);
801 vector<string> zgroup = tree[z].getGroup();
802 tree[z].setGroup(tree[i].getGroup());
803 tree[i].setGroup(zgroup);
805 string zname = tree[z].getName();
806 tree[z].setName(tree[i].getName());
807 tree[i].setName(zname);
809 map<string,int> gcount_hold = tree[z].pcount;
810 tree[z].pcount = (tree[i].pcount);
811 tree[i].pcount = (gcount_hold);
814 for (int k = 0; k < (tree[i].getGroup()).size(); k++) { groupNodeInfo[(tree[i].getGroup())[k]].push_back(i); }
815 for (int k = 0; k < (tree[z].getGroup()).size(); k++) { groupNodeInfo[(tree[z].getGroup())[k]].push_back(z); }
818 catch(exception& e) {
819 m->errorOut(e, "Tree", "randomLabels");
823 /**************************************************************************************************/
824 void Tree::randomBlengths() {
826 for(int i=numNodes-1;i>=0;i--){
827 int z = int((float)(i+1) * (float)(rand()) / ((float)RAND_MAX+1.0));
829 float bl_hold = tree[z].getBranchLength();
830 tree[z].setBranchLength(tree[i].getBranchLength());
831 tree[i].setBranchLength(bl_hold);
834 catch(exception& e) {
835 m->errorOut(e, "Tree", "randomBlengths");
839 /*************************************************************************************************/
840 void Tree::assembleRandomUnifracTree(vector<string> g) {
842 map<string, string> empty;
845 /*************************************************************************************************/
846 void Tree::assembleRandomUnifracTree(string groupA, string groupB) {
847 vector<string> temp; temp.push_back(groupA); temp.push_back(groupB);
849 map<string, string> empty;
853 /*************************************************************************************************/
854 //for now it's just random topology but may become random labels as well later that why this is such a simple function now...
855 void Tree::assembleRandomTree() {
857 map<string, string> empty;
860 /**************************************************************************************************/
862 void Tree::randomTopology() {
864 for(int i=0;i<numNodes;i++){
865 tree[i].setParent(-1);
867 for(int i=numLeaves;i<numNodes;i++){
868 tree[i].setChildren(-1, -1);
871 for(int i=numLeaves;i<numNodes;i++){
873 int rnd_index1, rnd_index2;
875 rnd_index1 = (int)(((double)rand() / (double) RAND_MAX)*i);
876 if(tree[rnd_index1].getParent() == -1){escape = 1;}
881 rnd_index2 = (int)(((double)rand() / (double) RAND_MAX)*i);
882 if(rnd_index2 != rnd_index1 && tree[rnd_index2].getParent() == -1){
887 tree[i].setChildren(rnd_index1,rnd_index2);
888 tree[i].setParent(-1);
889 tree[rnd_index1].setParent(i);
890 tree[rnd_index2].setParent(i);
893 catch(exception& e) {
894 m->errorOut(e, "Tree", "randomTopology");
898 /*****************************************************************/
899 void Tree::print(ostream& out) {
901 int root = findRoot();
902 printBranch(root, out, "branch");
905 catch(exception& e) {
906 m->errorOut(e, "Tree", "print");
910 /*****************************************************************/
911 void Tree::print(ostream& out, map<string, string> nameMap) {
913 int root = findRoot();
914 printBranch(root, out, nameMap);
917 catch(exception& e) {
918 m->errorOut(e, "Tree", "print");
922 /*****************************************************************/
923 void Tree::print(ostream& out, string mode) {
925 int root = findRoot();
926 printBranch(root, out, mode);
929 catch(exception& e) {
930 m->errorOut(e, "Tree", "print");
934 /*****************************************************************/
935 // This prints out the tree in Newick form.
936 void Tree::createNewickFile(string f) {
938 int root = findRoot();
942 m->openOutputFile(filename, out);
944 printBranch(root, out, "branch");
946 // you are at the end of the tree
950 catch(exception& e) {
951 m->errorOut(e, "Tree", "createNewickFile");
956 /*****************************************************************/
957 //This function finds the index of the root node.
959 int Tree::findRoot() {
961 for (int i = 0; i < numNodes; i++) {
963 if (tree[i].getParent() == -1) { return i; }
964 //cout << "i = " << i << endl;
965 //cout << "i's parent = " << tree[i].getParent() << endl;
969 catch(exception& e) {
970 m->errorOut(e, "Tree", "findRoot");
974 /*****************************************************************/
975 void Tree::printBranch(int node, ostream& out, map<string, string> names) {
978 // you are not a leaf
979 if (tree[node].getLChild() != -1) {
981 printBranch(tree[node].getLChild(), out, names);
983 printBranch(tree[node].getRChild(), out, names);
986 //if there is a branch length then print it
987 if (tree[node].getBranchLength() != -1) {
988 out << ":" << tree[node].getBranchLength();
991 }else { //you are a leaf
992 map<string, string>::iterator itNames = names.find(tree[node].getName());
994 string outputString = "";
995 if (itNames != names.end()) {
997 vector<string> dupNames;
998 m->splitAtComma((itNames->second), dupNames);
1000 if (dupNames.size() == 1) {
1001 outputString += tree[node].getName();
1002 if (tree[node].getBranchLength() != -1) {
1003 outputString += ":" + toString(tree[node].getBranchLength());
1006 outputString += "(";
1008 for (int u = 0; u < dupNames.size()-1; u++) {
1009 outputString += dupNames[u];
1011 if (tree[node].getBranchLength() != -1) {
1012 outputString += ":" + toString(0.0);
1014 outputString += ",";
1017 outputString += dupNames[dupNames.size()-1];
1018 if (tree[node].getBranchLength() != -1) {
1019 outputString += ":" + toString(0.0);
1022 outputString += ")";
1023 if (tree[node].getBranchLength() != -1) {
1024 outputString += ":" + toString(tree[node].getBranchLength());
1028 outputString = tree[node].getName();
1029 //if there is a branch length then print it
1030 if (tree[node].getBranchLength() != -1) {
1031 outputString += ":" + toString(tree[node].getBranchLength());
1034 m->mothurOut("[ERROR]: " + tree[node].getName() + " is not in your namefile, please correct."); m->mothurOutEndLine();
1037 out << outputString;
1041 catch(exception& e) {
1042 m->errorOut(e, "Tree", "printBranch");
1046 /*****************************************************************/
1047 void Tree::printBranch(int node, ostream& out, string mode) {
1050 // you are not a leaf
1051 if (tree[node].getLChild() != -1) {
1053 printBranch(tree[node].getLChild(), out, mode);
1055 printBranch(tree[node].getRChild(), out, mode);
1057 if (mode == "branch") {
1058 //if there is a branch length then print it
1059 if (tree[node].getBranchLength() != -1) {
1060 out << ":" << tree[node].getBranchLength();
1062 }else if (mode == "boot") {
1063 //if there is a label then print it
1064 if (tree[node].getLabel() != -1) {
1065 out << tree[node].getLabel();
1067 }else if (mode == "both") {
1068 if (tree[node].getLabel() != -1) {
1069 out << tree[node].getLabel();
1071 //if there is a branch length then print it
1072 if (tree[node].getBranchLength() != -1) {
1073 out << ":" << tree[node].getBranchLength();
1076 }else { //you are a leaf
1077 string leafGroup = tmap->getGroup(tree[node].getName());
1079 if (mode == "branch") {
1081 //if there is a branch length then print it
1082 if (tree[node].getBranchLength() != -1) {
1083 out << ":" << tree[node].getBranchLength();
1085 }else if (mode == "boot") {
1087 //if there is a label then print it
1088 if (tree[node].getLabel() != -1) {
1089 out << tree[node].getLabel();
1091 }else if (mode == "both") {
1092 out << tree[node].getName();
1093 if (tree[node].getLabel() != -1) {
1094 out << tree[node].getLabel();
1096 //if there is a branch length then print it
1097 if (tree[node].getBranchLength() != -1) {
1098 out << ":" << tree[node].getBranchLength();
1104 catch(exception& e) {
1105 m->errorOut(e, "Tree", "printBranch");
1109 /*****************************************************************/
1110 void Tree::printBranch(int node, ostream& out, string mode, vector<Node>& theseNodes) {
1113 // you are not a leaf
1114 if (theseNodes[node].getLChild() != -1) {
1116 printBranch(theseNodes[node].getLChild(), out, mode);
1118 printBranch(theseNodes[node].getRChild(), out, mode);
1120 if (mode == "branch") {
1121 //if there is a branch length then print it
1122 if (theseNodes[node].getBranchLength() != -1) {
1123 out << ":" << theseNodes[node].getBranchLength();
1125 }else if (mode == "boot") {
1126 //if there is a label then print it
1127 if (theseNodes[node].getLabel() != -1) {
1128 out << theseNodes[node].getLabel();
1130 }else if (mode == "both") {
1131 if (theseNodes[node].getLabel() != -1) {
1132 out << theseNodes[node].getLabel();
1134 //if there is a branch length then print it
1135 if (theseNodes[node].getBranchLength() != -1) {
1136 out << ":" << theseNodes[node].getBranchLength();
1139 }else { //you are a leaf
1140 string leafGroup = tmap->getGroup(theseNodes[node].getName());
1142 if (mode == "branch") {
1144 //if there is a branch length then print it
1145 if (theseNodes[node].getBranchLength() != -1) {
1146 out << ":" << theseNodes[node].getBranchLength();
1148 }else if (mode == "boot") {
1150 //if there is a label then print it
1151 if (theseNodes[node].getLabel() != -1) {
1152 out << theseNodes[node].getLabel();
1154 }else if (mode == "both") {
1155 out << theseNodes[node].getName();
1156 if (theseNodes[node].getLabel() != -1) {
1157 out << theseNodes[node].getLabel();
1159 //if there is a branch length then print it
1160 if (theseNodes[node].getBranchLength() != -1) {
1161 out << ":" << theseNodes[node].getBranchLength();
1167 catch(exception& e) {
1168 m->errorOut(e, "Tree", "printBranch");
1172 /*****************************************************************/
1174 void Tree::printTree() {
1176 for(int i=0;i<numNodes;i++){
1178 tree[i].printNode();
1183 /*****************************************************************/
1184 //this code is a mess and should be rethought...-slw
1185 void Tree::parseTreeFile() {
1187 //only takes names from the first tree and assumes that all trees use the same names.
1189 string filename = m->getTreeFile();
1190 ifstream filehandle;
1191 m->openInputFile(filename, filehandle);
1196 //ifyou are not a nexus file
1197 if((c = filehandle.peek()) != '#') {
1198 while((c = filehandle.peek()) != ';') {
1199 while ((c = filehandle.peek()) != ';') {
1200 // get past comments
1207 if((c == '(') && (comment != 1)){ break; }
1211 done = readTreeString(filehandle);
1212 if (done == 0) { break; }
1214 //ifyou are a nexus file
1215 }else if((c = filehandle.peek()) == '#') {
1218 // get past comments
1219 while(holder != "translate" && holder != "Translate"){
1220 if(holder == "[" || holder == "[!"){
1226 filehandle >> holder;
1228 //if there is no translate then you must read tree string otherwise use translate to get names
1229 if((holder == "tree") && (comment != 1)){
1230 //pass over the "tree rep.6878900 = "
1231 while (((c = filehandle.get()) != '(') && ((c = filehandle.peek()) != EOF)) {;}
1233 if(c == EOF) { break; }
1234 filehandle.putback(c); //put back first ( of tree.
1235 done = readTreeString(filehandle);
1240 if (done == 0) { break; }
1243 //use nexus translation rather than parsing tree to save time
1244 if((holder == "translate") || (holder == "Translate")) {
1246 string number, name, h;
1247 h = ""; // so it enters the loop the first time
1248 while((h != ";") && (number != ";")) {
1249 filehandle >> number;
1252 //c = , until done with translation then c = ;
1253 h = name.substr(name.length()-1, name.length());
1254 name.erase(name.end()-1); //erase the comma
1255 m->Treenames.push_back(number);
1257 if(number == ";") { m->Treenames.pop_back(); } //in case ';' from translation is on next line instead of next to last name
1262 //for (int i = 0; i < globaldata->Treenames.size(); i++) {
1263 //cout << globaldata->Treenames[i] << endl; }
1264 //cout << globaldata->Treenames.size() << endl;
1266 catch(exception& e) {
1267 m->errorOut(e, "Tree", "parseTreeFile");
1271 /*******************************************************/
1273 /*******************************************************/
1274 int Tree::readTreeString(ifstream& filehandle) {
1279 while((c = filehandle.peek()) != ';') {
1281 //cout << " at beginning of while " << k << endl;
1283 //to pass over labels in trees
1285 while((c!=',') && (c != -1) && (c!= ':') && (c!=';')){ c=filehandle.get(); }
1286 filehandle.putback(c);
1288 if(c == ';') { return 0; }
1289 if(c == -1) { return 0; }
1291 if((c != '(') && (c != ')') && (c != ',') && (c != ':') && (c != '\n') && (c != '\t') && (c != 32)) { //32 is space
1293 c = filehandle.get();
1295 //cout << k << endl;
1296 while ((c != '(') && (c != ')') && (c != ',') && (c != ':') && (c != '\n') && (c != 32) && (c != '\t')) {
1298 c = filehandle.get();
1300 //cout << " in name while " << k << endl;
1303 //cout << "name = " << name << endl;
1304 m->Treenames.push_back(name);
1305 filehandle.putback(c);
1307 //cout << " after putback" << k << endl;
1310 if(c == ':') { //read until you reach the end of the branch length
1311 while ((c != '(') && (c != ')') && (c != ',') && (c != ';') && (c != '\n') && (c != '\t') && (c != 32)) {
1312 c = filehandle.get();
1314 //cout << " in branch while " << k << endl;
1316 filehandle.putback(c);
1319 c = filehandle.get();
1321 //cout << " here after get " << k << endl;
1322 if(c == ';') { return 0; }
1323 if(c == ')') { filehandle.putback(c); }
1325 //cout << k << endl;
1330 catch(exception& e) {
1331 m->errorOut(e, "Tree", "readTreeString");
1336 /*******************************************************/
1338 /*******************************************************/