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 ParsimonyCommand::ParsimonyCommand() {
15 globaldata = GlobalData::getInstance();
17 //randomtree will tell us if user had their own treefile or if they just want the random distribution
18 randomtree = globaldata->getRandomTree();
20 //user has entered their own tree
21 if (randomtree == "") {
22 T = globaldata->gTree;
23 tmap = globaldata->gTreemap;
24 parsFile = globaldata->getTreeFile() + ".parsimony";
25 openOutputFile(parsFile, out);
26 sumFile = globaldata->getTreeFile() + ".psummary";
27 openOutputFile(sumFile, outSum);
28 distFile = globaldata->getTreeFile() + ".pdistrib";
29 openOutputFile(distFile, outDist);
30 //set users groups to analyze
32 }else { //user wants random distribution
33 savetmap = globaldata->gTreemap;
35 parsFile = randomtree + ".rd_parsimony";
36 openOutputFile(parsFile, out);
39 convert(globaldata->getIters(), iters); //how many random trees to generate
40 pars = new Parsimony(tmap);
44 cout << "Standard Error: " << e.what() << " has occurred in the ParsimonyCommand class Function ParsimonyCommand. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
48 cout << "An unknown error has occurred in the ParsimonyCommand class function ParsimonyCommand. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
52 /***********************************************************/
53 int ParsimonyCommand::execute() {
56 //get pscore for users tree
57 userData.resize(1,0); //data[0] = pscore.
58 randomData.resize(1,0); //data[0] = pscore.
61 outDist.setf(ios::fixed, ios::floatfield); outDist.setf(ios::showpoint);
62 outDist << "RandomTree#" << '\t' << "ParsScore" << endl;
64 if (randomtree == "") {
65 //get pscores for users trees
66 for (int i = 0; i < T.size(); i++) {
67 cout << "Processing tree " << i+1 << endl;
68 userData = pars->getValues(T[i]); //userData[0] = pscore
69 cout << "Tree " << i+1 << " parsimony score = " << userData[0] << endl;
71 it = uscoreFreq.find(userData[0]);
72 if (it == uscoreFreq.end()) {//new score
73 uscoreFreq[userData[0]] = 1;
74 }else{ uscoreFreq[userData[0]]++; }
76 //add users score to valid scores
77 validScores[userData[0]] = userData[0];
79 //save score for summary file
80 userTreeScores.push_back(userData[0]);
84 //get pscores for random trees
85 for (int j = 0; j < iters; j++) {
86 //create new tree with same num nodes and leaves as users
88 //create random relationships between nodes
89 randT->assembleRandomTree();
90 //get pscore of random tree
91 randomData = pars->getValues(randT);
93 //add trees pscore to map of scores
94 it2 = rscoreFreq.find(randomData[0]);
95 if (it2 != rscoreFreq.end()) {//already have that score
96 rscoreFreq[randomData[0]]++;
97 }else{//first time we have seen this score
98 rscoreFreq[randomData[0]] = 1;
101 //add randoms score to validscores
102 validScores[randomData[0]] = randomData[0];
104 //output info to pdistrib file
105 outDist << j+1 << '\t'<< '\t' << randomData[0] << endl;
110 //get pscores for random trees
111 for (int j = 0; j < iters; j++) {
112 //create new tree with same num nodes and leaves as users
114 //create random relationships between nodes
115 randT->assembleRandomTree();
116 //get pscore of random tree
117 randomData = pars->getValues(randT);
119 //add trees pscore to map of scores
120 it2 = rscoreFreq.find(randomData[0]);
121 if (it2 != rscoreFreq.end()) {//already have that score
122 rscoreFreq[randomData[0]]++;
123 }else{//first time we have seen this score
124 rscoreFreq[randomData[0]] = 1;
127 //add randoms score to validscores
128 validScores[randomData[0]] = randomData[0];
134 float rcumul = 0.0000;
135 float ucumul = 0.0000;
137 //this loop fills the cumulative maps and put 0.0000 in the score freq map to make it easier to print.
138 for (it = validScores.begin(); it != validScores.end(); it++) {
139 if (randomtree == "") {
140 it2 = uscoreFreq.find(it->first);
141 //user data has that score
142 if (it2 != uscoreFreq.end()) { uscoreFreq[it->first] /= T.size(); ucumul+= it2->second; }
143 else { uscoreFreq[it->first] = 0.0000; } //no user trees with that score
145 uCumul[it->first] = ucumul;
148 //make rscoreFreq map and rCumul
149 it2 = rscoreFreq.find(it->first);
150 //get percentage of random trees with that info
151 if (it2 != rscoreFreq.end()) { rscoreFreq[it->first] /= iters; rcumul+= it2->second; }
152 else { rscoreFreq[it->first] = 0.0000; } //no random trees with that score
153 rCumul[it->first] = rcumul;
156 //find the signifigance of each user trees score when compared to the random trees and save for printing the summary file
157 for (int h = 0; h < userTreeScores.size(); h++) {
158 UScoreSig.push_back(rCumul[userTreeScores[h]]);
161 printParsimonyFile();
164 //reset globaldata's treemap if you just did random distrib
165 if (randomtree != "") { globaldata->gTreemap = savetmap; }
167 //reset randomTree parameter to ""
168 globaldata->setRandomTree("");
169 //reset groups parameter
170 globaldata->Groups.clear();
175 catch(exception& e) {
176 cout << "Standard Error: " << e.what() << " has occurred in the ParsimonyCommand class Function execute. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
180 cout << "An unknown error has occurred in the ParsimonyCommand class function execute. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
185 /***********************************************************/
186 void ParsimonyCommand::printParsimonyFile() {
189 if (randomtree == "") {
190 out << "Score" << '\t' << "UserFreq" << '\t' << "UserCumul" << '\t' << "RandFreq" << '\t' << "RandCumul" << endl;
192 out << "Score" << '\t' << "RandFreq" << '\t' << "RandCumul" << endl;
196 out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
199 for (it = validScores.begin(); it != validScores.end(); it++) {
200 if (randomtree == "") {
201 out << setprecision(6) << it->first << '\t' << '\t' << uscoreFreq[it->first] << '\t' << uCumul[it->first] << '\t' << rscoreFreq[it->first] << '\t' << rCumul[it->first] << endl;
203 out << setprecision(6) << it->first << '\t' << '\t' << rscoreFreq[it->first] << '\t' << rCumul[it->first] << endl;
210 catch(exception& e) {
211 cout << "Standard Error: " << e.what() << " has occurred in the ParsimonyCommand class Function printParsimonyFile. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
215 cout << "An unknown error has occurred in the ParsimonyCommand class function printParsimonyFile. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
219 /***********************************************************/
220 void ParsimonyCommand::printUSummaryFile() {
223 outSum << "Tree#" << '\t' << "ParsScore" << '\t' << '\t' << "ParsSig" << endl;
226 outSum.setf(ios::fixed, ios::floatfield); outSum.setf(ios::showpoint);
229 for (int i = 0; i< T.size(); i++) {
230 outSum << setprecision(6) << i+1 << '\t' << '\t' << userTreeScores[i] << '\t' << UScoreSig[i] << endl;
235 catch(exception& e) {
236 cout << "Standard Error: " << e.what() << " has occurred in the ParsimonyCommand class Function printUSummaryFile. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
240 cout << "An unknown error has occurred in the ParsimonyCommand class function printUSummaryFile. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
245 /***********************************************************/
246 void ParsimonyCommand::getUserInput() {
250 tmap = new TreeMap();
252 cout << "Please enter the number of groups you would like to analyze: ";
257 numEachGroup.resize(numGroups, 0);
259 for (int i = 1; i <= numGroups; i++) {
260 cout << "Please enter the number of sequences in group " << i << ": ";
263 //set tmaps seqsPerGroup
264 tmap->seqsPerGroup[toString(i)] = num;
266 //set tmaps namesOfSeqs
267 for (int j = 0; j < num; j++) {
268 tmap->namesOfSeqs.push_back(toString(count));
269 tmap->treemap[toString(count)].groupname = toString(i);
274 //clears buffer so next command doesn't have error
278 //save tmap for later
279 globaldata->gTreemap = tmap;
282 catch(exception& e) {
283 cout << "Standard Error: " << e.what() << " has occurred in the ParsimonyCommand class Function getUserInput. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
287 cout << "An unknown error has occurred in the ParsimonyCommand class function getUserInput. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
291 /***********************************************************/
293 void ParsimonyCommand::setGroups() {
295 //if the user has not entered specific groups to analyze then do them all
296 if (globaldata->Groups.size() != 0) {
297 //check that groups are valid
298 for (int i = 0; i < globaldata->Groups.size(); i++) {
299 if (tmap->isValidGroup(globaldata->Groups[i]) != true) {
300 cout << globaldata->Groups[i] << " is not a valid group, and will be disregarded." << endl;
301 // erase the invalid group from globaldata->Groups
302 globaldata->Groups.erase (globaldata->Groups.begin()+i);
306 //if the user only entered invalid groups
307 if (globaldata->Groups.size() == 0) {
308 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;
309 for (int i = 0; i < tmap->namesOfGroups.size(); i++) {
310 globaldata->Groups.push_back(tmap->namesOfGroups[i]);
315 for (int i = 0; i < tmap->namesOfGroups.size(); i++) {
316 globaldata->Groups.push_back(tmap->namesOfGroups[i]);
320 catch(exception& e) {
321 cout << "Standard Error: " << e.what() << " has occurred in the ParsimonyCommand class Function setGroups. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
325 cout << "An unknown error has occurred in the ParsimonyCommand class function setGroups. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
330 /*****************************************************************/