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();
34 parseTreeFile(); m->runParse = false;
37 m->errorOut(e, "Tree", "Tree - just parse");
41 /*****************************************************************/
42 Tree::Tree(TreeMap* t) : tmap(t) {
44 m = MothurOut::getInstance();
46 if (m->runParse == true) { parseTreeFile(); m->runParse = false; }
47 //for(int i = 0; i < globaldata->Treenames.size(); i++) { cout << i << '\t' << globaldata->Treenames[i] << endl; }
48 numLeaves = m->Treenames.size();
49 numNodes = 2*numLeaves - 1;
51 tree.resize(numNodes);
53 //initialize groupNodeInfo
54 for (int i = 0; i < (tmap->getNamesOfGroups()).size(); i++) {
55 groupNodeInfo[(tmap->getNamesOfGroups())[i]].resize(0);
58 //initialize tree with correct number of nodes, name and group info.
59 for (int i = 0; i < numNodes; i++) {
60 //initialize leaf nodes
61 if (i <= (numLeaves-1)) {
62 tree[i].setName(m->Treenames[i]);
65 string group = tmap->getGroup(m->Treenames[i]);
67 vector<string> tempGroups; tempGroups.push_back(group);
68 tree[i].setGroup(tempGroups);
69 groupNodeInfo[group].push_back(i);
71 //set pcount and pGroup for groupname to 1.
72 tree[i].pcount[group] = 1;
73 tree[i].pGroups[group] = 1;
75 //Treemap knows name, group and index to speed up search
76 tmap->setIndex(m->Treenames[i], i);
78 //intialize non leaf nodes
79 }else if (i > (numLeaves-1)) {
81 vector<string> tempGroups;
82 tree[i].setGroup(tempGroups);
88 m->errorOut(e, "Tree", "Tree");
92 /*****************************************************************/
93 Tree::Tree(TreeMap* t, vector< vector<double> >& sims) : tmap(t) {
95 m = MothurOut::getInstance();
97 if (m->runParse == true) { parseTreeFile(); m->runParse = false; }
98 //for(int i = 0; i < globaldata->Treenames.size(); i++) { cout << i << '\t' << globaldata->Treenames[i] << endl; }
99 numLeaves = m->Treenames.size();
100 numNodes = 2*numLeaves - 1;
102 tree.resize(numNodes);
104 //initialize groupNodeInfo
105 for (int i = 0; i < (tmap->getNamesOfGroups()).size(); i++) {
106 groupNodeInfo[(tmap->getNamesOfGroups())[i]].resize(0);
109 //initialize tree with correct number of nodes, name and group info.
110 for (int i = 0; i < numNodes; i++) {
111 //initialize leaf nodes
112 if (i <= (numLeaves-1)) {
113 tree[i].setName(m->Treenames[i]);
116 string group = tmap->getGroup(m->Treenames[i]);
118 vector<string> tempGroups; tempGroups.push_back(group);
119 tree[i].setGroup(tempGroups);
120 groupNodeInfo[group].push_back(i);
122 //set pcount and pGroup for groupname to 1.
123 tree[i].pcount[group] = 1;
124 tree[i].pGroups[group] = 1;
126 //Treemap knows name, group and index to speed up search
127 tmap->setIndex(m->Treenames[i], i);
129 //intialize non leaf nodes
130 }else if (i > (numLeaves-1)) {
132 vector<string> tempGroups;
133 tree[i].setGroup(tempGroups);
137 //build tree from matrix
139 map<int, int> indexes; //maps row in simMatrix to vector index in the tree
140 int numGroups = (tmap->getNamesOfGroups()).size();
141 for (int g = 0; g < numGroups; g++) { indexes[g] = g; }
143 //do merges and create tree structure by setting parents and children
144 //there are numGroups - 1 merges to do
145 for (int i = 0; i < (numGroups - 1); i++) {
146 float largest = -1000.0;
148 if (m->control_pressed) { break; }
151 //find largest value in sims matrix by searching lower triangle
152 for (int j = 1; j < sims.size(); j++) {
153 for (int k = 0; k < j; k++) {
154 if (sims[j][k] > largest) { largest = sims[j][k]; row = j; column = k; }
158 //set non-leaf node info and update leaves to know their parents
160 tree[numGroups + i].setChildren(indexes[row], indexes[column]);
163 tree[indexes[row]].setParent(numGroups + i);
164 tree[indexes[column]].setParent(numGroups + i);
166 //blength = distance / 2;
167 float blength = ((1.0 - largest) / 2);
170 tree[indexes[row]].setBranchLength(blength - tree[indexes[row]].getLengthToLeaves());
171 tree[indexes[column]].setBranchLength(blength - tree[indexes[column]].getLengthToLeaves());
173 //set your length to leaves to your childs length plus branchlength
174 tree[numGroups + i].setLengthToLeaves(tree[indexes[row]].getLengthToLeaves() + tree[indexes[row]].getBranchLength());
178 indexes[row] = numGroups+i;
179 indexes[column] = numGroups+i;
181 //remove highest value that caused the merge.
182 sims[row][column] = -1000.0;
183 sims[column][row] = -1000.0;
185 //merge values in simsMatrix
186 for (int n = 0; n < sims.size(); n++) {
187 //row becomes merge of 2 groups
188 sims[row][n] = (sims[row][n] + sims[column][n]) / 2;
189 sims[n][row] = sims[row][n];
191 sims[column][n] = -1000.0;
192 sims[n][column] = -1000.0;
196 //adjust tree to make sure root to tip length is .5
197 int root = findRoot();
198 tree[root].setBranchLength((0.5 - tree[root].getLengthToLeaves()));
200 catch(exception& e) {
201 m->errorOut(e, "Tree", "Tree");
205 /*****************************************************************/
207 /*****************************************************************/
208 void Tree::addNamesToCounts(map<string, string> nameMap) {
210 //ex. seq1 seq2,seq3,se4
216 //before this function seq1.pcount = pasture -> 1
217 //after seq1.pcount = pasture -> 2, forest -> 1, ocean -> 1
219 //before this function seq1.pgroups = pasture -> 1
220 //after seq1.pgroups = pasture -> 1 since that is the dominant group
223 //go through each leaf and update its pcounts and pgroups
227 for (int i = 0; i < numLeaves; i++) {
229 string name = tree[i].getName();
231 map<string, string>::iterator itNames = nameMap.find(name);
233 if (itNames == nameMap.end()) { m->mothurOut(name + " is not in your name file, please correct."); m->mothurOutEndLine(); exit(1); }
235 vector<string> dupNames;
236 m->splitAtComma(nameMap[name], dupNames);
238 map<string, int>::iterator itCounts;
240 set<string> groupsAddedForThisNode;
241 for (int j = 0; j < dupNames.size(); j++) {
243 string group = tmap->getGroup(dupNames[j]);
245 if (dupNames[j] != name) {//you already added yourself in the constructor
247 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
250 itCounts = tree[i].pcount.find(group);
251 if (itCounts == tree[i].pcount.end()) { //new group, add it
252 tree[i].pcount[group] = 1;
254 tree[i].pcount[group]++;
258 itCounts = tree[i].pGroups.find(group);
259 if (itCounts == tree[i].pGroups.end()) { //new group, add it
260 tree[i].pGroups[group] = 1;
262 tree[i].pGroups[group]++;
266 if(tree[i].pGroups[group] > maxPars){
267 maxPars = tree[i].pGroups[group];
269 }else { groupsAddedForThisNode.insert(group); } //add it so you don't add it to groupNodeInfo again
272 if (maxPars > 1) { //then we have some more dominant groups
273 //erase all the groups that are less than maxPars because you found a more dominant group.
274 for(it=tree[i].pGroups.begin();it!=tree[i].pGroups.end();){
275 if(it->second < maxPars){
276 tree[i].pGroups.erase(it++);
279 //set one remaining groups to 1
280 for(it=tree[i].pGroups.begin();it!=tree[i].pGroups.end();it++){
281 tree[i].pGroups[it->first] = 1;
285 //update groups to reflect all the groups this node represents
286 vector<string> nodeGroups;
287 map<string, int>::iterator itGroups;
288 for (itGroups = tree[i].pcount.begin(); itGroups != tree[i].pcount.end(); itGroups++) {
289 nodeGroups.push_back(itGroups->first);
291 tree[i].setGroup(nodeGroups);
297 //cout << "addNamesToCounts\t" << (B - A) / CLOCKS_PER_SEC << endl;
300 catch(exception& e) {
301 m->errorOut(e, "Tree", "addNamesToCounts");
305 /*****************************************************************/
306 int Tree::getIndex(string searchName) {
308 //Treemap knows name, group and index to speed up search
309 // getIndex function will return the vector index or -1 if seq is not found.
310 int index = tmap->getIndex(searchName);
314 catch(exception& e) {
315 m->errorOut(e, "Tree", "getIndex");
319 /*****************************************************************/
321 void Tree::setIndex(string searchName, int index) {
323 //set index in treemap
324 tmap->setIndex(searchName, index);
326 catch(exception& e) {
327 m->errorOut(e, "Tree", "setIndex");
331 /*****************************************************************/
332 int Tree::assembleTree() {
336 //if user has given a names file we want to include that info in the pgroups and pcount info.
337 if(m->names.size() != 0) { addNamesToCounts(m->names); }
339 //build the pGroups in non leaf nodes to be used in the parsimony calcs.
340 for (int i = numLeaves; i < numNodes; i++) {
341 if (m->control_pressed) { return 1; }
343 tree[i].pGroups = (mergeGroups(i));
344 tree[i].pcount = (mergeGcounts(i));
347 //cout << "assembleTree\t" << (B-A) / CLOCKS_PER_SEC << endl;
350 catch(exception& e) {
351 m->errorOut(e, "Tree", "assembleTree");
355 /*****************************************************************/
356 int Tree::assembleTree(string n) {
359 //build the pGroups in non leaf nodes to be used in the parsimony calcs.
360 for (int i = numLeaves; i < numNodes; i++) {
361 if (m->control_pressed) { return 1; }
363 tree[i].pGroups = (mergeGroups(i));
364 tree[i].pcount = (mergeGcounts(i));
367 //cout << "assembleTree\t" << (B-A) / CLOCKS_PER_SEC << endl;
370 catch(exception& e) {
371 m->errorOut(e, "Tree", "assembleTree");
375 /*****************************************************************/
376 //assumes leaf node names are in groups and no names file - used by indicator command
377 void Tree::getSubTree(Tree* Ctree, vector<string> Groups) {
380 //copy Tree since we are going to destroy it
381 Tree* copy = new Tree(tmap);
382 copy->getCopy(Ctree);
383 copy->assembleTree("nonames");
385 //we want to select some of the leaf nodes to create the output tree
386 //go through the input Tree starting at parents of leaves
387 for (int i = 0; i < numNodes; i++) {
389 //initialize leaf nodes
390 if (i <= (numLeaves-1)) {
391 tree[i].setName(Groups[i]);
394 string group = tmap->getGroup(Groups[i]);
395 vector<string> tempGroups; tempGroups.push_back(group);
396 tree[i].setGroup(tempGroups);
397 groupNodeInfo[group].push_back(i);
399 //set pcount and pGroup for groupname to 1.
400 tree[i].pcount[group] = 1;
401 tree[i].pGroups[group] = 1;
403 //Treemap knows name, group and index to speed up search
404 tmap->setIndex(Groups[i], i);
406 //intialize non leaf nodes
407 }else if (i > (numLeaves-1)) {
409 vector<string> tempGroups;
410 tree[i].setGroup(tempGroups);
414 set<int> removedLeaves;
415 for (int i = 0; i < copy->getNumLeaves(); i++) {
417 if (removedLeaves.count(i) == 0) {
420 int parent = copy->tree[i].getParent();
424 if (m->inUsersGroups(copy->tree[i].getName(), Groups)) {
425 //find my siblings name
426 int parentRC = copy->tree[parent].getRChild();
427 int parentLC = copy->tree[parent].getLChild();
429 //if I am the right child, then my sib is the left child
430 int sibIndex = parentRC;
431 if (parentRC == i) { sibIndex = parentLC; }
433 string sibsName = copy->tree[sibIndex].getName();
435 //if yes, is my sibling
436 if ((m->inUsersGroups(sibsName, Groups)) || (sibsName == "")) {
437 //we both are okay no trimming required
439 //i am, my sib is not, so remove sib by setting my parent to my grandparent
440 int grandparent = copy->tree[parent].getParent();
441 int grandparentLC = copy->tree[grandparent].getLChild();
442 int grandparentRC = copy->tree[grandparent].getRChild();
444 //whichever of my granparents children was my parent now equals me
445 if (grandparentLC == parent) { grandparentLC = i; }
446 else { grandparentRC = i; }
448 copy->tree[i].setParent(grandparent);
449 copy->tree[i].setBranchLength((copy->tree[i].getBranchLength()+copy->tree[parent].getBranchLength()));
450 if (grandparent != -1) {
451 copy->tree[grandparent].setChildren(grandparentLC, grandparentRC);
453 removedLeaves.insert(sibIndex);
456 //find my siblings name
457 int parentRC = copy->tree[parent].getRChild();
458 int parentLC = copy->tree[parent].getLChild();
460 //if I am the right child, then my sib is the left child
461 int sibIndex = parentRC;
462 if (parentRC == i) { sibIndex = parentLC; }
464 string sibsName = copy->tree[sibIndex].getName();
466 //if no is my sibling
467 if ((m->inUsersGroups(sibsName, Groups)) || (sibsName == "")) {
468 //i am not, but my sib is
469 int grandparent = copy->tree[parent].getParent();
470 int grandparentLC = copy->tree[grandparent].getLChild();
471 int grandparentRC = copy->tree[grandparent].getRChild();
473 //whichever of my granparents children was my parent now equals my sib
474 if (grandparentLC == parent) { grandparentLC = sibIndex; }
475 else { grandparentRC = sibIndex; }
477 copy->tree[sibIndex].setParent(grandparent);
478 copy->tree[sibIndex].setBranchLength((copy->tree[sibIndex].getBranchLength()+copy->tree[parent].getBranchLength()));
479 if (grandparent != -1) {
480 copy->tree[grandparent].setChildren(grandparentLC, grandparentRC);
482 removedLeaves.insert(i);
484 //neither of us are, so we want to eliminate ourselves and our parent
485 //so set our parents sib to our great-grandparent
486 int parent = copy->tree[i].getParent();
487 int grandparent = copy->tree[parent].getParent();
489 if (grandparent != -1) {
490 int greatgrandparent = copy->tree[grandparent].getParent();
491 int greatgrandparentLC, greatgrandparentRC;
492 if (greatgrandparent != -1) {
493 greatgrandparentLC = copy->tree[greatgrandparent].getLChild();
494 greatgrandparentRC = copy->tree[greatgrandparent].getRChild();
497 int grandparentLC = copy->tree[grandparent].getLChild();
498 int grandparentRC = copy->tree[grandparent].getRChild();
500 parentsSibIndex = grandparentLC;
501 if (grandparentLC == parent) { parentsSibIndex = grandparentRC; }
503 //whichever of my greatgrandparents children was my grandparent
504 if (greatgrandparentLC == grandparent) { greatgrandparentLC = parentsSibIndex; }
505 else { greatgrandparentRC = parentsSibIndex; }
507 copy->tree[parentsSibIndex].setParent(greatgrandparent);
508 copy->tree[parentsSibIndex].setBranchLength((copy->tree[parentsSibIndex].getBranchLength()+copy->tree[grandparent].getBranchLength()));
509 if (greatgrandparent != -1) {
510 copy->tree[greatgrandparent].setChildren(greatgrandparentLC, greatgrandparentRC);
513 copy->tree[parent].setParent(-1);
514 //cout << "issues with making subtree" << endl;
516 removedLeaves.insert(sibIndex);
517 removedLeaves.insert(i);
525 for (int i = 0; i < copy->getNumNodes(); i++) {
527 if (copy->tree[i].getParent() == -1) { root = i; break; }
530 int nextSpot = numLeaves;
531 populateNewTree(copy->tree, root, nextSpot);
535 catch(exception& e) {
536 m->errorOut(e, "Tree", "getSubTree");
540 /*****************************************************************/
541 //assumes nameMap contains unique names as key or is empty.
542 //assumes numLeaves defined in tree constructor equals size of seqsToInclude and seqsToInclude only contains unique seqs.
543 int Tree::getSubTree(Tree* copy, vector<string> seqsToInclude, map<string, string> nameMap) {
546 if (numLeaves != seqsToInclude.size()) { m->mothurOut("[ERROR]: numLeaves does not equal numUniques, cannot create subtree.\n"); m->control_pressed = true; return 0; }
548 getSubTree(copy, seqsToInclude);
549 if (nameMap.size() != 0) { addNamesToCounts(nameMap); }
551 //build the pGroups in non leaf nodes to be used in the parsimony calcs.
552 for (int i = numLeaves; i < numNodes; i++) {
553 if (m->control_pressed) { return 1; }
555 tree[i].pGroups = (mergeGroups(i));
556 tree[i].pcount = (mergeGcounts(i));
561 catch(exception& e) {
562 m->errorOut(e, "Tree", "getSubTree");
566 /*****************************************************************/
567 int Tree::populateNewTree(vector<Node>& oldtree, int node, int& index) {
570 if (oldtree[node].getLChild() != -1) {
571 int rc = populateNewTree(oldtree, oldtree[node].getLChild(), index);
572 int lc = populateNewTree(oldtree, oldtree[node].getRChild(), index);
574 tree[index].setChildren(lc, rc);
575 tree[rc].setParent(index);
576 tree[lc].setParent(index);
578 tree[index].setBranchLength(oldtree[node].getBranchLength());
579 tree[rc].setBranchLength(oldtree[oldtree[node].getLChild()].getBranchLength());
580 tree[lc].setBranchLength(oldtree[oldtree[node].getRChild()].getBranchLength());
583 }else { //you are a leaf
584 int indexInNewTree = tmap->getIndex(oldtree[node].getName());
585 return indexInNewTree;
588 catch(exception& e) {
589 m->errorOut(e, "Tree", "populateNewTree");
593 /*****************************************************************/
594 void Tree::getCopy(Tree* copy) {
597 //for each node in the tree copy its info
598 for (int i = 0; i < numNodes; i++) {
600 tree[i].setName(copy->tree[i].getName());
603 tree[i].setGroup(copy->tree[i].getGroup());
606 tree[i].setBranchLength(copy->tree[i].getBranchLength());
609 tree[i].setParent(copy->tree[i].getParent());
612 tree[i].setChildren(copy->tree[i].getLChild(), copy->tree[i].getRChild());
614 //copy index in node and tmap
615 tree[i].setIndex(copy->tree[i].getIndex());
616 setIndex(copy->tree[i].getName(), getIndex(copy->tree[i].getName()));
619 tree[i].pGroups = copy->tree[i].pGroups;
622 tree[i].pcount = copy->tree[i].pcount;
625 groupNodeInfo = copy->groupNodeInfo;
628 catch(exception& e) {
629 m->errorOut(e, "Tree", "getCopy");
633 /*****************************************************************/
634 //returns a map with a groupname and the number of times that group was seen in the children
635 //for instance if your children are white and black then it would return a map with 2 entries
636 // p[white] = 1 and p[black] = 1. Now go up a level and merge that with a node who has p[white] = 1
637 //and you get p[white] = 2, p[black] = 1, but you erase the p[black] because you have a p value higher than 1.
639 map<string, int> Tree::mergeGroups(int i) {
641 int lc = tree[i].getLChild();
642 int rc = tree[i].getRChild();
644 //set parsimony groups to left child
645 map<string,int> parsimony = tree[lc].pGroups;
649 //look at right child groups and update maxPars if right child has something higher for that group.
650 for(it=tree[rc].pGroups.begin();it!=tree[rc].pGroups.end();it++){
651 it2 = parsimony.find(it->first);
652 if (it2 != parsimony.end()) {
653 parsimony[it->first]++;
655 parsimony[it->first] = 1;
658 if(parsimony[it->first] > maxPars){
659 maxPars = parsimony[it->first];
663 // this is true if right child had a greater parsimony for a certain group
665 //erase all the groups that are only 1 because you found something with 2.
666 for(it=parsimony.begin();it!=parsimony.end();){
668 parsimony.erase(it++);
671 //set one remaining groups to 1
672 //so with our above example p[white] = 2 would be left and it would become p[white] = 1
673 for(it=parsimony.begin();it!=parsimony.end();it++){
674 parsimony[it->first] = 1;
681 catch(exception& e) {
682 m->errorOut(e, "Tree", "mergeGroups");
686 /*****************************************************************/
687 //returns a map with a groupname and the number of times that group was seen in the children
688 //for instance if your children are white and black then it would return a map with 2 entries
689 // p[white] = 1 and p[black] = 1. Now go up a level and merge that with a node who has p[white] = 1
690 //and you get p[white] = 2, p[black] = 1, but you erase the p[black] because you have a p value higher than 1.
692 map<string, int> Tree::mergeUserGroups(int i, vector<string> g) {
695 int lc = tree[i].getLChild();
696 int rc = tree[i].getRChild();
698 //loop through nodes groups removing the ones the user doesn't want
699 for(it=tree[lc].pGroups.begin();it!=tree[lc].pGroups.end();){
700 if (m->inUsersGroups(it->first, g) != true) {
701 tree[lc].pGroups.erase(it++);
705 //loop through nodes groups removing the ones the user doesn't want
706 for(it=tree[rc].pGroups.begin();it!=tree[rc].pGroups.end();){
707 if (m->inUsersGroups(it->first, g) != true) {
708 tree[rc].pGroups.erase(it++);
712 //set parsimony groups to left child
713 map<string,int> parsimony = tree[lc].pGroups;
717 //look at right child groups and update maxPars if right child has something higher for that group.
718 for(it=tree[rc].pGroups.begin();it!=tree[rc].pGroups.end();it++){
719 it2 = parsimony.find(it->first);
720 if (it2 != parsimony.end()) {
721 parsimony[it->first]++;
723 parsimony[it->first] = 1;
726 if(parsimony[it->first] > maxPars){
727 maxPars = parsimony[it->first];
731 // this is true if right child had a greater parsimony for a certain group
733 //erase all the groups that are only 1 because you found something with 2.
734 for(it=parsimony.begin();it!=parsimony.end();){
736 parsimony.erase(it++);
740 for(it=parsimony.begin();it!=parsimony.end();it++){
741 parsimony[it->first] = 1;
747 catch(exception& e) {
748 m->errorOut(e, "Tree", "mergeUserGroups");
754 /**************************************************************************************************/
756 map<string,int> Tree::mergeGcounts(int position) {
758 map<string,int>::iterator pos;
760 int lc = tree[position].getLChild();
761 int rc = tree[position].getRChild();
763 map<string,int> sum = tree[lc].pcount;
765 for(it=tree[rc].pcount.begin();it!=tree[rc].pcount.end();it++){
766 sum[it->first] += it->second;
770 catch(exception& e) {
771 m->errorOut(e, "Tree", "mergeGcounts");
775 /**************************************************************************************************/
776 void Tree::randomLabels(vector<string> g) {
779 //initialize groupNodeInfo
780 for (int i = 0; i < (tmap->getNamesOfGroups()).size(); i++) {
781 groupNodeInfo[(tmap->getNamesOfGroups())[i]].resize(0);
784 for(int i = 0; i < numLeaves; i++){
786 //get random index to switch with
787 z = int((float)(i+1) * (float)(rand()) / ((float)RAND_MAX+1.0));
789 //you only want to randomize the nodes that are from a group the user wants analyzed, so
790 //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.
793 treez = m->inUsersGroups(tree[z].getGroup(), g);
794 treei = m->inUsersGroups(tree[i].getGroup(), g);
796 if ((treez == true) && (treei == true)) {
797 //switches node i and node z's info.
798 map<string,int> lib_hold = tree[z].pGroups;
799 tree[z].pGroups = (tree[i].pGroups);
800 tree[i].pGroups = (lib_hold);
802 vector<string> zgroup = tree[z].getGroup();
803 tree[z].setGroup(tree[i].getGroup());
804 tree[i].setGroup(zgroup);
806 string zname = tree[z].getName();
807 tree[z].setName(tree[i].getName());
808 tree[i].setName(zname);
810 map<string,int> gcount_hold = tree[z].pcount;
811 tree[z].pcount = (tree[i].pcount);
812 tree[i].pcount = (gcount_hold);
815 for (int k = 0; k < (tree[i].getGroup()).size(); k++) { groupNodeInfo[(tree[i].getGroup())[k]].push_back(i); }
816 for (int k = 0; k < (tree[z].getGroup()).size(); k++) { groupNodeInfo[(tree[z].getGroup())[k]].push_back(z); }
819 catch(exception& e) {
820 m->errorOut(e, "Tree", "randomLabels");
824 /**************************************************************************************************/
825 void Tree::randomBlengths() {
827 for(int i=numNodes-1;i>=0;i--){
828 int z = int((float)(i+1) * (float)(rand()) / ((float)RAND_MAX+1.0));
830 float bl_hold = tree[z].getBranchLength();
831 tree[z].setBranchLength(tree[i].getBranchLength());
832 tree[i].setBranchLength(bl_hold);
835 catch(exception& e) {
836 m->errorOut(e, "Tree", "randomBlengths");
840 /*************************************************************************************************/
841 void Tree::assembleRandomUnifracTree(vector<string> g) {
843 assembleTree("noNameCounts");
845 /*************************************************************************************************/
846 void Tree::assembleRandomUnifracTree(string groupA, string groupB) {
848 vector<string> temp; temp.push_back(groupA); temp.push_back(groupB);
850 assembleTree("noNameCounts");
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() {
859 /**************************************************************************************************/
861 void Tree::randomTopology() {
863 for(int i=0;i<numNodes;i++){
864 tree[i].setParent(-1);
866 for(int i=numLeaves;i<numNodes;i++){
867 tree[i].setChildren(-1, -1);
870 for(int i=numLeaves;i<numNodes;i++){
872 int rnd_index1, rnd_index2;
874 rnd_index1 = (int)(((double)rand() / (double) RAND_MAX)*i);
875 if(tree[rnd_index1].getParent() == -1){escape = 1;}
880 rnd_index2 = (int)(((double)rand() / (double) RAND_MAX)*i);
881 if(rnd_index2 != rnd_index1 && tree[rnd_index2].getParent() == -1){
886 tree[i].setChildren(rnd_index1,rnd_index2);
887 tree[i].setParent(-1);
888 tree[rnd_index1].setParent(i);
889 tree[rnd_index2].setParent(i);
892 catch(exception& e) {
893 m->errorOut(e, "Tree", "randomTopology");
897 /*****************************************************************/
898 void Tree::print(ostream& out) {
900 int root = findRoot();
901 printBranch(root, out, "branch");
904 catch(exception& e) {
905 m->errorOut(e, "Tree", "print");
909 /*****************************************************************/
910 void Tree::print(ostream& out, string mode) {
912 int root = findRoot();
913 printBranch(root, out, mode);
916 catch(exception& e) {
917 m->errorOut(e, "Tree", "print");
921 /*****************************************************************/
922 // This prints out the tree in Newick form.
923 void Tree::createNewickFile(string f) {
925 int root = findRoot();
929 m->openOutputFile(filename, out);
931 printBranch(root, out, "branch");
933 // you are at the end of the tree
937 catch(exception& e) {
938 m->errorOut(e, "Tree", "createNewickFile");
943 /*****************************************************************/
944 //This function finds the index of the root node.
946 int Tree::findRoot() {
948 for (int i = 0; i < numNodes; i++) {
950 if (tree[i].getParent() == -1) { return i; }
951 //cout << "i = " << i << endl;
952 //cout << "i's parent = " << tree[i].getParent() << endl;
956 catch(exception& e) {
957 m->errorOut(e, "Tree", "findRoot");
961 /*****************************************************************/
962 void Tree::printBranch(int node, ostream& out, string mode) {
965 // you are not a leaf
966 if (tree[node].getLChild() != -1) {
968 printBranch(tree[node].getLChild(), out, mode);
970 printBranch(tree[node].getRChild(), out, mode);
972 if (mode == "branch") {
973 //if there is a branch length then print it
974 if (tree[node].getBranchLength() != -1) {
975 out << ":" << tree[node].getBranchLength();
977 }else if (mode == "boot") {
978 //if there is a label then print it
979 if (tree[node].getLabel() != -1) {
980 out << tree[node].getLabel();
982 }else if (mode == "both") {
983 if (tree[node].getLabel() != -1) {
984 out << tree[node].getLabel();
986 //if there is a branch length then print it
987 if (tree[node].getBranchLength() != -1) {
988 out << ":" << tree[node].getBranchLength();
990 }else if (mode == "deunique") {
991 //if there is a branch length then print it
992 if (tree[node].getBranchLength() != -1) {
993 out << ":" << tree[node].getBranchLength();
996 }else { //you are a leaf
997 string leafGroup = tmap->getGroup(tree[node].getName());
999 if (mode == "branch") {
1001 //if there is a branch length then print it
1002 if (tree[node].getBranchLength() != -1) {
1003 out << ":" << tree[node].getBranchLength();
1005 }else if (mode == "boot") {
1007 //if there is a label then print it
1008 if (tree[node].getLabel() != -1) {
1009 out << tree[node].getLabel();
1011 }else if (mode == "both") {
1012 out << tree[node].getName();
1013 if (tree[node].getLabel() != -1) {
1014 out << tree[node].getLabel();
1016 //if there is a branch length then print it
1017 if (tree[node].getBranchLength() != -1) {
1018 out << ":" << tree[node].getBranchLength();
1020 }else if (mode == "deunique") {
1021 map<string, string>::iterator itNames = m->names.find(tree[node].getName());
1023 string outputString = "";
1024 if (itNames != m->names.end()) {
1026 vector<string> dupNames;
1027 m->splitAtComma((itNames->second), dupNames);
1029 if (dupNames.size() == 1) {
1030 outputString += tree[node].getName();
1031 if (tree[node].getBranchLength() != -1) {
1032 outputString += ":" + toString(tree[node].getBranchLength());
1035 outputString += "(";
1037 for (int u = 0; u < dupNames.size()-1; u++) {
1038 outputString += dupNames[u];
1040 if (tree[node].getBranchLength() != -1) {
1041 outputString += ":" + toString(0.0);
1043 outputString += ",";
1046 outputString += dupNames[dupNames.size()-1];
1047 if (tree[node].getBranchLength() != -1) {
1048 outputString += ":" + toString(0.0);
1051 outputString += ")";
1052 if (tree[node].getBranchLength() != -1) {
1053 outputString += ":" + toString(tree[node].getBranchLength());
1057 outputString = tree[node].getName();
1058 //if there is a branch length then print it
1059 if (tree[node].getBranchLength() != -1) {
1060 outputString += ":" + toString(tree[node].getBranchLength());
1063 m->mothurOut("[ERROR]: " + tree[node].getName() + " is not in your namefile, please correct."); m->mothurOutEndLine();
1066 out << outputString;
1071 catch(exception& e) {
1072 m->errorOut(e, "Tree", "printBranch");
1076 /*****************************************************************/
1077 void Tree::printBranch(int node, ostream& out, string mode, vector<Node>& theseNodes) {
1080 // you are not a leaf
1081 if (theseNodes[node].getLChild() != -1) {
1083 printBranch(theseNodes[node].getLChild(), out, mode);
1085 printBranch(theseNodes[node].getRChild(), out, mode);
1087 if (mode == "branch") {
1088 //if there is a branch length then print it
1089 if (theseNodes[node].getBranchLength() != -1) {
1090 out << ":" << theseNodes[node].getBranchLength();
1092 }else if (mode == "boot") {
1093 //if there is a label then print it
1094 if (theseNodes[node].getLabel() != -1) {
1095 out << theseNodes[node].getLabel();
1097 }else if (mode == "both") {
1098 if (theseNodes[node].getLabel() != -1) {
1099 out << theseNodes[node].getLabel();
1101 //if there is a branch length then print it
1102 if (theseNodes[node].getBranchLength() != -1) {
1103 out << ":" << theseNodes[node].getBranchLength();
1106 }else { //you are a leaf
1107 string leafGroup = tmap->getGroup(theseNodes[node].getName());
1109 if (mode == "branch") {
1111 //if there is a branch length then print it
1112 if (theseNodes[node].getBranchLength() != -1) {
1113 out << ":" << theseNodes[node].getBranchLength();
1115 }else if (mode == "boot") {
1117 //if there is a label then print it
1118 if (theseNodes[node].getLabel() != -1) {
1119 out << theseNodes[node].getLabel();
1121 }else if (mode == "both") {
1122 out << theseNodes[node].getName();
1123 if (theseNodes[node].getLabel() != -1) {
1124 out << theseNodes[node].getLabel();
1126 //if there is a branch length then print it
1127 if (theseNodes[node].getBranchLength() != -1) {
1128 out << ":" << theseNodes[node].getBranchLength();
1134 catch(exception& e) {
1135 m->errorOut(e, "Tree", "printBranch");
1139 /*****************************************************************/
1141 void Tree::printTree() {
1143 for(int i=0;i<numNodes;i++){
1145 tree[i].printNode();
1150 /*****************************************************************/
1151 //this code is a mess and should be rethought...-slw
1152 void Tree::parseTreeFile() {
1154 //only takes names from the first tree and assumes that all trees use the same names.
1156 string filename = m->getTreeFile();
1157 ifstream filehandle;
1158 m->openInputFile(filename, filehandle);
1163 //ifyou are not a nexus file
1164 if((c = filehandle.peek()) != '#') {
1165 while((c = filehandle.peek()) != ';') {
1166 while ((c = filehandle.peek()) != ';') {
1167 // get past comments
1174 if((c == '(') && (comment != 1)){ break; }
1178 done = readTreeString(filehandle);
1179 if (done == 0) { break; }
1181 //ifyou are a nexus file
1182 }else if((c = filehandle.peek()) == '#') {
1185 // get past comments
1186 while(holder != "translate" && holder != "Translate"){
1187 if(holder == "[" || holder == "[!"){
1193 filehandle >> holder;
1195 //if there is no translate then you must read tree string otherwise use translate to get names
1196 if((holder == "tree") && (comment != 1)){
1197 //pass over the "tree rep.6878900 = "
1198 while (((c = filehandle.get()) != '(') && ((c = filehandle.peek()) != EOF)) {;}
1200 if(c == EOF) { break; }
1201 filehandle.putback(c); //put back first ( of tree.
1202 done = readTreeString(filehandle);
1207 if (done == 0) { break; }
1210 //use nexus translation rather than parsing tree to save time
1211 if((holder == "translate") || (holder == "Translate")) {
1213 string number, name, h;
1214 h = ""; // so it enters the loop the first time
1215 while((h != ";") && (number != ";")) {
1216 filehandle >> number;
1219 //c = , until done with translation then c = ;
1220 h = name.substr(name.length()-1, name.length());
1221 name.erase(name.end()-1); //erase the comma
1222 m->Treenames.push_back(number);
1224 if(number == ";") { m->Treenames.pop_back(); } //in case ';' from translation is on next line instead of next to last name
1229 //for (int i = 0; i < globaldata->Treenames.size(); i++) {
1230 //cout << globaldata->Treenames[i] << endl; }
1231 //cout << globaldata->Treenames.size() << endl;
1233 catch(exception& e) {
1234 m->errorOut(e, "Tree", "parseTreeFile");
1238 /*******************************************************/
1240 /*******************************************************/
1241 int Tree::readTreeString(ifstream& filehandle) {
1246 while((c = filehandle.peek()) != ';') {
1248 //cout << " at beginning of while " << k << endl;
1250 //to pass over labels in trees
1252 while((c!=',') && (c != -1) && (c!= ':') && (c!=';')){ c=filehandle.get(); }
1253 filehandle.putback(c);
1255 if(c == ';') { return 0; }
1256 if(c == -1) { return 0; }
1258 if((c != '(') && (c != ')') && (c != ',') && (c != ':') && (c != '\n') && (c != '\t') && (c != 32)) { //32 is space
1260 c = filehandle.get();
1262 //cout << k << endl;
1263 while ((c != '(') && (c != ')') && (c != ',') && (c != ':') && (c != '\n') && (c != 32) && (c != '\t')) {
1265 c = filehandle.get();
1267 //cout << " in name while " << k << endl;
1270 //cout << "name = " << name << endl;
1271 m->Treenames.push_back(name);
1272 filehandle.putback(c);
1274 //cout << " after putback" << k << endl;
1277 if(c == ':') { //read until you reach the end of the branch length
1278 while ((c != '(') && (c != ')') && (c != ',') && (c != ';') && (c != '\n') && (c != '\t') && (c != 32)) {
1279 c = filehandle.get();
1281 //cout << " in branch while " << k << endl;
1283 filehandle.putback(c);
1286 c = filehandle.get();
1288 //cout << " here after get " << k << endl;
1289 if(c == ';') { return 0; }
1290 if(c == ')') { filehandle.putback(c); }
1292 //cout << k << endl;
1297 catch(exception& e) {
1298 m->errorOut(e, "Tree", "readTreeString");
1303 /*******************************************************/
1305 /*******************************************************/