*/
#include "unifracweightedcommand.h"
+#include "consensus.h"
+#include "subsample.h"
+#include "treereader.h"
//**********************************************************************************************************************
-vector<string> UnifracWeightedCommand::getValidParameters(){
+vector<string> UnifracWeightedCommand::setParameters(){
try {
- string Array[] = {"groups","iters","distance","random","processors","outputdir","inputdir"};
- vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
+ CommandParameter ptree("tree", "InputTypes", "", "", "none", "none", "none","weighted-wsummary",false,true,true); parameters.push_back(ptree);
+ CommandParameter pname("name", "InputTypes", "", "", "NameCount", "none", "none","",false,false,true); parameters.push_back(pname);
+ CommandParameter pcount("count", "InputTypes", "", "", "NameCount-CountGroup", "none", "none","",false,false,true); parameters.push_back(pcount);
+ CommandParameter pgroup("group", "InputTypes", "", "", "CountGroup", "none", "none","",false,false,true); parameters.push_back(pgroup);
+ 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,true); parameters.push_back(pprocessors);
+ CommandParameter psubsample("subsample", "String", "", "", "", "", "","",false,false); parameters.push_back(psubsample);
+ CommandParameter pconsensus("consensus", "Boolean", "", "F", "", "", "","tree",false,false); parameters.push_back(pconsensus);
+ CommandParameter prandom("random", "Boolean", "", "F", "", "", "","",false,false); parameters.push_back(prandom);
+ CommandParameter pdistance("distance", "Multiple", "column-lt-square-phylip", "column", "", "", "","phylip-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);
}
}
//**********************************************************************************************************************
-UnifracWeightedCommand::UnifracWeightedCommand(){
+string UnifracWeightedCommand::getHelpString(){
try {
- //initialize outputTypes
- vector<string> tempOutNames;
- outputTypes["weighted"] = tempOutNames;
- outputTypes["wsummary"] = tempOutNames;
- outputTypes["phylip"] = tempOutNames;
+ string helpString = "";
+ helpString += "The unifrac.weighted command parameters are tree, group, name, count, groups, iters, distance, processors, root, subsample, consensus 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 subsample parameter allows you to enter the size pergroup of the sample or you can set subsample=T and mothur will use the size of your smallest group. The subsample parameter may only be used with a group file.\n";
+ helpString += "The consensus parameter allows you to indicate you would like trees built from distance matrices created with the results, as well as a consensus tree built from these trees. Default=F.\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";
+ return helpString;
}
catch(exception& e) {
- m->errorOut(e, "UnifracWeightedCommand", "UnifracWeightedCommand");
+ m->errorOut(e, "UnifracWeightedCommand", "getHelpString");
exit(1);
}
}
//**********************************************************************************************************************
-vector<string> UnifracWeightedCommand::getRequiredParameters(){
- try {
- vector<string> myArray;
- return myArray;
- }
- catch(exception& e) {
- m->errorOut(e, "UnifracWeightedCommand", "getRequiredParameters");
- exit(1);
- }
+string UnifracWeightedCommand::getOutputPattern(string type) {
+ try {
+ string pattern = "";
+ if (type == "weighted") { pattern = "[filename],weighted-[filename],[tag],weighted"; }
+ else if (type == "wsummary") { pattern = "[filename],wsummary"; }
+ else if (type == "phylip") { pattern = "[filename],[tag],[tag2],dist"; }
+ else if (type == "column") { pattern = "[filename],[tag],[tag2],dist"; }
+ else if (type == "tree") { pattern = "[filename],[tag],[tag2],tre"; }
+ else { m->mothurOut("[ERROR]: No definition for type " + type + " output pattern.\n"); m->control_pressed = true; }
+
+ return pattern;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "UnifracWeightedCommand", "getOutputPattern");
+ exit(1);
+ }
}
//**********************************************************************************************************************
-vector<string> UnifracWeightedCommand::getRequiredFiles(){
+UnifracWeightedCommand::UnifracWeightedCommand(){
try {
- string Array[] = {"tree","group"};
- vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
-
- return myArray;
+ abort = true; calledHelp = true;
+ setParameters();
+ vector<string> tempOutNames;
+ outputTypes["weighted"] = tempOutNames;
+ outputTypes["wsummary"] = tempOutNames;
+ outputTypes["phylip"] = tempOutNames;
+ outputTypes["column"] = tempOutNames;
+ outputTypes["tree"] = tempOutNames;
}
catch(exception& e) {
- m->errorOut(e, "UnifracWeightedCommand", "getRequiredFiles");
+ m->errorOut(e, "UnifracWeightedCommand", "UnifracWeightedCommand");
exit(1);
}
}
+
/***********************************************************/
UnifracWeightedCommand::UnifracWeightedCommand(string option) {
try {
- globaldata = GlobalData::getInstance();
- abort = false;
- Groups.clear();
+ abort = false; calledHelp = false;
//allow user to run help
- if(option == "help") { help(); abort = true; }
+ if(option == "help") { help(); abort = true; calledHelp = true; }
+ else if(option == "citation") { citation(); abort = true; calledHelp = true;}
else {
- //valid paramters for this command
- string Array[] = {"groups","iters","distance","random","processors","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["weighted"] = tempOutNames;
outputTypes["wsummary"] = tempOutNames;
outputTypes["phylip"] = tempOutNames;
+ outputTypes["column"] = tempOutNames;
+ outputTypes["tree"] = 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; }
+ }
+
+ it = parameters.find("count");
+ //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["count"] = inputDir + it->second; }
+ }
}
+
+ //check for required parameters
+ treefile = validParameter.validFile(parameters, "tree", true);
+ if (treefile == "not open") { treefile = ""; 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; }
+ }else { m->setTreeFile(treefile); }
+
+ //check for required parameters
+ groupfile = validParameter.validFile(parameters, "group", true);
+ if (groupfile == "not open") { abort = true; }
+ else if (groupfile == "not found") { groupfile = ""; }
+ else { m->setGroupFile(groupfile); }
+
+ namefile = validParameter.validFile(parameters, "name", true);
+ if (namefile == "not open") { namefile = ""; abort = true; }
+ else if (namefile == "not found") { namefile = ""; }
+ else { m->setNameFile(namefile); }
+
+ countfile = validParameter.validFile(parameters, "count", true);
+ if (countfile == "not open") { countfile = ""; abort = true; }
+ else if (countfile == "not found") { countfile = ""; }
+ else { m->setCountTableFile(countfile); }
+
+ if ((namefile != "") && (countfile != "")) {
+ m->mothurOut("[ERROR]: you may only use one of the following: name or count."); m->mothurOutEndLine(); abort = true;
+ }
+
+ if ((groupfile != "") && (countfile != "")) {
+ m->mothurOut("[ERROR]: you may only use one of the following: group or count."); m->mothurOutEndLine(); abort=true;
+ }
+
+ outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){ outputDir = m->hasPath(treefile); }
+
//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->setGroups(Groups);
}
itersString = validParameter.validFile(parameters, "iters", false); if (itersString == "not found") { itersString = "1000"; }
- convert(itersString, iters);
+ m->mothurConvert(itersString, iters);
- string temp = validParameter.validFile(parameters, "distance", false); if (temp == "not found") { temp = "false"; }
- phylip = m->isTrue(temp);
-
- temp = validParameter.validFile(parameters, "random", false); if (temp == "not found") { temp = "F"; }
+ string temp = validParameter.validFile(parameters, "distance", false);
+ if (temp == "not found") { phylip = false; outputForm = ""; }
+ else{
+ if (temp=="phylip") { temp = "lt"; }
+ if ((temp == "lt") || (temp == "column") || (temp == "square")) { phylip = true; outputForm = temp; }
+ else { m->mothurOut("Options for distance are: lt, square, or column. Using lt."); m->mothurOutEndLine(); phylip = true; outputForm = "lt"; }
+ }
+
+ temp = validParameter.validFile(parameters, "random", false); if (temp == "not found") { temp = "F"; }
random = m->isTrue(temp);
- temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = "1"; }
- convert(temp, processors);
+ temp = validParameter.validFile(parameters, "root", false); if (temp == "not found") { temp = "F"; }
+ includeRoot = m->isTrue(temp);
- if (!random) { iters = 0; } //turn off random calcs
-
+ temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); }
+ m->setProcessors(temp);
+ m->mothurConvert(temp, processors);
+
+ temp = validParameter.validFile(parameters, "subsample", false); if (temp == "not found") { temp = "F"; }
+ if (m->isNumeric1(temp)) { m->mothurConvert(temp, subsampleSize); subsample = true; }
+ else {
+ if (m->isTrue(temp)) { subsample = true; subsampleSize = -1; } //we will set it to smallest group later
+ else { subsample = false; }
+ }
- 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);
-
- }
+ if (!subsample) { subsampleIters = 0; }
+ else { subsampleIters = iters; }
+
+ temp = validParameter.validFile(parameters, "consensus", false); if (temp == "not found") { temp = "F"; }
+ consensus = m->isTrue(temp);
+
+ if (subsample && random) { m->mothurOut("[ERROR]: random must be false, if subsample=t.\n"); abort=true; }
+ if (countfile == "") { if (subsample && (groupfile == "")) { m->mothurOut("[ERROR]: if subsample=t, a group file must be provided.\n"); abort=true; } }
+ else {
+ CountTable testCt;
+ if ((!testCt.testGroups(countfile)) && (subsample)) {
+ m->mothurOut("[ERROR]: if subsample=t, a count file with group info must be provided.\n"); abort=true;
+ }
+ }
+ if (subsample && (!phylip)) { phylip=true; outputForm = "lt"; }
+ if (consensus && (!subsample)) { m->mothurOut("[ERROR]: you cannot use consensus without subsample.\n"); abort=true; }
+
+ if (countfile=="") {
+ if (namefile == "") {
+ vector<string> files; files.push_back(treefile);
+ parser.getNameFile(files);
+ }
+ }
}
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 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 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) { return 0; }
+ if (abort == true) { if (calledHelp) { return 0; } return 2; }
- int start = time(NULL);
+ m->setTreeFile(treefile);
- //get weighted for users tree
- userData.resize(numComp,0); //data[0] = weightedscore AB, data[1] = weightedscore AC...
- randomData.resize(numComp,0); //data[0] = weightedscore AB, data[1] = weightedscore AC...
-
- //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; }
-
- 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");
- }
-
- 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; }
-
- //save users score
- for (int s=0; s<numComp; s++) {
- //add users score to vector of user scores
- uScores[s].push_back(userData[s]);
-
- //save users tree score for summary file
- utreeScores.push_back(userData[s]);
- }
-
- if (random) {
- vector<double> sums = weighted->getBranchLengthSums(T[i]);
-
- //calculate number of comparisons i.e. with groups A,B,C = AB, AC, BC = 3;
- 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]);
- namesOfGroupCombos.push_back(groups);
- }
- }
-
- lines.clear();
-
- #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
- if(processors != 1){
- int numPairs = namesOfGroupCombos.size();
- int numPairsPerProcessor = numPairs / processors;
-
- for (int i = 0; i < processors; i++) {
- int startPos = i * numPairsPerProcessor;
- if(i == processors - 1){
- numPairsPerProcessor = numPairs - i * numPairsPerProcessor;
- }
- lines.push_back(linePair(startPos, numPairsPerProcessor));
- }
- }
- #endif
-
-
- //get scores for random trees
- for (int j = 0; j < iters; j++) {
-
- #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
- if(processors == 1){
- driver(T[i], namesOfGroupCombos, 0, namesOfGroupCombos.size(), sums, rScores);
- }else{
- createProcesses(T[i], namesOfGroupCombos, sums, rScores);
- }
- #else
- driver(T[i], namesOfGroupCombos, 0, namesOfGroupCombos.size(), sums, 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; }
-
- //report progress
- m->mothurOut("Iter: " + toString(j+1)); m->mothurOutEndLine();
- }
- lines.clear();
-
- //find the signifigance of the score for summary file
- for (int f = 0; f < numComp; f++) {
- //sort random scores
- sort(rScores[f].begin(), rScores[f].end());
-
- //the index of the score higher than yours is returned
- //so if you have 1000 random trees the index returned is 100
- //then there are 900 trees with a score greater then you.
- //giving you a signifigance of 0.900
- int index = findIndex(userData[f], f); if (index == -1) { m->mothurOut("error in UnifracWeightedCommand"); m->mothurOutEndLine(); exit(1); } //error code
-
- //the signifigance is the number of trees with the users score or higher
- WScoreSig.push_back((iters-index)/(float)iters);
- }
-
- //out << "Tree# " << i << endl;
- calculateFreqsCumuls();
- printWeightedFile();
-
- delete output;
-
- }
-
- //clear data
- rScores.clear();
- uScores.clear();
- validScores.clear();
- }
+ TreeReader* reader;
+ if (countfile == "") { reader = new TreeReader(treefile, groupfile, namefile); }
+ else { reader = new TreeReader(treefile, countfile); }
+ T = reader->getTrees();
+ ct = T[0]->getCountTable();
+ delete reader;
+
+ if (m->control_pressed) { delete ct; for (int i = 0; i < T.size(); i++) { delete T[i]; } return 0; }
+
+ map<string, string> variables;
+ variables["[filename]"] = outputDir + m->getSimpleName(treefile);
+ sumFile = getOutputFileName("wsummary",variables);
+ m->openOutputFile(sumFile, outSum);
+ outputNames.push_back(sumFile); outputTypes["wsummary"].push_back(sumFile);
+ SharedUtil util;
+ string s; //to make work with setgroups
+ Groups = m->getGroups();
+ vector<string> nameGroups = ct->getNamesOfGroups();
+ if (nameGroups.size() < 2) { m->mothurOut("[ERROR]: You cannot run unifrac.weighted with less than 2 groups, aborting.\n"); delete ct; for (int i = 0; i < T.size(); i++) { delete T[i]; } return 0; }
+ util.setGroups(Groups, nameGroups, s, numGroups, "weighted"); //sets the groups the user wants to analyze
+ m->setGroups(Groups);
+
+ if (m->control_pressed) { delete ct; for (int i = 0; i < T.size(); i++) { delete T[i]; } return 0; }
+
+ Weighted weighted(includeRoot);
+
+ int start = time(NULL);
+
+ //set or check size
+ if (subsample) {
+ //user has not set size, set size = smallest samples size
+ if (subsampleSize == -1) {
+ vector<string> temp; temp.push_back(Groups[0]);
+ subsampleSize = ct->getGroupCount(Groups[0]); //num in first group
+ for (int i = 1; i < Groups.size(); i++) {
+ int thisSize = ct->getGroupCount(Groups[i]);
+ if (thisSize < subsampleSize) { subsampleSize = thisSize; }
+ }
+ m->mothurOut("\nSetting subsample size to " + toString(subsampleSize) + ".\n\n");
+ }else { //eliminate any too small groups
+ vector<string> newGroups = Groups;
+ Groups.clear();
+ for (int i = 0; i < newGroups.size(); i++) {
+ int thisSize = ct->getGroupCount(newGroups[i]);
+
+ if (thisSize >= subsampleSize) { Groups.push_back(newGroups[i]); }
+ else { m->mothurOut("You have selected a size that is larger than "+newGroups[i]+" number of sequences, removing "+newGroups[i]+".\n"); }
+ }
+ m->setGroups(Groups);
+ }
+ }
+
+ //here in case some groups are removed by subsample
+ util.getCombos(groupComb, Groups, numComp);
+
+ if (numComp < processors) { processors = numComp; }
+
+ if (consensus && (numComp < 2)) { m->mothurOut("consensus can only be used with numComparisions greater than 1, setting consensus=f.\n"); consensus=false; }
+
+ //get weighted scores for users trees
+ for (int i = 0; i < T.size(); i++) {
+
+ if (m->control_pressed) { delete ct; for (int i = 0; i < T.size(); i++) { delete T[i]; } outSum.close(); for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } 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...
+
+ vector<double> userData; userData.resize(numComp,0); //weighted score info for user tree. data[0] = weightedscore AB, data[1] = weightedscore AC...
+ vector<double> randomData; randomData.resize(numComp,0); //weighted score info for random trees. data[0] = weightedscore AB, data[1] = weightedscore AC...
+
+ if (random) {
+ variables["[filename]"] = outputDir + m->getSimpleName(treefile);
+ variables["[tag]"] = toString(i+1);
+ string wFileName = getOutputFileName("weighted", variables);
+ output = new ColumnFile(wFileName, itersString);
+ outputNames.push_back(wFileName); outputTypes["wweighted"].push_back(wFileName);
+ }
+
+ userData = weighted.getValues(T[i], processors, outputDir); //userData[0] = weightedscore
+ if (m->control_pressed) { delete ct; 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++) { m->mothurRemove(outputNames[i]); } return 0; }
+
+ //save users score
+ for (int s=0; s<numComp; s++) {
+ //add users score to vector of user scores
+ uScores[s].push_back(userData[s]);
+ //save users tree score for summary file
+ utreeScores.push_back(userData[s]);
+ }
+
+ if (random) { runRandomCalcs(T[i], userData); }
+
+ //clear data
+ rScores.clear();
+ uScores.clear();
+ validScores.clear();
+
+ //subsample loop
+ vector< vector<double> > calcDistsTotals; //each iter, each groupCombos dists. this will be used to make .dist files
+ for (int thisIter = 0; thisIter < subsampleIters; thisIter++) { //subsampleIters=0, if subsample=f.
+ if (m->control_pressed) { break; }
+
+ //copy to preserve old one - would do this in subsample but memory cleanup becomes messy.
+ CountTable* newCt = new CountTable();
+
+ int sampleTime = 0;
+ if (m->debug) { sampleTime = time(NULL); }
+
+ //uses method of setting groups to doNotIncludeMe
+ SubSample sample;
+ Tree* subSampleTree = sample.getSample(T[i], ct, newCt, subsampleSize);
+
+ if (m->debug) { m->mothurOut("[DEBUG]: iter " + toString(thisIter) + " took " + toString(time(NULL) - sampleTime) + " seconds to sample tree.\n"); }
+
+ //call new weighted function
+ vector<double> iterData; iterData.resize(numComp,0);
+ Weighted thisWeighted(includeRoot);
+ iterData = thisWeighted.getValues(subSampleTree, processors, outputDir); //userData[0] = weightedscore
+
+ //save data to make ave dist, std dist
+ calcDistsTotals.push_back(iterData);
+
+ delete newCt;
+ delete subSampleTree;
+
+ if((thisIter+1) % 100 == 0){ m->mothurOutJustToScreen(toString(thisIter+1)+"\n"); }
+ }
+
+ if (m->control_pressed) { delete ct; 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++) { m->mothurRemove(outputNames[i]); } return 0; }
+
+ if (subsample) { getAverageSTDMatrices(calcDistsTotals, i); }
+ if (consensus) { getConsensusTrees(calcDistsTotals, 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 ct; for (int i = 0; i < T.size(); i++) { delete T[i]; } outSum.close(); for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; }
+ if (phylip) { createPhylipFile(); }
+
printWSummaryFile();
- if (phylip) { createPhylipFile(); }
-
//clear out users groups
- globaldata->Groups.clear();
+ m->clearGroups();
+ delete ct;
+ for (int i = 0; i < T.size(); i++) { delete T[i]; }
+ if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; }
- if (m->control_pressed) {
- for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); }
- return 0;
+ m->mothurOut("It took " + toString(time(NULL) - start) + " secs to run unifrac.weighted."); m->mothurOutEndLine();
+
+ //set phylip file as new current phylipfile
+ string current = "";
+ itTypes = outputTypes.find("phylip");
+ if (itTypes != outputTypes.end()) {
+ if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setPhylipFile(current); }
}
- m->mothurOut("It took " + toString(time(NULL) - start) + " secs to run unifrac.weighted."); m->mothurOutEndLine();
+ //set column file as new current columnfile
+ itTypes = outputTypes.find("column");
+ if (itTypes != outputTypes.end()) {
+ if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setColumnFile(current); }
+ }
m->mothurOutEndLine();
m->mothurOut("Output File Names: "); m->mothurOutEndLine();
}
}
/**************************************************************************************************/
+int UnifracWeightedCommand::getAverageSTDMatrices(vector< vector<double> >& dists, int treeNum) {
+ try {
+ //we need to find the average distance and standard deviation for each groups distance
+
+ //finds sum
+ vector<double> averages = m->getAverages(dists);
+
+ //find standard deviation
+ vector<double> stdDev = m->getStandardDeviation(dists, averages);
+
+ //make matrix with scores in it
+ vector< vector<double> > avedists; //avedists.resize(m->getNumGroups());
+ for (int i = 0; i < m->getNumGroups(); i++) {
+ vector<double> temp;
+ for (int j = 0; j < m->getNumGroups(); j++) { temp.push_back(0.0); }
+ avedists.push_back(temp);
+ }
+
+ //make matrix with scores in it
+ vector< vector<double> > stddists; //stddists.resize(m->getNumGroups());
+ for (int i = 0; i < m->getNumGroups(); i++) {
+ vector<double> temp;
+ for (int j = 0; j < m->getNumGroups(); j++) { temp.push_back(0.0); }
+ //stddists[i].resize(m->getNumGroups(), 0.0);
+ stddists.push_back(temp);
+ }
+
+
+ //flip it so you can print it
+ int count = 0;
+ for (int r=0; r<m->getNumGroups(); r++) {
+ for (int l = 0; l < r; l++) {
+ avedists[r][l] = averages[count];
+ avedists[l][r] = averages[count];
+ stddists[r][l] = stdDev[count];
+ stddists[l][r] = stdDev[count];
+ count++;
+ }
+ }
+
+ map<string, string> variables;
+ variables["[filename]"] = outputDir + m->getSimpleName(treefile);
+ variables["[tag]"] = toString(treeNum+1);
+ variables["[tag2]"] = "weighted.ave";
+ string aveFileName = getOutputFileName("phylip",variables);
+ if (outputForm != "column") { outputNames.push_back(aveFileName); outputTypes["phylip"].push_back(aveFileName); }
+ else { outputNames.push_back(aveFileName); outputTypes["column"].push_back(aveFileName); }
+ ofstream out;
+ m->openOutputFile(aveFileName, out);
+
+ variables["[tag2]"] = "weighted.std";
+ string stdFileName = getOutputFileName("phylip",variables);
+ if (outputForm != "column") { outputNames.push_back(stdFileName); outputTypes["phylip"].push_back(stdFileName); }
+ else { outputNames.push_back(stdFileName); outputTypes["column"].push_back(stdFileName); }
+ ofstream outStd;
+ m->openOutputFile(stdFileName, outStd);
+
+ if ((outputForm == "lt") || (outputForm == "square")) {
+ //output numSeqs
+ out << m->getNumGroups() << endl;
+ outStd << m->getNumGroups() << endl;
+ }
+
+ //output to file
+ for (int r=0; r<m->getNumGroups(); r++) {
+ //output name
+ string name = (m->getGroups())[r];
+ if (name.length() < 10) { //pad with spaces to make compatible
+ while (name.length() < 10) { name += " "; }
+ }
+
+ if (outputForm == "lt") {
+ out << name << '\t';
+ outStd << name << '\t';
+
+ //output distances
+ for (int l = 0; l < r; l++) { out << avedists[r][l] << '\t'; outStd << stddists[r][l] << '\t';}
+ out << endl; outStd << endl;
+ }else if (outputForm == "square") {
+ out << name << '\t';
+ outStd << name << '\t';
+
+ //output distances
+ for (int l = 0; l < m->getNumGroups(); l++) { out << avedists[r][l] << '\t'; outStd << stddists[r][l] << '\t'; }
+ out << endl; outStd << endl;
+ }else{
+ //output distances
+ for (int l = 0; l < r; l++) {
+ string otherName = (m->getGroups())[l];
+ if (otherName.length() < 10) { //pad with spaces to make compatible
+ while (otherName.length() < 10) { otherName += " "; }
+ }
+
+ out << name << '\t' << otherName << avedists[r][l] << endl;
+ outStd << name << '\t' << otherName << stddists[r][l] << endl;
+ }
+ }
+ }
+ out.close();
+ outStd.close();
+
+ return 0;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "UnifracWeightedCommand", "getAverageSTDMatrices");
+ exit(1);
+ }
+}
+
+/**************************************************************************************************/
+int UnifracWeightedCommand::getConsensusTrees(vector< vector<double> >& dists, int treeNum) {
+ try {
+
+ //used in tree constructor
+ m->runParse = false;
+
+ ///create treemap class from groupmap for tree class to use
+ CountTable newCt;
+ set<string> nameMap;
+ map<string, string> groupMap;
+ set<string> gps;
+ for (int i = 0; i < m->getGroups().size(); i++) {
+ nameMap.insert(m->getGroups()[i]);
+ gps.insert(m->getGroups()[i]);
+ groupMap[m->getGroups()[i]] = m->getGroups()[i];
+ }
+ newCt.createTable(nameMap, groupMap, gps);
+
+ //clear old tree names if any
+ m->Treenames.clear();
+
+ //fills globaldatas tree names
+ m->Treenames = m->getGroups();
+
+ vector<Tree*> newTrees = buildTrees(dists, treeNum, newCt); //also creates .all.tre file containing the trees created
+
+ if (m->control_pressed) { return 0; }
+
+ Consensus con;
+ Tree* conTree = con.getTree(newTrees);
+
+ //create a new filename
+ map<string, string> variables;
+ variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(treefile));
+ variables["[tag]"] = toString(treeNum+1);
+ variables["[tag2]"] = "weighted.cons";
+ string conFile = getOutputFileName("tree",variables);
+ outputNames.push_back(conFile); outputTypes["tree"].push_back(conFile);
+ ofstream outTree;
+ m->openOutputFile(conFile, outTree);
+
+ if (conTree != NULL) { conTree->print(outTree, "boot"); delete conTree; }
+ outTree.close();
+
+ return 0;
+
+ }
+ catch(exception& e) {
+ m->errorOut(e, "UnifracWeightedCommand", "getConsensusTrees");
+ exit(1);
+ }
+}
+/**************************************************************************************************/
+
+vector<Tree*> UnifracWeightedCommand::buildTrees(vector< vector<double> >& dists, int treeNum, CountTable& myct) {
+ try {
+
+ vector<Tree*> trees;
+
+ //create a new filename
+ map<string, string> variables;
+ variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(treefile));
+ variables["[tag]"] = toString(treeNum+1);
+ variables["[tag2]"] = "weighted.all";
+ string outputFile = getOutputFileName("tree",variables);
+ outputNames.push_back(outputFile); outputTypes["tree"].push_back(outputFile);
+
+ ofstream outAll;
+ m->openOutputFile(outputFile, outAll);
+
+
+ for (int i = 0; i < dists.size(); i++) { //dists[0] are the dists for the first subsampled tree.
+
+ if (m->control_pressed) { break; }
+
+ //make matrix with scores in it
+ vector< vector<double> > sims; sims.resize(m->getNumGroups());
+ for (int j = 0; j < m->getNumGroups(); j++) {
+ sims[j].resize(m->getNumGroups(), 0.0);
+ }
+
+ int count = 0;
+ for (int r=0; r<m->getNumGroups(); r++) {
+ for (int l = 0; l < r; l++) {
+ double sim = -(dists[i][count]-1.0);
+ sims[r][l] = sim;
+ sims[l][r] = sim;
+ count++;
+ }
+ }
+
+ //create tree
+ Tree* tempTree = new Tree(&myct, sims);
+ tempTree->assembleTree();
+
+ trees.push_back(tempTree);
+
+ //print tree
+ tempTree->print(outAll);
+ }
+
+ outAll.close();
+
+ if (m->control_pressed) { for (int i = 0; i < trees.size(); i++) { delete trees[i]; trees[i] = NULL; } m->mothurRemove(outputFile); }
+
+ return trees;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "UnifracWeightedCommand", "buildTrees");
+ exit(1);
+ }
+}
+/**************************************************************************************************/
+
+int UnifracWeightedCommand::runRandomCalcs(Tree* thisTree, vector<double> usersScores) {
+ try {
+
+ //calculate number of comparisons i.e. with groups A,B,C = AB, AC, BC = 3;
+ vector< vector<string> > namesOfGroupCombos;
+ for (int a=0; a<numGroups; a++) {
+ for (int l = 0; l < a; l++) {
+ vector<string> groups; groups.push_back((m->getGroups())[a]); groups.push_back((m->getGroups())[l]);
+ namesOfGroupCombos.push_back(groups);
+ }
+ }
+
+ lines.clear();
+
+ //breakdown work between processors
+ int numPairs = namesOfGroupCombos.size();
+ int numPairsPerProcessor = numPairs / processors;
+
+ for (int i = 0; i < processors; i++) {
+ int startPos = i * numPairsPerProcessor;
+ if(i == processors - 1){ numPairsPerProcessor = numPairs - i * numPairsPerProcessor; }
+ lines.push_back(linePair(startPos, numPairsPerProcessor));
+ }
+
+
+ //get scores for random trees
+ for (int j = 0; j < iters; j++) {
+//#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
+ //if(processors == 1){
+ // driver(thisTree, namesOfGroupCombos, 0, namesOfGroupCombos.size(), rScores);
+ // }else{
+ createProcesses(thisTree, namesOfGroupCombos, rScores);
+ // }
+//#else
+ //driver(thisTree, namesOfGroupCombos, 0, namesOfGroupCombos.size(), rScores);
+//#endif
+
+
+ if (m->control_pressed) { delete ct; for (int i = 0; i < T.size(); i++) { delete T[i]; } delete output; outSum.close(); for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; }
+
+ }
+ lines.clear();
+
+ //find the signifigance of the score for summary file
+ for (int f = 0; f < numComp; f++) {
+ //sort random scores
+ sort(rScores[f].begin(), rScores[f].end());
+
+ //the index of the score higher than yours is returned
+ //so if you have 1000 random trees the index returned is 100
+ //then there are 900 trees with a score greater then you.
+ //giving you a signifigance of 0.900
+ int index = findIndex(usersScores[f], f); if (index == -1) { m->mothurOut("error in UnifracWeightedCommand"); m->mothurOutEndLine(); exit(1); } //error code
+
+ //the signifigance is the number of trees with the users score or higher
+ WScoreSig.push_back((iters-index)/(float)iters);
+ }
+
+ //out << "Tree# " << i << endl;
+ calculateFreqsCumuls();
+ printWeightedFile();
+
+ delete output;
+
+ return 0;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "UnifracWeightedCommand", "runRandomCalcs");
+ exit(1);
+ }
+}
+/**************************************************************************************************/
-int UnifracWeightedCommand::createProcesses(Tree* t, vector< vector<string> > namesOfGroupCombos, vector<double>& sums, vector< vector<double> >& scores) {
+int UnifracWeightedCommand::createProcesses(Tree* t, vector< vector<string> > namesOfGroupCombos, vector< vector<double> >& scores) {
try {
-#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
- int process = 1;
- int num = 0;
+ int process = 1;
vector<int> processIDS;
-
EstOutput results;
-
+
+#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
//loop through and create all the processes you want
while (process != processors) {
int pid = fork();
processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
process++;
}else if (pid == 0){
- driver(t, namesOfGroupCombos, lines[process].start, lines[process].num, sums, scores);
+ driver(t, namesOfGroupCombos, lines[process].start, lines[process].num, scores);
//pass numSeqs to parent
ofstream out;
out.close();
exit(0);
- }else { m->mothurOut("unable to spawn the necessary processes."); m->mothurOutEndLine(); exit(0); }
+ }else {
+ m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
+ for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
+ exit(0);
+ }
}
- driver(t, namesOfGroupCombos, lines[0].start, lines[0].num, sums, scores);
+ driver(t, namesOfGroupCombos, lines[0].start, lines[0].num, scores);
//force parent to wait until all the processes are done
for (int i=0;i<(processors-1);i++) {
double tempScore;
for (int j = lines[(i+1)].start; j < (lines[(i+1)].start + lines[(i+1)].num); j++) { in >> tempScore; scores[j].push_back(tempScore); }
in.close();
- remove(s.c_str());
+ m->mothurRemove(s);
}
+#else
+ //fill in functions
+ vector<weightedRandomData*> pDataArray;
+ DWORD dwThreadIdArray[processors-1];
+ HANDLE hThreadArray[processors-1];
+ vector<CountTable*> cts;
+ vector<Tree*> trees;
- return 0;
-#endif
+ //Create processor worker threads.
+ for( int i=1; i<processors; i++ ){
+ CountTable* copyCount = new CountTable();
+ copyCount->copy(ct);
+ Tree* copyTree = new Tree(copyCount);
+ copyTree->getCopy(t);
+
+ cts.push_back(copyCount);
+ trees.push_back(copyTree);
+
+ vector< vector<double> > copyScores = rScores;
+
+ weightedRandomData* tempweighted = new weightedRandomData(m, lines[i].start, lines[i].num, namesOfGroupCombos, copyTree, copyCount, includeRoot, copyScores);
+ pDataArray.push_back(tempweighted);
+ processIDS.push_back(i);
+
+ hThreadArray[i-1] = CreateThread(NULL, 0, MyWeightedRandomThreadFunction, pDataArray[i-1], 0, &dwThreadIdArray[i-1]);
+ }
+
+ driver(t, namesOfGroupCombos, lines[0].start, lines[0].num, scores);
+
+ //Wait until all threads have terminated.
+ WaitForMultipleObjects(processors-1, hThreadArray, TRUE, INFINITE);
+
+ //Close all thread handles and free memory allocations.
+ for(int i=0; i < pDataArray.size(); i++){
+ for (int j = pDataArray[i]->start; j < (pDataArray[i]->start+pDataArray[i]->num); j++) {
+ scores[j].push_back(pDataArray[i]->scores[j][pDataArray[i]->scores[j].size()-1]);
+ }
+ delete cts[i];
+ delete trees[i];
+ CloseHandle(hThreadArray[i]);
+ delete pDataArray[i];
+ }
+
+
+#endif
+
+ return 0;
}
catch(exception& e) {
m->errorOut(e, "UnifracWeightedCommand", "createProcesses");
}
/**************************************************************************************************/
-int UnifracWeightedCommand::driver(Tree* t, vector< vector<string> > namesOfGroupCombos, int start, int num, vector<double>& sums, vector< vector<double> >& scores) {
+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(ct);
+
+ Weighted weighted(includeRoot);
+
for (int h = start; h < (start+num); h++) {
if (m->control_pressed) { return 0; }
if (m->control_pressed) { delete randT; return 0; }
//get wscore of random tree
- EstOutput randomData = weighted->getValues(randT, groupA, groupB, sums);
+ EstOutput randomData = weighted.getValues(randT, groupA, groupB);
if (m->control_pressed) { delete randT; return 0; }
void UnifracWeightedCommand::printWSummaryFile() {
try {
//column headers
- outSum << "Tree#" << '\t' << "Groups" << '\t' << "WScore" << '\t' << "WSig" << endl;
- m->mothurOut("Tree#\tGroups\tWScore\tWSig"); m->mothurOutEndLine();
+ outSum << "Tree#" << '\t' << "Groups" << '\t' << "WScore" << '\t';
+ m->mothurOut("Tree#\tGroups\tWScore\t");
+ if (random) { outSum << "WSig"; m->mothurOut("WSig"); }
+ outSum << endl; m->mothurOutEndLine();
//format output
outSum.setf(ios::fixed, ios::floatfield); outSum.setf(ios::showpoint);
m->mothurOutJustToLog(toString(i+1) +"\t" + groupComb[j] +"\t" + toString(utreeScores[count]) +"\t<" + toString((1/float(iters))) + "\n");
}
}else{
- outSum << setprecision(6) << i+1 << '\t' << groupComb[j] << '\t' << utreeScores[count] << '\t' << "0.00" << endl;
- cout << setprecision(6) << i+1 << '\t' << groupComb[j] << '\t' << utreeScores[count] << '\t' << "0.00" << endl;
- m->mothurOutJustToLog(toString(i+1) +"\t" + groupComb[j] +"\t" + toString(utreeScores[count]) +"\t0.00\n");
+ outSum << setprecision(6) << i+1 << '\t' << groupComb[j] << '\t' << utreeScores[count] << endl;
+ cout << setprecision(6) << i+1 << '\t' << groupComb[j] << '\t' << utreeScores[count] << endl;
+ m->mothurOutJustToLog(toString(i+1) +"\t" + groupComb[j] +"\t" + toString(utreeScores[count]) +"\n");
}
count++;
}
//for each tree
for (int i = 0; i < T.size(); i++) {
- string phylipFileName = outputDir + m->getSimpleName(globaldata->getTreeFile()) + toString(i+1) + ".weighted.dist";
- outputNames.push_back(phylipFileName);
- outputTypes["phylip"].push_back(phylipFileName);
+ string phylipFileName;
+ map<string, string> variables;
+ variables["[filename]"] = outputDir + m->getSimpleName(treefile);
+ variables["[tag]"] = toString(i+1);
+ if ((outputForm == "lt") || (outputForm == "square")) {
+ variables["[tag2]"] = "weighted.phylip";
+ phylipFileName = getOutputFileName("phylip",variables);
+ outputNames.push_back(phylipFileName); outputTypes["phylip"].push_back(phylipFileName);
+ }else { //column
+ variables["[tag2]"] = "weighted.column";
+ phylipFileName = getOutputFileName("column",variables);
+ outputNames.push_back(phylipFileName); outputTypes["column"].push_back(phylipFileName);
+ }
+
+
ofstream out;
m->openOutputFile(phylipFileName, out);
- //output numSeqs
- out << globaldata->Groups.size() << endl;
-
+ if ((outputForm == "lt") || (outputForm == "square")) {
+ //output numSeqs
+ out << m->getNumGroups() << 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->getNumGroups());
+ for (int i = 0; i < m->getNumGroups(); i++) {
+ dists[i].resize(m->getNumGroups(), 0.0);
}
//flip it so you can print it
- for (int r=0; r<globaldata->Groups.size(); r++) {
- for (int l = r+1; l < globaldata->Groups.size(); l++) {
+ for (int r=0; r<m->getNumGroups(); r++) {
+ for (int l = 0; l < r; l++) {
dists[r][l] = utreeScores[count];
dists[l][r] = utreeScores[count];
count++;
}
//output to file
- for (int r=0; r<globaldata->Groups.size(); r++) {
+ for (int r=0; r<m->getNumGroups(); r++) {
//output name
- string name = globaldata->Groups[r];
+ string name = (m->getGroups())[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 < r; l++) { out << dists[r][l] << '\t'; }
- out << endl;
+ if (outputForm == "lt") {
+ out << name << '\t';
+
+ //output distances
+ for (int l = 0; l < r; l++) { out << dists[r][l] << '\t'; }
+ out << endl;
+ }else if (outputForm == "square") {
+ out << name << '\t';
+
+ //output distances
+ for (int l = 0; l < m->getNumGroups(); l++) { out << dists[r][l] << '\t'; }
+ out << endl;
+ }else{
+ //output distances
+ for (int l = 0; l < r; l++) {
+ string otherName = (m->getGroups())[l];
+ if (otherName.length() < 10) { //pad with spaces to make compatible
+ while (otherName.length() < 10) { otherName += " "; }
+ }
+
+ out << name << '\t' << otherName << dists[r][l] << endl;
+ }
+ }
}
out.close();
}
exit(1);
}
}
-
/***********************************************************/
projects / mothur.git / blobdiff