*/
#include "unifracunweightedcommand.h"
+#include "treereader.h"
+#include "subsample.h"
+#include "consensus.h"
//**********************************************************************************************************************
vector<string> UnifracUnweightedCommand::setParameters(){
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 psubsample("subsample", "String", "", "", "", "", "",false,false); parameters.push_back(psubsample);
+ CommandParameter pconsensus("consensus", "Boolean", "", "F", "", "", "",false,false); parameters.push_back(pconsensus);
+ 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);
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.unweighted command should be in the following format: unifrac.unweighted(groups=yourGroups, iters=yourIters).\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 of the subsampling, as well as a consensus tree built from these trees. Default=F.\n";
helpString += "Example unifrac.unweighted(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.unweighted command output two files: .unweighted and .uwsummary their descriptions are in the manual.\n";
exit(1);
}
}
+//**********************************************************************************************************************
+string UnifracUnweightedCommand::getOutputFileNameTag(string type, string inputName=""){
+ try {
+ string outputFileName = "";
+ map<string, vector<string> >::iterator it;
+
+ //is this a type this command creates
+ it = outputTypes.find(type);
+ if (it == outputTypes.end()) { m->mothurOut("[ERROR]: this command doesn't create a " + type + " output file.\n"); }
+ else {
+ if (type == "unweighted") { outputFileName = "unweighted"; }
+ else if (type == "uwsummary") { outputFileName = "uwsummary"; }
+ else if (type == "phylip") { outputFileName = "dist"; }
+ else if (type == "column") { outputFileName = "dist"; }
+ else if (type == "tree") { outputFileName = "tre"; }
+ else { m->mothurOut("[ERROR]: No definition for type " + type + " output file tag.\n"); m->control_pressed = true; }
+ }
+ return outputFileName;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "UnifracUnweightedCommand", "getOutputFileNameTag");
+ exit(1);
+ }
+}
+
//**********************************************************************************************************************
UnifracUnweightedCommand::UnifracUnweightedCommand(){
try {
outputTypes["uwsummary"] = tempOutNames;
outputTypes["phylip"] = tempOutNames;
outputTypes["column"] = tempOutNames;
+ outputTypes["tree"] = tempOutNames;
}
catch(exception& e) {
m->errorOut(e, "UnifracUnweightedCommand", "UnifracUnweightedCommand");
//allow user to run help
if(option == "help") { help(); abort = true; calledHelp = true; }
+ else if(option == "citation") { citation(); abort = true; calledHelp = true;}
else {
vector<string> myArray = setParameters();
outputTypes["uwsummary"] = tempOutNames;
outputTypes["phylip"] = tempOutNames;
outputTypes["column"] = tempOutNames;
+ outputTypes["tree"] = tempOutNames;
//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);
}
}
- m->runParse = true;
- m->Groups.clear();
- m->namesOfGroups.clear();
- m->Treenames.clear();
- m->names.clear();
-
- //check for required parameters
+ //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; }
- }
+ }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") { abort = true; }
+ if (namefile == "not open") { namefile = ""; abort = true; }
else if (namefile == "not found") { namefile = ""; }
+ else { m->setNameFile(namefile); }
- outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){ outputDir = ""; }
+ 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);
- m->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") { phylip = false; outputForm = ""; }
temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); }
m->setProcessors(temp);
- convert(temp, processors);
-
+ 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 (!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 (subsample && (groupfile == "")) { m->mothurOut("[ERROR]: if subsample=t, a group file 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 (!random) { iters = 0; } //turn off random calcs
//if user selects distance = true and no groups it won't calc the pairwise
if ((phylip) && (Groups.size() == 0)) {
groups = "all";
m->splitAtDash(groups, Groups);
- m->Groups = Groups;
+ m->setGroups(Groups);
+ }
+
+ if (namefile == "") {
+ vector<string> files; files.push_back(treefile);
+ parser.getNameFile(files);
}
}
m->setTreeFile(treefile);
- 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) + ".uwsummary";
+ TreeReader* reader = new TreeReader(treefile, groupfile, namefile);
+ T = reader->getTrees();
+ tmap = T[0]->getTreeMap();
+ map<string, string> nameMap = reader->getNames();
+ map<string, string> unique2Dup = reader->getNameMap();
+ delete reader;
+
+ sumFile = outputDir + m->getRootName(m->getSimpleName(treefile)) + getOutputFileNameTag("uwsummary");
outputNames.push_back(sumFile); outputTypes["uwsummary"].push_back(sumFile);
m->openOutputFile(sumFile, outSum);
- util = new SharedUtil();
- util->setGroups(m->Groups, tmap->namesOfGroups, allGroups, numGroups, "unweighted"); //sets the groups the user wants to analyze
- util->getCombos(groupComb, m->Groups, numComp);
- delete util;
-
- if (numGroups == 1) { numComp++; groupComb.push_back(allGroups); }
+ SharedUtil util;
+ Groups = m->getGroups();
+ vector<string> namesGroups = tmap->getNamesOfGroups();
+ util.setGroups(Groups, namesGroups, allGroups, numGroups, "unweighted"); //sets the groups the user wants to analyze
- unweighted = new Unweighted(tmap, includeRoot);
+ Unweighted unweighted(includeRoot);
int start = time(NULL);
-
- userData.resize(numComp,0); //data[0] = unweightedscore
- randomData.resize(numComp,0); //data[0] = unweightedscore
- //create new tree with same num nodes and leaves as users
-
+
+ //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 = (tmap->getNamesSeqs(temp)).size(); //num in first group
+ for (int i = 1; i < Groups.size(); i++) {
+ temp.clear(); temp.push_back(Groups[i]);
+ int thisSize = (tmap->getNamesSeqs(temp)).size();
+ 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++) {
+ vector<string> thisGroup; thisGroup.push_back(newGroups[i]);
+ vector<string> thisGroupsSeqs = tmap->getNamesSeqs(thisGroup);
+ int thisSize = thisGroupsSeqs.size();
+
+ 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);
+ }
+ }
+
+ util.getCombos(groupComb, Groups, numComp);
+ m->setGroups(Groups);
+
+ if (numGroups == 1) { numComp++; groupComb.push_back(allGroups); }
+
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; }
- outSum << "Tree#" << '\t' << "Groups" << '\t' << "UWScore" <<'\t' << "UWSig" << endl;
- m->mothurOut("Tree#\tGroups\tUWScore\tUWSig"); m->mothurOutEndLine();
+ outSum << "Tree#" << '\t' << "Groups" << '\t' << "UWScore" <<'\t';
+ m->mothurOut("Tree#\tGroups\tUWScore\t");
+ if (random) { outSum << "UWSig"; m->mothurOut("UWSig"); }
+ outSum << endl; m->mothurOutEndLine();
//get pscores for users trees
for (int i = 0; i < T.size(); i++) {
- if (m->control_pressed) {
- delete tmap; delete unweighted;
- 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;
- }
+ if (m->control_pressed) { delete tmap; 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;
+ counter = 0;
if (random) {
- output = new ColumnFile(outputDir + m->getSimpleName(treefile) + toString(i+1) + ".unweighted", itersString);
- outputNames.push_back(outputDir + m->getSimpleName(treefile) + toString(i+1) + ".unweighted");
- outputTypes["unweighted"].push_back(outputDir + m->getSimpleName(treefile) + toString(i+1) + ".unweighted");
+ output = new ColumnFile(outputDir + m->getSimpleName(treefile) + toString(i+1) + "." + getOutputFileNameTag("unweighted"), itersString);
+ outputNames.push_back(outputDir + m->getSimpleName(treefile) + toString(i+1) + "." + getOutputFileNameTag("unweighted"));
+ outputTypes["unweighted"].push_back(outputDir + m->getSimpleName(treefile) + toString(i+1) + "." + getOutputFileNameTag("unweighted"));
}
rCumul.resize(numComp);
utreeScores.resize(numComp);
UWScoreSig.resize(numComp);
+
+ vector<double> userData; userData.resize(numComp,0); //weighted score info for user tree. data[0] = weightedscore AB, data[1] = weightedscore AC...
- userData = unweighted->getValues(T[i], processors, outputDir); //userData[0] = unweightedscore
+ userData = unweighted.getValues(T[i], processors, outputDir); //userData[0] = unweightedscore
- if (m->control_pressed) { delete tmap; delete unweighted;
- 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; }
+ if (m->control_pressed) { delete tmap; 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; }
//output scores for each combination
for(int k = 0; k < numComp; k++) {
//add users score to validscores
validScores[userData[k]] = userData[k];
+
+ if (!random) { UWScoreSig[k].push_back(0.0); }
}
-
- //get unweighted scores for random trees - if random is false iters = 0
- for (int j = 0; j < iters; j++) {
-
- //we need a different getValues because when we swap the labels we only want to swap those in each pairwise comparison
- randomData = unweighted->getValues(T[i], "", "", processors, outputDir);
-
- if (m->control_pressed) { delete tmap; delete unweighted;
- 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; }
-
- for(int k = 0; k < numComp; k++) {
- //add trees unweighted score to map of scores
- map<float,float>::iterator it = rscoreFreq[k].find(randomData[k]);
- if (it != rscoreFreq[k].end()) {//already have that score
- rscoreFreq[k][randomData[k]]++;
- }else{//first time we have seen this score
- rscoreFreq[k][randomData[k]] = 1;
- }
-
- //add randoms score to validscores
- validScores[randomData[k]] = randomData[k];
- }
-
- //report progress
-// m->mothurOut("Iter: " + toString(j+1)); m->mothurOutEndLine();
- }
-
- for(int a = 0; a < numComp; a++) {
- float rcumul = 1.0000;
-
- if (random) {
- //this loop fills the cumulative maps and put 0.0000 in the score freq map to make it easier to print.
- for (map<float,float>::iterator it = validScores.begin(); it != validScores.end(); it++) {
- //make rscoreFreq map and rCumul
- map<float,float>::iterator it2 = rscoreFreq[a].find(it->first);
- rCumul[a][it->first] = rcumul;
- //get percentage of random trees with that info
- if (it2 != rscoreFreq[a].end()) { rscoreFreq[a][it->first] /= iters; rcumul-= it2->second; }
- else { rscoreFreq[a][it->first] = 0.0000; } //no random trees with that score
- }
- UWScoreSig[a].push_back(rCumul[a][userData[a]]);
- }else { UWScoreSig[a].push_back(0.0); }
-
- }
-
- if (m->control_pressed) { delete tmap; delete unweighted;
- 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; }
-
- //print output files
+
+ if (random) { runRandomCalcs(T[i], userData); }
+
+ if (m->control_pressed) { delete tmap; 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; }
+
+ int startSubsample = time(NULL);
+
+ //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.
+ TreeMap* newTmap = new TreeMap();
+ //newTmap->getCopy(*tmap);
+
+ //SubSample sample;
+ //Tree* subSampleTree = sample.getSample(T[i], newTmap, nameMap, subsampleSize);
+
+ //uses method of setting groups to doNotIncludeMe
+ SubSample sample;
+ Tree* subSampleTree = sample.getSample(T[i], tmap, newTmap, subsampleSize, unique2Dup);
+
+ //call new weighted function
+ vector<double> iterData; iterData.resize(numComp,0);
+ Unweighted thisUnweighted(includeRoot);
+ iterData = thisUnweighted.getValues(subSampleTree, processors, outputDir); //userData[0] = weightedscore
+
+ //save data to make ave dist, std dist
+ calcDistsTotals.push_back(iterData);
+
+ delete newTmap;
+ delete subSampleTree;
+
+ if((thisIter+1) % 100 == 0){ m->mothurOut(toString(thisIter+1)); m->mothurOutEndLine(); }
+ }
+ m->mothurOut("It took " + toString(time(NULL) - startSubsample) + " secs to run the subsampling."); m->mothurOutEndLine();
+
+ if (m->control_pressed) { delete tmap; 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); }
+
+ //print output files
printUWSummaryFile(i);
if (random) { printUnweightedFile(); delete output; }
if (phylip) { createPhylipFile(i); }
outSum.close();
- m->Groups.clear();
- delete tmap; delete unweighted;
+ delete tmap;
for (int i = 0; i < T.size(); i++) { delete T[i]; }
- if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
+ if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; }
m->mothurOut("It took " + toString(time(NULL) - start) + " secs to run unifrac.unweighted."); m->mothurOutEndLine();
exit(1);
}
}
+/**************************************************************************************************/
+int UnifracUnweightedCommand::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; averages.resize(numComp, 0);
+ for (int thisIter = 0; thisIter < subsampleIters; thisIter++) {
+ for (int i = 0; i < dists[thisIter].size(); i++) {
+ averages[i] += dists[thisIter][i];
+ }
+ }
+
+ //finds average.
+ for (int i = 0; i < averages.size(); i++) { averages[i] /= (float) subsampleIters; }
+
+ //find standard deviation
+ vector<double> stdDev; stdDev.resize(numComp, 0);
+
+ for (int thisIter = 0; thisIter < iters; thisIter++) { //compute the difference of each dist from the mean, and square the result of each
+ for (int j = 0; j < dists[thisIter].size(); j++) {
+ stdDev[j] += ((dists[thisIter][j] - averages[j]) * (dists[thisIter][j] - averages[j]));
+ }
+ }
+ for (int i = 0; i < stdDev.size(); i++) {
+ stdDev[i] /= (float) subsampleIters;
+ stdDev[i] = sqrt(stdDev[i]);
+ }
+
+ //make matrix with scores in it
+ vector< vector<double> > avedists; avedists.resize(m->getNumGroups());
+ for (int i = 0; i < m->getNumGroups(); i++) {
+ avedists[i].resize(m->getNumGroups(), 0.0);
+ }
+
+ //make matrix with scores in it
+ vector< vector<double> > stddists; stddists.resize(m->getNumGroups());
+ for (int i = 0; i < m->getNumGroups(); i++) {
+ stddists[i].resize(m->getNumGroups(), 0.0);
+ }
+
+ //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++;
+ }
+ }
+
+ string aveFileName = outputDir + m->getSimpleName(treefile) + toString(treeNum+1) + ".unweighted.ave." + getOutputFileNameTag("phylip");
+ 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);
+
+ string stdFileName = outputDir + m->getSimpleName(treefile) + toString(treeNum+1) + ".unweighted.std." + getOutputFileNameTag("phylip");
+ 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, "UnifracUnweightedCommand", "getAverageSTDMatrices");
+ exit(1);
+ }
+}
+
+/**************************************************************************************************/
+int UnifracUnweightedCommand::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
+ TreeMap newTmap;
+ newTmap.makeSim(m->getGroups());
+
+ //clear old tree names if any
+ m->Treenames.clear();
+
+ //fills globaldatas tree names
+ m->Treenames = m->getGroups();
+
+ vector<Tree*> newTrees = buildTrees(dists, treeNum, newTmap); //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
+ string conFile = outputDir + m->getRootName(m->getSimpleName(treefile)) + toString(treeNum+1) + ".unweighted.cons." + getOutputFileNameTag("tree");
+ 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, "UnifracUnweightedCommand", "getConsensusTrees");
+ exit(1);
+ }
+}
+/**************************************************************************************************/
+
+vector<Tree*> UnifracUnweightedCommand::buildTrees(vector< vector<double> >& dists, int treeNum, TreeMap& mytmap) {
+ try {
+
+ vector<Tree*> trees;
+
+ //create a new filename
+ string outputFile = outputDir + m->getRootName(m->getSimpleName(treefile)) + toString(treeNum+1) + ".unweighted.all." + getOutputFileNameTag("tree");
+ 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(&mytmap, sims);
+ map<string, string> empty;
+ tempTree->assembleTree(empty);
+
+ 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, "UnifracUnweightedCommand", "buildTrees");
+ exit(1);
+ }
+}
+/**************************************************************************************************/
+
+int UnifracUnweightedCommand::runRandomCalcs(Tree* thisTree, vector<double> usersScores) {
+ try {
+ vector<double> randomData; randomData.resize(numComp,0); //weighted score info for random trees. data[0] = weightedscore AB, data[1] = weightedscore AC...
+
+ Unweighted unweighted(includeRoot);
+
+ //get unweighted scores for random trees - if random is false iters = 0
+ for (int j = 0; j < iters; j++) {
+
+ //we need a different getValues because when we swap the labels we only want to swap those in each pairwise comparison
+ randomData = unweighted.getValues(thisTree, "", "", processors, outputDir);
+
+ if (m->control_pressed) { return 0; }
+
+ for(int k = 0; k < numComp; k++) {
+ //add trees unweighted score to map of scores
+ map<float,float>::iterator it = rscoreFreq[k].find(randomData[k]);
+ if (it != rscoreFreq[k].end()) {//already have that score
+ rscoreFreq[k][randomData[k]]++;
+ }else{//first time we have seen this score
+ rscoreFreq[k][randomData[k]] = 1;
+ }
+
+ //add randoms score to validscores
+ validScores[randomData[k]] = randomData[k];
+ }
+ }
+
+ for(int a = 0; a < numComp; a++) {
+ float rcumul = 1.0000;
+
+ //this loop fills the cumulative maps and put 0.0000 in the score freq map to make it easier to print.
+ for (map<float,float>::iterator it = validScores.begin(); it != validScores.end(); it++) {
+ //make rscoreFreq map and rCumul
+ map<float,float>::iterator it2 = rscoreFreq[a].find(it->first);
+ rCumul[a][it->first] = rcumul;
+ //get percentage of random trees with that info
+ if (it2 != rscoreFreq[a].end()) { rscoreFreq[a][it->first] /= iters; rcumul-= it2->second; }
+ else { rscoreFreq[a][it->first] = 0.0000; } //no random trees with that score
+ }
+ UWScoreSig[a].push_back(rCumul[a][usersScores[a]]);
+ }
+
+ return 0;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "UnifracUnweightedCommand", "runRandomCalcs");
+ exit(1);
+ }
+}
/***********************************************************/
void UnifracUnweightedCommand::printUnweightedFile() {
try {
m->mothurOutJustToLog(groupComb[a] + "\t" + toString(utreeScores[a][0]) + "\t<" + toString((1/float(iters))) + "\n");
}
}else{
- outSum << setprecision(6) << groupComb[a] << '\t' << utreeScores[a][0] << '\t' << "0.00" << endl;
- cout << setprecision(6) << groupComb[a] << '\t' << utreeScores[a][0] << '\t' << "0.00" << endl;
- m->mothurOutJustToLog(groupComb[a] + "\t" + toString(utreeScores[a][0]) + "\t0.00\n");
+ outSum << setprecision(6) << groupComb[a] << '\t' << utreeScores[a][0] << endl;
+ cout << setprecision(6) << groupComb[a] << '\t' << utreeScores[a][0] << endl;
+ m->mothurOutJustToLog(groupComb[a] + "\t" + toString(utreeScores[a][0]) + "\n");
}
}
try {
string phylipFileName;
if ((outputForm == "lt") || (outputForm == "square")) {
- phylipFileName = outputDir + m->getSimpleName(treefile) + toString(i+1) + ".unweighted.phylip.dist";
+ phylipFileName = outputDir + m->getSimpleName(treefile) + toString(i+1) + ".unweighted.phylip." + getOutputFileNameTag("phylip");
outputNames.push_back(phylipFileName); outputTypes["phylip"].push_back(phylipFileName);
}else { //column
- phylipFileName = outputDir + m->getSimpleName(treefile) + toString(i+1) + ".unweighted.column.dist";
+ phylipFileName = outputDir + m->getSimpleName(treefile) + toString(i+1) + ".unweighted.column." + getOutputFileNameTag("column");
outputNames.push_back(phylipFileName); outputTypes["column"].push_back(phylipFileName);
}
if ((outputForm == "lt") || (outputForm == "square")) {
//output numSeqs
- out << m->Groups.size() << endl;
+ out << m->getNumGroups() << endl;
}
//make matrix with scores in it
- 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);
+ 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
int count = 0;
- for (int r=0; r<m->Groups.size(); r++) {
+ for (int r=0; r<m->getNumGroups(); r++) {
for (int l = 0; l < r; l++) {
dists[r][l] = utreeScores[count][0];
dists[l][r] = utreeScores[count][0];
}
//output to file
- for (int r=0; r<m->Groups.size(); r++) {
+ for (int r=0; r<m->getNumGroups(); r++) {
//output name
- string name = m->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 < m->Groups.size(); l++) { out << dists[r][l] << '\t'; }
+ 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->Groups[l];
+ string otherName = (m->getGroups())[l];
if (otherName.length() < 10) { //pad with spaces to make compatible
while (otherName.length() < 10) { otherName += " "; }
}
m->errorOut(e, "UnifracUnweightedCommand", "createPhylipFile");
exit(1);
}
-}/*****************************************************************/
-int UnifracUnweightedCommand::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, "UnifracUnweightedCommand", "readNamesFile");
- exit(1);
- }
}
/***********************************************************/
projects / mothur.git / blobdiff