From b2dca66a02f8f82aa5528e531eace60fbbd2967b Mon Sep 17 00:00:00 2001 From: westcott Date: Wed, 11 Feb 2009 18:14:23 +0000 Subject: [PATCH] added errorchecking and help info on new unifrac and treeclimber code --- errorchecking.cpp | 9 +++- globaldata.cpp | 17 ++++++-- helpcommand.cpp | 38 ++++++++++++++++- treemap.cpp | 19 ++++++++- treemap.h | 3 +- unifracunweightedcommand.cpp | 27 +++++++++--- unifracweightedcommand.cpp | 34 ++++++++++++++-- unweighted.cpp | 23 +++++++++-- unweighted.h | 2 + validcommands.cpp | 2 +- validparameter.cpp | 2 +- weighted.cpp | 79 +++++++++++++++++++++++++++--------- weighted.h | 1 + 13 files changed, 213 insertions(+), 43 deletions(-) diff --git a/errorchecking.cpp b/errorchecking.cpp index 03844d5..067d246 100644 --- a/errorchecking.cpp +++ b/errorchecking.cpp @@ -90,7 +90,7 @@ bool ErrorCheck::checkInput(string input) { if (parameter == "name" ) { namefile = value; } if (parameter == "order" ) { orderfile = value; } if (parameter == "fasta" ) { fastafile = value; } - if (parameter == "treefile" ) { treefile = value; } + if (parameter == "tree" ) { treefile = value; } if (parameter == "group" ) { groupfile = value; } if (parameter == "cutoff" ) { cutoff = value; } if (parameter == "precision" ) { precision = value; } @@ -161,7 +161,7 @@ bool ErrorCheck::checkInput(string input) { if (parameter == "order" ) { orderfile = value; } if (parameter == "group" ) { groupfile = value; } if (parameter == "fasta" ) { fastafile = value; } - if (parameter == "treefile" ) { treefile = value; } + if (parameter == "tree" ) { treefile = value; } if (parameter == "cutoff" ) { cutoff = value; } if (parameter == "precision" ) { precision = value; } if (parameter == "iters" ) { iters = value; } @@ -254,6 +254,11 @@ bool ErrorCheck::checkInput(string input) { } } + if ((commandName == "unifrac.weighted") || (commandName == "unifrac.unweighted")) { + if (globaldata->gTree.size() == 0) {//no trees were read + cout << "You must execute the read.tree command, before you may execute the unifrac.weighted or unifrac.unweighted command." << endl; return false; } + } + //check for valid method if (commandName == "cluster") { if ((method == "furthest") || (method == "nearest") || (method == "average")) { } diff --git a/globaldata.cpp b/globaldata.cpp index bc836e8..d4eb998 100644 --- a/globaldata.cpp +++ b/globaldata.cpp @@ -103,7 +103,7 @@ void GlobalData::parseGlobalData(string commandString, string optionText){ if (key == "rabund" ) { rabundfile = value; inputFileName = value; fileroot = value; format = "rabund"; } if (key == "sabund" ) { sabundfile = value; inputFileName = value; fileroot = value; format = "sabund"; } if (key == "fasta" ) { fastafile = value; inputFileName = value; fileroot = value; format = "fasta"; } - if (key == "treefile" ) { treefile = value; inputFileName = value; fileroot = value; format = "tree"; } + if (key == "tree" ) { treefile = value; inputFileName = value; fileroot = value; format = "tree"; } if (key == "name" ) { namefile = value; } if (key == "order" ) { orderfile = value; } if (key == "group" ) { groupfile = value; } @@ -154,14 +154,14 @@ void GlobalData::parseGlobalData(string commandString, string optionText){ splitAtDash(value, lines); allLines = 0; } - if (key == "label") {//stores lines to be used in a set + if (key == "label") {//stores labels to be used in a set labels.clear(); label = value; line = ""; splitAtDash(value, labels); allLines = 0; } - if (key == "groups") {//stores lines to be used in a vector + if (key == "groups") {//stores groups to be used in a vector Groups.clear(); groups = value; splitAtDash(value, Groups); @@ -178,7 +178,7 @@ void GlobalData::parseGlobalData(string commandString, string optionText){ if (key == "rabund" ) { rabundfile = value; inputFileName = value; fileroot = value; format = "rabund"; } if (key == "sabund" ) { sabundfile = value; inputFileName = value; fileroot = value; format = "sabund"; } if (key == "fasta" ) { fastafile = value; inputFileName = value; fileroot = value; format = "fasta"; } - if (key == "treefile" ) { treefile = value; inputFileName = value; fileroot = value; format = "tree"; } + if (key == "tree" ) { treefile = value; inputFileName = value; fileroot = value; format = "tree"; } if (key == "name" ) { namefile = value; } if (key == "order" ) { orderfile = value; } if (key == "group" ) { groupfile = value; } @@ -190,6 +190,8 @@ void GlobalData::parseGlobalData(string commandString, string optionText){ if (key == "method" ) { method = value; } if (key == "fileroot" ) { fileroot = value; } if (key == "randomtree" ) { randomtree = value; } + if (key == "groups" ) { groups = value; } + if (key == "single") {//stores estimators in a vector singleEstimators.clear(); //clears out old values @@ -236,6 +238,12 @@ void GlobalData::parseGlobalData(string commandString, string optionText){ splitAtDash(value, labels); allLines = 0; } + if (key == "groups") {//stores groups to be used in a vector + Groups.clear(); + groups = value; + splitAtDash(value, Groups); + } + } //set format for shared @@ -345,6 +353,7 @@ void GlobalData::clear() { iters = "1000"; line = ""; label = ""; + groups = ""; jumble = "1"; //0 means don't jumble, 1 means jumble. randomtree = "0"; //0 means user will enter some user trees, 1 means they just want the random tree distribution. freq = "100"; diff --git a/helpcommand.cpp b/helpcommand.cpp index 6fe7931..bc17dbe 100644 --- a/helpcommand.cpp +++ b/helpcommand.cpp @@ -42,6 +42,12 @@ int HelpCommand::execute(){ cout << "The list parameter and group paramaters are required. When using the command the second way read.otu command parses the .list file" << "\n"; cout << "and separates it into groups. It outputs a .shared file containing the OTU information for each group. The read.otu command also outputs a .list file for each group. " << "\n"; cout << "Note: No spaces between parameter labels (i.e. list), '=' and parameters (i.e.yourListfile)." << "\n" << "\n"; + }else if (globaldata->helpRequest == "read.tree") { + cout << "The read.tree command must be run before you execute a unifrac.weighted, unifrac.unweighted. " << "\n"; + cout << "It also must be run before using the parsimony command, unless you are using the randomtree parameter." << "\n"; + cout << "The read.tree command should be in the following format: read.tree(tree=yourTreeFile, group=yourGroupFile)." << "\n"; + cout << "The tree and group parameters are both required." << "\n"; + cout << "Note: No spaces between parameter labels (i.e. list), '=' and parameters (i.e.yourListfile)." << "\n" << "\n"; }else if (globaldata->helpRequest == "cluster") { cout << "The cluster command can only be executed after a successful read.dist command." << "\n"; cout << "The cluster command parameter options are method, cuttoff and precision. No parameters are required." << "\n"; @@ -112,17 +118,45 @@ int HelpCommand::execute(){ cout << "The default value for jumble is 1 (meaning jumble, if it’s set to 0 then it will not jumble) and sharedsummary is sharedsobs-sharedChao-sharedAce-sharedJabund-sharedSorensonAbund-sharedJclass-sharedSorClass-sharedJest-sharedSorEst-SharedThetaYC-SharedThetaN" << "\n"; cout << "The label and line parameters are used to analyze specific lines in your input." << "\n"; cout << "Note: No spaces between parameter labels (i.e. list), '=' and parameters (i.e.yourListfile)." << "\n" << "\n"; + }else if (globaldata->helpRequest == "parsimony") { + cout << "The parsimony command can only be executed after a successful read.tree command, unless you use the randomtree parameter." << "\n"; + cout << "The parsimony command parameters are randomtree and iters. No parameters are required." << "\n"; + cout << "The parsimony command should be in the following format: parsimony(randomtree=yourRandomTreeValue, iters=yourIters)." << "\n"; + cout << "Example parsimony(randomtree=1, iters=500)." << "\n"; + cout << "The default value for randomTree is 0 (meaning you want to use the trees in your inputfile, randomtree=1 means you just want the random distribution of trees)," << "\n"; + cout << "and iters is 1000. The parsimony command output three files: .parsimony, .psummary and .pdistrib, their descriptions are in the manual." << "\n"; + cout << "Note: No spaces between parameter labels (i.e. list), '=' and parameters (i.e.yourListfile)." << "\n" << "\n"; + }else if (globaldata->helpRequest == "unifrac.weighted") { + cout << "The unifrac.weighted command can only be executed after a successful read.tree command." << "\n"; + cout << "The unifrac.weighted command parameters are groups and iters. No parameters are required." << "\n"; + cout << "The groups paramter allows you to specify which of the groups in your groupfile you would like analyzed. You must enter at least 2 valid groups." << "\n"; + cout << "The group names are separated by dashes. The iters parameter allows you to specify how many random trees you would like compared to your tree." << "\n"; + cout << "The unifrac.weighted command should be in the following format: unifrac.weighted(groups=yourGroups, iters=yourIters)." << "\n"; + cout << "Example unifrac.weighted(groups=A-B-C, iters=500)." << "\n"; + cout << "The default value for groups is all the groups in your groupfile, and iters is 1000." << "\n"; + cout << "The unifrac.weighted command output three files: .weighted, .wsummary and .wdistrib, their descriptions are in the manual." << "\n"; + cout << "Note: No spaces between parameter labels (i.e. list), '=' and parameters (i.e.yourListfile)." << "\n" << "\n"; + }else if (globaldata->helpRequest == "unifrac.unweighted") { + cout << "The unifrac.unweighted command can only be executed after a successful read.tree command." << "\n"; + cout << "The unifrac.unweighted command parameters are groups and iters. No parameters are required." << "\n"; + cout << "The groups paramter allows you to specify which of the groups in your groupfile you would like analyzed. You must enter at least 1 valid group." << "\n"; + cout << "The group names are separated by dashes. The iters parameter allows you to specify how many random trees you would like compared to your tree." << "\n"; + cout << "The unifrac.unweighted command should be in the following format: unifrac.unweighted(groups=yourGroups, iters=yourIters)." << "\n"; + cout << "Example unifrac.unweighted(groups=A-B-C, iters=500)." << "\n"; + cout << "The default value for groups is all the groups in your groupfile, and iters is 1000." << "\n"; + cout << "The unifrac.unweighted command output three files: .unweighted, .uwsummary and .uwdistrib, their descriptions are in the manual." << "\n"; + cout << "Note: No spaces between parameter labels (i.e. list), '=' and parameters (i.e.yourListfile)." << "\n" << "\n"; }else if (globaldata->helpRequest == "quit") { cout << "The quit command will terminate Dotur and should be in the following format: " << "\n"; cout << "quit()" << "\n" << "\n"; }else if (globaldata->helpRequest == "") { - cout << "Valid commands are read.dist(), read.otu(), cluster(), deconvolute(), collect.single(), rarefaction.single(), summary.single(), collect.shared(), rarefaction.shared(), summary.shared(), quit(), help()." << "\n"; + cout << "Valid commands are read.dist(), read.otu(), read.tree(), cluster(), deconvolute(), collect.single(), rarefaction.single(), summary.single(), collect.shared(), rarefaction.shared(), summary.shared(), parsimony(), unifrac.weighted(), unifrac.unweighted(), quit(), help()." << "\n"; cout << "For more information about a specific command type 'help(commandName)' i.e. 'help(read.dist)'" << endl; }else { cout << globaldata->helpRequest << " is not a valid command" << endl; } - cout << endl << "For further assistance please refer to the Mothur manual, or contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + cout << endl << "For further assistance please refer to the Mothur manual on our wiki at http://schloss.micro.umass.edu/mothur/, or contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; return 0; } diff --git a/treemap.cpp b/treemap.cpp index bee0d11..6314eb0 100644 --- a/treemap.cpp +++ b/treemap.cpp @@ -100,7 +100,24 @@ void TreeMap::setNamesOfGroups(string seqGroup) { namesOfGroups.push_back(seqGroup); //new group } } - +/************************************************************/ +bool TreeMap::isValidGroup(string groupname) { + try { + for (int i = 0; i < namesOfGroups.size(); i++) { + if (groupname == namesOfGroups[i]) { return true; } + } + + return false; + } + catch(exception& e) { + cout << "Standard Error: " << e.what() << " has occurred in the TreeMap class Function isValidGroup. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } + catch(...) { + cout << "An unknown error has occurred in the TreeMap class function isValidGroup. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } +} /***********************************************************************/ void TreeMap::print(ostream& output){ diff --git a/treemap.h b/treemap.h index 60f75b3..60fb924 100644 --- a/treemap.h +++ b/treemap.h @@ -34,11 +34,12 @@ public: int getNumSeqs(); void setIndex(string, int); //sequencename, index int getIndex(string); //returns vector index of sequence + bool isValidGroup(string); //return true if string is a valid group string getGroup(string); vector namesOfGroups; vector namesOfSeqs; map seqsPerGroup; //groupname, number of seqs in that group. - map treemap; //sequence name and groupname + map treemap; //sequence name and void print(ostream&); private: diff --git a/unifracunweightedcommand.cpp b/unifracunweightedcommand.cpp index 8adff2d..7706944 100644 --- a/unifracunweightedcommand.cpp +++ b/unifracunweightedcommand.cpp @@ -23,6 +23,23 @@ UnifracUnweightedCommand::UnifracUnweightedCommand() { distFile = globaldata->getTreeFile() + ".uwdistrib"; openOutputFile(distFile, outDist); + //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 (tmap->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; + } + } + convert(globaldata->getIters(), iters); //how many random trees to generate unweighted = new Unweighted(tmap); @@ -93,6 +110,8 @@ int UnifracUnweightedCommand::execute() { outDist << i+1 << '\t' << '\t'<< j+1 << '\t' << '\t' << randomData[0] << endl; } + saveRandomScores(); //save all random scores for unweighted file + //find the signifigance of the score float rcumul = 0.0000; for (it = rscoreFreq.begin(); it != rscoreFreq.end(); it++) { @@ -105,8 +124,7 @@ int UnifracUnweightedCommand::execute() { //save the signifigance of the users score for printing later UWScoreSig.push_back(rCumul[userData[0]]); - saveRandomScores(); //save all random scores for unweighted file - + //clear random data rscoreFreq.clear(); //you clear this because in the summary file you want the unweighted signifinance to be relative to these 1000 trees. rCumul.clear(); @@ -207,15 +225,14 @@ void UnifracUnweightedCommand::printUWSummaryFile() { /***********************************************************/ void UnifracUnweightedCommand::saveRandomScores() { try { - //update total map with new random scores for (it = rscoreFreq.begin(); it != rscoreFreq.end(); it++) { //does this score already exist in the total map it2 = totalrscoreFreq.find(it->first); //if yes then add them if (it2 != totalrscoreFreq.end()) { - it2->second += it->second; + totalrscoreFreq[it->first] += rscoreFreq[it->first]; }else{ //its a new score - totalrscoreFreq[it->first] = 1; + totalrscoreFreq[it->first] = rscoreFreq[it->first]; } } } diff --git a/unifracweightedcommand.cpp b/unifracweightedcommand.cpp index ff2c245..4a8325a 100644 --- a/unifracweightedcommand.cpp +++ b/unifracweightedcommand.cpp @@ -22,8 +22,30 @@ UnifracWeightedCommand::UnifracWeightedCommand() { openOutputFile(sumFile, outSum); distFile = globaldata->getTreeFile() + ".wdistrib"; openOutputFile(distFile, outDist); - - numGroups = tmap->getNumGroups(); + + //if the user has not entered specific groups to analyze then do them all + if (globaldata->Groups.size() == 0) { + numGroups = tmap->getNumGroups(); + }else { + //check that groups are valid + for (int i = 0; i < globaldata->Groups.size(); i++) { + if (tmap->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) { + numGroups = tmap->getNumGroups(); + cout << "When using the groups parameter you must have at least 2 valid groups. I will run the command using all the groups in your groupfile." << endl; + }else if (globaldata->Groups.size() == 1) { + cout << "When using the groups parameter you must have at least 2 valid groups. I will run the command using all the groups in your groupfile." << endl; + numGroups = tmap->getNumGroups(); + globaldata->Groups.clear(); + }else { numGroups = globaldata->Groups.size(); } + } //calculate number of comparisons i.e. with groups A,B,C = AB, AC, BC = 3; numComp = 0; @@ -32,7 +54,11 @@ UnifracWeightedCommand::UnifracWeightedCommand() { numComp += i; for (int l = n; l < numGroups; l++) { //set group comparison labels - groupComb.push_back(tmap->namesOfGroups[i-1]+tmap->namesOfGroups[l]); + if (globaldata->Groups.size() != 0) { + groupComb.push_back(globaldata->Groups[i-1]+globaldata->Groups[l]); + }else { + groupComb.push_back(tmap->namesOfGroups[i-1]+tmap->namesOfGroups[l]); + } } n++; } @@ -174,6 +200,8 @@ int UnifracWeightedCommand::execute() { //reset randomTree parameter to 0 globaldata->setRandomTree("0"); + //clear out users groups + globaldata->Groups.clear(); delete randT; diff --git a/unweighted.cpp b/unweighted.cpp index a7cf9c9..4416e79 100644 --- a/unweighted.cpp +++ b/unweighted.cpp @@ -13,7 +13,8 @@ EstOutput Unweighted::getValues(Tree* t) { try { - + globaldata = GlobalData::getInstance(); + //clear out old values data.resize(1,0); penalty.resize(t->getNumLeaves(), 0); @@ -60,7 +61,7 @@ EstOutput Unweighted::getValues(Tree* t) { map::iterator pos; for(pos=unique.begin();pos!=unique.end();pos++){ - if(pos->first!="xxx"){ + if((pos->first!="xxx") && (inUsersGroups(pos->first))){ UW += unique[pos->first]; } } @@ -85,4 +86,20 @@ EstOutput Unweighted::getValues(Tree* t) { } -/**************************************************************************************************/ \ No newline at end of file +/**************************************************************************************************/ +bool Unweighted::inUsersGroups(string groupname) { + try { + for (int i = 0; i < globaldata->Groups.size(); i++) { + if (groupname == globaldata->Groups[i]) { return true; } + } + return false; + } + catch(exception& e) { + cout << "Standard Error: " << e.what() << " has occurred in the Unweighted class Function inUsersGroups. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } + catch(...) { + cout << "An unknown error has occurred in the Unweighted class function inUsersGroups. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } +} \ No newline at end of file diff --git a/unweighted.h b/unweighted.h index d39b5e7..e3bf3a8 100644 --- a/unweighted.h +++ b/unweighted.h @@ -24,9 +24,11 @@ class Unweighted : public TreeCalculator { EstOutput getValues(Tree*); private: + GlobalData* globaldata; EstOutput data; vector penalty; TreeMap* tmap; + bool inUsersGroups(string); }; diff --git a/validcommands.cpp b/validcommands.cpp index 8f06fef..53ee743 100644 --- a/validcommands.cpp +++ b/validcommands.cpp @@ -55,7 +55,7 @@ bool ValidCommands::isValidCommand(string command) { if ((commands.find(command)) != (commands.end())) { return true; }else{ - cout << command << " is not a valid command in Mothur. Valid commands are read.dist(), read.otu(), cluster(), collect.single(), collect.shared(), rarefaction.single(), rarefaction.shared(), summary.single(), summary.shared(), quit(), help()." << endl; + cout << command << " is not a valid command in Mothur. Valid commands are read.dist(), read.otu(), read.tree(), cluster(), deconvolute(), collect.single(), collect.shared(), rarefaction.single(), rarefaction.shared(), summary.single(), summary.shared(), parsimony(), unifrac.weighted(), unifrac.unweighted(), quit(), help()." << endl; return false; } diff --git a/validparameter.cpp b/validparameter.cpp index 7072825..fbb368e 100644 --- a/validparameter.cpp +++ b/validparameter.cpp @@ -23,7 +23,7 @@ ValidParameters::ValidParameters() { parameters["group"] = "group"; parameters["order"] = "order"; parameters["fasta"] = "fasta"; - parameters["treefile"] = "treefile"; + parameters["tree"] = "tree"; parameters["fileroot"] = "fileroot"; parameters["cutoff"] = "cutoff"; parameters["method"] = "method"; diff --git a/weighted.cpp b/weighted.cpp index 481b34e..bcdee11 100644 --- a/weighted.cpp +++ b/weighted.cpp @@ -13,15 +13,26 @@ EstOutput Weighted::getValues(Tree* t) { try { - - int numGroups = tmap->getNumGroups(); + globaldata = GlobalData::getInstance(); + int numGroups; + + //if the user has not entered specific groups to analyze then do them all + if (globaldata->Groups.size() == 0) { + numGroups = tmap->getNumGroups(); + }else { + numGroups = globaldata->Groups.size(); + } //calculate number of comparisons i.e. with groups A,B,C = AB, AC, BC = 3; int n = 1; for (int i=1; inamesOfGroups[i-1]+tmap->namesOfGroups[l]] = 0.0; + if (globaldata->Groups.size() == 0) { + WScore[tmap->namesOfGroups[i-1]+tmap->namesOfGroups[l]] = 0.0; + }else { + WScore[globaldata->Groups[i-1]+globaldata->Groups[l]] = 0.0; + } } } @@ -70,25 +81,48 @@ EstOutput Weighted::getValues(Tree* t) { for (int b=1; btree[i].pcount.find(tmap->namesOfGroups[b-1]); - //if it does u = # of its descendants with a certain group / total number in tree with a certain group - if (it != t->tree[i].pcount.end()) { - u = (double) t->tree[i].pcount[tmap->namesOfGroups[b-1]] / (double) tmap->seqsPerGroup[tmap->namesOfGroups[b-1]]; - }else { u = 0.00; } + //the user has not entered specific groups + if (globaldata->Groups.size() == 0) { + //does this node have descendants from group b-1 + it = t->tree[i].pcount.find(tmap->namesOfGroups[b-1]); + //if it does u = # of its descendants with a certain group / total number in tree with a certain group + if (it != t->tree[i].pcount.end()) { + u = (double) t->tree[i].pcount[tmap->namesOfGroups[b-1]] / (double) tmap->seqsPerGroup[tmap->namesOfGroups[b-1]]; + }else { u = 0.00; } - //does this node have descendants from group l - it = t->tree[i].pcount.find(tmap->namesOfGroups[l]); - //if it does subtract their percentage from u - if (it != t->tree[i].pcount.end()) { - u -= (double) t->tree[i].pcount[tmap->namesOfGroups[l]] / (double) tmap->seqsPerGroup[tmap->namesOfGroups[l]]; - } + //does this node have descendants from group l + it = t->tree[i].pcount.find(tmap->namesOfGroups[l]); + //if it does subtract their percentage from u + if (it != t->tree[i].pcount.end()) { + u -= (double) t->tree[i].pcount[tmap->namesOfGroups[l]] / (double) tmap->seqsPerGroup[tmap->namesOfGroups[l]]; + } - u = abs(u) * t->tree[i].getBranchLength(); + u = abs(u) * t->tree[i].getBranchLength(); - //save groupcombs u value - WScore[tmap->namesOfGroups[b-1]+tmap->namesOfGroups[l]] += u; - + //save groupcombs u value + WScore[tmap->namesOfGroups[b-1]+tmap->namesOfGroups[l]] += u; + + //the user has entered specific groups + }else { + //does this node have descendants from group b-1 + it = t->tree[i].pcount.find(globaldata->Groups[b-1]); + //if it does u = # of its descendants with a certain group / total number in tree with a certain group + if (it != t->tree[i].pcount.end()) { + u = (double) t->tree[i].pcount[globaldata->Groups[b-1]] / (double) tmap->seqsPerGroup[globaldata->Groups[b-1]]; + }else { u = 0.00; } + + //does this node have descendants from group l + it = t->tree[i].pcount.find(globaldata->Groups[l]); + //if it does subtract their percentage from u + if (it != t->tree[i].pcount.end()) { + u -= (double) t->tree[i].pcount[globaldata->Groups[l]] / (double) tmap->seqsPerGroup[globaldata->Groups[l]]; + } + + u = abs(u) * t->tree[i].getBranchLength(); + + //save groupcombs u value + WScore[globaldata->Groups[b-1]+globaldata->Groups[l]] += u; + } } n++; } @@ -99,7 +133,12 @@ EstOutput Weighted::getValues(Tree* t) { n = 1; for (int i=1; inamesOfGroups[i-1]+tmap->namesOfGroups[l]] / D); + //the user has not entered specific groups + if (globaldata->Groups.size() == 0) { + UN = (WScore[tmap->namesOfGroups[i-1]+tmap->namesOfGroups[l]] / D); + }else {//they have entered specific groups + UN = (WScore[globaldata->Groups[i-1]+globaldata->Groups[l]] / D); + } if (isnan(UN) || isinf(UN)) { UN = 0; } data.push_back(UN); } diff --git a/weighted.h b/weighted.h index b0d4433..7e4be7c 100644 --- a/weighted.h +++ b/weighted.h @@ -24,6 +24,7 @@ class Weighted : public TreeCalculator { EstOutput getValues(Tree*); private: + GlobalData* globaldata; EstOutput data; TreeMap* tmap; map::iterator it; -- 2.39.2