X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=tree.cpp;h=1e115bddad1a451c181c4c26b90b758788df1155;hb=dbc9352256135d7b3e2b17dddc5c506770e2df92;hp=2ad46f044f71556a9726c93a96d4fcd5ae42a6b6;hpb=d97b619c4297b1274c754d73a64792ba656b0a79;p=mothur.git diff --git a/tree.cpp b/tree.cpp index 2ad46f0..1e115bd 100644 --- a/tree.cpp +++ b/tree.cpp @@ -26,8 +26,6 @@ Tree::Tree() { 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); //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; @@ -38,8 +36,6 @@ Tree::Tree() { }else if (i > (numLeaves-1)) { tree[i].setName(""); tree[i].setGroup(""); - //the node knows its index - tree[i].setIndex(i); } } } @@ -54,7 +50,6 @@ Tree::Tree() { } /*****************************************************************/ - int Tree::getIndex(string searchName) { try { //Treemap knows name, group and index to speed up search @@ -155,10 +150,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 +160,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 +178,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 +200,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,49 +294,35 @@ map Tree::mergeGcounts(int position) { } /**************************************************************************************************/ -void Tree::randomLabels() { +void Tree::randomLabels(vector g) { try { - //set up the groups the user wants to include - setGroups(); - - 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; - - while(escape == 1){ - //get random index to switch with - z = int((float)(i+1) * (float)(rand()) / ((float)RAND_MAX+1.0)); - - if(tree[z].pGroups.size() != 0){ - escape = 0; - } - } + //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; - - //leaves have only one group so you can just set it to begin() - it = tree[z].pGroups.begin(); - treez = inUsersGroups(it->first, globaldata->Groups); - - it = tree[i].pGroups.begin(); - treei = inUsersGroups(it->first, globaldata->Groups); + + treez = inUsersGroups(tree[z].getGroup(), g); + treei = inUsersGroups(tree[i].getGroup(), g); 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); - - tree[z].setGroup(tree[z].pGroups.begin()->first); - tree[i].setGroup(tree[i].pGroups.begin()->first); - + + 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); @@ -289,6 +340,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--){ @@ -309,8 +393,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(); } @@ -365,16 +454,18 @@ void Tree::randomTopology() { /*****************************************************************/ // 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); // 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"; @@ -419,7 +510,7 @@ void Tree::printBranch(int node) { printBranch(tree[node].getRChild()); out << ")"; }else { //you are a leaf - tree[node].printNode(out); //prints out name and branch length + out << tree[node].getGroup() << ":" << tree[node].getBranchLength(); } } @@ -434,41 +525,16 @@ void Tree::printBranch(int node) { } /*****************************************************************/ -void Tree::setGroups() { - try { - //if the user has not entered specific groups to analyze then do them all - if (globaldata->Groups.size() != 0) { - //check that groups are valid - for (int i = 0; i < globaldata->Groups.size(); i++) { - if (globaldata->gTreemap->isValidGroup(globaldata->Groups[i]) != true) { - cout << globaldata->Groups[i] << " is not a valid group, and will be disregarded." << endl; - // erase the invalid group from globaldata->Groups - globaldata->Groups.erase (globaldata->Groups.begin()+i); - } - } - - //if the user only entered invalid groups - if (globaldata->Groups.size() == 0) { - cout << "When using the groups parameter you must have at least 1 valid group. I will run the command using all the groups in your groupfile." << endl; - for (int i = 0; i < globaldata->gTreemap->namesOfGroups.size(); i++) { - globaldata->Groups.push_back(globaldata->gTreemap->namesOfGroups[i]); - } - } - - }else { - for (int i = 0; i < globaldata->gTreemap->namesOfGroups.size(); i++) { - globaldata->Groups.push_back(globaldata->gTreemap->namesOfGroups[i]); - } - } - } - catch(exception& e) { - cout << "Standard Error: " << e.what() << " has occurred in the Tree class Function setGroups. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } - catch(...) { - cout << "An unknown error has occurred in the Tree class function setGroups. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; - exit(1); - } +void Tree::printTree() { + + for(int i=0;i