X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=tree.cpp;h=e138121661f19f77d1c8603036793f86a95dfe67;hb=477e76a8a79b60f6cd4253324dd830bdea25e3e9;hp=2689b037f6a494474ce1e008ba5efb1d250f0791;hpb=5a3592c6478d5d786ec20e4bee71854ad92fdb8c;p=mothur.git diff --git a/tree.cpp b/tree.cpp index 2689b03..e138121 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 @@ -157,7 +154,7 @@ 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; @@ -183,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 @@ -211,7 +208,7 @@ map Tree::mergeGroups(int i) { // 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) { +map Tree::mergeUserGroups(int i, vector g) { try { int lc = tree[i].getLChild(); @@ -219,12 +216,12 @@ map Tree::mergeUserGroups(int i) { //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, globaldata->Groups) != true) { tree[lc].pGroups.erase(it->first); } + 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, globaldata->Groups) != true) { tree[rc].pGroups.erase(it->first); } + if (inUsersGroups(it->first, g) != true) { tree[rc].pGroups.erase(it->first); } } //set parsimony groups to left child @@ -245,24 +242,20 @@ map Tree::mergeUserGroups(int i) { 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); - it--; } } - //set one remaining groups to 1 - //so with our above example p[white] = 2 would be left and it would become p[white] = 1 for(it=parsimony.begin();it!=parsimony.end();it++){ parsimony[it->first] = 1; } + } - } - return parsimony; } catch(exception& e) { @@ -303,12 +296,9 @@ 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 = 0; i < numLeaves; i++){ int z; //get random index to switch with @@ -318,8 +308,8 @@ void Tree::randomLabels() { //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; - treez = inUsersGroups(tree[z].getGroup(), globaldata->Groups); - treei = inUsersGroups(tree[i].getGroup(), 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. @@ -354,33 +344,27 @@ void Tree::randomLabels() { void Tree::randomLabels(string groupA, string groupB) { try { - 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. - if (((tree[z].getGroup() == groupA) || (tree[z].getGroup() == groupB)) && ((tree[i].getGroup() == groupA) || (tree[i].getGroup() == groupB))) { - //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); + int numSeqsA = globaldata->gTreemap->seqsPerGroup[groupA]; + int numSeqsB = globaldata->gTreemap->seqsPerGroup[groupB]; + + vector randomGroups(numSeqsA+numSeqsB, groupA); + for(int i=numSeqsA;i gcount_hold = tree[z].pcount; - tree[z].pcount = (tree[i].pcount); - tree[i].pcount = (gcount_hold); + int randomCounter = 0; + for(int i=0;i g) { + randomLabels(g); assembleTree(); } /*************************************************************************************************/ @@ -469,7 +453,22 @@ void Tree::randomTopology() { exit(1); } } - +/*****************************************************************/ +void Tree::print(ostream& out) { + try { + int root = findRoot(); + printBranch(root, out); + 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); + } +} /*****************************************************************/ // This prints out the tree in Newick form. void Tree::createNewickFile(string f) { @@ -479,7 +478,7 @@ void Tree::createNewickFile(string f) { filename = f; openOutputFile(filename, out); - printBranch(root); + printBranch(root, out); // you are at the end of the tree out << ";" << endl; @@ -502,7 +501,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; } @@ -517,18 +516,26 @@ int Tree::findRoot() { } /*****************************************************************/ -void Tree::printBranch(int node) { +void Tree::printBranch(int node, ostream& out) { try { // you are not a leaf if (tree[node].getLChild() != -1) { out << "("; - printBranch(tree[node].getLChild()); + printBranch(tree[node].getLChild(), out); out << ","; - printBranch(tree[node].getRChild()); + printBranch(tree[node].getRChild(), out); out << ")"; + //if there is a branch length then print it + if (tree[node].getBranchLength() != -1) { + out << ":" << tree[node].getBranchLength(); + } }else { //you are a leaf - out << tree[node].getName() << ":" << tree[node].getBranchLength(); + out << tree[node].getGroup(); + //if there is a branch length then print it + if (tree[node].getBranchLength() != -1) { + out << ":" << tree[node].getBranchLength(); + } } } @@ -544,46 +551,6 @@ 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