From 692e0c1f69a78b568dc85cbdcea9fb6c189e2e6c Mon Sep 17 00:00:00 2001 From: westcott Date: Wed, 18 Feb 2009 18:28:27 +0000 Subject: [PATCH] parsimony command now uses groups, fixed bug with unweighted groups --- helpcommand.cpp | 10 ++++--- parsimony.cpp | 56 +++++++++++++++++++++++++++++++++++- parsimony.h | 4 +++ parsimonycommand.cpp | 26 +++++++++++++++++ unifracunweightedcommand.cpp | 14 +++++++-- unifracweightedcommand.cpp | 1 + 6 files changed, 104 insertions(+), 7 deletions(-) diff --git a/helpcommand.cpp b/helpcommand.cpp index e5247f9..b627c14 100644 --- a/helpcommand.cpp +++ b/helpcommand.cpp @@ -128,8 +128,10 @@ int HelpCommand::execute(){ cout << "Note: No spaces between parameter labels (i.e. line), '=' and parameters (i.e.yourLines)." << "\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 random parameter." << "\n"; - cout << "The parsimony command parameters are random and iters. No parameters are required." << "\n"; - cout << "The parsimony command should be in the following format: parsimony(random=yourOutputFilename, iters=yourIters)." << "\n"; + cout << "The parsimony command parameters are random, groups and iters. No parameters are required." << "\n"; + cout << "The groups parameter 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 parsimony command should be in the following format: parsimony(random=yourOutputFilename, groups=yourGroups, iters=yourIters)." << "\n"; cout << "Example parsimony(random=out, iters=500)." << "\n"; cout << "The default value for random is "" (meaning you want to use the trees in your inputfile, randomtree=out means you just want the random distribution of trees outputted to out.rd_parsimony)," << "\n"; cout << "and iters is 1000. The parsimony command output three files: .parsimony, .psummary and .pdistrib, their descriptions are in the manual." << "\n"; @@ -137,7 +139,7 @@ int HelpCommand::execute(){ }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 groups parameter 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"; @@ -147,7 +149,7 @@ int HelpCommand::execute(){ }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 groups parameter 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"; diff --git a/parsimony.cpp b/parsimony.cpp index 6c18470..2a673b0 100644 --- a/parsimony.cpp +++ b/parsimony.cpp @@ -12,6 +12,8 @@ /**************************************************************************************************/ EstOutput Parsimony::getValues(Tree* t) { try { + globaldata = GlobalData::getInstance(); + data.resize(1,0); int score = 0; @@ -20,8 +22,42 @@ EstOutput Parsimony::getValues(Tree* t) { int lc = t->tree[i].getLChild(); int rc = t->tree[i].getRChild(); - if(t->tree[i].pGroups.size() > t->tree[rc].pGroups.size() || t->tree[i].pGroups.size() > t->tree[lc].pGroups.size()){ + int iSize = 0; + int rcSize = 0; + int lcSize = 0; + + //add in all the groups the users wanted + for (it = t->tree[i].pGroups.begin(); it != t->tree[i].pGroups.end(); it++) { + if (inUsersGroups(it->first) == true) { iSize++; } + } + + //if that leaves no groups give it 1 so it will cause no change to parent + if (iSize == 0) { iSize++; } + + //add in all the groups the users wanted + for (it = t->tree[rc].pGroups.begin(); it != t->tree[rc].pGroups.end(); it++) { + + if (inUsersGroups(it->first) == true) { rcSize++; } + } + + //if that leaves no groups give it 1 so it will cause no change to parent + if (rcSize == 0) { rcSize++; } + + + //add in all the groups the users wanted + for (it = t->tree[lc].pGroups.begin(); it != t->tree[lc].pGroups.end(); it++) { + + if (inUsersGroups(it->first) == true) { lcSize++; } + } + + //if that leaves no groups give it 1 so it will cause no change to parent + if (lcSize == 0) { lcSize++; } + + + //if you have more groups than either of your kids then theres been a change. + if(iSize > rcSize || iSize > lcSize){ score++; + } } @@ -40,3 +76,21 @@ EstOutput Parsimony::getValues(Tree* t) { } /**************************************************************************************************/ + +bool Parsimony::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 Parsimony class Function inUsersGroups. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } + catch(...) { + cout << "An unknown error has occurred in the Parsimony class function inUsersGroups. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n"; + exit(1); + } +} +/**************************************************************************************************/ diff --git a/parsimony.h b/parsimony.h index b5d7e20..8955ef7 100644 --- a/parsimony.h +++ b/parsimony.h @@ -13,6 +13,7 @@ #include "treecalculator.h" #include "treemap.h" +#include "globaldata.hpp" /***********************************************************************/ @@ -24,8 +25,11 @@ class Parsimony : public TreeCalculator { EstOutput getValues(Tree*); private: + GlobalData* globaldata; EstOutput data; TreeMap* tmap; + bool inUsersGroups(string); + map::iterator it; }; /***********************************************************************/ diff --git a/parsimonycommand.cpp b/parsimonycommand.cpp index 989ea12..534ff0b 100644 --- a/parsimonycommand.cpp +++ b/parsimonycommand.cpp @@ -27,6 +27,30 @@ ParsimonyCommand::ParsimonyCommand() { openOutputFile(sumFile, outSum); distFile = globaldata->getTreeFile() + ".pdistrib"; 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; + for (int i = 0; i < tmap->namesOfGroups.size(); i++) { + globaldata->Groups.push_back(tmap->namesOfGroups[i]); + } + } + }else { + for (int i = 0; i < tmap->namesOfGroups.size(); i++) { + globaldata->Groups.push_back(tmap->namesOfGroups[i]); + } + } }else { //user wants random distribution savetmap = globaldata->gTreemap; @@ -164,6 +188,8 @@ int ParsimonyCommand::execute() { //reset randomTree parameter to "" globaldata->setRandomTree(""); + //reset groups parameter + globaldata->Groups.clear(); return 0; diff --git a/unifracunweightedcommand.cpp b/unifracunweightedcommand.cpp index 3eedb2e..bf2129c 100644 --- a/unifracunweightedcommand.cpp +++ b/unifracunweightedcommand.cpp @@ -37,9 +37,16 @@ UnifracUnweightedCommand::UnifracUnweightedCommand() { //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 < tmap->namesOfGroups.size(); i++) { + globaldata->Groups.push_back(tmap->namesOfGroups[i]); + } + } + }else { + for (int i = 0; i < tmap->namesOfGroups.size(); i++) { + globaldata->Groups.push_back(tmap->namesOfGroups[i]); + } } - + convert(globaldata->getIters(), iters); //how many random trees to generate unweighted = new Unweighted(tmap); @@ -154,6 +161,9 @@ int UnifracUnweightedCommand::execute() { printUnweightedFile(); printUWSummaryFile(); + //reset groups parameter + globaldata->Groups.clear(); + delete randT; return 0; diff --git a/unifracweightedcommand.cpp b/unifracweightedcommand.cpp index 47b23d6..fef08ff 100644 --- a/unifracweightedcommand.cpp +++ b/unifracweightedcommand.cpp @@ -200,6 +200,7 @@ int UnifracWeightedCommand::execute() { //reset randomTree parameter to 0 globaldata->setRandomTree("0"); + //clear out users groups globaldata->Groups.clear(); -- 2.39.2