From 7762fd5ec7582a78a31bb690a5c05b76e5b899c2 Mon Sep 17 00:00:00 2001 From: westcott Date: Wed, 23 Feb 2011 21:56:08 +0000 Subject: [PATCH] fixed bug in unifrac commands with nexus translation if files don't match --- readtree.cpp | 19 +++++++++++++++---- tree.cpp | 24 ++++++++++++++++-------- unifracunweightedcommand.cpp | 2 +- unweighted.cpp | 1 - weighted.cpp | 16 +++++++--------- 5 files changed, 39 insertions(+), 23 deletions(-) diff --git a/readtree.cpp b/readtree.cpp index 8158a57..e8aeacc 100644 --- a/readtree.cpp +++ b/readtree.cpp @@ -129,7 +129,10 @@ int ReadNewickTree::read() { } //if you are a nexus file }else if ((c = filehandle.peek()) == '#') { - holder = nexusTranslation(); //reads file through the translation and updates treemap + //get right number of seqs from nexus file. + Tree* temp = new Tree(); delete temp; + + nexusTranslation(); //reads file through the translation and updates treemap while((c = filehandle.peek()) != EOF) { // get past comments while ((c = filehandle.peek()) != EOF) { @@ -179,7 +182,7 @@ string ReadNewickTree::nexusTranslation() { try { holder = ""; - //int numSeqs = globaldata->gTreemap->getNumSeqs(); //must save this some when we clear old names we can still know how many sequences there were + int numSeqs = globaldata->Treenames.size(); //must save this some when we clear old names we can still know how many sequences there were int comment = 0; // get past comments @@ -197,7 +200,7 @@ string ReadNewickTree::nexusTranslation() { //update treemap globaldata->gTreemap->namesOfSeqs.clear(); - char c; + /*char c; string number, name; while ((c = filehandle.peek()) != EOF) { @@ -213,7 +216,15 @@ string ReadNewickTree::nexusTranslation() { if ((name == "tree") || (name == ";") ) { break; } if (lastChar == ',') { name.erase(name.end()-1); } //erase the comma - + */ + cout << "numseqs = " << numSeqs << endl; + string number, name; + for(int i=0;i> number; + filehandle >> name; + name.erase(name.end()-1); //erase the comma + //insert new one with new name string group = globaldata->gTreemap->getGroup(name); globaldata->gTreemap->treemap[toString(number)].groupname = group; diff --git a/tree.cpp b/tree.cpp index 06f92cf..02f4004 100644 --- a/tree.cpp +++ b/tree.cpp @@ -330,7 +330,9 @@ void Tree::getSubTree(Tree* copy, vector Groups) { copy->tree[i].setParent(grandparent); copy->tree[i].setBranchLength((copy->tree[i].getBranchLength()+copy->tree[parent].getBranchLength())); - copy->tree[grandparent].setChildren(grandparentLC, grandparentRC); + if (grandparent != -1) { + copy->tree[grandparent].setChildren(grandparentLC, grandparentRC); + } removedLeaves.insert(sibIndex); } }else{ @@ -357,7 +359,9 @@ void Tree::getSubTree(Tree* copy, vector Groups) { copy->tree[sibIndex].setParent(grandparent); copy->tree[sibIndex].setBranchLength((copy->tree[sibIndex].getBranchLength()+copy->tree[parent].getBranchLength())); - copy->tree[grandparent].setChildren(grandparentLC, grandparentRC); + if (grandparent != -1) { + copy->tree[grandparent].setChildren(grandparentLC, grandparentRC); + } removedLeaves.insert(i); }else{ //neither of us are, so we want to eliminate ourselves and our parent @@ -367,8 +371,11 @@ void Tree::getSubTree(Tree* copy, vector Groups) { int parentsSibIndex; if (grandparent != -1) { int greatgrandparent = copy->tree[grandparent].getParent(); - int greatgrandparentLC = copy->tree[greatgrandparent].getLChild(); - int greatgrandparentRC = copy->tree[greatgrandparent].getRChild(); + int greatgrandparentLC, greatgrandparentRC; + if (greatgrandparent != -1) { + greatgrandparentLC = copy->tree[greatgrandparent].getLChild(); + greatgrandparentRC = copy->tree[greatgrandparent].getRChild(); + } int grandparentLC = copy->tree[grandparent].getLChild(); int grandparentRC = copy->tree[grandparent].getRChild(); @@ -382,10 +389,12 @@ void Tree::getSubTree(Tree* copy, vector Groups) { copy->tree[parentsSibIndex].setParent(greatgrandparent); copy->tree[parentsSibIndex].setBranchLength((copy->tree[parentsSibIndex].getBranchLength()+copy->tree[grandparent].getBranchLength())); - copy->tree[greatgrandparent].setChildren(greatgrandparentLC, greatgrandparentRC); + if (greatgrandparent != -1) { + copy->tree[greatgrandparent].setChildren(greatgrandparentLC, greatgrandparentRC); + } }else{ - copy->tree[parent].setChildren(-1, -1); - cout << "issues with making subtree" << endl; + copy->tree[parent].setParent(-1); + //cout << "issues with making subtree" << endl; } removedLeaves.insert(sibIndex); removedLeaves.insert(i); @@ -403,7 +412,6 @@ void Tree::getSubTree(Tree* copy, vector Groups) { int nextSpot = numLeaves; populateNewTree(copy->tree, root, nextSpot); - } catch(exception& e) { m->errorOut(e, "Tree", "getCopy"); diff --git a/unifracunweightedcommand.cpp b/unifracunweightedcommand.cpp index fbbc4f1..63487e0 100644 --- a/unifracunweightedcommand.cpp +++ b/unifracunweightedcommand.cpp @@ -432,7 +432,7 @@ void UnifracUnweightedCommand::createPhylipFile(int i) { out << name << '\t'; //output distances - for (int l = 0; l < globaldata->Groups.size(); l++) { out << setprecision(11) << dists[r][l] << '\t'; } + for (int l = 0; l < globaldata->Groups.size(); l++) { out << dists[r][l] << '\t'; } out << endl; }else{ //output distances diff --git a/unweighted.cpp b/unweighted.cpp index d0eea3d..f0c7c40 100644 --- a/unweighted.cpp +++ b/unweighted.cpp @@ -260,7 +260,6 @@ EstOutput Unweighted::getValues(Tree* t, string groupA, string groupB, int p, st processors = p; outputDir = o; - //if the users enters no groups then give them the score of all groups int numGroups = globaldata->Groups.size(); diff --git a/weighted.cpp b/weighted.cpp index a0d593c..b158815 100644 --- a/weighted.cpp +++ b/weighted.cpp @@ -183,7 +183,7 @@ EstOutput Weighted::driver(Tree* t, vector< vector > namesOfGroupCombos, D[count] += weightedSum; } - + //adding the wieghted sums from group l for (int j = 0; j < t->groupNodeInfo[groupB].size(); j++) { //the leaf nodes that have seqs from group l map::iterator it = t->tree[t->groupNodeInfo[groupB][j]].pcount.find(groupB); @@ -194,15 +194,13 @@ EstOutput Weighted::driver(Tree* t, vector< vector > namesOfGroupCombos, D[count] += weightedSum; } - count++; } - + //calculate u for the group comb - for (int h = start; h < (start+num); h++) { //report progress - m->mothurOut("Processing combo: " + toString(h)); m->mothurOutEndLine(); + //m->mothurOut("Processing combo: " + toString(h)); m->mothurOutEndLine(); string groupA = namesOfGroupCombos[h][0]; string groupB = namesOfGroupCombos[h][1]; @@ -224,6 +222,7 @@ EstOutput Weighted::driver(Tree* t, vector< vector > namesOfGroupCombos, //does this node have descendants from group l it = t->tree[i].pcount.find(groupB); + //if it does subtract their percentage from u if (it != t->tree[i].pcount.end()) { u -= (double) t->tree[i].pcount[groupB] / (double) tmap->seqsPerGroup[groupB]; @@ -236,9 +235,8 @@ EstOutput Weighted::driver(Tree* t, vector< vector > namesOfGroupCombos, WScore[(groupA+groupB)] += u; } } - - } + } /********************************************************/ @@ -372,7 +370,7 @@ double Weighted::getLengthToRoot(Tree* t, int v, string groupA, string groupB) { int sib = lc; if (lc == index) { sib = rc; } - + map::iterator itGroup; int pcountSize = 0; itGroup = t->tree[sib].pcount.find(groupA); @@ -402,12 +400,12 @@ double Weighted::getLengthToRoot(Tree* t, int v, string groupA, string groupB) { //get all nodes above the root to add so we don't add their u values above index = *(rootForGrouping[grouping].begin()); + while(t->tree[index].getParent() != -1){ int parent = t->tree[index].getParent(); rootForGrouping[grouping].insert(parent); index = parent; } - return sum; } catch(exception& e) { -- 2.39.2