fixed bug in unifrac commands with nexus translation if files don't match

diff --git a/readtree.cpp b/readtree.cpp
index 8158a577f478b71139ce5c6f91eeb9dedc4b0d62..e8aeacc4431688ca0dfc8f94bf546c1fea2a91dd 100644 (file)
--- 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<numSeqs;i++){
+                       
+                       filehandle >> 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 06f92cf46e22cbb28f4ff28596977d5b62fb7e02..02f400407cfaa5f8c9283fa4f96d8bbc2fb0c3d8 100644 (file)
--- a/tree.cpp
+++ b/tree.cpp
@@ -330,7 +330,9 @@ void Tree::getSubTree(Tree* copy, vector<string> 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<string> 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<string> 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<string> 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<string> 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 fbbc4f1de67cebeb1b9e57a4dce47f499e272f63..63487e05cfc77ecfd6baa705b3a6a756a009cac4 100644 (file)
--- 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 d0eea3dc2b9f1f507213627e6d18df9a4e6afe68..f0c7c4051a7f7d83214cb921ae40570e7f4f8dd7 100644 (file)
--- 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 a0d593c803f691821afdb680e79be80f179b0afa..b158815f90f266e07d64d4efa173fb7c4edaa0a0 100644 (file)
--- a/weighted.cpp
+++ b/weighted.cpp
@@ -183,7 +183,7 @@ EstOutput Weighted::driver(Tree* t, vector< vector<string> > 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<string, int>::iterator it = t->tree[t->groupNodeInfo[groupB][j]].pcount.find(groupB);
@@ -194,15 +194,13 @@ EstOutput Weighted::driver(Tree* t, vector< vector<string> > 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<string> > 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<string> > 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<string, int>::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) {