From 03dca3b32a903c3f29fbcf5b410b19d6ab6dae63 Mon Sep 17 00:00:00 2001 From: Sarah Westcott Date: Wed, 11 Apr 2012 08:07:05 -0400 Subject: [PATCH] added subsample and consensus parameters to unifrac.weighted command --- indicatorcommand.cpp | 1 - preclustercommand.cpp | 28 +- subsample.cpp | 110 +++++- subsample.h | 9 +- tree.cpp | 199 ++++++++--- tree.h | 6 +- treegroupscommand.cpp | 10 +- treemap.cpp | 54 +++ treemap.h | 16 +- trimseqscommand.cpp | 11 + unifracweightedcommand.cpp | 684 +++++++++++++++++++++++++++---------- unifracweightedcommand.h | 12 +- weighted.cpp | 4 +- 13 files changed, 897 insertions(+), 247 deletions(-) diff --git a/indicatorcommand.cpp b/indicatorcommand.cpp index 01d8c2e..fe99c46 100644 --- a/indicatorcommand.cpp +++ b/indicatorcommand.cpp @@ -334,7 +334,6 @@ int IndicatorCommand::execute(){ //no longer need original tree, we have output tree to use and label for (int i = 0; i < T.size(); i++) { delete T[i]; } - if (m->control_pressed) { if (designfile != "") { delete designMap; } if (sharedfile != "") { for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } } diff --git a/preclustercommand.cpp b/preclustercommand.cpp index f2fbc80..bcff0fc 100644 --- a/preclustercommand.cpp +++ b/preclustercommand.cpp @@ -314,7 +314,17 @@ int PreClusterCommand::createProcessesGroups(SequenceParser* parser, string newF 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){ + outputNames.clear(); num = driverGroups(parser, newFName + toString(getpid()) + ".temp", newNName + toString(getpid()) + ".temp", newMFile, lines[process].start, lines[process].end, groups); + + string tempFile = toString(getpid()) + ".outputNames.temp"; + ofstream outTemp; + m->openOutputFile(tempFile, outTemp); + + outTemp << outputNames.size(); + for (int i = 0; i < outputNames.size(); i++) { outTemp << outputNames[i] << endl; } + outTemp.close(); + exit(0); }else { m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine(); @@ -331,7 +341,23 @@ int PreClusterCommand::createProcessesGroups(SequenceParser* parser, string newF int temp = processIDS[i]; wait(&temp); } - + + for (int i = 0; i < processIDS.size(); i++) { + string tempFile = toString(processIDS[i]) + ".outputNames.temp"; + ifstream intemp; + m->openInputFile(tempFile, intemp); + + int num; + intemp >> num; + for (int k = 0; k < num; k++) { + string name = ""; + intemp >> name; m->gobble(intemp); + + outputNames.push_back(name); outputTypes["map"].push_back(name); + } + intemp.close(); + m->mothurRemove(tempFile); + } #else ////////////////////////////////////////////////////////////////////////////////////////////////////// diff --git a/subsample.cpp b/subsample.cpp index e6dd845..b1e78a4 100644 --- a/subsample.cpp +++ b/subsample.cpp @@ -8,6 +8,114 @@ #include "subsample.h" +//********************************************************************************************************************** +Tree* SubSample::getSample(Tree* T, TreeMap* tmap, map whole, int size) { + try { + Tree* newTree = NULL; + + vector subsampledSeqs = getSample(tmap, size); + map sampledNameMap = deconvolute(whole, subsampledSeqs); + + //remove seqs not in sample from treemap + for (int i = 0; i < tmap->namesOfSeqs.size(); i++) { + //is that name in the subsample? + int count = 0; + for (int j = 0; j < subsampledSeqs.size(); j++) { + if (tmap->namesOfSeqs[i] == subsampledSeqs[j]) { break; } //found it + count++; + } + + if (m->control_pressed) { return newTree; } + + //if you didnt find it, remove it + if (count == subsampledSeqs.size()) { + tmap->removeSeq(tmap->namesOfSeqs[i]); + i--; //need this because removeSeq removes name from namesOfSeqs + } + } + + //create new tree + int numUniques = sampledNameMap.size(); + if (sampledNameMap.size() == 0) { numUniques = subsampledSeqs.size(); } + + newTree = new Tree(numUniques, tmap); //numNodes, treemap + newTree->getSubTree(T, subsampledSeqs, sampledNameMap); + + return newTree; + } + catch(exception& e) { + m->errorOut(e, "SubSample", "getSample-Tree"); + exit(1); + } +} +//********************************************************************************************************************** +//assumes whole maps dupName -> uniqueName +map SubSample::deconvolute(map whole, vector& wanted) { + try { + map nameMap; + + //whole will be empty if user gave no name file, so we don't need to make a new one + if (whole.size() == 0) { return nameMap; } + + vector newWanted; + for (int i = 0; i < wanted.size(); i++) { + + if (m->control_pressed) { break; } + + string dupName = wanted[i]; + + map::iterator itWhole = whole.find(dupName); + if (itWhole != whole.end()) { + string repName = itWhole->second; + + //do we already have this rep? + map::iterator itName = nameMap.find(repName); + if (itName != nameMap.end()) { //add this seqs to dups list + (itName->second) += "," + dupName; + }else { //first sighting of this seq + nameMap[repName] = dupName; + newWanted.push_back(repName); + } + }else { m->mothurOut("[ERROR]: "+dupName+" is not in your name file, please correct.\n"); m->control_pressed = true; } + } + + wanted = newWanted; + return nameMap; + } + catch(exception& e) { + m->errorOut(e, "SubSample", "deconvolute"); + exit(1); + } +} +//********************************************************************************************************************** +vector SubSample::getSample(TreeMap* tMap, int size) { + try { + vector sample; + + vector Groups = tMap->getNamesOfGroups(); + for (int i = 0; i < Groups.size(); i++) { + + if (m->control_pressed) { break; } + + vector thisGroup; thisGroup.push_back(Groups[i]); + vector thisGroupsSeqs = tMap->getNamesSeqs(thisGroup); + int thisSize = thisGroupsSeqs.size(); + + if (thisSize >= size) { + + random_shuffle(thisGroupsSeqs.begin(), thisGroupsSeqs.end()); + + for (int j = 0; j < size; j++) { sample.push_back(thisGroupsSeqs[j]); } + }else { m->mothurOut("[ERROR]: You have selected a size that is larger than "+Groups[i]+" number of sequences.\n"); m->control_pressed = true; } + } + + return sample; + } + catch(exception& e) { + m->errorOut(e, "SubSample", "getSample-TreeMap"); + exit(1); + } +} //********************************************************************************************************************** vector SubSample::getSample(vector& thislookup, int size) { try { @@ -64,7 +172,7 @@ vector SubSample::getSample(vector& thislookup, int } catch(exception& e) { - m->errorOut(e, "SubSample", "getSample"); + m->errorOut(e, "SubSample", "getSample-shared"); exit(1); } } diff --git a/subsample.h b/subsample.h index 09c7dcd..aaf5244 100644 --- a/subsample.h +++ b/subsample.h @@ -11,6 +11,8 @@ #include "mothurout.h" #include "sharedrabundvector.h" +#include "treemap.h" +#include "tree.h" //subsampling overwrites the sharedRabunds. If you need to reuse the original use the getSamplePreserve function. @@ -22,12 +24,17 @@ class SubSample { ~SubSample() {} vector getSample(vector&, int); //returns the bin labels for the subsample, mothurOuts binlabels are preserved so you can run this multiple times. Overwrites original vector passed in, if you need to preserve it deep copy first. - + + Tree* getSample(Tree*, TreeMap*, map, int); //creates new subsampled tree, destroys treemap so copy if needed. private: MothurOut* m; int eliminateZeroOTUS(vector&); + + vector getSample(TreeMap*, int); //returns map contains names of seqs in subsample -> group. + map deconvolute(map wholeSet, vector& subsampleWanted); //returns new nameMap containing only subsampled names, and removes redundants from subsampled wanted because it makes the new nameMap. + }; diff --git a/tree.cpp b/tree.cpp index d9d71ae..432811e 100644 --- a/tree.cpp +++ b/tree.cpp @@ -89,11 +89,123 @@ Tree::Tree(TreeMap* t) : tmap(t) { exit(1); } } - +/*****************************************************************/ +Tree::Tree(TreeMap* t, vector< vector >& sims) : tmap(t) { + try { + m = MothurOut::getInstance(); + + if (m->runParse == true) { parseTreeFile(); m->runParse = false; } + //for(int i = 0; i < globaldata->Treenames.size(); i++) { cout << i << '\t' << globaldata->Treenames[i] << endl; } + numLeaves = m->Treenames.size(); + numNodes = 2*numLeaves - 1; + + tree.resize(numNodes); + + //initialize groupNodeInfo + for (int i = 0; i < (tmap->getNamesOfGroups()).size(); i++) { + groupNodeInfo[(tmap->getNamesOfGroups())[i]].resize(0); + } + + //initialize tree with correct number of nodes, name and group info. + for (int i = 0; i < numNodes; i++) { + //initialize leaf nodes + if (i <= (numLeaves-1)) { + tree[i].setName(m->Treenames[i]); + + //save group info + string group = tmap->getGroup(m->Treenames[i]); + + vector tempGroups; tempGroups.push_back(group); + tree[i].setGroup(tempGroups); + groupNodeInfo[group].push_back(i); + + //set pcount and pGroup for groupname to 1. + tree[i].pcount[group] = 1; + tree[i].pGroups[group] = 1; + + //Treemap knows name, group and index to speed up search + tmap->setIndex(m->Treenames[i], i); + + //intialize non leaf nodes + }else if (i > (numLeaves-1)) { + tree[i].setName(""); + vector tempGroups; + tree[i].setGroup(tempGroups); + } + } + + //build tree from matrix + //initialize indexes + map indexes; //maps row in simMatrix to vector index in the tree + int numGroups = (tmap->getNamesOfGroups()).size(); + for (int g = 0; g < numGroups; g++) { indexes[g] = g; } + + //do merges and create tree structure by setting parents and children + //there are numGroups - 1 merges to do + for (int i = 0; i < (numGroups - 1); i++) { + float largest = -1000.0; + + if (m->control_pressed) { break; } + + int row, column; + //find largest value in sims matrix by searching lower triangle + for (int j = 1; j < sims.size(); j++) { + for (int k = 0; k < j; k++) { + if (sims[j][k] > largest) { largest = sims[j][k]; row = j; column = k; } + } + } + + //set non-leaf node info and update leaves to know their parents + //non-leaf + tree[numGroups + i].setChildren(indexes[row], indexes[column]); + + //parents + tree[indexes[row]].setParent(numGroups + i); + tree[indexes[column]].setParent(numGroups + i); + + //blength = distance / 2; + float blength = ((1.0 - largest) / 2); + + //branchlengths + tree[indexes[row]].setBranchLength(blength - tree[indexes[row]].getLengthToLeaves()); + tree[indexes[column]].setBranchLength(blength - tree[indexes[column]].getLengthToLeaves()); + + //set your length to leaves to your childs length plus branchlength + tree[numGroups + i].setLengthToLeaves(tree[indexes[row]].getLengthToLeaves() + tree[indexes[row]].getBranchLength()); + + + //update index + indexes[row] = numGroups+i; + indexes[column] = numGroups+i; + + //remove highest value that caused the merge. + sims[row][column] = -1000.0; + sims[column][row] = -1000.0; + + //merge values in simsMatrix + for (int n = 0; n < sims.size(); n++) { + //row becomes merge of 2 groups + sims[row][n] = (sims[row][n] + sims[column][n]) / 2; + sims[n][row] = sims[row][n]; + //delete column + sims[column][n] = -1000.0; + sims[n][column] = -1000.0; + } + } + + //adjust tree to make sure root to tip length is .5 + int root = findRoot(); + tree[root].setBranchLength((0.5 - tree[root].getLengthToLeaves())); + } + catch(exception& e) { + m->errorOut(e, "Tree", "Tree"); + exit(1); + } +} /*****************************************************************/ Tree::~Tree() {} /*****************************************************************/ -void Tree::addNamesToCounts() { +void Tree::addNamesToCounts(map nameMap) { try { //ex. seq1 seq2,seq3,se4 // seq1 = pasture @@ -116,12 +228,12 @@ void Tree::addNamesToCounts() { string name = tree[i].getName(); - map::iterator itNames = m->names.find(name); + map::iterator itNames = nameMap.find(name); - if (itNames == m->names.end()) { m->mothurOut(name + " is not in your name file, please correct."); m->mothurOutEndLine(); exit(1); } + if (itNames == nameMap.end()) { m->mothurOut(name + " is not in your name file, please correct."); m->mothurOutEndLine(); exit(1); } else { vector dupNames; - m->splitAtComma(m->names[name], dupNames); + m->splitAtComma(nameMap[name], dupNames); map::iterator itCounts; int maxPars = 1; @@ -222,7 +334,7 @@ int Tree::assembleTree() { //float A = clock(); //if user has given a names file we want to include that info in the pgroups and pcount info. - if(m->names.size() != 0) { addNamesToCounts(); } + if(m->names.size() != 0) { addNamesToCounts(m->names); } //build the pGroups in non leaf nodes to be used in the parsimony calcs. for (int i = numLeaves; i < numNodes; i++) { @@ -261,9 +373,15 @@ int Tree::assembleTree(string n) { } } /*****************************************************************/ -void Tree::getSubTree(Tree* copy, vector Groups) { +//assumes leaf node names are in groups and no names file - used by indicator command +void Tree::getSubTree(Tree* Ctree, vector Groups) { try { - + + //copy Tree since we are going to destroy it + Tree* copy = new Tree(tmap); + copy->getCopy(Ctree); + copy->assembleTree("nonames"); + //we want to select some of the leaf nodes to create the output tree //go through the input Tree starting at parents of leaves for (int i = 0; i < numNodes; i++) { @@ -408,12 +526,40 @@ void Tree::getSubTree(Tree* copy, vector Groups) { //you found the root if (copy->tree[i].getParent() == -1) { root = i; break; } } - + int nextSpot = numLeaves; populateNewTree(copy->tree, root, nextSpot); + + delete copy; } catch(exception& e) { - m->errorOut(e, "Tree", "getCopy"); + m->errorOut(e, "Tree", "getSubTree"); + exit(1); + } +} +/*****************************************************************/ +//assumes nameMap contains unique names as key or is empty. +//assumes numLeaves defined in tree constructor equals size of seqsToInclude and seqsToInclude only contains unique seqs. +int Tree::getSubTree(Tree* copy, vector seqsToInclude, map nameMap) { + try { + + if (numLeaves != seqsToInclude.size()) { m->mothurOut("[ERROR]: numLeaves does not equal numUniques, cannot create subtree.\n"); m->control_pressed = true; return 0; } + + getSubTree(copy, seqsToInclude); + if (nameMap.size() != 0) { addNamesToCounts(nameMap); } + + //build the pGroups in non leaf nodes to be used in the parsimony calcs. + for (int i = numLeaves; i < numNodes; i++) { + if (m->control_pressed) { return 1; } + + tree[i].pGroups = (mergeGroups(i)); + tree[i].pcount = (mergeGcounts(i)); + } + + return 0; + } + catch(exception& e) { + m->errorOut(e, "Tree", "getSubTree"); exit(1); } } @@ -627,7 +773,6 @@ map Tree::mergeGcounts(int position) { } } /**************************************************************************************************/ - void Tree::randomLabels(vector g) { try { @@ -676,37 +821,7 @@ void Tree::randomLabels(vector g) { exit(1); } } -/************************************************************************************************** - -void Tree::randomLabels(string groupA, string groupB) { - try { - int numSeqsA = globaldata->gTreemap->seqsPerGroup[groupA]; - int numSeqsB = globaldata->gTreemap->seqsPerGroup[groupB]; - - vector randomGroups(numSeqsA+numSeqsB, groupA); - for(int i=numSeqsA;ierrorOut(e, "Tree", "randomLabels"); - exit(1); - } -} -**************************************************************************************************/ +/**************************************************************************************************/ void Tree::randomBlengths() { try { for(int i=numNodes-1;i>=0;i--){ diff --git a/tree.h b/tree.h index 2d9d4f8..0b61c6e 100644 --- a/tree.h +++ b/tree.h @@ -19,13 +19,17 @@ public: Tree(string); //do not use tree generated by this constructor its just to extract the treenames, its a chicken before the egg thing that needs to be revisited. Tree(int, TreeMap*); Tree(TreeMap*); //to generate a tree from a file + Tree(TreeMap*, vector< vector >&); //create tree from sim matrix ~Tree(); void getCopy(Tree*); //makes tree a copy of the one passed in. void getSubTree(Tree*, vector); //makes tree a that contains only the names passed in. + int getSubTree(Tree* originalToCopy, vector seqToInclude, map nameMap); //used with (int, TreeMap) constructor. SeqsToInclude contains subsample wanted - assumes these are unique seqs and size of vector=numLeaves passed into constructor. nameMap is unique -> redundantList can be empty if no namesfile was provided. + void assembleRandomTree(); void assembleRandomUnifracTree(vector); void assembleRandomUnifracTree(string, string); + void createNewickFile(string); int getIndex(string); void setIndex(string, int); @@ -54,7 +58,7 @@ private: map mergeGroups(int); //returns a map with a groupname and the number of times that group was seen in the children map mergeGcounts(int); - void addNamesToCounts(); + void addNamesToCounts(map); void randomTopology(); void randomBlengths(); void randomLabels(vector); diff --git a/treegroupscommand.cpp b/treegroupscommand.cpp index 0150a7a..4a77211 100644 --- a/treegroupscommand.cpp +++ b/treegroupscommand.cpp @@ -482,9 +482,9 @@ int TreeGroupCommand::execute(){ Tree* TreeGroupCommand::createTree(vector< vector >& simMatrix){ try { //create tree - t = new Tree(tmap); + t = new Tree(tmap, simMatrix); - //initialize index + /* //initialize index map index; //maps row in simMatrix to vector index in the tree for (int g = 0; g < numGroups; g++) { index[g] = g; } @@ -544,7 +544,7 @@ Tree* TreeGroupCommand::createTree(vector< vector >& simMatrix){ //adjust tree to make sure root to tip length is .5 int root = t->findRoot(); t->tree[root].setBranchLength((0.5 - t->tree[root].getLengthToLeaves())); - + */ //assemble tree t->assembleTree(); @@ -1004,11 +1004,13 @@ int TreeGroupCommand::process(vector thisLookup) { if (m->control_pressed) { for (int k = 0; k < trees.size(); k++) { delete trees[k]; } } Consensus consensus; + //clear old tree names if any + m->Treenames.clear(); m->Treenames = m->getGroups(); //may have changed if subsample eliminated groups Tree* conTree = consensus.getTree(trees, tmap); //create a new filename string conFile = outputDir + m->getRootName(m->getSimpleName(inputfile)) + treeCalculators[i]->getName() + "." + thisLookup[0]->getLabel() + ".cons.tre"; - outputNames.push_back(conFile); outputTypes["tree"].push_back(outputFile); + outputNames.push_back(conFile); outputTypes["tree"].push_back(conFile); ofstream outTree; m->openOutputFile(conFile, outTree); diff --git a/treemap.cpp b/treemap.cpp index 1fc5c01..450d8fb 100644 --- a/treemap.cpp +++ b/treemap.cpp @@ -229,6 +229,60 @@ void TreeMap::makeSim(ListVector* list) { exit(1); } } +/************************************************************/ +int TreeMap::getCopy(TreeMap* copy){ + try { + + namesOfGroups = copy->getNamesOfGroups(); + numGroups = copy->getNumGroups(); + namesOfSeqs = copy->namesOfSeqs; + seqsPerGroup = copy->seqsPerGroup; + treemap = copy->treemap; + + return 0; + } + catch(exception& e) { + m->errorOut(e, "TreeMap", "getCopy"); + exit(1); + } +} +/************************************************************/ +vector TreeMap::getNamesSeqs(){ + try { + + vector names; + + for(it = treemap.begin(); it != treemap.end(); it++){ + names.push_back(it->first); + } + + return names; + } + catch(exception& e) { + m->errorOut(e, "TreeMap", "getNamesSeqs"); + exit(1); + } +} +/************************************************************/ +vector TreeMap::getNamesSeqs(vector picked){ + try { + + vector names; + + for(it = treemap.begin(); it != treemap.end(); it++){ + //if you are belong to one the the groups in the picked vector add you + if (m->inUsersGroups(it->second.groupname, picked)) { + names.push_back(it->first); + } + } + + return names; + } + catch(exception& e) { + m->errorOut(e, "TreeMap", "getNamesSeqs"); + exit(1); + } +} /************************************************************/ diff --git a/treemap.h b/treemap.h index 7ed8d04..fc9c369 100644 --- a/treemap.h +++ b/treemap.h @@ -42,13 +42,19 @@ public: sort(namesOfGroups.begin(), namesOfGroups.end()); return namesOfGroups; } - vector namesOfSeqs; - map seqsPerGroup; //groupname, number of seqs in that group. - map treemap; //sequence name and - void print(ostream&); + + void print(ostream&); void makeSim(vector); //takes groupmap info and fills treemap for use by tree.shared command. void makeSim(ListVector*); //takes listvector info and fills treemap for use by tree.shared command. - + vector getNamesSeqs(); + vector getNamesSeqs(vector); //get names of seqs belonging to a group or set of groups + int getCopy(TreeMap*); + + vector namesOfSeqs; + map seqsPerGroup; //groupname, number of seqs in that group. + map treemap; //sequence name and + + private: vector namesOfGroups; ifstream fileHandle; diff --git a/trimseqscommand.cpp b/trimseqscommand.cpp index f00743c..b480702 100644 --- a/trimseqscommand.cpp +++ b/trimseqscommand.cpp @@ -437,6 +437,17 @@ int TrimSeqsCommand::execute(){ Sequence currSeq(in); m->gobble(in); out << currSeq.getName() << '\t' << it->second << endl; + + if (nameFile != "") { + map::iterator itName = nameMap.find(currSeq.getName()); + if (itName != nameMap.end()) { + vector thisSeqsNames; + m->splitAtChar(itName->second, thisSeqsNames, ','); + for (int k = 1; k < thisSeqsNames.size(); k++) { //start at 1 to skip self + out << thisSeqsNames[k] << '\t' << it->second << endl; + } + }else { m->mothurOut("[ERROR]: " + currSeq.getName() + " is not in your namefile, please correct."); m->mothurOutEndLine(); } + } } in.close(); out.close(); diff --git a/unifracweightedcommand.cpp b/unifracweightedcommand.cpp index b3a54c9..e596db2 100644 --- a/unifracweightedcommand.cpp +++ b/unifracweightedcommand.cpp @@ -8,6 +8,8 @@ */ #include "unifracweightedcommand.h" +#include "consensus.h" +#include "subsample.h" //********************************************************************************************************************** vector UnifracWeightedCommand::setParameters(){ @@ -18,7 +20,9 @@ vector UnifracWeightedCommand::setParameters(){ 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 psubsample("subsample", "String", "", "", "", "", "",false,false); parameters.push_back(psubsample); + CommandParameter pconsensus("consensus", "Boolean", "", "F", "", "", "",false,false); parameters.push_back(pconsensus); + 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); @@ -37,14 +41,16 @@ vector UnifracWeightedCommand::setParameters(){ 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 unifrac.weighted command parameters are tree, group, name, 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 unifrac.weighted command should be in the following format: unifrac.weighted(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, 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"; @@ -66,6 +72,7 @@ UnifracWeightedCommand::UnifracWeightedCommand(){ outputTypes["wsummary"] = tempOutNames; outputTypes["phylip"] = tempOutNames; outputTypes["column"] = tempOutNames; + outputTypes["tree"] = tempOutNames; } catch(exception& e) { m->errorOut(e, "UnifracWeightedCommand", "UnifracWeightedCommand"); @@ -102,6 +109,7 @@ UnifracWeightedCommand::UnifracWeightedCommand(string option) { outputTypes["wsummary"] = 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); @@ -159,7 +167,7 @@ UnifracWeightedCommand::UnifracWeightedCommand(string option) { 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 @@ -190,9 +198,25 @@ UnifracWeightedCommand::UnifracWeightedCommand(string option) { 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 (!random) { iters = 0; } //turn off random calcs - + 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 (namefile == "") { vector files; files.push_back(treefile); parser.getNameFile(files); @@ -214,7 +238,389 @@ int UnifracWeightedCommand::execute() { m->setTreeFile(treefile); - if (groupfile != "") { + readTrees(); if (m->control_pressed) { delete tmap; for (int i = 0; i < T.size(); i++) { delete T[i]; } return 0; } + + sumFile = outputDir + m->getSimpleName(treefile) + ".wsummary"; + 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 nameGroups = tmap->getNamesOfGroups(); + util.setGroups(Groups, nameGroups, s, numGroups, "weighted"); //sets the groups the user wants to analyze + m->setGroups(Groups); + + if (m->control_pressed) { delete tmap; for (int i = 0; i < T.size(); i++) { delete T[i]; } return 0; } + + Weighted weighted(tmap, 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 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 newGroups = Groups; + Groups.clear(); + for (int i = 0; i < newGroups.size(); i++) { + vector thisGroup; thisGroup.push_back(newGroups[i]); + vector 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); + } + } + + //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 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; + rScores.resize(numComp); //data[0] = weightedscore AB, data[1] = weightedscore AC... + uScores.resize(numComp); //data[0] = weightedscore AB, data[1] = weightedscore AC... + + vector userData; userData.resize(numComp,0); //weighted score info for user tree. data[0] = weightedscore AB, data[1] = weightedscore AC... + vector randomData; randomData.resize(numComp,0); //weighted score info for random trees. data[0] = weightedscore AB, data[1] = weightedscore AC... + + if (random) { + 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) { 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; } + + //save users score + for (int s=0; s > 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 tree needs it and memory cleanup becomes messy. + TreeMap* newTmap = new TreeMap(); + newTmap->getCopy(tmap); + + SubSample sample; + Tree* subSampleTree = sample.getSample(T[i], newTmap, nameMap, subsampleSize); + + //call new weighted function + vector iterData; iterData.resize(numComp,0); + Weighted thisWeighted(newTmap, includeRoot); + iterData = thisWeighted.getValues(subSampleTree, processors, outputDir); //userData[0] = weightedscore + + //save data to make ave dist, std dist + calcDistsTotals.push_back(iterData); + + delete newTmap; + delete subSampleTree; + } + + 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); } + } + + + 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; } + + if (phylip) { createPhylipFile(); } + + printWSummaryFile(); + + //clear out users groups + m->clearGroups(); + 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++) { m->mothurRemove(outputNames[i]); } 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); } + } + + //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(); + for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); } + m->mothurOutEndLine(); + + return 0; + + } + catch(exception& e) { + m->errorOut(e, "UnifracWeightedCommand", "execute"); + exit(1); + } +} +/**************************************************************************************************/ +int UnifracWeightedCommand::getAverageSTDMatrices(vector< vector >& dists, int treeNum) { + try { + //we need to find the average distance and standard deviation for each groups distance + + //finds sum + vector 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 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 > 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 > 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; rgetNumGroups(); 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) + ".weighted.ave.dist"; + outputNames.push_back(aveFileName); outputTypes["phylip"].push_back(aveFileName); + + ofstream out; + m->openOutputFile(aveFileName, out); + + string stdFileName = outputDir + m->getSimpleName(treefile) + toString(treeNum+1) + ".weighted.std.dist"; + outputNames.push_back(stdFileName); outputTypes["phylip"].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; rgetNumGroups(); 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 >& dists, int treeNum) { + try { + + //used in tree constructor + m->runParse = false; + + //create treemap class from groupmap for tree class to use + TreeMap* newTmap = new TreeMap(); + newTmap->makeSim(m->getGroups()); + + //clear old tree names if any + m->Treenames.clear(); + + //fills globaldatas tree names + m->Treenames = m->getGroups(); + + vector newTrees = buildTrees(dists, treeNum, newTmap); //also creates .all.tre file containing the trees created + + if (m->control_pressed) { delete newTmap; return 0; } + + Consensus con; + Tree* conTree = con.getTree(newTrees, newTmap); + + //create a new filename + string conFile = outputDir + m->getRootName(m->getSimpleName(treefile)) + toString(treeNum+1) + ".weighted.cons.tre"; + 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(); + delete newTmap; + + return 0; + + } + catch(exception& e) { + m->errorOut(e, "UnifracWeightedCommand", "getConsensusTrees"); + exit(1); + } +} +/**************************************************************************************************/ + +vector UnifracWeightedCommand::buildTrees(vector< vector >& dists, int treeNum, TreeMap* mytmap) { + try { + + vector trees; + + //create a new filename + string outputFile = outputDir + m->getRootName(m->getSimpleName(treefile)) + toString(treeNum+1) + ".weighted.all.tre"; + 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 > 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; rgetNumGroups(); 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); + 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::readTrees() { + try { + + if (groupfile != "") { //read in group map info. tmap = new TreeMap(groupfile); tmap->readMap(); @@ -275,183 +681,91 @@ int UnifracWeightedCommand::execute() { } } } - - 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 - Groups = m->getGroups(); - vector nameGroups = tmap->getNamesOfGroups(); - util->setGroups(Groups, nameGroups, s, numGroups, "weighted"); //sets the groups the user wants to analyze - util->getCombos(groupComb, Groups, numComp); - m->setGroups(Groups); - delete util; - - weighted = new Weighted(tmap, includeRoot); - - int start = time(NULL); - - //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... - - if (numComp < processors) { processors = numComp; } - - //get weighted scores for users trees - for (int i = 0; i < T.size(); i++) { - - 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++) { 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... - - if (random) { - 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) { 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++) { m->mothurRemove(outputNames[i]); } return 0; } - - //save users score - for (int s=0; s > namesOfGroupCombos; - for (int a=0; a groups; groups.push_back((m->getGroups())[a]); groups.push_back((m->getGroups())[l]); - namesOfGroupCombos.push_back(groups); - } - } - - lines.clear(); - - #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix) - 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) || (__linux__) || (__unix__) || (__unix) - if(processors == 1){ - driver(T[i], namesOfGroupCombos, 0, namesOfGroupCombos.size(), rScores); - }else{ - createProcesses(T[i], namesOfGroupCombos, rScores); - } - #else - driver(T[i], namesOfGroupCombos, 0, namesOfGroupCombos.size(), rScores); - #endif - - 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++) { m->mothurRemove(outputNames[i]); } 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(); - } - - - 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++) { m->mothurRemove(outputNames[i]); } return 0; } - - printWSummaryFile(); - - if (phylip) { createPhylipFile(); } - - //clear out users groups - m->clearGroups(); - delete tmap; delete weighted; - 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; - } - - 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); } - } - - //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(); - for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); } - m->mothurOutEndLine(); - - return 0; - + return 0; + } + catch(exception& e) { + m->errorOut(e, "UnifracWeightedCommand", "readTrees"); + exit(1); } +} +/**************************************************************************************************/ + +int UnifracWeightedCommand::runRandomCalcs(Tree* thisTree, vector usersScores) { + try { + + //calculate number of comparisons i.e. with groups A,B,C = AB, AC, BC = 3; + vector< vector > namesOfGroupCombos; + for (int a=0; a groups; groups.push_back((m->getGroups())[a]); groups.push_back((m->getGroups())[l]); + namesOfGroupCombos.push_back(groups); + } + } + + lines.clear(); + +#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix) + 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) || (__linux__) || (__unix__) || (__unix) + if(processors == 1){ + driver(thisTree, namesOfGroupCombos, 0, namesOfGroupCombos.size(), rScores); + }else{ + createProcesses(thisTree, namesOfGroupCombos, rScores); + } +#else + driver(T[i], namesOfGroupCombos, 0, namesOfGroupCombos.size(), rScores); +#endif + + if (m->control_pressed) { delete tmap; 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; } + + //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(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", "execute"); + m->errorOut(e, "UnifracWeightedCommand", "runRandomCalcs"); exit(1); } } @@ -524,7 +838,9 @@ int UnifracWeightedCommand::createProcesses(Tree* t, vector< vector > na int UnifracWeightedCommand::driver(Tree* t, vector< vector > namesOfGroupCombos, int start, int num, vector< vector >& scores) { try { Tree* randT = new Tree(tmap); - + + Weighted weighted(tmap, includeRoot); + for (int h = start; h < (start+num); h++) { if (m->control_pressed) { return 0; } @@ -542,7 +858,7 @@ int UnifracWeightedCommand::driver(Tree* t, vector< vector > namesOfGrou if (m->control_pressed) { delete randT; return 0; } //get wscore of random tree - EstOutput randomData = weighted->getValues(randT, groupA, groupB); + EstOutput randomData = weighted.getValues(randT, groupA, groupB); if (m->control_pressed) { delete randT; return 0; } @@ -789,7 +1105,7 @@ int UnifracWeightedCommand::readNamesFile() { m->splitAtComma(second, dupNames); for (int i = 0; i < dupNames.size(); i++) { - nameMap[dupNames[i]] = dupNames[i]; + nameMap[dupNames[i]] = first; 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; } diff --git a/unifracweightedcommand.h b/unifracweightedcommand.h index b1db317..e1ebd16 100644 --- a/unifracweightedcommand.h +++ b/unifracweightedcommand.h @@ -51,21 +51,18 @@ class UnifracWeightedCommand : public Command { vector WScoreSig; //tree weighted score signifigance when compared to random trees - percentage of random trees with that score or lower. vector groupComb; // AB. AC, BC... TreeMap* tmap; - Weighted* weighted; string sumFile, outputDir; int iters, numGroups, numComp, counter; - EstOutput userData; //weighted score info for user tree - EstOutput randomData; //weighted score info for random trees vector< vector > rScores; //vector each group comb has an entry vector< vector > uScores; //vector each group comb has an entry vector< map > rScoreFreq; //map -vector entry for each combination. vector< map > rCumul; //map -vector entry for each c map validScores; //map contains scores from random - bool abort, phylip, random, includeRoot; + bool abort, phylip, random, includeRoot, subsample, consensus; string groups, itersString, outputForm, treefile, groupfile, namefile; vector Groups, outputNames; //holds groups to be used - int processors, numUniquesInName; + int processors, numUniquesInName, subsampleSize, subsampleIters; ofstream outSum; map nameMap; @@ -78,6 +75,11 @@ class UnifracWeightedCommand : public Command { int createProcesses(Tree*, vector< vector >, vector< vector >&); int driver(Tree*, vector< vector >, int, int, vector< vector >&); int readNamesFile(); + int runRandomCalcs(Tree*, vector); + int readTrees(); + vector buildTrees(vector< vector >&, int, TreeMap*); + int getConsensusTrees(vector< vector >&, int); + int getAverageSTDMatrices(vector< vector >&, int); }; diff --git a/weighted.cpp b/weighted.cpp index 7a31da4..b0a642b 100644 --- a/weighted.cpp +++ b/weighted.cpp @@ -87,7 +87,7 @@ EstOutput Weighted::createProcesses(Tree* t, vector< vector > namesOfGro EstOutput Myresults; Myresults = driver(t, namesOfGroupCombos, lines[process].start, lines[process].num); - m->mothurOut("Merging results."); m->mothurOutEndLine(); + //m->mothurOut("Merging results."); m->mothurOutEndLine(); //pass numSeqs to parent ofstream out; @@ -142,7 +142,7 @@ EstOutput Weighted::createProcesses(Tree* t, vector< vector > namesOfGro m->mothurRemove(s); } - m->mothurOut("DONE."); m->mothurOutEndLine(); m->mothurOutEndLine(); + //m->mothurOut("DONE."); m->mothurOutEndLine(); m->mothurOutEndLine(); return results; #endif -- 2.39.2