From: westcott Date: Tue, 7 Dec 2010 13:03:17 +0000 (+0000) Subject: fixed indicator grouping problem X-Git-Url: https://git.donarmstrong.com/?p=mothur.git;a=commitdiff_plain;h=3247d888e7aafc4a65ec9062a94dfd166c2c5b1d fixed indicator grouping problem --- diff --git a/indicatorcommand.cpp b/indicatorcommand.cpp index 9f96f35..01627c0 100644 --- a/indicatorcommand.cpp +++ b/indicatorcommand.cpp @@ -318,7 +318,7 @@ int IndicatorCommand::GetIndicatorSpecies(Tree*& T){ //you need the distances to leaf to decide grouping below //this will also set branch lengths if the tree does not include them - map distToLeaf = getLengthToLeaf(T); + map distToRoot = getDistToRoot(T); //for each node for (int i = T->getNumLeaves(); i < T->getNumNodes(); i++) { @@ -336,9 +336,8 @@ int IndicatorCommand::GetIndicatorSpecies(Tree*& T){ //get nodes that will be a valid grouping //you are valid if you are not one of my descendants - //AND your distToLeaf is <= mine - //AND your distToLeaf is >= my smallest childs - //AND you were not added as part of a larger groupings + //AND your distToRoot is >= mine + //AND you were not added as part of a larger grouping. Largest nodes are added first. set groupsAlreadyAdded; //create a grouping with my grouping @@ -358,7 +357,9 @@ int IndicatorCommand::GetIndicatorSpecies(Tree*& T){ for (int j = (T->getNumNodes()-1); j >= 0; j--) { - if ((descendantNodes[i].count(j) == 0) && (distToLeaf[j] <= distToLeaf[i]) && ((distToLeaf[j] >= distToLeaf[T->tree[i].getLChild()]) || (distToLeaf[j] >= distToLeaf[T->tree[i].getRChild()]))) { + + + if ((descendantNodes[i].count(j) == 0) && (distToRoot[j] >= distToRoot[i])) { vector subset; int count = 0; int doneCount = nodeToDescendants[j].size(); @@ -377,10 +378,10 @@ int IndicatorCommand::GetIndicatorSpecies(Tree*& T){ } } - if (groupsAlreadyAdded.size() != lookup.size()) { cout << i << '\t' << groupsAlreadyAdded.size() << '\t' << lookup.size() << endl; m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); } - for (int k = 0; k < lookup.size(); k++) { - if (groupsAlreadyAdded.count(lookup[k]->getGroup()) == 0) { cout << lookup[k]->getGroup() << endl; } - } + if (groupsAlreadyAdded.size() != lookup.size()) { m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); } + //for (int k = 0; k < lookup.size(); k++) { + // if (groupsAlreadyAdded.count(lookup[k]->getGroup()) == 0) { cout << lookup[k]->getGroup() << endl; } + //} indicatorValues = getValues(groupings); @@ -389,9 +390,8 @@ int IndicatorCommand::GetIndicatorSpecies(Tree*& T){ //get nodes that will be a valid grouping //you are valid if you are not one of my descendants - //AND your distToLeaf is <= mine - //AND your distToLeaf is >= my smallest childs - //AND you were not added as part of a larger grouping + //AND your distToRoot is >= mine + //AND you were not added as part of a larger grouping. Largest nodes are added first. set groupsAlreadyAdded; //create a grouping with my grouping @@ -410,7 +410,7 @@ int IndicatorCommand::GetIndicatorSpecies(Tree*& T){ if (subset.size() != 0) { groupings.push_back(subset); } for (int j = (T->getNumNodes()-1); j >= 0; j--) { - if ((descendantNodes[i].count(j) == 0) && (distToLeaf[j] <= distToLeaf[i]) && ((distToLeaf[j] >= distToLeaf[T->tree[i].getLChild()]) || (distToLeaf[j] >= distToLeaf[T->tree[i].getRChild()]))) { + if ((descendantNodes[i].count(j) == 0) && (distToRoot[j] >= distToRoot[i])) { vector subset; int count = 0; int doneCount = nodeToDescendants[j].size(); @@ -576,9 +576,9 @@ vector IndicatorCommand::getValues(vector< vector >& } } //********************************************************************************************************************** -//you need the distances to leaf to decide groupings +//you need the distances to root to decide groupings //this will also set branch lengths if the tree does not include them -map IndicatorCommand::getLengthToLeaf(Tree*& T){ +map IndicatorCommand::getDistToRoot(Tree*& T){ try { map dists; @@ -587,13 +587,13 @@ map IndicatorCommand::getLengthToLeaf(Tree*& T){ if (T->tree[i].getBranchLength() > 0.0) { hasBranchLengths = true; break; } } - for (int i = 0; i < T->getNumNodes(); i++) { - - int lc = T->tree[i].getLChild(); - int rc = T->tree[i].getRChild(); - - //if you have no branch length, set it then calc - if (!hasBranchLengths) { + //set branchlengths if needed + if (!hasBranchLengths) { + for (int i = 0; i < T->getNumNodes(); i++) { + + int lc = T->tree[i].getLChild(); + int rc = T->tree[i].getRChild(); + if (lc == -1) { // you are a leaf //if you are a leaf set you priliminary length to 1.0, this may adjust later T->tree[i].setBranchLength(1.0); @@ -604,28 +604,32 @@ map IndicatorCommand::getLengthToLeaf(Tree*& T){ float rdist = dists[rc]; float greater = ldist; - if (rdist > greater) { greater = rdist; dists[i] = ldist + 1.0; } + if (rdist > greater) { greater = rdist; dists[i] = ldist + 1.0;} else { dists[i] = rdist + 1.0; } + //branch length = difference + 1 T->tree[lc].setBranchLength((abs(ldist-greater) + 1.0)); T->tree[rc].setBranchLength((abs(rdist-greater) + 1.0)); } - - }else{ - if (lc == -1) { dists[i] = T->tree[i].getBranchLength(); } - else { //smaller of my two children distances plus my branch length - //look at your children's length to leaf - float ldist = dists[lc]; - float rdist = dists[rc]; - - float smaller = ldist; - if (rdist < smaller) { smaller = rdist; } - - dists[i] = smaller + T->tree[i].getBranchLength(); + } + } + + dists.clear(); + + for (int i = 0; i < T->getNumNodes(); i++) { + + double sum = 0.0; + int index = i; + + while(T->tree[index].getParent() != -1){ + if (T->tree[index].getBranchLength() != -1) { + sum += abs(T->tree[index].getBranchLength()); } + index = T->tree[index].getParent(); } + dists[i] = sum; } return dists; diff --git a/indicatorcommand.h b/indicatorcommand.h index 1a41ffd..0ec100a 100644 --- a/indicatorcommand.h +++ b/indicatorcommand.h @@ -47,7 +47,7 @@ private: set getDescendantList(Tree*&, int, map >, map >&); vector getValues(vector< vector >&); vector getValues(vector< vector >&); - map getLengthToLeaf(Tree*&); + map getDistToRoot(Tree*&); }; diff --git a/mothur b/mothur index a9d33c7..835aa5c 100755 Binary files a/mothur and b/mothur differ diff --git a/trimseqscommand.cpp b/trimseqscommand.cpp index fce4763..3c0d378 100644 --- a/trimseqscommand.cpp +++ b/trimseqscommand.cpp @@ -360,12 +360,12 @@ int TrimSeqsCommand::execute(){ #endif if (m->control_pressed) { return 0; } - cout << "done with driver " << endl; + for(int i=0;iisBlank(fastaFileNames[i])) { cout << fastaFileNames[i] << " was blank" << endl; remove(fastaFileNames[i].c_str()); } - else if (filesToRemove.count(fastaFileNames[i]) > 0) { cout << fastaFileNames[i] << " was on the remove list" << endl; remove(fastaFileNames[i].c_str()); } + if (m->isBlank(fastaFileNames[i])) { remove(fastaFileNames[i].c_str()); } + else if (filesToRemove.count(fastaFileNames[i]) > 0) { remove(fastaFileNames[i].c_str()); } else { ifstream inFASTA; string seqName; @@ -382,7 +382,7 @@ int TrimSeqsCommand::execute(){ if(itCombo->second == i){ thisGroup = itCombo->first; combos.erase(itCombo); break; } } }else{ thisGroup = groupVector[i]; } - cout << thisGroup << '\t' << i << '\t' << comboStarts << endl; + while(!inFASTA.eof()){ if(inFASTA.get() == '>'){ inFASTA >> seqName; @@ -394,12 +394,12 @@ int TrimSeqsCommand::execute(){ inFASTA.close(); } } - cout << "done with fastaFileNames " << endl; + if(qFileName != ""){ for(int i=0;iisBlank(qualFileNames[i])) { cout << qualFileNames[i] << " was blank" << endl; remove(qualFileNames[i].c_str()); } - else if (filesToRemove.count(qualFileNames[i]) > 0) { cout << qualFileNames[i] << " was on the remove list" << endl; remove(qualFileNames[i].c_str()); } + if (m->isBlank(qualFileNames[i])) { remove(qualFileNames[i].c_str()); } + else if (filesToRemove.count(qualFileNames[i]) > 0) { remove(qualFileNames[i].c_str()); } else { ifstream inQual; string seqName; @@ -419,7 +419,7 @@ int TrimSeqsCommand::execute(){ } } } - cout << "done with qualFileNames " << endl; + if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } @@ -460,8 +460,6 @@ int TrimSeqsCommand::driverCreateTrim(string filename, string qFileName, string } ofstream outGroups; - //vector fastaFileNames; - //vector qualFileNames; if (oligoFile != "") { m->openOutputFile(groupFile, outGroups); @@ -753,7 +751,7 @@ int TrimSeqsCommand::createProcessesCreateTrim(string filename, string qFileName } } - if (allFiles) { m->mothurOut("Done with allfile"); m->mothurOutEndLine(); } + if (allFiles) { m->mothurOut("Done with allfiles"); m->mothurOutEndLine(); } } return exitCommand;