2 * unifracunweightedcommand.cpp
5 * Created by Sarah Westcott on 2/9/09.
6 * Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
10 #include "unifracunweightedcommand.h"
12 //**********************************************************************************************************************
13 vector<string> UnifracUnweightedCommand::setParameters(){
15 CommandParameter ptree("tree", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(ptree);
16 CommandParameter pgroup("group", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(pgroup);
17 CommandParameter pname("name", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(pname);
18 CommandParameter pgroups("groups", "String", "", "", "", "", "",false,false); parameters.push_back(pgroups);
19 CommandParameter piters("iters", "Number", "", "1000", "", "", "",false,false); parameters.push_back(piters);
20 CommandParameter pprocessors("processors", "Number", "", "1", "", "", "",false,false); parameters.push_back(pprocessors);
21 CommandParameter prandom("random", "Boolean", "", "F", "", "", "",false,false); parameters.push_back(prandom);
22 CommandParameter pdistance("distance", "Multiple", "column-lt-square", "column", "", "", "",false,false); parameters.push_back(pdistance);
23 CommandParameter proot("root", "Boolean", "F", "", "", "", "",false,false); parameters.push_back(proot);
24 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
25 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
27 vector<string> myArray;
28 for (int i = 0; i < parameters.size(); i++) { myArray.push_back(parameters[i].name); }
32 m->errorOut(e, "UnifracUnweightedCommand", "setParameters");
36 //**********************************************************************************************************************
37 string UnifracUnweightedCommand::getHelpString(){
39 string helpString = "";
40 helpString += "The unifrac.unweighted command parameters are tree, group, name, groups, iters, distance, processors, root and random. tree parameter is required unless you have valid current tree file.\n";
41 helpString += "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";
42 helpString += "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";
43 helpString += "The distance parameter allows you to create a distance file from the results. The default is false. You may set distance to lt, square or column.\n";
44 helpString += "The random parameter allows you to shut off the comparison to random trees. The default is false, meaning compare don't your trees with randomly generated trees.\n";
45 helpString += "The root parameter allows you to include the entire root in your calculations. The default is false, meaning stop at the root for this comparision instead of the root of the entire tree.\n";
46 helpString += "The processors parameter allows you to specify the number of processors to use. The default is 1.\n";
47 helpString += "The unifrac.unweighted command should be in the following format: unifrac.unweighted(groups=yourGroups, iters=yourIters).\n";
48 helpString += "Example unifrac.unweighted(groups=A-B-C, iters=500).\n";
49 helpString += "The default value for groups is all the groups in your groupfile, and iters is 1000.\n";
50 helpString += "The unifrac.unweighted command output two files: .unweighted and .uwsummary their descriptions are in the manual.\n";
51 helpString += "Note: No spaces between parameter labels (i.e. groups), '=' and parameters (i.e.yourGroups).\n";
55 m->errorOut(e, "UnifracUnweightedCommand", "getHelpString");
59 //**********************************************************************************************************************
60 UnifracUnweightedCommand::UnifracUnweightedCommand(){
62 abort = true; calledHelp = true;
64 vector<string> tempOutNames;
65 outputTypes["unweighted"] = tempOutNames;
66 outputTypes["uwsummary"] = tempOutNames;
67 outputTypes["phylip"] = tempOutNames;
68 outputTypes["column"] = tempOutNames;
71 m->errorOut(e, "UnifracUnweightedCommand", "UnifracUnweightedCommand");
75 /***********************************************************/
76 UnifracUnweightedCommand::UnifracUnweightedCommand(string option) {
78 abort = false; calledHelp = false;
81 //allow user to run help
82 if(option == "help") { help(); abort = true; calledHelp = true; }
85 vector<string> myArray = setParameters();
87 OptionParser parser(option);
88 map<string,string> parameters = parser.getParameters();
89 map<string,string>::iterator it;
91 ValidParameters validParameter;
93 //check to make sure all parameters are valid for command
94 for (map<string,string>::iterator it = parameters.begin(); it != parameters.end(); it++) {
95 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
98 //initialize outputTypes
99 vector<string> tempOutNames;
100 outputTypes["unweighted"] = tempOutNames;
101 outputTypes["uwsummary"] = tempOutNames;
102 outputTypes["phylip"] = tempOutNames;
103 outputTypes["column"] = tempOutNames;
105 //if the user changes the input directory command factory will send this info to us in the output parameter
106 string inputDir = validParameter.validFile(parameters, "inputdir", false);
107 if (inputDir == "not found"){ inputDir = ""; }
110 it = parameters.find("tree");
111 //user has given a template file
112 if(it != parameters.end()){
113 path = m->hasPath(it->second);
114 //if the user has not given a path then, add inputdir. else leave path alone.
115 if (path == "") { parameters["tree"] = inputDir + it->second; }
118 it = parameters.find("group");
119 //user has given a template file
120 if(it != parameters.end()){
121 path = m->hasPath(it->second);
122 //if the user has not given a path then, add inputdir. else leave path alone.
123 if (path == "") { parameters["group"] = inputDir + it->second; }
126 it = parameters.find("name");
127 //user has given a template file
128 if(it != parameters.end()){
129 path = m->hasPath(it->second);
130 //if the user has not given a path then, add inputdir. else leave path alone.
131 if (path == "") { parameters["name"] = inputDir + it->second; }
137 //check for required parameters
138 treefile = validParameter.validFile(parameters, "tree", true);
139 if (treefile == "not open") { abort = true; }
140 else if (treefile == "not found") { //if there is a current design file, use it
141 treefile = m->getTreeFile();
142 if (treefile != "") { m->mothurOut("Using " + treefile + " as input file for the tree parameter."); m->mothurOutEndLine(); }
143 else { m->mothurOut("You have no current tree file and the tree parameter is required."); m->mothurOutEndLine(); abort = true; }
146 //check for required parameters
147 groupfile = validParameter.validFile(parameters, "group", true);
148 if (groupfile == "not open") { abort = true; }
149 else if (groupfile == "not found") { groupfile = ""; }
151 namefile = validParameter.validFile(parameters, "name", true);
152 if (namefile == "not open") { abort = true; }
153 else if (namefile == "not found") { namefile = ""; }
155 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){ outputDir = ""; }
157 //check for optional parameter and set defaults
158 // ...at some point should added some additional type checking...
159 groups = validParameter.validFile(parameters, "groups", false);
160 if (groups == "not found") { groups = ""; }
162 m->splitAtDash(groups, Groups);
166 itersString = validParameter.validFile(parameters, "iters", false); if (itersString == "not found") { itersString = "1000"; }
167 convert(itersString, iters);
169 string temp = validParameter.validFile(parameters, "distance", false);
170 if (temp == "not found") { phylip = false; outputForm = ""; }
172 if ((temp == "lt") || (temp == "column") || (temp == "square")) { phylip = true; outputForm = temp; }
173 else { m->mothurOut("Options for distance are: lt, square, or column. Using lt."); m->mothurOutEndLine(); phylip = true; outputForm = "lt"; }
176 temp = validParameter.validFile(parameters, "random", false); if (temp == "not found") { temp = "f"; }
177 random = m->isTrue(temp);
179 temp = validParameter.validFile(parameters, "root", false); if (temp == "not found") { temp = "F"; }
180 includeRoot = m->isTrue(temp);
182 temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); }
183 m->setProcessors(temp);
184 convert(temp, processors);
186 if (!random) { iters = 0; } //turn off random calcs
188 //if user selects distance = true and no groups it won't calc the pairwise
189 if ((phylip) && (Groups.size() == 0)) {
191 m->splitAtDash(groups, Groups);
197 catch(exception& e) {
198 m->errorOut(e, "UnifracUnweightedCommand", "UnifracUnweightedCommand");
203 /***********************************************************/
204 int UnifracUnweightedCommand::execute() {
207 if (abort == true) { if (calledHelp) { return 0; } return 2; }
209 if (groupfile != "") {
210 //read in group map info.
211 tmap = new TreeMap(groupfile);
213 }else{ //fake out by putting everyone in one group
214 Tree* tree = new Tree(treefile); delete tree; //extracts names from tree to make faked out groupmap
215 tmap = new TreeMap();
217 for (int i = 0; i < m->Treenames.size(); i++) { tmap->addSeq(m->Treenames[i], "Group1"); }
220 if (namefile != "") { readNamesFile(); }
222 read = new ReadNewickTree(treefile);
223 int readOk = read->read(tmap);
225 if (readOk != 0) { m->mothurOut("Read Terminated."); m->mothurOutEndLine(); delete tmap; delete read; return 0; }
227 read->AssembleTrees();
228 T = read->getTrees();
231 //make sure all files match
232 //if you provide a namefile we will use the numNames in the namefile as long as the number of unique match the tree names size.
234 if (namefile != "") {
235 if (numUniquesInName == m->Treenames.size()) { numNamesInTree = nameMap.size(); }
236 else { numNamesInTree = m->Treenames.size(); }
237 }else { numNamesInTree = m->Treenames.size(); }
240 //output any names that are in group file but not in tree
241 if (numNamesInTree < tmap->getNumSeqs()) {
242 for (int i = 0; i < tmap->namesOfSeqs.size(); i++) {
243 //is that name in the tree?
245 for (int j = 0; j < m->Treenames.size(); j++) {
246 if (tmap->namesOfSeqs[i] == m->Treenames[j]) { break; } //found it
250 if (m->control_pressed) {
251 delete tmap; for (int i = 0; i < T.size(); i++) { delete T[i]; }
252 for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } outputTypes.clear();
257 //then you did not find it so report it
258 if (count == m->Treenames.size()) {
259 //if it is in your namefile then don't remove
260 map<string, string>::iterator it = nameMap.find(tmap->namesOfSeqs[i]);
262 if (it == nameMap.end()) {
263 m->mothurOut(tmap->namesOfSeqs[i] + " is in your groupfile and not in your tree. It will be disregarded."); m->mothurOutEndLine();
264 tmap->removeSeq(tmap->namesOfSeqs[i]);
265 i--; //need this because removeSeq removes name from namesOfSeqs
271 sumFile = outputDir + m->getSimpleName(treefile) + ".uwsummary";
272 outputNames.push_back(sumFile); outputTypes["uwsummary"].push_back(sumFile);
273 m->openOutputFile(sumFile, outSum);
275 util = new SharedUtil();
276 util->setGroups(m->Groups, tmap->namesOfGroups, allGroups, numGroups, "unweighted"); //sets the groups the user wants to analyze
277 util->getCombos(groupComb, m->Groups, numComp);
280 if (numGroups == 1) { numComp++; groupComb.push_back(allGroups); }
282 unweighted = new Unweighted(tmap, includeRoot);
284 int start = time(NULL);
286 userData.resize(numComp,0); //data[0] = unweightedscore
287 randomData.resize(numComp,0); //data[0] = unweightedscore
288 //create new tree with same num nodes and leaves as users
290 if (numComp < processors) { processors = numComp; }
292 outSum << "Tree#" << '\t' << "Groups" << '\t' << "UWScore" <<'\t' << "UWSig" << endl;
293 m->mothurOut("Tree#\tGroups\tUWScore\tUWSig"); m->mothurOutEndLine();
295 //get pscores for users trees
296 for (int i = 0; i < T.size(); i++) {
297 if (m->control_pressed) {
298 delete tmap; delete unweighted;
299 for (int i = 0; i < T.size(); i++) { delete T[i]; }
301 for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); }
308 output = new ColumnFile(outputDir + m->getSimpleName(treefile) + toString(i+1) + ".unweighted", itersString);
309 outputNames.push_back(outputDir + m->getSimpleName(treefile) + toString(i+1) + ".unweighted");
310 outputTypes["unweighted"].push_back(outputDir + m->getSimpleName(treefile) + toString(i+1) + ".unweighted");
314 //get unweighted for users tree
315 rscoreFreq.resize(numComp);
316 rCumul.resize(numComp);
317 utreeScores.resize(numComp);
318 UWScoreSig.resize(numComp);
320 userData = unweighted->getValues(T[i], processors, outputDir); //userData[0] = unweightedscore
322 if (m->control_pressed) { delete tmap; delete unweighted;
323 for (int i = 0; i < T.size(); i++) { delete T[i]; }if (random) { delete output; } outSum.close(); for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); }return 0; }
325 //output scores for each combination
326 for(int k = 0; k < numComp; k++) {
328 utreeScores[k].push_back(userData[k]);
330 //add users score to validscores
331 validScores[userData[k]] = userData[k];
334 //get unweighted scores for random trees - if random is false iters = 0
335 for (int j = 0; j < iters; j++) {
337 //we need a different getValues because when we swap the labels we only want to swap those in each pairwise comparison
338 randomData = unweighted->getValues(T[i], "", "", processors, outputDir);
340 if (m->control_pressed) { delete tmap; delete unweighted;
341 for (int i = 0; i < T.size(); i++) { delete T[i]; }if (random) { delete output; } outSum.close(); for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
343 for(int k = 0; k < numComp; k++) {
344 //add trees unweighted score to map of scores
345 map<float,float>::iterator it = rscoreFreq[k].find(randomData[k]);
346 if (it != rscoreFreq[k].end()) {//already have that score
347 rscoreFreq[k][randomData[k]]++;
348 }else{//first time we have seen this score
349 rscoreFreq[k][randomData[k]] = 1;
352 //add randoms score to validscores
353 validScores[randomData[k]] = randomData[k];
357 // m->mothurOut("Iter: " + toString(j+1)); m->mothurOutEndLine();
360 for(int a = 0; a < numComp; a++) {
361 float rcumul = 1.0000;
364 //this loop fills the cumulative maps and put 0.0000 in the score freq map to make it easier to print.
365 for (map<float,float>::iterator it = validScores.begin(); it != validScores.end(); it++) {
366 //make rscoreFreq map and rCumul
367 map<float,float>::iterator it2 = rscoreFreq[a].find(it->first);
368 rCumul[a][it->first] = rcumul;
369 //get percentage of random trees with that info
370 if (it2 != rscoreFreq[a].end()) { rscoreFreq[a][it->first] /= iters; rcumul-= it2->second; }
371 else { rscoreFreq[a][it->first] = 0.0000; } //no random trees with that score
373 UWScoreSig[a].push_back(rCumul[a][userData[a]]);
374 }else { UWScoreSig[a].push_back(0.0); }
378 if (m->control_pressed) { delete tmap; delete unweighted;
379 for (int i = 0; i < T.size(); i++) { delete T[i]; }if (random) { delete output; } outSum.close(); for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
382 printUWSummaryFile(i);
383 if (random) { printUnweightedFile(); delete output; }
384 if (phylip) { createPhylipFile(i); }
396 delete tmap; delete unweighted;
397 for (int i = 0; i < T.size(); i++) { delete T[i]; }
399 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
401 m->mothurOut("It took " + toString(time(NULL) - start) + " secs to run unifrac.unweighted."); m->mothurOutEndLine();
403 //set phylip file as new current phylipfile
405 itTypes = outputTypes.find("phylip");
406 if (itTypes != outputTypes.end()) {
407 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setPhylipFile(current); }
410 //set column file as new current columnfile
411 itTypes = outputTypes.find("column");
412 if (itTypes != outputTypes.end()) {
413 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setColumnFile(current); }
416 m->mothurOutEndLine();
417 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
418 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
419 m->mothurOutEndLine();
424 catch(exception& e) {
425 m->errorOut(e, "UnifracUnweightedCommand", "execute");
429 /***********************************************************/
430 void UnifracUnweightedCommand::printUnweightedFile() {
435 tags.push_back("Score");
436 tags.push_back("RandFreq"); tags.push_back("RandCumul");
438 for(int a = 0; a < numComp; a++) {
439 output->initFile(groupComb[a], tags);
441 for (map<float,float>::iterator it = validScores.begin(); it != validScores.end(); it++) {
442 data.push_back(it->first); data.push_back(rscoreFreq[a][it->first]); data.push_back(rCumul[a][it->first]);
443 output->output(data);
449 catch(exception& e) {
450 m->errorOut(e, "UnifracUnweightedCommand", "printUnweightedFile");
455 /***********************************************************/
456 void UnifracUnweightedCommand::printUWSummaryFile(int i) {
460 outSum.setf(ios::fixed, ios::floatfield); outSum.setf(ios::showpoint);
464 for(int a = 0; a < numComp; a++) {
465 outSum << i+1 << '\t';
466 m->mothurOut(toString(i+1) + "\t");
469 if (UWScoreSig[a][0] > (1/(float)iters)) {
470 outSum << setprecision(6) << groupComb[a] << '\t' << utreeScores[a][0] << '\t' << setprecision(itersString.length()) << UWScoreSig[a][0] << endl;
471 cout << setprecision(6) << groupComb[a] << '\t' << utreeScores[a][0] << '\t' << setprecision(itersString.length()) << UWScoreSig[a][0] << endl;
472 m->mothurOutJustToLog(groupComb[a] + "\t" + toString(utreeScores[a][0]) + "\t" + toString(UWScoreSig[a][0])+ "\n");
474 outSum << setprecision(6) << groupComb[a] << '\t' << utreeScores[a][0] << '\t' << setprecision(itersString.length()) << "<" << (1/float(iters)) << endl;
475 cout << setprecision(6) << groupComb[a] << '\t' << utreeScores[a][0] << '\t' << setprecision(itersString.length()) << "<" << (1/float(iters)) << endl;
476 m->mothurOutJustToLog(groupComb[a] + "\t" + toString(utreeScores[a][0]) + "\t<" + toString((1/float(iters))) + "\n");
479 outSum << setprecision(6) << groupComb[a] << '\t' << utreeScores[a][0] << '\t' << "0.00" << endl;
480 cout << setprecision(6) << groupComb[a] << '\t' << utreeScores[a][0] << '\t' << "0.00" << endl;
481 m->mothurOutJustToLog(groupComb[a] + "\t" + toString(utreeScores[a][0]) + "\t0.00\n");
486 catch(exception& e) {
487 m->errorOut(e, "UnifracUnweightedCommand", "printUWSummaryFile");
491 /***********************************************************/
492 void UnifracUnweightedCommand::createPhylipFile(int i) {
494 string phylipFileName;
495 if ((outputForm == "lt") || (outputForm == "square")) {
496 phylipFileName = outputDir + m->getSimpleName(treefile) + toString(i+1) + ".unweighted.phylip.dist";
497 outputNames.push_back(phylipFileName); outputTypes["phylip"].push_back(phylipFileName);
499 phylipFileName = outputDir + m->getSimpleName(treefile) + toString(i+1) + ".unweighted.column.dist";
500 outputNames.push_back(phylipFileName); outputTypes["column"].push_back(phylipFileName);
504 m->openOutputFile(phylipFileName, out);
506 if ((outputForm == "lt") || (outputForm == "square")) {
508 out << m->Groups.size() << endl;
511 //make matrix with scores in it
512 vector< vector<float> > dists; dists.resize(m->Groups.size());
513 for (int i = 0; i < m->Groups.size(); i++) {
514 dists[i].resize(m->Groups.size(), 0.0);
517 //flip it so you can print it
519 for (int r=0; r<m->Groups.size(); r++) {
520 for (int l = 0; l < r; l++) {
521 dists[r][l] = utreeScores[count][0];
522 dists[l][r] = utreeScores[count][0];
528 for (int r=0; r<m->Groups.size(); r++) {
530 string name = m->Groups[r];
531 if (name.length() < 10) { //pad with spaces to make compatible
532 while (name.length() < 10) { name += " "; }
535 if (outputForm == "lt") {
539 for (int l = 0; l < r; l++) { out << dists[r][l] << '\t'; }
541 }else if (outputForm == "square") {
545 for (int l = 0; l < m->Groups.size(); l++) { out << dists[r][l] << '\t'; }
549 for (int l = 0; l < r; l++) {
550 string otherName = m->Groups[l];
551 if (otherName.length() < 10) { //pad with spaces to make compatible
552 while (otherName.length() < 10) { otherName += " "; }
555 out << name << '\t' << otherName << dists[r][l] << endl;
561 catch(exception& e) {
562 m->errorOut(e, "UnifracUnweightedCommand", "createPhylipFile");
565 }/*****************************************************************/
566 int UnifracUnweightedCommand::readNamesFile() {
569 numUniquesInName = 0;
572 m->openInputFile(namefile, in);
574 string first, second;
575 map<string, string>::iterator itNames;
578 in >> first >> second; m->gobble(in);
582 itNames = m->names.find(first);
583 if (itNames == m->names.end()) {
584 m->names[first] = second;
586 //we need a list of names in your namefile to use above when removing extra seqs above so we don't remove them
587 vector<string> dupNames;
588 m->splitAtComma(second, dupNames);
590 for (int i = 0; i < dupNames.size(); i++) {
591 nameMap[dupNames[i]] = dupNames[i];
592 if ((groupfile == "") && (i != 0)) { tmap->addSeq(dupNames[i], "Group1"); }
594 }else { m->mothurOut(first + " has already been seen in namefile, disregarding names file."); m->mothurOutEndLine(); in.close(); m->names.clear(); namefile = ""; return 1; }
600 catch(exception& e) {
601 m->errorOut(e, "UnifracUnweightedCommand", "readNamesFile");
605 /***********************************************************/