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\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 //valid paramters for this command
86 string Array[] = {"groups","iters","distance","random","root", "processors","outputdir","inputdir"};
87 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
89 OptionParser parser(option);
90 map<string,string> parameters = parser.getParameters();
91 map<string,string>::iterator it;
93 ValidParameters validParameter;
95 //check to make sure all parameters are valid for command
96 for (map<string,string>::iterator it = parameters.begin(); it != parameters.end(); it++) {
97 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
100 //initialize outputTypes
101 vector<string> tempOutNames;
102 outputTypes["unweighted"] = tempOutNames;
103 outputTypes["uwsummary"] = tempOutNames;
104 outputTypes["phylip"] = tempOutNames;
105 outputTypes["column"] = tempOutNames;
107 //if the user changes the input directory command factory will send this info to us in the output parameter
108 string inputDir = validParameter.validFile(parameters, "inputdir", false);
109 if (inputDir == "not found"){ inputDir = ""; }
112 it = parameters.find("tree");
113 //user has given a template file
114 if(it != parameters.end()){
115 path = m->hasPath(it->second);
116 //if the user has not given a path then, add inputdir. else leave path alone.
117 if (path == "") { parameters["tree"] = inputDir + it->second; }
120 it = parameters.find("group");
121 //user has given a template file
122 if(it != parameters.end()){
123 path = m->hasPath(it->second);
124 //if the user has not given a path then, add inputdir. else leave path alone.
125 if (path == "") { parameters["group"] = inputDir + it->second; }
128 it = parameters.find("name");
129 //user has given a template file
130 if(it != parameters.end()){
131 path = m->hasPath(it->second);
132 //if the user has not given a path then, add inputdir. else leave path alone.
133 if (path == "") { parameters["name"] = inputDir + it->second; }
139 //check for required parameters
140 treefile = validParameter.validFile(parameters, "tree", true);
141 if (treefile == "not open") { abort = true; }
142 else if (treefile == "not found") { //if there is a current design file, use it
143 treefile = m->getTreeFile();
144 if (treefile != "") { m->mothurOut("Using " + treefile + " as input file for the tree parameter."); m->mothurOutEndLine(); }
145 else { m->mothurOut("You have no current tree file and the tree parameter is required."); m->mothurOutEndLine(); abort = true; }
148 //check for required parameters
149 groupfile = validParameter.validFile(parameters, "group", true);
150 if (groupfile == "not open") { abort = true; }
151 else if (groupfile == "not found") { groupfile = ""; }
153 namefile = validParameter.validFile(parameters, "name", true);
154 if (namefile == "not open") { abort = true; }
155 else if (namefile == "not found") { namefile = ""; }
157 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){ outputDir = ""; }
159 //check for optional parameter and set defaults
160 // ...at some point should added some additional type checking...
161 groups = validParameter.validFile(parameters, "groups", false);
162 if (groups == "not found") { groups = ""; }
164 m->splitAtDash(groups, Groups);
168 itersString = validParameter.validFile(parameters, "iters", false); if (itersString == "not found") { itersString = "1000"; }
169 convert(itersString, iters);
171 string temp = validParameter.validFile(parameters, "distance", false);
172 if (temp == "not found") { phylip = false; outputForm = ""; }
174 if ((temp == "lt") || (temp == "column") || (temp == "square")) { phylip = true; outputForm = temp; }
175 else { m->mothurOut("Options for distance are: lt, square, or column. Using lt."); m->mothurOutEndLine(); phylip = true; outputForm = "lt"; }
178 temp = validParameter.validFile(parameters, "random", false); if (temp == "not found") { temp = "f"; }
179 random = m->isTrue(temp);
181 temp = validParameter.validFile(parameters, "root", false); if (temp == "not found") { temp = "F"; }
182 includeRoot = m->isTrue(temp);
184 temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); }
185 m->setProcessors(temp);
186 convert(temp, processors);
188 if (!random) { iters = 0; } //turn off random calcs
190 //if user selects distance = true and no groups it won't calc the pairwise
191 if ((phylip) && (Groups.size() == 0)) {
193 m->splitAtDash(groups, Groups);
199 catch(exception& e) {
200 m->errorOut(e, "UnifracUnweightedCommand", "UnifracUnweightedCommand");
205 /***********************************************************/
206 int UnifracUnweightedCommand::execute() {
209 if (abort == true) { if (calledHelp) { return 0; } return 2; }
211 if (groupfile != "") {
212 //read in group map info.
213 tmap = new TreeMap(groupfile);
215 }else{ //fake out by putting everyone in one group
216 Tree* tree = new Tree(treefile); delete tree; //extracts names from tree to make faked out groupmap
217 tmap = new TreeMap();
219 for (int i = 0; i < m->Treenames.size(); i++) { tmap->addSeq(m->Treenames[i], "Group1"); }
222 if (namefile != "") { readNamesFile(); }
224 read = new ReadNewickTree(treefile);
225 int readOk = read->read(tmap);
227 if (readOk != 0) { m->mothurOut("Read Terminated."); m->mothurOutEndLine(); delete tmap; delete read; return 0; }
229 read->AssembleTrees();
230 T = read->getTrees();
233 //make sure all files match
234 //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.
236 if (namefile != "") {
237 if (numUniquesInName == m->Treenames.size()) { numNamesInTree = nameMap.size(); }
238 else { numNamesInTree = m->Treenames.size(); }
239 }else { numNamesInTree = m->Treenames.size(); }
242 //output any names that are in group file but not in tree
243 if (numNamesInTree < tmap->getNumSeqs()) {
244 for (int i = 0; i < tmap->namesOfSeqs.size(); i++) {
245 //is that name in the tree?
247 for (int j = 0; j < m->Treenames.size(); j++) {
248 if (tmap->namesOfSeqs[i] == m->Treenames[j]) { break; } //found it
252 if (m->control_pressed) {
253 delete tmap; for (int i = 0; i < T.size(); i++) { delete T[i]; }
254 for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } outputTypes.clear();
259 //then you did not find it so report it
260 if (count == m->Treenames.size()) {
261 //if it is in your namefile then don't remove
262 map<string, string>::iterator it = nameMap.find(tmap->namesOfSeqs[i]);
264 if (it == nameMap.end()) {
265 m->mothurOut(tmap->namesOfSeqs[i] + " is in your groupfile and not in your tree. It will be disregarded."); m->mothurOutEndLine();
266 tmap->removeSeq(tmap->namesOfSeqs[i]);
267 i--; //need this because removeSeq removes name from namesOfSeqs
273 sumFile = outputDir + m->getSimpleName(treefile) + ".uwsummary";
274 outputNames.push_back(sumFile); outputTypes["uwsummary"].push_back(sumFile);
275 m->openOutputFile(sumFile, outSum);
277 util = new SharedUtil();
278 util->setGroups(m->Groups, tmap->namesOfGroups, allGroups, numGroups, "unweighted"); //sets the groups the user wants to analyze
279 util->getCombos(groupComb, m->Groups, numComp);
282 if (numGroups == 1) { numComp++; groupComb.push_back(allGroups); }
284 unweighted = new Unweighted(tmap, includeRoot);
286 int start = time(NULL);
288 userData.resize(numComp,0); //data[0] = unweightedscore
289 randomData.resize(numComp,0); //data[0] = unweightedscore
290 //create new tree with same num nodes and leaves as users
292 if (numComp < processors) { processors = numComp; }
294 outSum << "Tree#" << '\t' << "Groups" << '\t' << "UWScore" <<'\t' << "UWSig" << endl;
295 m->mothurOut("Tree#\tGroups\tUWScore\tUWSig"); m->mothurOutEndLine();
297 //get pscores for users trees
298 for (int i = 0; i < T.size(); i++) {
299 if (m->control_pressed) {
300 delete tmap; delete unweighted;
301 for (int i = 0; i < T.size(); i++) { delete T[i]; }
303 for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); }
310 output = new ColumnFile(outputDir + m->getSimpleName(treefile) + toString(i+1) + ".unweighted", itersString);
311 outputNames.push_back(outputDir + m->getSimpleName(treefile) + toString(i+1) + ".unweighted");
312 outputTypes["unweighted"].push_back(outputDir + m->getSimpleName(treefile) + toString(i+1) + ".unweighted");
316 //get unweighted for users tree
317 rscoreFreq.resize(numComp);
318 rCumul.resize(numComp);
319 utreeScores.resize(numComp);
320 UWScoreSig.resize(numComp);
322 userData = unweighted->getValues(T[i], processors, outputDir); //userData[0] = unweightedscore
324 if (m->control_pressed) { delete tmap; delete unweighted;
325 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; }
327 //output scores for each combination
328 for(int k = 0; k < numComp; k++) {
330 utreeScores[k].push_back(userData[k]);
332 //add users score to validscores
333 validScores[userData[k]] = userData[k];
336 //get unweighted scores for random trees - if random is false iters = 0
337 for (int j = 0; j < iters; j++) {
339 //we need a different getValues because when we swap the labels we only want to swap those in each pairwise comparison
340 randomData = unweighted->getValues(T[i], "", "", processors, outputDir);
342 if (m->control_pressed) { delete tmap; delete unweighted;
343 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; }
345 for(int k = 0; k < numComp; k++) {
346 //add trees unweighted score to map of scores
347 map<float,float>::iterator it = rscoreFreq[k].find(randomData[k]);
348 if (it != rscoreFreq[k].end()) {//already have that score
349 rscoreFreq[k][randomData[k]]++;
350 }else{//first time we have seen this score
351 rscoreFreq[k][randomData[k]] = 1;
354 //add randoms score to validscores
355 validScores[randomData[k]] = randomData[k];
359 // m->mothurOut("Iter: " + toString(j+1)); m->mothurOutEndLine();
362 for(int a = 0; a < numComp; a++) {
363 float rcumul = 1.0000;
366 //this loop fills the cumulative maps and put 0.0000 in the score freq map to make it easier to print.
367 for (map<float,float>::iterator it = validScores.begin(); it != validScores.end(); it++) {
368 //make rscoreFreq map and rCumul
369 map<float,float>::iterator it2 = rscoreFreq[a].find(it->first);
370 rCumul[a][it->first] = rcumul;
371 //get percentage of random trees with that info
372 if (it2 != rscoreFreq[a].end()) { rscoreFreq[a][it->first] /= iters; rcumul-= it2->second; }
373 else { rscoreFreq[a][it->first] = 0.0000; } //no random trees with that score
375 UWScoreSig[a].push_back(rCumul[a][userData[a]]);
376 }else { UWScoreSig[a].push_back(0.0); }
380 if (m->control_pressed) { delete tmap; delete unweighted;
381 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; }
384 printUWSummaryFile(i);
385 if (random) { printUnweightedFile(); delete output; }
386 if (phylip) { createPhylipFile(i); }
398 delete tmap; delete unweighted;
399 for (int i = 0; i < T.size(); i++) { delete T[i]; }
401 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
403 m->mothurOut("It took " + toString(time(NULL) - start) + " secs to run unifrac.unweighted."); m->mothurOutEndLine();
405 //set phylip file as new current phylipfile
407 itTypes = outputTypes.find("phylip");
408 if (itTypes != outputTypes.end()) {
409 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setPhylipFile(current); }
412 //set column file as new current columnfile
413 itTypes = outputTypes.find("column");
414 if (itTypes != outputTypes.end()) {
415 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setColumnFile(current); }
418 m->mothurOutEndLine();
419 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
420 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
421 m->mothurOutEndLine();
426 catch(exception& e) {
427 m->errorOut(e, "UnifracUnweightedCommand", "execute");
431 /***********************************************************/
432 void UnifracUnweightedCommand::printUnweightedFile() {
437 tags.push_back("Score");
438 tags.push_back("RandFreq"); tags.push_back("RandCumul");
440 for(int a = 0; a < numComp; a++) {
441 output->initFile(groupComb[a], tags);
443 for (map<float,float>::iterator it = validScores.begin(); it != validScores.end(); it++) {
444 data.push_back(it->first); data.push_back(rscoreFreq[a][it->first]); data.push_back(rCumul[a][it->first]);
445 output->output(data);
451 catch(exception& e) {
452 m->errorOut(e, "UnifracUnweightedCommand", "printUnweightedFile");
457 /***********************************************************/
458 void UnifracUnweightedCommand::printUWSummaryFile(int i) {
462 outSum.setf(ios::fixed, ios::floatfield); outSum.setf(ios::showpoint);
466 for(int a = 0; a < numComp; a++) {
467 outSum << i+1 << '\t';
468 m->mothurOut(toString(i+1) + "\t");
471 if (UWScoreSig[a][0] > (1/(float)iters)) {
472 outSum << setprecision(6) << groupComb[a] << '\t' << utreeScores[a][0] << '\t' << setprecision(itersString.length()) << UWScoreSig[a][0] << endl;
473 cout << setprecision(6) << groupComb[a] << '\t' << utreeScores[a][0] << '\t' << setprecision(itersString.length()) << UWScoreSig[a][0] << endl;
474 m->mothurOutJustToLog(groupComb[a] + "\t" + toString(utreeScores[a][0]) + "\t" + toString(UWScoreSig[a][0])+ "\n");
476 outSum << setprecision(6) << groupComb[a] << '\t' << utreeScores[a][0] << '\t' << setprecision(itersString.length()) << "<" << (1/float(iters)) << endl;
477 cout << setprecision(6) << groupComb[a] << '\t' << utreeScores[a][0] << '\t' << setprecision(itersString.length()) << "<" << (1/float(iters)) << endl;
478 m->mothurOutJustToLog(groupComb[a] + "\t" + toString(utreeScores[a][0]) + "\t<" + toString((1/float(iters))) + "\n");
481 outSum << setprecision(6) << groupComb[a] << '\t' << utreeScores[a][0] << '\t' << "0.00" << endl;
482 cout << setprecision(6) << groupComb[a] << '\t' << utreeScores[a][0] << '\t' << "0.00" << endl;
483 m->mothurOutJustToLog(groupComb[a] + "\t" + toString(utreeScores[a][0]) + "\t0.00\n");
488 catch(exception& e) {
489 m->errorOut(e, "UnifracUnweightedCommand", "printUWSummaryFile");
493 /***********************************************************/
494 void UnifracUnweightedCommand::createPhylipFile(int i) {
496 string phylipFileName;
497 if ((outputForm == "lt") || (outputForm == "square")) {
498 phylipFileName = outputDir + m->getSimpleName(treefile) + toString(i+1) + ".unweighted.phylip.dist";
499 outputNames.push_back(phylipFileName); outputTypes["phylip"].push_back(phylipFileName);
501 phylipFileName = outputDir + m->getSimpleName(treefile) + toString(i+1) + ".unweighted.column.dist";
502 outputNames.push_back(phylipFileName); outputTypes["column"].push_back(phylipFileName);
506 m->openOutputFile(phylipFileName, out);
508 if ((outputForm == "lt") || (outputForm == "square")) {
510 out << m->Groups.size() << endl;
513 //make matrix with scores in it
514 vector< vector<float> > dists; dists.resize(m->Groups.size());
515 for (int i = 0; i < m->Groups.size(); i++) {
516 dists[i].resize(m->Groups.size(), 0.0);
519 //flip it so you can print it
521 for (int r=0; r<m->Groups.size(); r++) {
522 for (int l = 0; l < r; l++) {
523 dists[r][l] = utreeScores[count][0];
524 dists[l][r] = utreeScores[count][0];
530 for (int r=0; r<m->Groups.size(); r++) {
532 string name = m->Groups[r];
533 if (name.length() < 10) { //pad with spaces to make compatible
534 while (name.length() < 10) { name += " "; }
537 if (outputForm == "lt") {
541 for (int l = 0; l < r; l++) { out << dists[r][l] << '\t'; }
543 }else if (outputForm == "square") {
547 for (int l = 0; l < m->Groups.size(); l++) { out << dists[r][l] << '\t'; }
551 for (int l = 0; l < r; l++) {
552 string otherName = m->Groups[l];
553 if (otherName.length() < 10) { //pad with spaces to make compatible
554 while (otherName.length() < 10) { otherName += " "; }
557 out << name << '\t' << otherName << dists[r][l] << endl;
563 catch(exception& e) {
564 m->errorOut(e, "UnifracUnweightedCommand", "createPhylipFile");
567 }/*****************************************************************/
568 int UnifracUnweightedCommand::readNamesFile() {
571 numUniquesInName = 0;
574 m->openInputFile(namefile, in);
576 string first, second;
577 map<string, string>::iterator itNames;
580 in >> first >> second; m->gobble(in);
584 itNames = m->names.find(first);
585 if (itNames == m->names.end()) {
586 m->names[first] = second;
588 //we need a list of names in your namefile to use above when removing extra seqs above so we don't remove them
589 vector<string> dupNames;
590 m->splitAtComma(second, dupNames);
592 for (int i = 0; i < dupNames.size(); i++) {
593 nameMap[dupNames[i]] = dupNames[i];
594 if ((groupfile == "") && (i != 0)) { tmap->addSeq(dupNames[i], "Group1"); }
596 }else { m->mothurOut(first + " has already been seen in namefile, disregarding names file."); m->mothurOutEndLine(); in.close(); m->names.clear(); namefile = ""; return 1; }
602 catch(exception& e) {
603 m->errorOut(e, "UnifracUnweightedCommand", "readNamesFile");
607 /***********************************************************/