X-Git-Url: https://git.donarmstrong.com/?p=mothur.git;a=blobdiff_plain;f=phylodiversitycommand.cpp;h=339649c7dcc814ac409e5b8ec1dfd432b544e43d;hp=d15899c2f2f853ba9529c9b65f74ce67885c8a63;hb=a8e2df1b96a57f5f29576b08361b86a96a8eff4f;hpb=0e051b4cfda410b0d441da6ff2f96d4bbe1d9e5a diff --git a/phylodiversitycommand.cpp b/phylodiversitycommand.cpp index d15899c..339649c 100644 --- a/phylodiversitycommand.cpp +++ b/phylodiversitycommand.cpp @@ -8,24 +8,26 @@ */ #include "phylodiversitycommand.h" +#include "treereader.h" //********************************************************************************************************************** vector PhyloDiversityCommand::setParameters(){ try { - CommandParameter ptree("tree", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(ptree); - CommandParameter pgroup("group", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(pgroup); - CommandParameter pname("name", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(pname); - CommandParameter pgroups("groups", "String", "", "", "", "", "",false,false); parameters.push_back(pgroups); - CommandParameter piters("iters", "Number", "", "1000", "", "", "",false,false); parameters.push_back(piters); - CommandParameter pfreq("freq", "Number", "", "100", "", "", "",false,false); parameters.push_back(pfreq); - CommandParameter pprocessors("processors", "Number", "", "1", "", "", "",false,false); parameters.push_back(pprocessors); - CommandParameter prarefy("rarefy", "Boolean", "", "F", "", "", "",false,false); parameters.push_back(prarefy); - CommandParameter psummary("summary", "Boolean", "", "T", "", "", "",false,false); parameters.push_back(psummary); - CommandParameter pcollect("collect", "Boolean", "", "F", "", "", "",false,false); parameters.push_back(pcollect); - CommandParameter pscale("scale", "Boolean", "", "F", "", "", "",false,false); parameters.push_back(pscale); - CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir); - CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir); + CommandParameter ptree("tree", "InputTypes", "", "", "none", "none", "none","phylodiv",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 pfreq("freq", "Number", "", "100", "", "", "","",false,false); parameters.push_back(pfreq); + CommandParameter pprocessors("processors", "Number", "", "1", "", "", "","",false,false,true); parameters.push_back(pprocessors); + CommandParameter prarefy("rarefy", "Boolean", "", "F", "", "", "","rarefy",false,false); parameters.push_back(prarefy); + CommandParameter psummary("summary", "Boolean", "", "T", "", "", "","summary",false,false); parameters.push_back(psummary); + CommandParameter pcollect("collect", "Boolean", "", "F", "", "", "","",false,false); parameters.push_back(pcollect); + CommandParameter pscale("scale", "Boolean", "", "F", "", "", "","",false,false); parameters.push_back(pscale); + CommandParameter pinputdir("inputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(pinputdir); + CommandParameter poutputdir("outputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(poutputdir); vector myArray; for (int i = 0; i < parameters.size(); i++) { myArray.push_back(parameters[i].name); } @@ -40,7 +42,7 @@ vector PhyloDiversityCommand::setParameters(){ string PhyloDiversityCommand::getHelpString(){ try { string helpString = ""; - helpString += "The phylo.diversity command parameters are tree, group, name, groups, iters, freq, processors, scale, rarefy, collect and summary. tree and group are required, unless you have valid current files.\n"; + helpString += "The phylo.diversity command parameters are tree, group, name, count, groups, iters, freq, processors, scale, rarefy, collect and summary. tree and group are required, unless you have valid current files.\n"; helpString += "The groups parameter allows you to specify which of the groups in your groupfile you would like analyzed. The group names are separated by dashes. By default all groups are used.\n"; helpString += "The iters parameter allows you to specify the number of randomizations to preform, by default iters=1000, if you set rarefy to true.\n"; helpString += "The freq parameter is used indicate when to output your data, by default it is set to 100. But you can set it to a percentage of the number of sequence. For example freq=0.10, means 10%. \n"; @@ -60,7 +62,23 @@ string PhyloDiversityCommand::getHelpString(){ exit(1); } } - +//********************************************************************************************************************** +string PhyloDiversityCommand::getOutputPattern(string type) { + try { + string pattern = ""; + + if (type == "phylodiv") { pattern = "[filename],[tag],phylodiv"; } + else if (type == "rarefy") { pattern = "[filename],[tag],phylodiv.rarefaction"; } + else if (type == "summary") { pattern = "[filename],[tag],phylodiv.summary"; } + else { m->mothurOut("[ERROR]: No definition for type " + type + " output pattern.\n"); m->control_pressed = true; } + + return pattern; + } + catch(exception& e) { + m->errorOut(e, "PhyloDiversityCommand", "getOutputPattern"); + exit(1); + } +} //********************************************************************************************************************** PhyloDiversityCommand::PhyloDiversityCommand(){ @@ -84,6 +102,7 @@ PhyloDiversityCommand::PhyloDiversityCommand(string option) { //allow user to run help if(option == "help") { help(); abort = true; calledHelp = true; } + else if(option == "citation") { citation(); abort = true; calledHelp = true;} else { vector myArray = setParameters();; @@ -133,46 +152,58 @@ PhyloDiversityCommand::PhyloDiversityCommand(string option) { //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; } + } } - m->runParse = true; - m->Groups.clear(); - m->namesOfGroups.clear(); - m->Treenames.clear(); - m->names.clear(); - //check for required parameters treefile = validParameter.validFile(parameters, "tree", true); - if (treefile == "not open") { abort = 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") { - //if there is a current design file, use it - groupfile = m->getGroupFile(); - if (groupfile != "") { m->mothurOut("Using " + groupfile + " as input file for the group parameter."); m->mothurOutEndLine(); } - else { m->mothurOut("You have no current group file and the group parameter is required."); m->mothurOutEndLine(); abort = true; } - } + if (groupfile == "not open") { groupfile = ""; 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); } + + 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); } string temp; temp = validParameter.validFile(parameters, "freq", false); if (temp == "not found") { temp = "100"; } - convert(temp, freq); + m->mothurConvert(temp, freq); temp = validParameter.validFile(parameters, "iters", false); if (temp == "not found") { temp = "1000"; } - convert(temp, iters); + m->mothurConvert(temp, iters); temp = validParameter.validFile(parameters, "rarefy", false); if (temp == "not found") { temp = "F"; } rarefy = m->isTrue(temp); @@ -189,16 +220,23 @@ PhyloDiversityCommand::PhyloDiversityCommand(string option) { temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); } m->setProcessors(temp); - convert(temp, processors); + m->mothurConvert(temp, processors); groups = validParameter.validFile(parameters, "groups", false); if (groups == "not found") { groups = ""; } else { m->splitAtDash(groups, Groups); - m->Groups = Groups; + m->setGroups(Groups); } if ((!collect) && (!rarefy) && (!summary)) { m->mothurOut("No outputs selected. You must set either collect, rarefy or summary to true, summary=T by default."); m->mothurOutEndLine(); abort=true; } + + if (countfile=="") { + if (namefile == "") { + vector files; files.push_back(treefile); + parser.getNameFile(files); + } + } } } @@ -214,81 +252,39 @@ int PhyloDiversityCommand::execute(){ if (abort == true) { if (calledHelp) { return 0; } return 2; } + int start = time(NULL); + m->setTreeFile(treefile); - - //read in group map info. - tmap = new TreeMap(groupfile); - tmap->readMap(); - - 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(); - vector trees = 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 < trees.size(); i++) { delete trees[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::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 - } - } - } - } - - SharedUtil* util = new SharedUtil(); - util->setGroups(m->Groups, tmap->namesOfGroups, "treegroup"); //sets the groups the user wants to analyze - delete util; + TreeReader* reader; + if (countfile == "") { reader = new TreeReader(treefile, groupfile, namefile); } + else { reader = new TreeReader(treefile, countfile); } + vector trees = reader->getTrees(); + ct = trees[0]->getCountTable(); + delete reader; + + SharedUtil util; + vector mGroups = m->getGroups(); + vector tGroups = ct->getNamesOfGroups(); + util.setGroups(mGroups, tGroups, "phylo.diversity"); //sets the groups the user wants to analyze //incase the user had some mismatches between the tree and group files we don't want group xxx to be analyzed - for (int i = 0; i < m->Groups.size(); i++) { if (m->Groups[i] == "xxx") { m->Groups.erase(m->Groups.begin()+i); break; } } + for (int i = 0; i < mGroups.size(); i++) { if (mGroups[i] == "xxx") { mGroups.erase(mGroups.begin()+i); break; } } + m->setGroups(mGroups); vector outputNames; //for each of the users trees for(int i = 0; i < trees.size(); i++) { - if (m->control_pressed) { delete tmap; for (int j = 0; j < trees.size(); j++) { delete trees[j]; } for (int j = 0; j < outputNames.size(); j++) { remove(outputNames[j].c_str()); } return 0; } + if (m->control_pressed) { delete ct; for (int j = 0; j < trees.size(); j++) { delete trees[j]; } for (int j = 0; j < outputNames.size(); j++) { m->mothurRemove(outputNames[j]); } return 0; } ofstream outSum, outRare, outCollect; - string outSumFile = outputDir + m->getRootName(m->getSimpleName(treefile)) + toString(i+1) + ".phylodiv.summary"; - string outRareFile = outputDir + m->getRootName(m->getSimpleName(treefile)) + toString(i+1) + ".phylodiv.rarefaction"; - string outCollectFile = outputDir + m->getRootName(m->getSimpleName(treefile)) + toString(i+1) + ".phylodiv"; + map variables; + variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(treefile)); + variables["[tag]"] = toString(i+1); + string outSumFile = getOutputFileName("summary",variables); + string outRareFile = getOutputFileName("rarefy",variables); + string outCollectFile = getOutputFileName("phylodiv",variables); if (summary) { m->openOutputFile(outSumFile, outSum); outputNames.push_back(outSumFile); outputTypes["summary"].push_back(outSumFile); } if (rarefy) { m->openOutputFile(outRareFile, outRare); outputNames.push_back(outRareFile); outputTypes["rarefy"].push_back(outRareFile); } @@ -299,10 +295,35 @@ int PhyloDiversityCommand::execute(){ //create a vector containing indexes of leaf nodes, randomize it, select nodes to send to calculator vector randomLeaf; for (int j = 0; j < numLeafNodes; j++) { - if (m->inUsersGroups(trees[i]->tree[j].getGroup(), m->Groups) == true) { //is this a node from the group the user selected. + if (m->inUsersGroups(trees[i]->tree[j].getGroup(), mGroups) == true) { //is this a node from the group the user selected. randomLeaf.push_back(j); } } + + /* float sum = 0; + vector sums; sums.resize(m->getGroups().size(), 0); + vector sumsAboveRoot; sumsAboveRoot.resize(m->getGroups().size(), 0); + for (int j = 0; j < trees[i]->getNumNodes(); j++) { + if (trees[i]->tree[j].getBranchLength() < 0) { cout << j << '\t' << trees[i]->tree[j].getName() << '\t' << trees[i]->tree[j].getBranchLength() << endl; } + + sum += abs(trees[i]->tree[j].getBranchLength()); + for (int k = 0; k < m->getGroups().size(); k++) { + map::iterator itGroup = trees[i]->tree[j].pcount.find(m->getGroups()[k]); + if (itGroup != trees[i]->tree[j].pcount.end()) { //this branch belongs to a group we care about + if (j < rootForGroup[m->getGroups()[k]]) { + sums[k] += abs(trees[i]->tree[j].getBranchLength()); + }else { + sumsAboveRoot[k] += abs(trees[i]->tree[j].getBranchLength()); + } + } + } + } + cout << sum << endl; //exit(1); + + for (int k = 0; k < m->getGroups().size(); k++) { + cout << m->getGroups()[k] << "root node = " << rootForGroup[m->getGroups()[k]] << "sum below root = " << sums[k] << "sum above root = " << sumsAboveRoot[k] << endl; + } + exit(1); */ numLeafNodes = randomLeaf.size(); //reset the number of leaf nodes you are using @@ -314,15 +335,16 @@ int PhyloDiversityCommand::execute(){ //find largest group total int largestGroup = 0; - for (int j = 0; j < m->Groups.size(); j++) { - if (tmap->seqsPerGroup[m->Groups[j]] > largestGroup) { largestGroup = tmap->seqsPerGroup[m->Groups[j]]; } + for (int j = 0; j < mGroups.size(); j++) { + int numSeqsThisGroup = ct->getGroupCount(mGroups[j]); + if (numSeqsThisGroup > largestGroup) { largestGroup = numSeqsThisGroup; } //initialize diversity - diversity[m->Groups[j]].resize(tmap->seqsPerGroup[m->Groups[j]]+1, 0.0); //numSampled + diversity[mGroups[j]].resize(numSeqsThisGroup+1, 0.0); //numSampled //groupA 0.0 0.0 //initialize sumDiversity - sumDiversity[m->Groups[j]].resize(tmap->seqsPerGroup[m->Groups[j]]+1, 0.0); + sumDiversity[mGroups[j]].resize(numSeqsThisGroup+1, 0.0); } //convert freq percentage to number @@ -336,44 +358,35 @@ int PhyloDiversityCommand::execute(){ if(largestGroup % increment != 0){ numSampledList.insert(largestGroup); } //add other groups ending points - for (int j = 0; j < m->Groups.size(); j++) { - if (numSampledList.count(diversity[m->Groups[j]].size()-1) == 0) { numSampledList.insert(diversity[m->Groups[j]].size()-1); } + for (int j = 0; j < mGroups.size(); j++) { + if (numSampledList.count(diversity[mGroups[j]].size()-1) == 0) { numSampledList.insert(diversity[mGroups[j]].size()-1); } } - #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) - if(processors == 1){ - driver(trees[i], diversity, sumDiversity, iters, increment, randomLeaf, numSampledList, outCollect, outSum, true); - }else{ - if (rarefy) { - vector procIters; - - int numItersPerProcessor = iters / processors; - - //divide iters between processes - for (int h = 0; h < processors; h++) { - if(h == processors - 1){ - numItersPerProcessor = iters - h * numItersPerProcessor; - } - procIters.push_back(numItersPerProcessor); - } - - createProcesses(procIters, trees[i], diversity, sumDiversity, iters, increment, randomLeaf, numSampledList, outCollect, outSum); - - }else{ //no need to paralellize if you dont want to rarefy - driver(trees[i], diversity, sumDiversity, iters, increment, randomLeaf, numSampledList, outCollect, outSum, true); - } - } - - #else - driver(trees[i], diversity, sumDiversity, iters, increment, randomLeaf, numSampledList, outCollect, outSum, true); - #endif - + if (rarefy) { + vector procIters; + int numItersPerProcessor = iters / processors; + + //divide iters between processes + for (int h = 0; h < processors; h++) { + if(h == processors - 1){ numItersPerProcessor = iters - h * numItersPerProcessor; } + procIters.push_back(numItersPerProcessor); + } + + createProcesses(procIters, trees[i], diversity, sumDiversity, iters, increment, randomLeaf, numSampledList, outCollect, outSum); + + }else{ //no need to paralellize if you dont want to rarefy + driver(trees[i], diversity, sumDiversity, iters, increment, randomLeaf, numSampledList, outCollect, outSum, true); + } + if (rarefy) { printData(numSampledList, sumDiversity, outRare, iters); } } - 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 phylo.diversity."); m->mothurOutEndLine(); + + m->mothurOutEndLine(); m->mothurOut("Output File Names: "); m->mothurOutEndLine(); for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); } @@ -390,12 +403,13 @@ int PhyloDiversityCommand::execute(){ //********************************************************************************************************************** int PhyloDiversityCommand::createProcesses(vector& procIters, Tree* t, map< string, vector >& div, map >& sumDiv, int numIters, int increment, vector& randomLeaf, set& numSampledList, ofstream& outCollect, ofstream& outSum){ try { - #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) - int process = 1; + int process = 1; vector processIDS; map< string, vector >::iterator itSum; - + + #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix) + //loop through and create all the processes you want while (process != processors) { int pid = fork(); @@ -462,7 +476,57 @@ int PhyloDiversityCommand::createProcesses(vector& procIters, Tree* t, map< } in.close(); - remove(inTemp.c_str()); + m->mothurRemove(inTemp); + } +#else + + //fill in functions + vector pDataArray; + DWORD dwThreadIdArray[processors-1]; + HANDLE hThreadArray[processors-1]; + vector cts; + vector trees; + map rootForGroup = getRootForGroups(t); + + //Create processor worker threads. + for( int i=1; icopy(ct); + Tree* copyTree = new Tree(copyCount); + copyTree->getCopy(t); + + cts.push_back(copyCount); + trees.push_back(copyTree); + + map > copydiv = div; + map > copysumDiv = sumDiv; + vector copyrandomLeaf = randomLeaf; + set copynumSampledList = numSampledList; + map copyRootForGrouping = rootForGroup; + + phylodivData* temp = new phylodivData(m, procIters[i], copydiv, copysumDiv, copyTree, copyCount, increment, copyrandomLeaf, copynumSampledList, copyRootForGrouping); + pDataArray.push_back(temp); + processIDS.push_back(i); + + hThreadArray[i-1] = CreateThread(NULL, 0, MyPhyloDivThreadFunction, pDataArray[i-1], 0, &dwThreadIdArray[i-1]); + } + + driver(t, div, sumDiv, procIters[0], increment, randomLeaf, numSampledList, outCollect, outSum, true); + + //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 (itSum = pDataArray[i]->sumDiv.begin(); itSum != pDataArray[i]->sumDiv.end(); itSum++) { + for (int k = 0; k < (itSum->second).size(); k++) { + sumDiv[itSum->first][k] += (itSum->second)[k]; + } + } + delete cts[i]; + delete trees[i]; + CloseHandle(hThreadArray[i]); + delete pDataArray[i]; } #endif @@ -479,47 +543,59 @@ int PhyloDiversityCommand::createProcesses(vector& procIters, Tree* t, map< int PhyloDiversityCommand::driver(Tree* t, map< string, vector >& div, map >& sumDiv, int numIters, int increment, vector& randomLeaf, set& numSampledList, ofstream& outCollect, ofstream& outSum, bool doSumCollect){ try { int numLeafNodes = randomLeaf.size(); - + vector mGroups = m->getGroups(); + + map rootForGroup = getRootForGroups(t); //maps groupName to root node in tree. "root" for group may not be the trees root and we don't want to include the extra branches. + for (int l = 0; l < numIters; l++) { random_shuffle(randomLeaf.begin(), randomLeaf.end()); - + //initialize counts map counts; - map< string, set > countedBranch; - for (int j = 0; j < m->Groups.size(); j++) { counts[m->Groups[j]] = 0; countedBranch[m->Groups[j]].insert(-2); } //add dummy index to initialize countedBranch sets + vector< map > countedBranch; + for (int i = 0; i < t->getNumNodes(); i++) { + map temp; + for (int j = 0; j < mGroups.size(); j++) { temp[mGroups[j]] = false; } + countedBranch.push_back(temp); + } + + for (int j = 0; j < mGroups.size(); j++) { counts[mGroups[j]] = 0; } for(int k = 0; k < numLeafNodes; k++){ if (m->control_pressed) { return 0; } //calc branch length of randomLeaf k - vector br = calcBranchLength(t, randomLeaf[k], countedBranch); + vector br = calcBranchLength(t, randomLeaf[k], countedBranch, rootForGroup); //for each group in the groups update the total branch length accounting for the names file vector groups = t->tree[randomLeaf[k]].getGroup(); for (int j = 0; j < groups.size(); j++) { - int numSeqsInGroupJ = 0; - map::iterator it; - it = t->tree[randomLeaf[k]].pcount.find(groups[j]); - if (it != t->tree[randomLeaf[k]].pcount.end()) { //this leaf node contains seqs from group j - numSeqsInGroupJ = it->second; - } - - if (numSeqsInGroupJ != 0) { div[groups[j]][(counts[groups[j]]+1)] = div[groups[j]][counts[groups[j]]] + br[j]; } - - for (int s = (counts[groups[j]]+2); s <= (counts[groups[j]]+numSeqsInGroupJ); s++) { - div[groups[j]][s] = div[groups[j]][s-1]; //update counts, but don't add in redundant branch lengths - } - counts[groups[j]] += numSeqsInGroupJ; + + if (m->inUsersGroups(groups[j], mGroups)) { + int numSeqsInGroupJ = 0; + map::iterator it; + it = t->tree[randomLeaf[k]].pcount.find(groups[j]); + if (it != t->tree[randomLeaf[k]].pcount.end()) { //this leaf node contains seqs from group j + numSeqsInGroupJ = it->second; + } + + if (numSeqsInGroupJ != 0) { div[groups[j]][(counts[groups[j]]+1)] = div[groups[j]][counts[groups[j]]] + br[j]; } + + for (int s = (counts[groups[j]]+2); s <= (counts[groups[j]]+numSeqsInGroupJ); s++) { + div[groups[j]][s] = div[groups[j]][s-1]; //update counts, but don't add in redundant branch lengths + } + counts[groups[j]] += numSeqsInGroupJ; + } } } if (rarefy) { //add this diversity to the sum - for (int j = 0; j < m->Groups.size(); j++) { - for (int g = 0; g < div[m->Groups[j]].size(); g++) { - sumDiv[m->Groups[j]][g] += div[m->Groups[j]][g]; + for (int j = 0; j < mGroups.size(); j++) { + for (int g = 0; g < div[mGroups[j]].size(); g++) { + sumDiv[mGroups[j]][g] += div[mGroups[j]][g]; } } } @@ -545,17 +621,19 @@ void PhyloDiversityCommand::printSumData(map< string, vector >& div, ofst out << "Groups\tnumSampled\tphyloDiversity" << endl; out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint); - - for (int j = 0; j < m->Groups.size(); j++) { - int numSampled = (div[m->Groups[j]].size()-1); - out << m->Groups[j] << '\t' << numSampled << '\t'; + + vector mGroups = m->getGroups(); + for (int j = 0; j < mGroups.size(); j++) { + int numSampled = (div[mGroups[j]].size()-1); + out << mGroups[j] << '\t' << numSampled << '\t'; float score; - if (scale) { score = (div[m->Groups[j]][numSampled] / (float)numIters) / (float)numSampled; } - else { score = div[m->Groups[j]][numSampled] / (float)numIters; } + if (scale) { score = (div[mGroups[j]][numSampled] / (float)numIters) / (float)numSampled; } + else { score = div[mGroups[j]][numSampled] / (float)numIters; } out << setprecision(4) << score << endl; + //cout << mGroups[j] << '\t' << numSampled << '\t'<< setprecision(4) << score << endl; } out.close(); @@ -572,7 +650,8 @@ void PhyloDiversityCommand::printData(set& num, map< string, vector try { out << "numSampled\t"; - for (int i = 0; i < m->Groups.size(); i++) { out << m->Groups[i] << '\t'; } + vector mGroups = m->getGroups(); + for (int i = 0; i < mGroups.size(); i++) { out << mGroups[i] << '\t'; } out << endl; out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint); @@ -581,12 +660,12 @@ void PhyloDiversityCommand::printData(set& num, map< string, vector int numSampled = *it; out << numSampled << '\t'; - - for (int j = 0; j < m->Groups.size(); j++) { - if (numSampled < div[m->Groups[j]].size()) { + + for (int j = 0; j < mGroups.size(); j++) { + if (numSampled < div[mGroups[j]].size()) { float score; - if (scale) { score = (div[m->Groups[j]][numSampled] / (float)numIters) / (float)numSampled; } - else { score = div[m->Groups[j]][numSampled] / (float)numIters; } + if (scale) { score = (div[mGroups[j]][numSampled] / (float)numIters) / (float)numSampled; } + else { score = div[mGroups[j]][numSampled] / (float)numIters; } out << setprecision(4) << score << '\t'; }else { out << "NA" << '\t'; } @@ -604,9 +683,9 @@ void PhyloDiversityCommand::printData(set& num, map< string, vector } //********************************************************************************************************************** //need a vector of floats one branch length for every group the node represents. -vector PhyloDiversityCommand::calcBranchLength(Tree* t, int leaf, map< string, set >& counted){ +vector PhyloDiversityCommand::calcBranchLength(Tree* t, int leaf, vector< map >& counted, map roots){ try { - + //calc the branch length //while you aren't at root vector sums; @@ -615,127 +694,103 @@ vector PhyloDiversityCommand::calcBranchLength(Tree* t, int leaf, map< st vector groups = t->tree[leaf].getGroup(); sums.resize(groups.size(), 0.0); - map > tempTotals; //maps node to total Branch Length - map< string, set > tempCounted; - set::iterator it; - - //you are a leaf + + //you are a leaf if(t->tree[index].getBranchLength() != -1){ for (int k = 0; k < groups.size(); k++) { - sums[k] += abs(t->tree[index].getBranchLength()); - counted[groups[k]].insert(index); + sums[k] += abs(t->tree[index].getBranchLength()); } } - - for (int k = 0; k < groups.size(); k++) { - tempTotals[groups[k]][index] = 0.0; - } - - index = t->tree[index].getParent(); - + + + index = t->tree[index].getParent(); + //while you aren't at root while(t->tree[index].getParent() != -1){ - + if (m->control_pressed) { return sums; } - int pcountSize = 0; for (int k = 0; k < groups.size(); k++) { - map::iterator itGroup = t->tree[index].pcount.find(groups[k]); - if (itGroup != t->tree[index].pcount.end()) { pcountSize++; } - - //do both your chidren have have descendants from the users groups? - int lc = t->tree[index].getLChild(); - int rc = t->tree[index].getRChild(); - - int LpcountSize = 0; - itGroup = t->tree[lc].pcount.find(groups[k]); - if (itGroup != t->tree[lc].pcount.end()) { LpcountSize++; } - - int RpcountSize = 0; - itGroup = t->tree[rc].pcount.find(groups[k]); - if (itGroup != t->tree[rc].pcount.end()) { RpcountSize++; } - - //if yes, add your childrens tempTotals - if ((LpcountSize != 0) && (RpcountSize != 0)) { - sums[k] += tempTotals[groups[k]][lc] + tempTotals[groups[k]][rc]; - - for (it = tempCounted[groups[k]].begin(); it != tempCounted[groups[k]].end(); it++) { counted[groups[k]].insert(*it); } - - //cout << "added to total " << tempTotals[lc] << '\t' << tempTotals[rc] << endl; - if (t->tree[index].getBranchLength() != -1) { - if (counted[groups[k]].count(index) == 0) { - tempTotals[groups[k]][index] = abs(t->tree[index].getBranchLength()); - tempCounted[groups[k]].insert(index); - }else{ - tempTotals[groups[k]][index] = 0.0; - } - }else { - tempTotals[groups[k]][index] = 0.0; - } - }else { //if no, your tempTotal is your childrens temp totals + your branch length - tempTotals[groups[k]][index] = tempTotals[groups[k]][lc] + tempTotals[groups[k]][rc]; - - if (counted[groups[k]].count(index) == 0) { - tempTotals[groups[k]][index] += abs(t->tree[index].getBranchLength()); - tempCounted[groups[k]].insert(index); - } - - } - //cout << "temptotal = "<< tempTotals[i] << endl; - } - - index = t->tree[index].getParent(); - } - + + if (index >= roots[groups[k]]) { counted[index][groups[k]] = true; } //if you are at this groups "root", then say we are done + + if (!counted[index][groups[k]]){ //if counted[index][groups[k] is true this groups has already added all br from here to root, so quit early + if (t->tree[index].getBranchLength() != -1) { + sums[k] += abs(t->tree[index].getBranchLength()); + } + counted[index][groups[k]] = true; + } + } + index = t->tree[index].getParent(); + } + return sums; - + } catch(exception& e) { m->errorOut(e, "PhyloDiversityCommand", "calcBranchLength"); exit(1); } } -/*****************************************************************/ -int PhyloDiversityCommand::readNamesFile() { +//********************************************************************************************************************** +map PhyloDiversityCommand::getRootForGroups(Tree* t){ try { - m->names.clear(); - numUniquesInName = 0; - - ifstream in; - m->openInputFile(namefile, in); - - string first, second; - map::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 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; + map roots; //maps group to root for group, may not be root of tree + map done; + + //initialize root for all groups to -1 + for (int k = 0; k < (t->getCountTable())->getNamesOfGroups().size(); k++) { done[(t->getCountTable())->getNamesOfGroups()[k]] = false; } + + for (int i = 0; i < t->getNumLeaves(); i++) { + + vector groups = t->tree[i].getGroup(); + + int index = t->tree[i].getParent(); + + for (int j = 0; j < groups.size(); j++) { + + if (done[groups[j]] == false) { //we haven't found the root for this group yet, initialize it + done[groups[j]] = true; + roots[groups[j]] = i; //set root to self to start + } + + //while you aren't at root + while(t->tree[index].getParent() != -1){ + + if (m->control_pressed) { return roots; } + + //do both your chidren have have descendants from the users groups? + int lc = t->tree[index].getLChild(); + int rc = t->tree[index].getRChild(); + + int LpcountSize = 0; + map:: iterator itGroup = t->tree[lc].pcount.find(groups[j]); + if (itGroup != t->tree[lc].pcount.end()) { LpcountSize++; } + + int RpcountSize = 0; + itGroup = t->tree[rc].pcount.find(groups[j]); + if (itGroup != t->tree[rc].pcount.end()) { RpcountSize++; } + + if ((LpcountSize != 0) && (RpcountSize != 0)) { //possible root + if (index > roots[groups[j]]) { roots[groups[j]] = index; } + }else { ;} + + index = t->tree[index].getParent(); + } + //} + } + } + + + + return roots; + } catch(exception& e) { - m->errorOut(e, "PhyloDiversityCommand", "readNamesFile"); + m->errorOut(e, "PhyloDiversityCommand", "getRootForGroups"); exit(1); } } - //**********************************************************************************************************************