5 * Created by Sarah Westcott on 1/26/09.
6 * Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
10 #include "parsimonycommand.h"
12 //**********************************************************************************************************************
13 vector<string> ParsimonyCommand::getValidParameters(){
15 string Array[] = {"random","groups","iters","outputdir","inputdir"};
16 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
20 m->errorOut(e, "ParsimonyCommand", "getValidParameters");
24 //**********************************************************************************************************************
25 ParsimonyCommand::ParsimonyCommand(){
28 //initialize outputTypes
29 vector<string> tempOutNames;
30 outputTypes["parsimony"] = tempOutNames;
31 outputTypes["psummary"] = tempOutNames;
34 m->errorOut(e, "ParsimonyCommand", "ParsimonyCommand");
38 //**********************************************************************************************************************
39 vector<string> ParsimonyCommand::getRequiredParameters(){
41 vector<string> myArray;
45 m->errorOut(e, "ParsimonyCommand", "getRequiredParameters");
49 //**********************************************************************************************************************
50 vector<string> ParsimonyCommand::getRequiredFiles(){
52 vector<string> myArray;
56 m->errorOut(e, "ParsimonyCommand", "getRequiredFiles");
60 /***********************************************************/
61 ParsimonyCommand::ParsimonyCommand(string option) {
63 globaldata = GlobalData::getInstance();
67 //allow user to run help
68 if(option == "help") { help(); abort = true; }
71 //valid paramters for this command
72 string Array[] = {"random","groups","iters","outputdir","inputdir"};
73 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
75 OptionParser parser(option);
76 map<string, string> parameters = parser.getParameters();
78 ValidParameters validParameter;
80 //check to make sure all parameters are valid for command
81 for (map<string,string>::iterator it = parameters.begin(); it != parameters.end(); it++) {
82 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
85 //initialize outputTypes
86 vector<string> tempOutNames;
87 outputTypes["parsimony"] = tempOutNames;
88 outputTypes["psummary"] = tempOutNames;
90 randomtree = validParameter.validFile(parameters, "random", false); if (randomtree == "not found") { randomtree = ""; }
92 //are you trying to use parsimony without reading a tree or saying you want random distribution
93 if (randomtree == "") {
94 if (globaldata->gTree.size() == 0) {
95 m->mothurOut("You must read a treefile and a groupfile or set the randomtree parameter to the output filename you wish, before you may execute the parsimony command."); m->mothurOutEndLine(); abort = true; }
98 //if the user changes the output directory command factory will send this info to us in the output parameter
99 string outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){ outputDir = ""; }
101 //check for optional parameter and set defaults
102 // ...at some point should added some additional type checking...
103 groups = validParameter.validFile(parameters, "groups", false);
104 if (groups == "not found") { groups = ""; globaldata->Groups.clear(); }
106 m->splitAtDash(groups, Groups);
107 globaldata->Groups = Groups;
110 itersString = validParameter.validFile(parameters, "iters", false); if (itersString == "not found") { itersString = "1000"; }
111 convert(itersString, iters);
113 if (abort == false) {
114 //randomtree will tell us if user had their own treefile or if they just want the random distribution
115 //user has entered their own tree
116 if (randomtree == "") {
117 T = globaldata->gTree;
118 tmap = globaldata->gTreemap;
120 if(outputDir == "") { outputDir += m->hasPath(globaldata->getTreeFile()); }
121 output = new ColumnFile(outputDir + m->getSimpleName(globaldata->getTreeFile()) + ".parsimony", itersString);
122 outputNames.push_back(outputDir + m->getSimpleName(globaldata->getTreeFile()) + ".parsimony");
123 outputTypes["parsimony"].push_back(outputDir + m->getSimpleName(globaldata->getTreeFile()) + ".parsimony");
125 sumFile = outputDir + m->getSimpleName(globaldata->getTreeFile()) + ".psummary";
126 m->openOutputFile(sumFile, outSum);
127 outputNames.push_back(sumFile);
128 outputTypes["psummary"].push_back(sumFile);
129 }else { //user wants random distribution
130 savetmap = globaldata->gTreemap;
133 if(outputDir == "") { outputDir += m->hasPath(randomtree); }
134 output = new ColumnFile(outputDir+ m->getSimpleName(randomtree), itersString);
135 outputNames.push_back(outputDir+ m->getSimpleName(randomtree));
136 outputTypes["parsimony"].push_back(outputDir+ m->getSimpleName(randomtree));
139 //set users groups to analyze
140 util = new SharedUtil();
141 util->setGroups(globaldata->Groups, tmap->namesOfGroups, allGroups, numGroups, "parsimony"); //sets the groups the user wants to analyze
142 util->getCombos(groupComb, globaldata->Groups, numComp);
144 if (numGroups == 1) { numComp++; groupComb.push_back(allGroups); }
146 pars = new Parsimony(tmap);
154 catch(exception& e) {
155 m->errorOut(e, "ParsimonyCommand", "ParsimonyCommand");
160 //**********************************************************************************************************************
162 void ParsimonyCommand::help(){
164 m->mothurOut("The parsimony command can only be executed after a successful read.tree command, unless you use the random parameter.\n");
165 m->mothurOut("The parsimony command parameters are random, groups and iters. No parameters are required.\n");
166 m->mothurOut("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");
167 m->mothurOut("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");
168 m->mothurOut("The parsimony command should be in the following format: parsimony(random=yourOutputFilename, groups=yourGroups, iters=yourIters).\n");
169 m->mothurOut("Example parsimony(random=out, iters=500).\n");
170 m->mothurOut("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");
171 m->mothurOut("and iters is 1000. The parsimony command output two files: .parsimony and .psummary their descriptions are in the manual.\n");
172 m->mothurOut("Note: No spaces between parameter labels (i.e. random), '=' and parameters (i.e.yourOutputFilename).\n\n");
174 catch(exception& e) {
175 m->errorOut(e, "ParsimonyCommand", "help");
181 /***********************************************************/
182 int ParsimonyCommand::execute() {
185 if (abort == true) { return 0; }
188 reading = new Progress("Comparing to random:", iters);
190 if (m->control_pressed) {
191 delete reading; delete pars; delete util; delete output;
192 if (randomtree == "") { outSum.close(); }
193 for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } outputTypes.clear();
194 globaldata->Groups.clear();
199 //get pscore for users tree
200 userData.resize(numComp,0); //data = AB, AC, BC, ABC.
201 randomData.resize(numComp,0); //data = AB, AC, BC, ABC.
202 rscoreFreq.resize(numComp);
203 uscoreFreq.resize(numComp);
204 rCumul.resize(numComp);
205 uCumul.resize(numComp);
206 userTreeScores.resize(numComp);
207 UScoreSig.resize(numComp);
209 if (randomtree == "") {
210 //get pscores for users trees
211 for (int i = 0; i < T.size(); i++) {
212 userData = pars->getValues(T[i]); //data = AB, AC, BC, ABC.
214 if (m->control_pressed) {
215 delete reading; delete pars; delete util; delete output;
216 if (randomtree == "") { outSum.close(); }
217 for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } outputTypes.clear();
218 globaldata->Groups.clear();
223 //output scores for each combination
224 for(int k = 0; k < numComp; k++) {
227 map<int,double>::iterator it = uscoreFreq[k].find(userData[k]);
228 if (it == uscoreFreq[k].end()) {//new score
229 uscoreFreq[k][userData[k]] = 1;
230 }else{ uscoreFreq[k][userData[k]]++; }
232 //add users score to valid scores
233 validScores[userData[k]] = userData[k];
235 //save score for summary file
236 userTreeScores[k].push_back(userData[k]);
240 //get pscores for random trees
241 for (int j = 0; j < iters; j++) {
243 //create new tree with same num nodes and leaves as users
246 //create random relationships between nodes
247 randT->assembleRandomTree();
249 //get pscore of random tree
250 randomData = pars->getValues(randT);
252 if (m->control_pressed) {
253 delete reading; delete pars; delete util; delete output; delete randT;
254 if (randomtree == "") { outSum.close(); }
255 for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } outputTypes.clear();
256 globaldata->Groups.clear();
260 for(int r = 0; r < numComp; r++) {
261 //add trees pscore to map of scores
262 map<int,double>::iterator it = rscoreFreq[r].find(randomData[r]);
263 if (it != rscoreFreq[r].end()) {//already have that score
264 rscoreFreq[r][randomData[r]]++;
265 }else{//first time we have seen this score
266 rscoreFreq[r][randomData[r]] = 1;
269 //add randoms score to validscores
270 validScores[randomData[r]] = randomData[r];
273 //update progress bar
280 //get pscores for random trees
281 for (int j = 0; j < iters; j++) {
283 //create new tree with same num nodes and leaves as users
285 //create random relationships between nodes
287 randT->assembleRandomTree();
289 if (m->control_pressed) {
290 delete reading; delete pars; delete util; delete output; delete randT;
291 globaldata->gTreemap = savetmap;
292 for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } outputTypes.clear();
293 globaldata->Groups.clear();
298 //get pscore of random tree
299 randomData = pars->getValues(randT);
301 if (m->control_pressed) {
302 delete reading; delete pars; delete util; delete output; delete randT;
303 globaldata->gTreemap = savetmap;
304 for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } outputTypes.clear();
305 globaldata->Groups.clear();
309 for(int r = 0; r < numComp; r++) {
310 //add trees pscore to map of scores
311 map<int,double>::iterator it = rscoreFreq[r].find(randomData[r]);
312 if (it != rscoreFreq[r].end()) {//already have that score
313 rscoreFreq[r][randomData[r]]++;
314 }else{//first time we have seen this score
315 rscoreFreq[r][randomData[r]] = 1;
318 //add randoms score to validscores
319 validScores[randomData[r]] = randomData[r];
322 //update progress bar
329 for(int a = 0; a < numComp; a++) {
330 float rcumul = 0.0000;
331 float ucumul = 0.0000;
332 //this loop fills the cumulative maps and put 0.0000 in the score freq map to make it easier to print.
333 for (map<int,double>::iterator it = validScores.begin(); it != validScores.end(); it++) {
334 if (randomtree == "") {
335 map<int,double>::iterator it2 = uscoreFreq[a].find(it->first);
336 //user data has that score
337 if (it2 != uscoreFreq[a].end()) { uscoreFreq[a][it->first] /= T.size(); ucumul+= it2->second; }
338 else { uscoreFreq[a][it->first] = 0.0000; } //no user trees with that score
340 uCumul[a][it->first] = ucumul;
343 //make rscoreFreq map and rCumul
344 map<int,double>::iterator it2 = rscoreFreq[a].find(it->first);
345 //get percentage of random trees with that info
346 if (it2 != rscoreFreq[a].end()) { rscoreFreq[a][it->first] /= iters; rcumul+= it2->second; }
347 else { rscoreFreq[a][it->first] = 0.0000; } //no random trees with that score
348 rCumul[a][it->first] = rcumul;
351 //find the signifigance of each user trees score when compared to the random trees and save for printing the summary file
352 for (int h = 0; h < userTreeScores[a].size(); h++) {
353 UScoreSig[a].push_back(rCumul[a][userTreeScores[a][h]]);
357 if (m->control_pressed) {
358 delete reading; delete pars; delete util; delete output;
359 if (randomtree == "") { outSum.close(); }
360 else { globaldata->gTreemap = savetmap; }
361 for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } outputTypes.clear();
362 globaldata->Groups.clear();
366 //finish progress bar
371 printParsimonyFile();
372 if (randomtree == "") { printUSummaryFile(); }
374 //reset globaldata's treemap if you just did random distrib
375 if (randomtree != "") {
376 //memory leak prevention
377 //if (globaldata->gTreemap != NULL) { delete globaldata->gTreemap; }
378 globaldata->gTreemap = savetmap;
381 //reset groups parameter
382 globaldata->Groups.clear();
384 if (m->control_pressed) {
385 delete pars; delete util; delete output;
386 for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } outputTypes.clear();
390 m->mothurOutEndLine();
391 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
392 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
393 m->mothurOutEndLine();
399 catch(exception& e) {
400 m->errorOut(e, "ParsimonyCommand", "execute");
405 /***********************************************************/
406 void ParsimonyCommand::printParsimonyFile() {
411 if (randomtree == "") {
412 tags.push_back("Score"); tags.push_back("UserFreq"); tags.push_back("UserCumul"); tags.push_back("RandFreq"); tags.push_back("RandCumul");
414 tags.push_back("Score"); tags.push_back("RandFreq"); tags.push_back("RandCumul");
417 for(int a = 0; a < numComp; a++) {
418 output->initFile(groupComb[a], tags);
420 for (map<int,double>::iterator it = validScores.begin(); it != validScores.end(); it++) {
421 if (randomtree == "") {
422 data.push_back(it->first); data.push_back(uscoreFreq[a][it->first]); data.push_back(uCumul[a][it->first]); data.push_back(rscoreFreq[a][it->first]); data.push_back(rCumul[a][it->first]);
424 data.push_back(it->first); data.push_back(rscoreFreq[a][it->first]); data.push_back(rCumul[a][it->first]);
426 output->output(data);
432 catch(exception& e) {
433 m->errorOut(e, "ParsimonyCommand", "printParsimonyFile");
437 /***********************************************************/
438 int ParsimonyCommand::printUSummaryFile() {
441 outSum << "Tree#" << '\t' << "Groups" << '\t' << "ParsScore" << '\t' << "ParsSig" << endl;
442 m->mothurOut("Tree#\tGroups\tParsScore\tParsSig"); m->mothurOutEndLine();
445 outSum.setf(ios::fixed, ios::floatfield); outSum.setf(ios::showpoint);
449 for (int i = 0; i< T.size(); i++) {
450 for(int a = 0; a < numComp; a++) {
451 if (m->control_pressed) { outSum.close(); return 0; }
452 if (UScoreSig[a][i] > (1/(float)iters)) {
453 outSum << setprecision(6) << i+1 << '\t' << groupComb[a] << '\t' << userTreeScores[a][i] << setprecision(itersString.length()) << '\t' << UScoreSig[a][i] << endl;
454 cout << setprecision(6) << i+1 << '\t' << groupComb[a] << '\t' << userTreeScores[a][i] << setprecision(itersString.length()) << '\t' << UScoreSig[a][i] << endl;
455 m->mothurOutJustToLog(toString(i+1) + "\t" + groupComb[a] + "\t" + toString(userTreeScores[a][i]) + "\t" + toString(UScoreSig[a][i])); m->mothurOutEndLine();
457 outSum << setprecision(6) << i+1 << '\t' << groupComb[a] << '\t' << userTreeScores[a][i] << setprecision(itersString.length()) << '\t' << "<" << (1/float(iters)) << endl;
458 cout << setprecision(6) << i+1 << '\t' << groupComb[a] << '\t' << userTreeScores[a][i] << setprecision(itersString.length()) << '\t' << "<" << (1/float(iters)) << endl;
459 m->mothurOutJustToLog(toString(i+1) + "\t" + groupComb[a] + "\t" + toString(userTreeScores[a][i]) + "\t" + toString((1/float(iters)))); m->mothurOutEndLine();
467 catch(exception& e) {
468 m->errorOut(e, "ParsimonyCommand", "printUSummaryFile");
473 /***********************************************************/
474 void ParsimonyCommand::getUserInput() {
478 tmap = new TreeMap();
480 m->mothurOut("Please enter the number of groups you would like to analyze: ");
482 m->mothurOutJustToLog(toString(numGroups)); m->mothurOutEndLine();
486 numEachGroup.resize(numGroups, 0);
488 for (int i = 1; i <= numGroups; i++) {
489 m->mothurOut("Please enter the number of sequences in group " + toString(i) + ": ");
491 m->mothurOutJustToLog(toString(num)); m->mothurOutEndLine();
493 //set tmaps seqsPerGroup
494 tmap->seqsPerGroup[toString(i)] = num;
495 tmap->namesOfGroups.push_back(toString(i));
497 //set tmaps namesOfSeqs
498 for (int j = 0; j < num; j++) {
499 tmap->namesOfSeqs.push_back(toString(count));
500 tmap->treemap[toString(count)].groupname = toString(i);
505 //clears buffer so next command doesn't have error
509 //save tmap for later
510 //memory leak prevention
511 //if (globaldata->gTreemap != NULL) { delete globaldata->gTreemap; }
512 globaldata->gTreemap = tmap;
513 globaldata->Treenames = tmap->namesOfSeqs;
516 catch(exception& e) {
517 m->errorOut(e, "ParsimonyCommand", "getUserInput");
522 /***********************************************************/