X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=tree.cpp;h=a6c9c07a0714fb850498a3d6f4439d4247602468;hb=c76b65598a0b85183126bd764c83f354e9f490d1;hp=94ab326d8f6ba32ededf4f5b2240dc5852e4816c;hpb=d037597badc8d18e235c59f0c1114180edb7f98f;p=mothur.git diff --git a/tree.cpp b/tree.cpp index 94ab326..a6c9c07 100644 --- a/tree.cpp +++ b/tree.cpp @@ -13,9 +13,9 @@ /*****************************************************************/ Tree::Tree() { try { - globaldata = GlobalData::getInstance(); - numLeaves = globaldata->gTreemap->getNumSeqs(); + + numLeaves = globaldata->Treenames.size(); numNodes = 2*numLeaves - 1; tree.resize(numNodes); @@ -24,22 +24,18 @@ Tree::Tree() { for (int i = 0; i < numNodes; i++) { //initialize leaf nodes if (i <= (numLeaves-1)) { - tree[i].setName(globaldata->gTreemap->namesOfSeqs[i]); - tree[i].setGroup(globaldata->gTreemap->getGroup(globaldata->gTreemap->namesOfSeqs[i])); - //the node knows its index - tree[i].setIndex(i); + tree[i].setName(globaldata->Treenames[i]); + tree[i].setGroup(globaldata->gTreemap->getGroup(globaldata->Treenames[i])); //set pcount and pGroup for groupname to 1. - tree[i].pcount[globaldata->gTreemap->getGroup(globaldata->gTreemap->namesOfSeqs[i])] = 1; - tree[i].pGroups[globaldata->gTreemap->getGroup(globaldata->gTreemap->namesOfSeqs[i])] = 1; + tree[i].pcount[globaldata->gTreemap->getGroup(globaldata->Treenames[i])] = 1; + tree[i].pGroups[globaldata->gTreemap->getGroup(globaldata->Treenames[i])] = 1; //Treemap knows name, group and index to speed up search - globaldata->gTreemap->setIndex(globaldata->gTreemap->namesOfSeqs[i], i); + globaldata->gTreemap->setIndex(globaldata->Treenames[i], i); //intialize non leaf nodes }else if (i > (numLeaves-1)) { tree[i].setName(""); tree[i].setGroup(""); - //the node knows its index - tree[i].setIndex(i); } } } @@ -54,7 +50,8 @@ Tree::Tree() { } /*****************************************************************/ - +Tree::~Tree() {} +/*****************************************************************/ int Tree::getIndex(string searchName) { try { //Treemap knows name, group and index to speed up search @@ -155,10 +152,9 @@ void Tree::getCopy(Tree* copy) { map Tree::mergeGroups(int i) { try { - int lc = tree[i].getLChild(); int rc = tree[i].getRChild(); - + //set parsimony groups to left child map parsimony = tree[lc].pGroups; @@ -166,7 +162,12 @@ map Tree::mergeGroups(int i) { //look at right child groups and update maxPars if right child has something higher for that group. for(it=tree[rc].pGroups.begin();it!=tree[rc].pGroups.end();it++){ - parsimony[it->first]++; + it2 = parsimony.find(it->first); + if (it2 != parsimony.end()) { + parsimony[it->first]++; + }else { + parsimony[it->first] = 1; + } if(parsimony[it->first] > maxPars){ maxPars = parsimony[it->first]; @@ -179,7 +180,7 @@ map Tree::mergeGroups(int i) { for(it=parsimony.begin();it!=parsimony.end();it++){ if(it->second == 1){ parsimony.erase(it->first); - it--; +// it--; } } //set one remaining groups to 1 @@ -201,6 +202,72 @@ map Tree::mergeGroups(int i) { exit(1); } } +/*****************************************************************/ +//returns a map with a groupname and the number of times that group was seen in the children +//for instance if your children are white and black then it would return a map with 2 entries +// p[white] = 1 and p[black] = 1. Now go up a level and merge that with a node who has p[white] = 1 +//and you get p[white] = 2, p[black] = 1, but you erase the p[black] because you have a p value higher than 1. + +map Tree::mergeUserGroups(int i, vector g) { + try { + + int lc = tree[i].getLChild(); + int rc = tree[i].getRChild(); + + //loop through nodes groups removing the ones the user doesn't want + for (it = tree[lc].pGroups.begin(); it != tree[lc].pGroups.end(); it++) { + if (inUsersGroups(it->first, g) != true) { tree[lc].pGroups.erase(it->first); } + } + + //loop through nodes groups removing the ones the user doesn't want + for (it = tree[rc].pGroups.begin(); it != tree[rc].pGroups.end(); it++) { + if (inUsersGroups(it->first, g) != true) { tree[rc].pGroups.erase(it->first); } + } + + //set parsimony groups to left child + map parsimony = tree[lc].pGroups; + + int maxPars = 1; + + //look at right child groups and update maxPars if right child has something higher for that group. + for(it=tree[rc].pGroups.begin();it!=tree[rc].pGroups.end();it++){ + it2 = parsimony.find(it->first); + if (it2 != parsimony.end()) { + parsimony[it->first]++; + }else { + parsimony[it->first] = 1; + } + + if(parsimony[it->first] > maxPars){ + maxPars = parsimony[it->first]; + } + } + + // this is true if right child had a greater parsimony for a certain group + if(maxPars > 1){ + //erase all the groups that are only 1 because you found something with 2. + for(it=parsimony.begin();it!=parsimony.end();it++){ + if(it->second == 1){ + parsimony.erase(it->first); + } + } + for(it=parsimony.begin();it!=parsimony.end();it++){ + parsimony[it->first] = 1; + } + } + + return parsimony; + } + catch(exception& e) { + cout << "Standard Error: " << e.what() << " has occurred in the Tree class Function mergeGroups. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } + catch(...) { + cout << "An unknown error has occurred in the Tree class function mergeGroups. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } +} + /**************************************************************************************************/ @@ -229,35 +296,39 @@ map Tree::mergeGcounts(int position) { } /**************************************************************************************************/ -void Tree::randomLabels() { +void Tree::randomLabels(vector g) { try { - for(int i=numLeaves-1;i>=0;i--){ - if(tree[i].pGroups.size() == 0){ - continue; - } - int escape = 1; + for(int i = 0; i < numLeaves; i++){ int z; + //get random index to switch with + z = int((float)(i+1) * (float)(rand()) / ((float)RAND_MAX+1.0)); + + //you only want to randomize the nodes that are from a group the user wants analyzed, so + //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. + bool treez, treei; - while(escape == 1){ - z = int((float)(i+1) * (float)(rand()) / ((float)RAND_MAX+1.0)); + treez = inUsersGroups(tree[z].getGroup(), g); + treei = inUsersGroups(tree[i].getGroup(), g); - if(tree[z].pGroups.size() != 0){ - escape = 0; - } + if ((treez == true) && (treei == true)) { + //switches node i and node z's info. + map lib_hold = tree[z].pGroups; + tree[z].pGroups = (tree[i].pGroups); + tree[i].pGroups = (lib_hold); + + string zgroup = tree[z].getGroup(); + tree[z].setGroup(tree[i].getGroup()); + tree[i].setGroup(zgroup); + + string zname = tree[z].getName(); + tree[z].setName(tree[i].getName()); + tree[i].setName(zname); + + map gcount_hold = tree[z].pcount; + tree[z].pcount = (tree[i].pcount); + tree[i].pcount = (gcount_hold); } - - - map lib_hold = tree[z].pGroups; - tree[z].pGroups = (tree[i].pGroups); - tree[i].pGroups = (lib_hold); - - tree[z].setGroup(tree[z].pGroups.begin()->first); - tree[i].setGroup(tree[i].pGroups.begin()->first); - - map gcount_hold = tree[z].pcount; - tree[z].pcount = (tree[i].pcount); - tree[i].pcount = (gcount_hold); } } catch(exception& e) { @@ -271,6 +342,39 @@ void Tree::randomLabels() { } /**************************************************************************************************/ +void Tree::randomLabels(string groupA, string groupB) { + try { + int numSeqsA = globaldata->gTreemap->seqsPerGroup[groupA]; + int numSeqsB = globaldata->gTreemap->seqsPerGroup[groupB]; + + vector randomGroups(numSeqsA+numSeqsB, groupA); + for(int i=numSeqsA;i=0;i--){ @@ -291,8 +395,13 @@ void Tree::randomBlengths() { } } /*************************************************************************************************/ -void Tree::assembleRandomUnifracTree() { - randomLabels(); +void Tree::assembleRandomUnifracTree(vector g) { + randomLabels(g); + assembleTree(); +} +/*************************************************************************************************/ +void Tree::assembleRandomUnifracTree(string groupA, string groupB) { + randomLabels(groupA, groupB); assembleTree(); } @@ -344,19 +453,53 @@ void Tree::randomTopology() { exit(1); } } +/*****************************************************************/ +void Tree::print(ostream& out) { + try { + int root = findRoot(); + printBranch(root, out, "branch"); + out << ";" << endl; + } + catch(exception& e) { + cout << "Standard Error: " << e.what() << " has occurred in the Tree class Function print. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } + catch(...) { + cout << "An unknown error has occurred in the Tree class function print. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } +} +/*****************************************************************/ +void Tree::printForBoot(ostream& out) { + try { + int root = findRoot(); + printBranch(root, out, "boot"); + out << ";" << endl; + } + catch(exception& e) { + cout << "Standard Error: " << e.what() << " has occurred in the Tree class Function printForBoot. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } + catch(...) { + cout << "An unknown error has occurred in the Tree class function printForBoot. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } +} /*****************************************************************/ // This prints out the tree in Newick form. -void Tree::createNewickFile() { +void Tree::createNewickFile(string f) { try { int root = findRoot(); - filename = getRootName(globaldata->getTreeFile()) + "newick"; + //filename = getRootName(globaldata->getTreeFile()) + "newick"; + filename = f; openOutputFile(filename, out); - printBranch(root); + printBranch(root, out, "branch"); // you are at the end of the tree out << ";" << endl; + out.close(); } catch(exception& e) { cout << "Standard Error: " << e.what() << " has occurred in the Tree class Function createNewickFile. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; @@ -375,7 +518,7 @@ int Tree::findRoot() { try { for (int i = 0; i < numNodes; i++) { //you found the root - if (tree[i].getParent() == -1) { return i; } + if (tree[i].getParent() == -1) { return i; } } return -1; } @@ -390,18 +533,40 @@ int Tree::findRoot() { } /*****************************************************************/ -void Tree::printBranch(int node) { +void Tree::printBranch(int node, ostream& out, string mode) { try { // you are not a leaf if (tree[node].getLChild() != -1) { out << "("; - printBranch(tree[node].getLChild()); + printBranch(tree[node].getLChild(), out, mode); out << ","; - printBranch(tree[node].getRChild()); + printBranch(tree[node].getRChild(), out, mode); out << ")"; + if (mode == "branch") { + //if there is a branch length then print it + if (tree[node].getBranchLength() != -1) { + out << ":" << tree[node].getBranchLength(); + } + }else if (mode == "boot") { + //if there is a label then print it + if (tree[node].getLabel() != -1) { + out << tree[node].getLabel(); + } + } }else { //you are a leaf - tree[node].printNode(out); //prints out name and branch length + out << tree[node].getGroup(); + if (mode == "branch") { + //if there is a branch length then print it + if (tree[node].getBranchLength() != -1) { + out << ":" << tree[node].getBranchLength(); + } + }else if (mode == "boot") { + //if there is a label then print it + if (tree[node].getLabel() != -1) { + out << tree[node].getLabel(); + } + } } } @@ -417,5 +582,15 @@ void Tree::printBranch(int node) { /*****************************************************************/ +void Tree::printTree() { + + for(int i=0;i