#include "unifracweightedcommand.h"
//**********************************************************************************************************************
-vector<string> UnifracWeightedCommand::getValidParameters(){
+vector<string> UnifracWeightedCommand::setParameters(){
try {
- string Array[] = {"groups","iters","distance","random","processors","root","outputdir","inputdir"};
- vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
+ CommandParameter ptree("tree", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(ptree);
+ CommandParameter pgroup("group", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(pgroup);
+ CommandParameter pname("name", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(pname);
+ CommandParameter pgroups("groups", "String", "", "", "", "", "",false,false); parameters.push_back(pgroups);
+ CommandParameter piters("iters", "Number", "", "1000", "", "", "",false,false); parameters.push_back(piters);
+ CommandParameter pprocessors("processors", "Number", "", "1", "", "", "",false,false); parameters.push_back(pprocessors);
+ CommandParameter prandom("random", "Boolean", "", "F", "", "", "",false,false); parameters.push_back(prandom);
+ CommandParameter pdistance("distance", "Multiple", "column-lt-square", "column", "", "", "",false,false); parameters.push_back(pdistance);
+ CommandParameter proot("root", "Boolean", "F", "", "", "", "",false,false); parameters.push_back(proot);
+ CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
+ CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
+
+ vector<string> myArray;
+ for (int i = 0; i < parameters.size(); i++) { myArray.push_back(parameters[i].name); }
return myArray;
}
catch(exception& e) {
- m->errorOut(e, "UnifracWeightedCommand", "getValidParameters");
+ m->errorOut(e, "UnifracWeightedCommand", "setParameters");
+ exit(1);
+ }
+}
+//**********************************************************************************************************************
+string UnifracWeightedCommand::getHelpString(){
+ try {
+ string helpString = "";
+ helpString += "The unifrac.weighted 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";
+ helpString += "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";
+ 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";
+ helpString += "The distance parameter allows you to create a distance file from the results. The default is false.\n";
+ helpString += "The random parameter allows you to shut off the comparison to random trees. The default is false, meaning don't compare your trees with randomly generated trees.\n";
+ 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";
+ helpString += "The processors parameter allows you to specify the number of processors to use. The default is 1.\n";
+ helpString += "The unifrac.weighted command should be in the following format: unifrac.weighted(groups=yourGroups, iters=yourIters).\n";
+ helpString += "Example unifrac.weighted(groups=A-B-C, iters=500).\n";
+ helpString += "The default value for groups is all the groups in your groupfile, and iters is 1000.\n";
+ helpString += "The unifrac.weighted command output two files: .weighted and .wsummary their descriptions are in the manual.\n";
+ helpString += "Note: No spaces between parameter labels (i.e. groups), '=' and parameters (i.e.yourGroups).\n\n";
+ return helpString;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "UnifracWeightedCommand", "getHelpString");
exit(1);
}
}
UnifracWeightedCommand::UnifracWeightedCommand(){
try {
abort = true; calledHelp = true;
+ setParameters();
vector<string> tempOutNames;
outputTypes["weighted"] = tempOutNames;
outputTypes["wsummary"] = tempOutNames;
exit(1);
}
}
-//**********************************************************************************************************************
-vector<string> UnifracWeightedCommand::getRequiredParameters(){
- try {
- vector<string> myArray;
- return myArray;
- }
- catch(exception& e) {
- m->errorOut(e, "UnifracWeightedCommand", "getRequiredParameters");
- exit(1);
- }
-}
-//**********************************************************************************************************************
-vector<string> UnifracWeightedCommand::getRequiredFiles(){
- try {
- string Array[] = {"tree","group"};
- vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
- return myArray;
- }
- catch(exception& e) {
- m->errorOut(e, "UnifracWeightedCommand", "getRequiredFiles");
- exit(1);
- }
-}
/***********************************************************/
UnifracWeightedCommand::UnifracWeightedCommand(string option) {
try {
- globaldata = GlobalData::getInstance();
abort = false; calledHelp = false;
- Groups.clear();
//allow user to run help
if(option == "help") { help(); abort = true; calledHelp = true; }
else {
- //valid paramters for this command
- string Array[] = {"groups","iters","distance","random","processors","root","outputdir","inputdir"};
- vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
+ vector<string> myArray = setParameters();
OptionParser parser(option);
map<string,string> parameters=parser.getParameters();
+ map<string,string>::iterator it;
ValidParameters validParameter;
outputTypes["phylip"] = tempOutNames;
outputTypes["column"] = tempOutNames;
- if (globaldata->gTree.size() == 0) {//no trees were read
- m->mothurOut("You must execute the read.tree command, before you may execute the unifrac.weighted command."); m->mothurOutEndLine(); abort = true; }
-
- //if the user changes the output directory command factory will send this info to us in the output parameter
- outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){
- outputDir = "";
- outputDir += m->hasPath(globaldata->inputFileName); //if user entered a file with a path then preserve it
+ //if the user changes the input directory command factory will send this info to us in the output parameter
+ string inputDir = validParameter.validFile(parameters, "inputdir", false);
+ if (inputDir == "not found"){ inputDir = ""; }
+ else {
+ string path;
+ it = parameters.find("tree");
+ //user has given a template file
+ if(it != parameters.end()){
+ path = m->hasPath(it->second);
+ //if the user has not given a path then, add inputdir. else leave path alone.
+ if (path == "") { parameters["tree"] = inputDir + it->second; }
+ }
+
+ it = parameters.find("group");
+ //user has given a template file
+ if(it != parameters.end()){
+ path = m->hasPath(it->second);
+ //if the user has not given a path then, add inputdir. else leave path alone.
+ if (path == "") { parameters["group"] = inputDir + it->second; }
+ }
+
+ it = parameters.find("name");
+ //user has given a template file
+ if(it != parameters.end()){
+ path = m->hasPath(it->second);
+ //if the user has not given a path then, add inputdir. else leave path alone.
+ if (path == "") { parameters["name"] = inputDir + it->second; }
+ }
}
+
+ m->runParse = true;
+
+ //check for required parameters
+ treefile = validParameter.validFile(parameters, "tree", true);
+ if (treefile == "not open") { abort = true; }
+ else if (treefile == "not found") { //if there is a current design file, use it
+ treefile = m->getTreeFile();
+ if (treefile != "") { m->mothurOut("Using " + treefile + " as input file for the tree parameter."); m->mothurOutEndLine(); }
+ else { m->mothurOut("You have no current tree file and the tree parameter is required."); m->mothurOutEndLine(); abort = true; }
+ }
+
+ //check for required parameters
+ groupfile = validParameter.validFile(parameters, "group", true);
+ if (groupfile == "not open") { abort = true; }
+ else if (groupfile == "not found") { groupfile = ""; }
+
+ namefile = validParameter.validFile(parameters, "name", true);
+ if (namefile == "not open") { abort = true; }
+ else if (namefile == "not found") { namefile = ""; }
+
+ outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){ outputDir = ""; }
+
//check for optional parameter and set defaults
// ...at some point should added some additional type checking...
if (groups == "not found") { groups = ""; }
else {
m->splitAtDash(groups, Groups);
- globaldata->Groups = Groups;
+ m->Groups = Groups;
}
itersString = validParameter.validFile(parameters, "iters", false); if (itersString == "not found") { itersString = "1000"; }
temp = validParameter.validFile(parameters, "root", false); if (temp == "not found") { temp = "F"; }
includeRoot = m->isTrue(temp);
- temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = "1"; }
- convert(temp, processors);
+ temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); }
+ m->setProcessors(temp);
+ convert(temp, processors);
if (!random) { iters = 0; } //turn off random calcs
-
-
- if (abort == false) {
- T = globaldata->gTree;
- tmap = globaldata->gTreemap;
- sumFile = outputDir + m->getSimpleName(globaldata->getTreeFile()) + ".wsummary";
- m->openOutputFile(sumFile, outSum);
- outputNames.push_back(sumFile); outputTypes["wsummary"].push_back(sumFile);
-
- util = new SharedUtil();
- string s; //to make work with setgroups
- util->setGroups(globaldata->Groups, tmap->namesOfGroups, s, numGroups, "weighted"); //sets the groups the user wants to analyze
- util->getCombos(groupComb, globaldata->Groups, numComp);
-
- weighted = new Weighted(tmap, includeRoot);
-
- }
}
exit(1);
}
}
-//**********************************************************************************************************************
-
-void UnifracWeightedCommand::help(){
- try {
- m->mothurOut("The unifrac.weighted command can only be executed after a successful read.tree command.\n");
- m->mothurOut("The unifrac.weighted command parameters are groups, iters, distance, processors, root and random. No parameters are required.\n");
- 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 2 valid groups.\n");
- 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");
- m->mothurOut("The distance parameter allows you to create a distance file from the results. The default is false.\n");
- m->mothurOut("The random parameter allows you to shut off the comparison to random trees. The default is false, meaning don't compare your trees with randomly generated trees.\n");
- m->mothurOut("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");
- m->mothurOut("The processors parameter allows you to specify the number of processors to use. The default is 1.\n");
- m->mothurOut("The unifrac.weighted command should be in the following format: unifrac.weighted(groups=yourGroups, iters=yourIters).\n");
- m->mothurOut("Example unifrac.weighted(groups=A-B-C, iters=500).\n");
- m->mothurOut("The default value for groups is all the groups in your groupfile, and iters is 1000.\n");
- m->mothurOut("The unifrac.weighted command output two files: .weighted and .wsummary their descriptions are in the manual.\n");
- m->mothurOut("Note: No spaces between parameter labels (i.e. groups), '=' and parameters (i.e.yourGroups).\n\n");
- }
- catch(exception& e) {
- m->errorOut(e, "UnifracWeightedCommand", "help");
- exit(1);
- }
-}
-
/***********************************************************/
int UnifracWeightedCommand::execute() {
try {
if (abort == true) { if (calledHelp) { return 0; } return 2; }
+ if (groupfile != "") {
+ //read in group map info.
+ tmap = new TreeMap(groupfile);
+ tmap->readMap();
+ }else{ //fake out by putting everyone in one group
+ Tree* tree = new Tree(treefile); delete tree; //extracts names from tree to make faked out groupmap
+ tmap = new TreeMap();
+
+ for (int i = 0; i < m->Treenames.size(); i++) { tmap->addSeq(m->Treenames[i], "Group1"); }
+ }
+
+ if (namefile != "") { readNamesFile(); }
+
+ read = new ReadNewickTree(treefile);
+ int readOk = read->read(tmap);
+
+ if (readOk != 0) { m->mothurOut("Read Terminated."); m->mothurOutEndLine(); delete tmap; delete read; return 0; }
+
+ read->AssembleTrees();
+ T = read->getTrees();
+ delete read;
+
+ //make sure all files match
+ //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.
+ int numNamesInTree;
+ if (namefile != "") {
+ if (numUniquesInName == m->Treenames.size()) { numNamesInTree = nameMap.size(); }
+ else { numNamesInTree = m->Treenames.size(); }
+ }else { numNamesInTree = m->Treenames.size(); }
+
+
+ //output any names that are in group file but not in tree
+ if (numNamesInTree < tmap->getNumSeqs()) {
+ for (int i = 0; i < tmap->namesOfSeqs.size(); i++) {
+ //is that name in the tree?
+ int count = 0;
+ for (int j = 0; j < m->Treenames.size(); j++) {
+ if (tmap->namesOfSeqs[i] == m->Treenames[j]) { break; } //found it
+ count++;
+ }
+
+ if (m->control_pressed) {
+ delete tmap; for (int i = 0; i < T.size(); i++) { delete T[i]; }
+ for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } outputTypes.clear();
+ m->Groups.clear();
+ return 0;
+ }
+
+ //then you did not find it so report it
+ if (count == m->Treenames.size()) {
+ //if it is in your namefile then don't remove
+ map<string, string>::iterator it = nameMap.find(tmap->namesOfSeqs[i]);
+
+ if (it == nameMap.end()) {
+ m->mothurOut(tmap->namesOfSeqs[i] + " is in your groupfile and not in your tree. It will be disregarded."); m->mothurOutEndLine();
+ tmap->removeSeq(tmap->namesOfSeqs[i]);
+ i--; //need this because removeSeq removes name from namesOfSeqs
+ }
+ }
+ }
+ }
+
+ sumFile = outputDir + m->getSimpleName(treefile) + ".wsummary";
+ m->openOutputFile(sumFile, outSum);
+ outputNames.push_back(sumFile); outputTypes["wsummary"].push_back(sumFile);
+
+ util = new SharedUtil();
+ string s; //to make work with setgroups
+ util->setGroups(m->Groups, tmap->namesOfGroups, s, numGroups, "weighted"); //sets the groups the user wants to analyze
+ util->getCombos(groupComb, m->Groups, numComp);
+ delete util;
+
+ weighted = new Weighted(tmap, includeRoot);
+
int start = time(NULL);
//get weighted for users tree
//get weighted scores for users trees
for (int i = 0; i < T.size(); i++) {
- if (m->control_pressed) { outSum.close(); for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
+ if (m->control_pressed) { delete tmap; delete weighted;
+ for (int i = 0; i < T.size(); i++) { delete T[i]; } outSum.close(); for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
counter = 0;
rScores.resize(numComp); //data[0] = weightedscore AB, data[1] = weightedscore AC...
uScores.resize(numComp); //data[0] = weightedscore AB, data[1] = weightedscore AC...
if (random) {
- output = new ColumnFile(outputDir + m->getSimpleName(globaldata->getTreeFile()) + toString(i+1) + ".weighted", itersString);
- outputNames.push_back(outputDir + m->getSimpleName(globaldata->getTreeFile()) + toString(i+1) + ".weighted");
- outputTypes["weighted"].push_back(outputDir + m->getSimpleName(globaldata->getTreeFile()) + toString(i+1) + ".weighted");
+ output = new ColumnFile(outputDir + m->getSimpleName(treefile) + toString(i+1) + ".weighted", itersString);
+ outputNames.push_back(outputDir + m->getSimpleName(treefile) + toString(i+1) + ".weighted");
+ outputTypes["weighted"].push_back(outputDir + m->getSimpleName(treefile) + toString(i+1) + ".weighted");
}
userData = weighted->getValues(T[i], processors, outputDir); //userData[0] = weightedscore
- if (m->control_pressed) { if (random) { delete output; } outSum.close(); for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
+ if (m->control_pressed) { delete tmap; delete weighted;
+ 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; }
//save users score
for (int s=0; s<numComp; s++) {
vector< vector<string> > namesOfGroupCombos;
for (int a=0; a<numGroups; a++) {
for (int l = 0; l < a; l++) {
- vector<string> groups; groups.push_back(globaldata->Groups[a]); groups.push_back(globaldata->Groups[l]);
+ vector<string> groups; groups.push_back(m->Groups[a]); groups.push_back(m->Groups[l]);
namesOfGroupCombos.push_back(groups);
}
}
driver(T[i], namesOfGroupCombos, 0, namesOfGroupCombos.size(), rScores);
#endif
- if (m->control_pressed) { delete output; outSum.close(); for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
+ if (m->control_pressed) { delete tmap; delete weighted;
+ for (int i = 0; i < T.size(); i++) { delete T[i]; } delete output; outSum.close(); for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
//report progress
// m->mothurOut("Iter: " + toString(j+1)); m->mothurOutEndLine();
}
- if (m->control_pressed) { outSum.close(); for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
+ if (m->control_pressed) { delete tmap; delete weighted;
+ for (int i = 0; i < T.size(); i++) { delete T[i]; } outSum.close(); for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
printWSummaryFile();
if (phylip) { createPhylipFile(); }
//clear out users groups
- globaldata->Groups.clear();
+ m->Groups.clear();
+ delete tmap; delete weighted;
+ for (int i = 0; i < T.size(); i++) { delete T[i]; }
if (m->control_pressed) {
/**************************************************************************************************/
int UnifracWeightedCommand::driver(Tree* t, vector< vector<string> > namesOfGroupCombos, int start, int num, vector< vector<double> >& scores) {
try {
- Tree* randT = new Tree();
+ Tree* randT = new Tree(tmap);
for (int h = start; h < (start+num); h++) {
string phylipFileName;
if ((outputForm == "lt") || (outputForm == "square")) {
- phylipFileName = outputDir + m->getSimpleName(globaldata->getTreeFile()) + toString(i+1) + ".weighted.phylip.dist";
+ phylipFileName = outputDir + m->getSimpleName(treefile) + toString(i+1) + ".weighted.phylip.dist";
outputNames.push_back(phylipFileName); outputTypes["phylip"].push_back(phylipFileName);
}else { //column
- phylipFileName = outputDir + m->getSimpleName(globaldata->getTreeFile()) + toString(i+1) + ".weighted.column.dist";
+ phylipFileName = outputDir + m->getSimpleName(treefile) + toString(i+1) + ".weighted.column.dist";
outputNames.push_back(phylipFileName); outputTypes["column"].push_back(phylipFileName);
}
if ((outputForm == "lt") || (outputForm == "square")) {
//output numSeqs
- out << globaldata->Groups.size() << endl;
+ out << m->Groups.size() << endl;
}
//make matrix with scores in it
- vector< vector<float> > dists; dists.resize(globaldata->Groups.size());
- for (int i = 0; i < globaldata->Groups.size(); i++) {
- dists[i].resize(globaldata->Groups.size(), 0.0);
+ vector< vector<float> > dists; dists.resize(m->Groups.size());
+ for (int i = 0; i < m->Groups.size(); i++) {
+ dists[i].resize(m->Groups.size(), 0.0);
}
//flip it so you can print it
- for (int r=0; r<globaldata->Groups.size(); r++) {
+ for (int r=0; r<m->Groups.size(); r++) {
for (int l = 0; l < r; l++) {
dists[r][l] = utreeScores[count];
dists[l][r] = utreeScores[count];
}
//output to file
- for (int r=0; r<globaldata->Groups.size(); r++) {
+ for (int r=0; r<m->Groups.size(); r++) {
//output name
- string name = globaldata->Groups[r];
+ string name = m->Groups[r];
if (name.length() < 10) { //pad with spaces to make compatible
while (name.length() < 10) { name += " "; }
}
out << name << '\t';
//output distances
- for (int l = 0; l < globaldata->Groups.size(); l++) { out << dists[r][l] << '\t'; }
+ for (int l = 0; l < m->Groups.size(); l++) { out << dists[r][l] << '\t'; }
out << endl;
}else{
//output distances
for (int l = 0; l < r; l++) {
- string otherName = globaldata->Groups[l];
+ string otherName = m->Groups[l];
if (otherName.length() < 10) { //pad with spaces to make compatible
while (otherName.length() < 10) { otherName += " "; }
}
exit(1);
}
}
-
+/*****************************************************************/
+int UnifracWeightedCommand::readNamesFile() {
+ try {
+ m->names.clear();
+ numUniquesInName = 0;
+
+ ifstream in;
+ m->openInputFile(namefile, in);
+
+ string first, second;
+ map<string, string>::iterator itNames;
+
+ while(!in.eof()) {
+ in >> first >> second; m->gobble(in);
+
+ numUniquesInName++;
+
+ itNames = m->names.find(first);
+ if (itNames == m->names.end()) {
+ m->names[first] = second;
+
+ //we need a list of names in your namefile to use above when removing extra seqs above so we don't remove them
+ vector<string> dupNames;
+ m->splitAtComma(second, dupNames);
+
+ for (int i = 0; i < dupNames.size(); i++) {
+ nameMap[dupNames[i]] = dupNames[i];
+ if ((groupfile == "") && (i != 0)) { tmap->addSeq(dupNames[i], "Group1"); }
+ }
+ }else { m->mothurOut(first + " has already been seen in namefile, disregarding names file."); m->mothurOutEndLine(); in.close(); m->names.clear(); namefile = ""; return 1; }
+ }
+ in.close();
+
+ return 0;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "UnifracWeightedCommand", "readNamesFile");
+ exit(1);
+ }
+}
/***********************************************************/
projects / mothur.git / blobdiff