X-Git-Url: https://git.donarmstrong.com/?p=mothur.git;a=blobdiff_plain;f=phylodiversitycommand.cpp;h=339649c7dcc814ac409e5b8ec1dfd432b544e43d;hp=29e2ef0ba39ba4fbe73bf88949bf7a4756d50578;hb=cf9987b67aa49777a4c91c2d21f96e58bf17aa82;hpb=f6b5fea6e0fcc395235dfbfb466b9b6a031b540a diff --git a/phylodiversitycommand.cpp b/phylodiversitycommand.cpp index 29e2ef0..339649c 100644 --- a/phylodiversitycommand.cpp +++ b/phylodiversitycommand.cpp @@ -14,19 +14,20 @@ 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); } @@ -41,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"; @@ -61,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(){ @@ -135,6 +152,14 @@ 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; } + } } //check for required parameters @@ -158,6 +183,19 @@ PhyloDiversityCommand::PhyloDiversityCommand(string option) { 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; @@ -193,10 +231,12 @@ PhyloDiversityCommand::PhyloDiversityCommand(string option) { 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 (namefile == "") { - vector files; files.push_back(treefile); - parser.getNameFile(files); - } + if (countfile=="") { + if (namefile == "") { + vector files; files.push_back(treefile); + parser.getNameFile(files); + } + } } } @@ -215,14 +255,16 @@ int PhyloDiversityCommand::execute(){ int start = time(NULL); m->setTreeFile(treefile); - TreeReader* reader = new TreeReader(treefile, groupfile, namefile); + TreeReader* reader; + if (countfile == "") { reader = new TreeReader(treefile, groupfile, namefile); } + else { reader = new TreeReader(treefile, countfile); } vector trees = reader->getTrees(); - tmap = trees[0]->getTreeMap(); + ct = trees[0]->getCountTable(); delete reader; SharedUtil util; vector mGroups = m->getGroups(); - vector tGroups = tmap->getNamesOfGroups(); + 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 @@ -234,12 +276,15 @@ int PhyloDiversityCommand::execute(){ //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++) { m->mothurRemove(outputNames[j]); } 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); } @@ -254,6 +299,31 @@ int PhyloDiversityCommand::execute(){ 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 @@ -265,15 +335,16 @@ int PhyloDiversityCommand::execute(){ //find largest group total int largestGroup = 0; - for (int j = 0; j < mGroups.size(); j++) { - if (tmap->seqsPerGroup[mGroups[j]] > largestGroup) { largestGroup = tmap->seqsPerGroup[mGroups[j]]; } + for (int j = 0; j < mGroups.size(); j++) { + int numSeqsThisGroup = ct->getGroupCount(mGroups[j]); + if (numSeqsThisGroup > largestGroup) { largestGroup = numSeqsThisGroup; } //initialize diversity - diversity[mGroups[j]].resize(tmap->seqsPerGroup[mGroups[j]]+1, 0.0); //numSampled + diversity[mGroups[j]].resize(numSeqsThisGroup+1, 0.0); //numSampled //groupA 0.0 0.0 //initialize sumDiversity - sumDiversity[mGroups[j]].resize(tmap->seqsPerGroup[mGroups[j]]+1, 0.0); + sumDiversity[mGroups[j]].resize(numSeqsThisGroup+1, 0.0); } //convert freq percentage to number @@ -291,41 +362,29 @@ int PhyloDiversityCommand::execute(){ if (numSampledList.count(diversity[mGroups[j]].size()-1) == 0) { numSampledList.insert(diversity[mGroups[j]].size()-1); } } - #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix) - 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++) { m->mothurRemove(outputNames[i]); } return 0; } - m->mothurOut("It took " + toString(time(NULL) - start) + " secs to run unifrac.unweighted."); m->mothurOutEndLine(); + m->mothurOut("It took " + toString(time(NULL) - start) + " secs to run phylo.diversity."); m->mothurOutEndLine(); m->mothurOutEndLine(); @@ -344,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) || (__linux__) || (__unix__) || (__unix) - 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(); @@ -418,6 +478,56 @@ int PhyloDiversityCommand::createProcesses(vector& procIters, Tree* t, map< in.close(); 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 @@ -439,7 +549,7 @@ int PhyloDiversityCommand::driver(Tree* t, map< string, vector >& div, ma for (int l = 0; l < numIters; l++) { random_shuffle(randomLeaf.begin(), randomLeaf.end()); - cout << l << endl; + //initialize counts map counts; vector< map > countedBranch; @@ -523,6 +633,7 @@ void PhyloDiversityCommand::printSumData(map< string, vector >& div, ofst 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(); @@ -587,7 +698,7 @@ vector PhyloDiversityCommand::calcBranchLength(Tree* t, int leaf, vector< //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()); + sums[k] += abs(t->tree[index].getBranchLength()); } } @@ -628,7 +739,7 @@ map PhyloDiversityCommand::getRootForGroups(Tree* t){ map done; //initialize root for all groups to -1 - for (int k = 0; k < (t->getTreeMap())->getNamesOfGroups().size(); k++) { done[(t->getTreeMap())->getNamesOfGroups()[k]] = false; } + for (int k = 0; k < (t->getCountTable())->getNamesOfGroups().size(); k++) { done[(t->getCountTable())->getNamesOfGroups()[k]] = false; } for (int i = 0; i < t->getNumLeaves(); i++) { @@ -638,10 +749,10 @@ map PhyloDiversityCommand::getRootForGroups(Tree* t){ for (int j = 0; j < groups.size(); j++) { - if (done[groups[j]] == false) { //we haven't found the root for this group yet - - done[groups[j]] = true; - roots[groups[j]] = i; //set root to self to start + 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){ @@ -661,15 +772,17 @@ map PhyloDiversityCommand::getRootForGroups(Tree* t){ if (itGroup != t->tree[rc].pcount.end()) { RpcountSize++; } if ((LpcountSize != 0) && (RpcountSize != 0)) { //possible root - roots[groups[j]] = index; + if (index > roots[groups[j]]) { roots[groups[j]] = index; } }else { ;} index = t->tree[index].getParent(); } - } + //} } } + + return roots; }