From 7b287636ea67fb2254b119c21b4057a177d3ce19 Mon Sep 17 00:00:00 2001 From: pschloss Date: Mon, 21 Feb 2011 21:14:58 +0000 Subject: [PATCH] reworked amova / homova / anosim --- Mothur.xcodeproj/project.pbxproj | 4 +- amovacommand.cpp | 754 ++---------------------- amovacommand.h | 25 +- anosimcommand.cpp | 950 ++++++++----------------------- anosimcommand.h | 39 +- homovacommand.cpp | 922 +++++++----------------------- homovacommand.h | 48 +- readphylipvector.cpp | 2 + venn.cpp | 28 +- 9 files changed, 547 insertions(+), 2225 deletions(-) diff --git a/Mothur.xcodeproj/project.pbxproj b/Mothur.xcodeproj/project.pbxproj index e479725..902e927 100644 --- a/Mothur.xcodeproj/project.pbxproj +++ b/Mothur.xcodeproj/project.pbxproj @@ -1625,7 +1625,7 @@ attributes = { ORGANIZATIONNAME = "Schloss Lab"; }; - buildConfigurationList = 1DEB928908733DD80010E9CD /* Build configuration list for PBXProject "mothur" */; + buildConfigurationList = 1DEB928908733DD80010E9CD /* Build configuration list for PBXProject "Mothur" */; compatibilityVersion = "Xcode 3.1"; developmentRegion = English; hasScannedForEncodings = 1; @@ -2048,7 +2048,7 @@ defaultConfigurationIsVisible = 0; defaultConfigurationName = Release; }; - 1DEB928908733DD80010E9CD /* Build configuration list for PBXProject "mothur" */ = { + 1DEB928908733DD80010E9CD /* Build configuration list for PBXProject "Mothur" */ = { isa = XCConfigurationList; buildConfigurations = ( 1DEB928A08733DD80010E9CD /* Debug */, diff --git a/amovacommand.cpp b/amovacommand.cpp index 4dfe43b..b9b5dee 100644 --- a/amovacommand.cpp +++ b/amovacommand.cpp @@ -8,55 +8,13 @@ */ #include "amovacommand.h" -#include "sharedutilities.h" -//#include "sharedsobscollectsummary.h" -//#include "sharedchao1.h" -//#include "sharedace.h" -//#include "sharednseqs.h" -//#include "sharedjabund.h" -//#include "sharedsorabund.h" -//#include "sharedjclass.h" -//#include "sharedsorclass.h" -//#include "sharedjest.h" -//#include "sharedsorest.h" -//#include "sharedthetayc.h" -//#include "sharedthetan.h" -//#include "sharedkstest.h" -//#include "whittaker.h" -//#include "sharedochiai.h" -//#include "sharedanderbergs.h" -//#include "sharedkulczynski.h" -//#include "sharedkulczynskicody.h" -//#include "sharedlennon.h" -//#include "sharedmorisitahorn.h" -//#include "sharedbraycurtis.h" -//#include "sharedjackknife.h" -//#include "whittaker.h" -//#include "odum.h" -//#include "canberra.h" -//#include "structeuclidean.h" -//#include "structchord.h" -//#include "hellinger.h" -//#include "manhattan.h" -//#include "structpearson.h" -//#include "soergel.h" -//#include "spearman.h" -//#include "structkulczynski.h" -//#include "structchi2.h" -//#include "speciesprofile.h" -//#include "hamming.h" -//#include "gower.h" -//#include "memchi2.h" -//#include "memchord.h" -//#include "memeuclidean.h" -//#include "mempearson.h" +#include "readphylipvector.h" +#include "groupmap.h" //********************************************************************************************************************** vector AmovaCommand::getValidParameters(){ try { - string Array[] = {"groups","label","outputdir","iters","phylip","design","alpha", - //"sets","processors" - "inputdir"}; + string Array[] = {"outputdir","iters","phylip","design","alpha", "inputdir"}; vector myArray (Array, Array+(sizeof(Array)/sizeof(string))); return myArray; } @@ -105,19 +63,14 @@ vector AmovaCommand::getRequiredFiles(){ AmovaCommand::AmovaCommand(string option) { try { - globaldata = GlobalData::getInstance(); abort = false; calledHelp = false; -// allLines = 1; -// labels.clear(); //allow user to run help if(option == "help") { help(); abort = true; calledHelp = true; } else { //valid paramters for this command - string AlignArray[] = {"groups","label","outputdir","iters","phylip","alpha", -// "design","sets","processors", - "inputdir"}; + string AlignArray[] = {"design","outputdir","iters","phylip","alpha", "inputdir"}; vector myArray (AlignArray, AlignArray+(sizeof(AlignArray)/sizeof(string))); OptionParser parser(option); @@ -163,10 +116,12 @@ AmovaCommand::AmovaCommand(string option) { phylipFileName = validParameter.validFile(parameters, "phylip", true); if (phylipFileName == "not open") { phylipFileName = ""; abort = true; } else if (phylipFileName == "not found") { phylipFileName = ""; } - else { - globaldata->newRead(); - globaldata->setPhylipFile(phylipFileName); - } + else if (designFileName == "not found") { + designFileName = ""; + m->mothurOut("You must provide an phylip file."); + m->mothurOutEndLine(); + abort = true; + } //check for required parameters designFileName = validParameter.validFile(parameters, "design", true); @@ -178,151 +133,14 @@ AmovaCommand::AmovaCommand(string option) { abort = true; } - string temp = validParameter.validFile(parameters, "iters", false); if (temp == "not found") { temp = "1000"; } + string temp = validParameter.validFile(parameters, "iters", false); + if (temp == "not found") { temp = "1000"; } convert(temp, iters); - temp = validParameter.validFile(parameters, "alpha", false); if (temp == "not found") { temp = "0.05"; } + temp = validParameter.validFile(parameters, "alpha", false); + if (temp == "not found") { temp = "0.05"; } convert(temp, experimentwiseAlpha); - -// make sure the user has already run the read.otu command -// if ((globaldata->getSharedFile() == "")) { -// if ((phylipfile == "")) { -// m->mothurOut("You must read a list and a group, a shared file or provide a distance matrix before you can use the amova command."); -// m->mothurOutEndLine(); -// abort = true; -// } -// }else { sharedfile = globaldata->getSharedFile(); } -// -// //use distance matrix if one is provided -// if ((sharedfile != "") && (phylipfile != "")) { sharedfile = ""; } -// -// //check for optional parameter and set defaults -// // ...at some point should added some additional type checking... -// label = validParameter.validFile(parameters, "label", false); -// if (label == "not found") { label = ""; } -// else { -// if(label != "all") { m->splitAtDash(label, labels); allLines = 0; } -// else { allLines = 1; } -// } -// -// //if the user has not specified any labels use the ones from read.otu -// if (label == "") { -// allLines = globaldata->allLines; -// labels = globaldata->labels; -// } -// -// groups = validParameter.validFile(parameters, "groups", false); -// if (groups == "not found") { groups = ""; } -// else { -// m->splitAtDash(groups, Groups); -// globaldata->Groups = Groups; -// } -// -// sets = validParameter.validFile(parameters, "sets", false); -// if (sets == "not found") { sets = ""; pickedGroups = false; } -// else { -// pickedGroups = true; -// m->splitAtDash(sets, Sets); -// } -// -// temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = "1"; } -// convert(temp, processors); -// -// vector Estimators; -// calc = validParameter.validFile(parameters, "calc", false); -// if (calc == "not found") { calc = "morisitahorn"; } -// m->splitAtDash(calc, Estimators); -// -// if (abort == false) { -// -// ValidCalculators* validCalculator = new ValidCalculators(); -// -// for (int i=0; iisValidCalculator("treegroup", Estimators[i]) == true) { -// if (Estimators[i] == "sharedsobs") { -// calculators.push_back(new SharedSobsCS()); -// }else if (Estimators[i] == "sharedchao") { -// calculators.push_back(new SharedChao1()); -// }else if (Estimators[i] == "sharedace") { -// calculators.push_back(new SharedAce()); -// }else if (Estimators[i] == "jabund") { -// calculators.push_back(new JAbund()); -// }else if (Estimators[i] == "sorabund") { -// calculators.push_back(new SorAbund()); -// }else if (Estimators[i] == "jclass") { -// calculators.push_back(new Jclass()); -// }else if (Estimators[i] == "sorclass") { -// calculators.push_back(new SorClass()); -// }else if (Estimators[i] == "jest") { -// calculators.push_back(new Jest()); -// }else if (Estimators[i] == "sorest") { -// calculators.push_back(new SorEst()); -// }else if (Estimators[i] == "thetayc") { -// calculators.push_back(new ThetaYC()); -// }else if (Estimators[i] == "thetan") { -// calculators.push_back(new ThetaN()); -// }else if (Estimators[i] == "kstest") { -// calculators.push_back(new KSTest()); -// }else if (Estimators[i] == "sharednseqs") { -// calculators.push_back(new SharedNSeqs()); -// }else if (Estimators[i] == "ochiai") { -// calculators.push_back(new Ochiai()); -// }else if (Estimators[i] == "anderberg") { -// calculators.push_back(new Anderberg()); -// }else if (Estimators[i] == "kulczynski") { -// calculators.push_back(new Kulczynski()); -// }else if (Estimators[i] == "kulczynskicody") { -// calculators.push_back(new KulczynskiCody()); -// }else if (Estimators[i] == "lennon") { -// calculators.push_back(new Lennon()); -// }else if (Estimators[i] == "morisitahorn") { -// calculators.push_back(new MorHorn()); -// }else if (Estimators[i] == "braycurtis") { -// calculators.push_back(new BrayCurtis()); -// }else if (Estimators[i] == "whittaker") { -// calculators.push_back(new Whittaker()); -// }else if (Estimators[i] == "odum") { -// calculators.push_back(new Odum()); -// }else if (Estimators[i] == "canberra") { -// calculators.push_back(new Canberra()); -// }else if (Estimators[i] == "structeuclidean") { -// calculators.push_back(new StructEuclidean()); -// }else if (Estimators[i] == "structchord") { -// calculators.push_back(new StructChord()); -// }else if (Estimators[i] == "hellinger") { -// calculators.push_back(new Hellinger()); -// }else if (Estimators[i] == "manhattan") { -// calculators.push_back(new Manhattan()); -// }else if (Estimators[i] == "structpearson") { -// calculators.push_back(new StructPearson()); -// }else if (Estimators[i] == "soergel") { -// calculators.push_back(new Soergel()); -// }else if (Estimators[i] == "spearman") { -// calculators.push_back(new Spearman()); -// }else if (Estimators[i] == "structkulczynski") { -// calculators.push_back(new StructKulczynski()); -// }else if (Estimators[i] == "speciesprofile") { -// calculators.push_back(new SpeciesProfile()); -// }else if (Estimators[i] == "hamming") { -// calculators.push_back(new Hamming()); -// }else if (Estimators[i] == "structchi2") { -// calculators.push_back(new StructChi2()); -// }else if (Estimators[i] == "gower") { -// calculators.push_back(new Gower()); -// }else if (Estimators[i] == "memchi2") { -// calculators.push_back(new MemChi2()); -// }else if (Estimators[i] == "memchord") { -// calculators.push_back(new MemChord()); -// }else if (Estimators[i] == "memeuclidean") { -// calculators.push_back(new MemEuclidean()); -// }else if (Estimators[i] == "mempearson") { -// calculators.push_back(new MemPearson()); -// } -// } -// } -// } } - } catch(exception& e) { m->errorOut(e, "AmovaCommand", "AmovaCommand"); @@ -335,19 +153,13 @@ AmovaCommand::AmovaCommand(string option) { void AmovaCommand::help(){ try { m->mothurOut("Referenced: Anderson MJ (2001). A new method for non-parametric multivariate analysis of variance. Austral Ecol 26: 32-46.\n"); -// m->mothurOut("The amova command can only be executed after a successful read.otu command of a list and group or shared file, or by providing a phylip formatted distance matrix.\n"); m->mothurOut("The amova command outputs a .amova file. \n"); m->mothurOut("The amova command parameters are phylip, iters, and alpha. The phylip and design parameters are required.\n"); m->mothurOut("The design parameter allows you to assign your samples to groups when you are running amova. It is required. \n"); m->mothurOut("The design file looks like the group file. It is a 2 column tab delimited file, where the first column is the sample name and the second column is the group the sample belongs to.\n"); -// m->mothurOut("The sets parameter allows you to specify which of the sets in your designfile you would like to analyze. The set names are separated by dashes. THe default is all sets in the designfile. To run the pairwise comparisons of all sets use sets=all.\n"); m->mothurOut("The iters parameter allows you to set number of randomization for the P value. The default is 1000. \n"); -// m->mothurOut("The groups parameter allows you to specify which of the groups you would like included. The group names are separated by dashes. groups=all will find all pairwise comparisons. \n"); -// m->mothurOut("The label parameter allows you to select what distance levels you would like, and are also separated by dashes.\n"); -// m->mothurOut("The processors parameter allows you to specify how many processors you would like to use. The default is 1. \n"); m->mothurOut("The amova command should be in the following format: amova(phylip=file.dist, design=file.design).\n"); m->mothurOut("Note: No spaces between parameter labels (i.e. iters), '=' and parameters (i.e. 1000).\n\n"); - } catch(exception& e) { m->errorOut(e, "AmovaCommand", "help"); @@ -375,7 +187,6 @@ int AmovaCommand::execute(){ //read in distance matrix and square it ReadPhylipVector readMatrix(phylipFileName); vector sampleNames = readMatrix.read(distanceMatrix); - int numSamples = sampleNames.size(); for(int i=0;i > origGroupSampleMap; - for(int i=0;igetGroup(sampleNames[i])].push_back(i); } int numGroups = origGroupSampleMap.size(); @@ -395,8 +206,6 @@ int AmovaCommand::execute(){ string AMOVAFileName = outputDir + m->getRootName(m->getSimpleName(phylipFileName)) + "amova"; m->openOutputFile(AMOVAFileName, AMOVAFile); outputNames.push_back(AMOVAFileName); outputTypes["amova"].push_back(AMOVAFileName); - - experimentwiseAlpha = 0.05; double fullANOVAPValue = runAMOVA(AMOVAFile, origGroupSampleMap, experimentwiseAlpha); if(fullANOVAPValue <= experimentwiseAlpha && numGroups > 2){ @@ -418,232 +227,17 @@ int AmovaCommand::execute(){ runAMOVA(AMOVAFile, pairwiseGroupSampleMap, pairwiseAlpha); } } - cout << "Experiment-wise error rate: " << experimentwiseAlpha << endl; - cout << "Pair-wise error rate (Bonferroni): " << pairwiseAlpha << endl; + m->mothurOut("Experiment-wise error rate: " + toString(experimentwiseAlpha) + '\n'); + m->mothurOut("Pair-wise error rate (Bonferroni): " + toString(pairwiseAlpha) + '\n'); } else{ - cout << "Experiment-wise error rate: " << experimentwiseAlpha << endl; + m->mothurOut("Experiment-wise error rate: " + toString(experimentwiseAlpha) + '\n'); } - - cout << "If you have borderline P-values, you should try increasing the number of iterations" << endl; - - -// //make sure sets are all in designMap -// SharedUtil* util = new SharedUtil(); -// util->setGroups(Sets, designMap->namesOfGroups); -// delete util; -// -// //if the user pickedGroups or set sets=all, then we want to figure out how to split the pairwise comparison between processors -// if (pickedGroups) { -// for (int a=0; a groups; groups.push_back(Sets[a]); groups.push_back(Sets[l]); -// namesOfGroupCombos.push_back(groups); -// } -// } -// }else { //one giant compare -// vector groups; -// for (int a=0; amothurOut("Not enough sets, I need at least 2 valid sets. Unable to complete command."); -// m->mothurOutEndLine(); m->control_pressed = true; -// } -// -// #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) -// if(processors != 1){ -// int numPairs = namesOfGroupCombos.size(); -// int numPairsPerProcessor = numPairs / processors; -// -// for (int i = 0; i < processors; i++) { -// int startPos = i * numPairsPerProcessor; -// if(i == processors - 1){ -// numPairsPerProcessor = numPairs - i * numPairsPerProcessor; -// } -// lines.push_back(linePair(startPos, numPairsPerProcessor)); -// } -// } -// #endif -// -// if (sharedfile != "") { //create distance matrix for each label -// -// if (outputDir == "") { outputDir = m->hasPath(sharedfile); } -// -// //for each calculator -// for(int i = 0 ; i < calculators.size(); i++) { -// -// //create a new filename -// ofstream out; -// string outputFile = outputDir + m->getRootName(m->getSimpleName(sharedfile)) + calculators[i]->getName() + ".amova"; -// m->openOutputFile(outputFile, out); -// outputNames.push_back(outputFile); outputTypes["amova"].push_back(outputFile); -// -// //print headers -// out << "label\tgroupsCompared\tMeanSquaredWithin\tMeanSquaredAmong\tFstatistic\tpValue" << endl; -// out.close(); -// m->mothurOut("label\tgroupsCompared\tMeanSquaredWithin\tMeanSquaredAmong\tFstatistic\tpValue\n"); -// } -// -// InputData input(sharedfile, "sharedfile"); -// vector lookup = input.getSharedRAbundVectors(); -// string lastLabel = lookup[0]->getLabel(); -// -// //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label. -// set processedLabels; -// set userLabels = labels; -// -// //sanity check between shared file groups and design file groups -// for (int i = 0; i < lookup.size(); i++) { -// string group = designMap->getGroup(lookup[i]->getGroup()); -// -// if (group == "not found") { m->control_pressed = true; m->mothurOut("[ERROR]: " + lookup[i]->getGroup() + " is not in your design file, please correct."); m->mothurOutEndLine(); } -// } -// -// //as long as you are not at the end of the file or done wih the lines you want -// while((lookup[0] != NULL) && ((allLines == 1) || (userLabels.size() != 0))) { -// -// if (m->control_pressed) { for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } globaldata->Groups.clear(); delete designMap; for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; } -// -// if(allLines == 1 || labels.count(lookup[0]->getLabel()) == 1){ -// -// m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine(); -// process(lookup); -// -// processedLabels.insert(lookup[0]->getLabel()); -// userLabels.erase(lookup[0]->getLabel()); -// } -// -// if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) { -// string saveLabel = lookup[0]->getLabel(); -// -// for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } -// lookup = input.getSharedRAbundVectors(lastLabel); -// m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine(); -// -// process(lookup); -// -// processedLabels.insert(lookup[0]->getLabel()); -// userLabels.erase(lookup[0]->getLabel()); -// -// //restore real lastlabel to save below -// lookup[0]->setLabel(saveLabel); -// } -// -// lastLabel = lookup[0]->getLabel(); -// //prevent memory leak -// for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; lookup[i] = NULL; } -// -// if (m->control_pressed) { globaldata->Groups.clear(); delete designMap; for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; } -// -// //get next line to process -// lookup = input.getSharedRAbundVectors(); -// } -// -// if (m->control_pressed) { globaldata->Groups.clear(); delete designMap; for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; } -// -// //output error messages about any remaining user labels -// set::iterator it; -// bool needToRun = false; -// for (it = userLabels.begin(); it != userLabels.end(); it++) { -// m->mothurOut("Your file does not include the label " + *it); -// if (processedLabels.count(lastLabel) != 1) { -// m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine(); -// needToRun = true; -// }else { -// m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine(); -// } -// } -// -// //run last label if you need to -// if (needToRun == true) { -// for (int i = 0; i < lookup.size(); i++) { if (lookup[i] != NULL) { delete lookup[i]; } } -// lookup = input.getSharedRAbundVectors(lastLabel); -// -// m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine(); -// -// process(lookup); -// -// for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } -// } -// -// //reset groups parameter -// globaldata->Groups.clear(); -// -// -// }else { //user provided distance matrix -// -// if (outputDir == "") { outputDir = m->hasPath(phylipFile); } -// -// //create a new filename -// ofstream out; -// string outputFile = outputDir + m->getRootName(m->getSimpleName(phylipFile)) + "amova"; -// m->openOutputFile(outputFile, out); -// outputNames.push_back(outputFile); outputTypes["amova"].push_back(outputFile); -// -// //print headers -// out << "groupsCompared\tMeanSquaredWithin\tMeanSquaredAmong\tFstatistic\tpValue" << endl; -// out.close(); -// m->mothurOut("groupsCompared\tMeanSquaredWithin\tMeanSquaredAmong\tFstatistic\tpValue\n"); -// -// ReadPhylipVector readMatrix(phylipFile); -// vector< vector > completeMatrix; -// vector names = readMatrix.read(completeMatrix); -// -// #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) -// if(processors == 1){ -// driver(0, namesOfGroupCombos.size(), names, "", completeMatrix); -// }else{ -// int process = 1; -// vector processIDS; -// -// //loop through and create all the processes you want -// while (process != processors) { -// int pid = fork(); -// -// if (pid > 0) { -// processIDS.push_back(pid); -// process++; -// }else if (pid == 0){ -// driver(lines[process].start, lines[process].num, names, toString(getpid()), completeMatrix); -// exit(0); -// }else { -// m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine(); -// for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); } -// exit(0); -// } -// } -// -// //do my part -// driver(lines[0].start, lines[0].num, names, "", completeMatrix); -// -// //force parent to wait until all the processes are done -// for (int i=0;i<(processors-1);i++) { -// int temp = processIDS[i]; -// wait(&temp); -// } -// -// //append files -// string outputFile = outputDir + m->getRootName(m->getSimpleName(phylipFile)) + "amova"; -// for (int j = 0; j < processIDS.size(); j++) { -// m->appendFiles((outputFile + toString(processIDS[j])), outputFile); -// remove((outputFile + toString(processIDS[j])).c_str()); -// } -// -// } -// #else -// driver(0, namesOfGroupCombos.size(), names, "", completeMatrix); -// #endif -// -// } + m->mothurOut("If you have borderline P-values, you should try increasing the number of iterations\n"); + AMOVAFile.close(); delete designMap; - if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; } - m->mothurOutEndLine(); m->mothurOut("Output File Names: "); m->mothurOutEndLine(); for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); } @@ -664,10 +258,10 @@ double AmovaCommand::runAMOVA(ofstream& AMOVAFile, map > gro map >::iterator it; int numGroups = groupSampleMap.size(); - int numSamples = 0; + int totalNumSamples = 0; for(it = groupSampleMap.begin();it!=groupSampleMap.end();it++){ - numSamples += it->second.size(); + totalNumSamples += it->second.size(); } double ssTotalOrig = calcSSTotal(groupSampleMap); @@ -689,29 +283,47 @@ double AmovaCommand::runAMOVA(ofstream& AMOVAFile, map > gro //print anova table it = groupSampleMap.begin(); - cout << it->first; + AMOVAFile << it->first; + m->mothurOut(it->first); it++; for(it;it!=groupSampleMap.end();it++){ - cout << " vs. " << it->first; - }cout << endl; - cout << "ANOVA\tAmong\tWithin\tTotal" << endl; - cout << "SS\t" << ssAmongOrig << '\t' << ssWithinOrig << '\t' << ssTotalOrig << endl; + AMOVAFile << '-' << it->first; + m->mothurOut('-' + it->first); + } + + AMOVAFile << "\tAmong\tWithin\tTotal" << endl; + m->mothurOut("\tAmong\tWithin\tTotal\n"); + + AMOVAFile << "SS\t" << ssAmongOrig << '\t' << ssWithinOrig << '\t' << ssTotalOrig << endl; + m->mothurOut("SS\t" + toString(ssAmongOrig) + '\t' + toString(ssWithinOrig) + '\t' + toString(ssTotalOrig) + '\n'); + + int dfAmong = numGroups - 1; double MSAmong = ssAmongOrig / (double) dfAmong; + int dfWithin = totalNumSamples - numGroups; double MSWithin = ssWithinOrig / (double) dfWithin; + int dfTotal = totalNumSamples - 1; double Fs = MSAmong / MSWithin; - int dfAmong = numGroups - 1; double MSAmong = ssAmongOrig / (double) dfAmong; - int dfWithin = numSamples - numGroups; double MSWithin = ssWithinOrig / (double) dfWithin; - int dfTotal = numSamples - 1; double Fs = MSAmong / MSWithin; + AMOVAFile << "df\t" << dfAmong << '\t' << dfWithin << '\t' << dfTotal << endl; + m->mothurOut("df\t" + toString(dfAmong) + '\t' + toString(dfWithin) + '\t' + toString(dfTotal) + '\n'); + + AMOVAFile << "MS\t" << MSAmong << '\t' << MSWithin << endl << endl; + m->mothurOut("MS\t" + toString(MSAmong) + '\t' + toString(MSWithin) + "\n\n"); + + AMOVAFile << "Fs:\t" << Fs << endl; + m->mothurOut("Fs:\t" + toString(Fs) + '\n'); - cout << "df\t" << dfAmong << '\t' << dfWithin << '\t' << dfTotal << endl; - cout << "MS\t" << MSAmong << '\t' << MSWithin << endl << endl; - cout << "Fs:\t" << Fs << endl; - cout << "p-value: " << pString; - if(pValue < alpha){ cout << " ***"; } - cout << endl << endl; + AMOVAFile << "p-value: " << pString; + m->mothurOut("p-value: " + pString); + + if(pValue < alpha){ + AMOVAFile << "*"; + m->mothurOut("*"); + } + AMOVAFile << endl << endl; + m->mothurOutEndLine();m->mothurOutEndLine(); return pValue; } catch(exception& e) { - m->errorOut(e, "AmovaCommand", "calcAmova"); + m->errorOut(e, "AmovaCommand", "runAMOVA"); exit(1); } } @@ -743,7 +355,7 @@ map > AmovaCommand::getRandomizedGroups(maperrorOut(e, "AmovaCommand", "randomizeGroups"); + m->errorOut(e, "AmovaCommand", "getRandomizedGroups"); exit(1); } } @@ -813,261 +425,9 @@ double AmovaCommand::calcSSWithin(map >& groupSampleMap) { return ssWithin; } catch(exception& e) { - m->errorOut(e, "AmovaCommand", "calcWithin"); + m->errorOut(e, "AmovaCommand", "calcSSWithin"); exit(1); } } //********************************************************************************************************************** - -//int AmovaCommand::process(vector thisLookup) { -// try{ -// -//#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) -// if(processors == 1){ -// driver(0, namesOfGroupCombos.size(), thisLookup, ""); -// }else{ -// int process = 1; -// vector processIDS; -// -// //loop through and create all the processes you want -// while (process != processors) { -// int pid = fork(); -// -// if (pid > 0) { -// processIDS.push_back(pid); -// process++; -// }else if (pid == 0){ -// driver(lines[process].start, lines[process].num, thisLookup, toString(getpid())); -// exit(0); -// }else { -// m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine(); -// for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); } -// exit(0); -// } -// } -// -// //do my part -// driver(lines[0].start, lines[0].num, thisLookup, ""); -// -// //force parent to wait until all the processes are done -// for (int i=0;i<(processors-1);i++) { -// int temp = processIDS[i]; -// wait(&temp); -// } -// -// //append files -// for(int i = 0 ; i < calculators.size(); i++) { -// string outputFile = outputDir + m->getRootName(m->getSimpleName(sharedfile)) + calculators[i]->getName() + ".amova"; -// -// for (int j = 0; j < processIDS.size(); j++) { -// m->appendFiles((outputFile + toString(processIDS[j])), outputFile); -// remove((outputFile + toString(processIDS[j])).c_str()); -// } -// } -// } -//#else -// driver(0, namesOfGroupCombos.size(), thisLookUp, ""); -//#endif -// -// return 0; -// } -// catch(exception& e) { -// m->errorOut(e, "AmovaCommand", "process"); -// exit(1); -// } -//} - -//********************************************************************************************************************** - -//int AmovaCommand::driver(int start, int num, vector thisLookup, string pidValue) { -// try { -// vector subset; -// EstOutput data; -// -// //for each calculator -// for(int i = 0 ; i < calculators.size(); i++) { -// -// //create a new filename -// ofstream out; -// string outputFile = outputDir + m->getRootName(m->getSimpleName(sharedfile)) + calculators[i]->getName() + ".amova" + pidValue; -// m->openOutputFileAppend(outputFile, out); -// out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint); -// -// //for each combo -// for (int c = start; c < (start+num); c++) { -// -// if (m->control_pressed) { out.close(); return 0; } -// -// //get set names -// vector setNames; -// for (int j = 0; j < namesOfGroupCombos[c].size(); j++) { setNames.push_back(namesOfGroupCombos[c][j]); } -// -// vector thisCombosLookup; -// vector thisCombosLookupSets; //what set each sharedRabund is from to be used when calculating SSWithin -// for (int k = 0; k < thisLookup.size(); k++) { -// string thisGroup = thisLookup[k]->getGroup(); -// -// //is this group for a set we want to compare?? -// if (m->inUsersGroups(designMap->getGroup(thisGroup), setNames)) { -// thisCombosLookup.push_back(thisLookup[k]); -// thisCombosLookupSets.push_back(designMap->getGroup(thisGroup)); -// } -// -// } -// -// int numGroups = thisCombosLookup.size(); -// -// //calc the distance matrix -// matrix.clear(); -// matrix.resize(numGroups); -// for (int k = 0; k < matrix.size(); k++) { for (int j = 0; j < matrix.size(); j++) { matrix[k].push_back(1.0); } } -// -// if (thisCombosLookup.size() == 0) { -// m->mothurOut("[ERROR]: Missing shared info for sets. Skipping comparison."); m->mothurOutEndLine(); -// }else{ -// -// m->mothurOut(thisLookup[0]->getLabel() + '\t'); -// out << thisLookup[0]->getLabel() << '\t'; -// if (setNames.size() == 2) { -// m->mothurOut(setNames[0] + '-' + setNames[1] + '\t'); -// out << setNames[0] << "-" << setNames[1] << '\t'; -// } -// else { -// out << "all" << '\t'; -// m->mothurOut("all\t"); -// } -// -// for (int k = 0; k < thisCombosLookup.size(); k++) { -// for (int l = k; l < thisCombosLookup.size(); l++) { -// -// if (m->control_pressed) { out.close(); return 0; } -// -// if (k != l) { //we dont need to similiarity of a groups to itself -// //get estimated similarity between 2 groups -// subset.clear(); //clear out old pair of sharedrabunds -// //add new pair of sharedrabunds -// subset.push_back(thisCombosLookup[k]); subset.push_back(thisCombosLookup[l]); -// -// //if this calc needs all groups to calculate the pair load all groups -// if (calculators[i]->getNeedsAll()) { -// //load subset with rest of lookup for those calcs that need everyone to calc for a pair -// for (int w = 0; w < thisCombosLookup.size(); w++) { -// if ((w != k) && (w != l)) { subset.push_back(thisCombosLookup[w]); } -// } -// } -// -// data = calculators[i]->getValues(subset); //saves the calculator outputs -// -// //save values in similarity matrix -// matrix[k][l] = 1.0 - data[0]; -// matrix[l][k] = 1.0 - data[0]; -// } -// } -// } -// -// //calc amova -// calcAmova(out, setNames.size(), thisCombosLookupSets); -// } -// } -// -// out.close(); -// } -// -// return 0; -// -// } -// catch(exception& e) { -// m->errorOut(e, "AmovaCommand", "driver"); -// exit(1); -// } -//} -// -//********************************************************************************************************************** -// -//int AmovaCommand::driver(int start, int num, vector names, string pidValue, vector< vector >& completeMatrix) { -// try { -// -// //create a new filename -// ofstream out; -// string outputFile = outputDir + m->getRootName(m->getSimpleName(phylipFile)) + "amova" + pidValue; -// m->openOutputFileAppend(outputFile, out); -// out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint); -// -// //for each combo -// for (int c = start; c < (start+num); c++) { -// -// if (m->control_pressed) { out.close(); return 0; } -// -// //get set names -// vector setNames; -// for (int j = 0; j < namesOfGroupCombos[c].size(); j++) { setNames.push_back(namesOfGroupCombos[c][j]); } -// -// vector thisCombosSets; //what set each line in the distance matrix is from to be used when calculating SSWithin -// set indexes; -// for (int k = 0; k < names.size(); k++) { -// //is this group for a set we want to compare?? -// if (m->inUsersGroups(designMap->getGroup(names[k]), setNames)) { -// thisCombosSets.push_back(designMap->getGroup(names[k])); -// indexes.insert(k); //save indexes of valid rows in matrix for submatrix -// } -// } -// -// int numGroups = thisCombosSets.size(); -// -// //calc the distance matrix -// matrix.clear(); -// matrix.resize(numGroups); -// -// for (int k = 0; k < matrix.size(); k++) { for (int j = 0; j < matrix.size(); j++) { matrix[k].push_back(1.0); } } -// -// if (thisCombosSets.size() == 0) { -// m->mothurOut("[ERROR]: Missing distance info for sets. Skipping comparison."); m->mothurOutEndLine(); -// }else{ -// -// if (setNames.size() == 2) { -// out << setNames[0] << "-" << setNames[1] << '\t'; -// m->mothurOut(setNames[0] + '-' + setNames[1] + '\t'); -// } -// else { -// out << "all" << '\t'; -// m->mothurOut("all\t"); -// } -// -// //fill submatrix -// int rowCount = 0; -// for (int j = 0; j < completeMatrix.size(); j++) { -// -// if (indexes.count(j) != 0) { //we want at least part of this row -// int columnCount = 0; -// for (int k = 0; k < completeMatrix[j].size(); k++) { -// -// if (indexes.count(k) != 0) { //we want this distance -// matrix[rowCount][columnCount] = completeMatrix[j][k]; -// matrix[columnCount][rowCount] = completeMatrix[j][k]; -// columnCount++; -// } -// } -// rowCount++; -// } -// } -// -// //calc amova -// calcAmova(out, setNames.size(), thisCombosSets); -// } -// } -// -// out.close(); -// -// -// return 0; -// -// } -// catch(exception& e) { -// m->errorOut(e, "AmovaCommand", "driver"); -// exit(1); -// } -//} -// -//********************************************************************************************************************** - diff --git a/amovacommand.h b/amovacommand.h index d9fe1e6..51a4200 100644 --- a/amovacommand.h +++ b/amovacommand.h @@ -11,12 +11,9 @@ */ #include "command.hpp" -#include "inputdata.h" -#include "sharedrabundvector.h" -#include "validcalculator.h" -#include "readphylipvector.h" -class GlobalData; +//class GlobalData; +class GroupMap; class AmovaCommand : public Command { @@ -39,7 +36,6 @@ private: bool abort; - GlobalData* globaldata; map > outputTypes; vector outputNames; @@ -48,23 +44,6 @@ private: vector< vector > distanceMatrix; int iters; double experimentwiseAlpha; - -// struct linePair { -// int start; -// int num; -// linePair(int i, int j) : start(i), num(j) {} -// }; -// vector lines; -// -// vector< vector > namesOfGroupCombos; -// vector Groups, outputNames, Sets; -// int processors; -// string groups, sets, calc, sharedfile, label, allLines, pickedGroups; -// vector calculators; -// set labels; //holds labels to be used -// int driver(int, int, vector, string); -// int driver(int, int, vector, string, vector< vector >&); -// int process(vector); }; #endif diff --git a/anosimcommand.cpp b/anosimcommand.cpp index 9f537b8..bb0f495 100644 --- a/anosimcommand.cpp +++ b/anosimcommand.cpp @@ -8,53 +8,13 @@ */ #include "anosimcommand.h" -#include "sharedutilities.h" -#include "sharedsobscollectsummary.h" -#include "sharedchao1.h" -#include "sharedace.h" -#include "sharednseqs.h" -#include "sharedjabund.h" -#include "sharedsorabund.h" -#include "sharedjclass.h" -#include "sharedsorclass.h" -#include "sharedjest.h" -#include "sharedsorest.h" -#include "sharedthetayc.h" -#include "sharedthetan.h" -#include "sharedkstest.h" -#include "whittaker.h" -#include "sharedochiai.h" -#include "sharedanderbergs.h" -#include "sharedkulczynski.h" -#include "sharedkulczynskicody.h" -#include "sharedlennon.h" -#include "sharedmorisitahorn.h" -#include "sharedbraycurtis.h" -#include "sharedjackknife.h" -#include "whittaker.h" -#include "odum.h" -#include "canberra.h" -#include "structeuclidean.h" -#include "structchord.h" -#include "hellinger.h" -#include "manhattan.h" -#include "structpearson.h" -#include "soergel.h" -#include "spearman.h" -#include "structkulczynski.h" -#include "structchi2.h" -#include "speciesprofile.h" -#include "hamming.h" -#include "gower.h" -#include "memchi2.h" -#include "memchord.h" -#include "memeuclidean.h" -#include "mempearson.h" +#include "inputdata.h" +#include "readphylipvector.h" //********************************************************************************************************************** vector AnosimCommand::getValidParameters(){ try { - string Array[] = {"groups","label","outputdir","iters","phylip","design","sets","processors","inputdir"}; + string Array[] = {"outputdir","iters","phylip","design", "alpha","inputdir"}; vector myArray (Array, Array+(sizeof(Array)/sizeof(string))); return myArray; } @@ -103,17 +63,14 @@ vector AnosimCommand::getRequiredFiles(){ AnosimCommand::AnosimCommand(string option) { try { - globaldata = GlobalData::getInstance(); abort = false; calledHelp = false; - allLines = 1; - labels.clear(); //allow user to run help if(option == "help") { help(); abort = true; calledHelp = true; } else { //valid paramters for this command - string AlignArray[] = {"groups","label","outputdir","iters","phylip","design","sets","processors","inputdir"}; + string AlignArray[] = {"outputdir","iters","phylip","design", "alpha","inputdir"}; vector myArray (AlignArray, AlignArray+(sizeof(AlignArray)/sizeof(string))); OptionParser parser(option); @@ -156,155 +113,33 @@ AnosimCommand::AnosimCommand(string option) { } } - phylipfile = validParameter.validFile(parameters, "phylip", true); - if (phylipfile == "not open") { phylipfile = ""; abort = true; } - else if (phylipfile == "not found") { phylipfile = ""; } - else { globaldata->newRead(); format = "phylip"; globaldata->setPhylipFile(phylipfile); } + phylipFileName = validParameter.validFile(parameters, "phylip", true); + if (phylipFileName == "not open") { phylipFileName = ""; abort = true; } + else if (phylipFileName == "not found") { phylipFileName = ""; } + else if (designFileName == "not found") { + designFileName = ""; + m->mothurOut("You must provide an phylip file."); + m->mothurOutEndLine(); + abort = true; + } //check for required parameters - designfile = validParameter.validFile(parameters, "design", true); - if (designfile == "not open") { abort = true; } - else if (designfile == "not found") { designfile = ""; m->mothurOut("You must provide an design file."); m->mothurOutEndLine(); abort = true; } - - //make sure the user has already run the read.otu command - if ((globaldata->getSharedFile() == "")) { - if ((phylipfile == "")) { - m->mothurOut("You must read a list and a group, a shared file or provide a distance matrix before you can use the anosim command."); m->mothurOutEndLine(); abort = true; - } - }else { sharedfile = globaldata->getSharedFile(); } - - //use distance matrix if one is provided - if ((sharedfile != "") && (phylipfile != "")) { sharedfile = ""; } - - //check for optional parameter and set defaults - // ...at some point should added some additional type checking... - label = validParameter.validFile(parameters, "label", false); - if (label == "not found") { label = ""; } - else { - if(label != "all") { m->splitAtDash(label, labels); allLines = 0; } - else { allLines = 1; } - } - - //if the user has not specified any labels use the ones from read.otu - if (label == "") { - allLines = globaldata->allLines; - labels = globaldata->labels; - } - - groups = validParameter.validFile(parameters, "groups", false); - if (groups == "not found") { groups = ""; } - else { - m->splitAtDash(groups, Groups); - globaldata->Groups = Groups; - } - - sets = validParameter.validFile(parameters, "sets", false); - if (sets == "not found") { sets = ""; pickedGroups = false; } - else { - pickedGroups = true; - m->splitAtDash(sets, Sets); - } - - - string temp = validParameter.validFile(parameters, "iters", false); if (temp == "not found") { temp = "1000"; } + designFileName = validParameter.validFile(parameters, "design", true); + if (designFileName == "not open") { abort = true; } + else if (designFileName == "not found") { + designFileName = ""; + m->mothurOut("You must provide an design file."); + m->mothurOutEndLine(); + abort = true; + } + + string temp = validParameter.validFile(parameters, "iters", false); + if (temp == "not found") { temp = "1000"; } convert(temp, iters); - temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = "1"; } - convert(temp, processors); - - vector Estimators; - calc = validParameter.validFile(parameters, "calc", false); - if (calc == "not found") { calc = "morisitahorn"; } - m->splitAtDash(calc, Estimators); - - if (abort == false) { - - ValidCalculators* validCalculator = new ValidCalculators(); - - for (int i=0; iisValidCalculator("treegroup", Estimators[i]) == true) { - if (Estimators[i] == "sharedsobs") { - calculators.push_back(new SharedSobsCS()); - }else if (Estimators[i] == "sharedchao") { - calculators.push_back(new SharedChao1()); - }else if (Estimators[i] == "sharedace") { - calculators.push_back(new SharedAce()); - }else if (Estimators[i] == "jabund") { - calculators.push_back(new JAbund()); - }else if (Estimators[i] == "sorabund") { - calculators.push_back(new SorAbund()); - }else if (Estimators[i] == "jclass") { - calculators.push_back(new Jclass()); - }else if (Estimators[i] == "sorclass") { - calculators.push_back(new SorClass()); - }else if (Estimators[i] == "jest") { - calculators.push_back(new Jest()); - }else if (Estimators[i] == "sorest") { - calculators.push_back(new SorEst()); - }else if (Estimators[i] == "thetayc") { - calculators.push_back(new ThetaYC()); - }else if (Estimators[i] == "thetan") { - calculators.push_back(new ThetaN()); - }else if (Estimators[i] == "kstest") { - calculators.push_back(new KSTest()); - }else if (Estimators[i] == "sharednseqs") { - calculators.push_back(new SharedNSeqs()); - }else if (Estimators[i] == "ochiai") { - calculators.push_back(new Ochiai()); - }else if (Estimators[i] == "anderberg") { - calculators.push_back(new Anderberg()); - }else if (Estimators[i] == "kulczynski") { - calculators.push_back(new Kulczynski()); - }else if (Estimators[i] == "kulczynskicody") { - calculators.push_back(new KulczynskiCody()); - }else if (Estimators[i] == "lennon") { - calculators.push_back(new Lennon()); - }else if (Estimators[i] == "morisitahorn") { - calculators.push_back(new MorHorn()); - }else if (Estimators[i] == "braycurtis") { - calculators.push_back(new BrayCurtis()); - }else if (Estimators[i] == "whittaker") { - calculators.push_back(new Whittaker()); - }else if (Estimators[i] == "odum") { - calculators.push_back(new Odum()); - }else if (Estimators[i] == "canberra") { - calculators.push_back(new Canberra()); - }else if (Estimators[i] == "structeuclidean") { - calculators.push_back(new StructEuclidean()); - }else if (Estimators[i] == "structchord") { - calculators.push_back(new StructChord()); - }else if (Estimators[i] == "hellinger") { - calculators.push_back(new Hellinger()); - }else if (Estimators[i] == "manhattan") { - calculators.push_back(new Manhattan()); - }else if (Estimators[i] == "structpearson") { - calculators.push_back(new StructPearson()); - }else if (Estimators[i] == "soergel") { - calculators.push_back(new Soergel()); - }else if (Estimators[i] == "spearman") { - calculators.push_back(new Spearman()); - }else if (Estimators[i] == "structkulczynski") { - calculators.push_back(new StructKulczynski()); - }else if (Estimators[i] == "speciesprofile") { - calculators.push_back(new SpeciesProfile()); - }else if (Estimators[i] == "hamming") { - calculators.push_back(new Hamming()); - }else if (Estimators[i] == "structchi2") { - calculators.push_back(new StructChi2()); - }else if (Estimators[i] == "gower") { - calculators.push_back(new Gower()); - }else if (Estimators[i] == "memchi2") { - calculators.push_back(new MemChi2()); - }else if (Estimators[i] == "memchord") { - calculators.push_back(new MemChord()); - }else if (Estimators[i] == "memeuclidean") { - calculators.push_back(new MemEuclidean()); - }else if (Estimators[i] == "mempearson") { - calculators.push_back(new MemPearson()); - } - } - } - } + temp = validParameter.validFile(parameters, "alpha", false); + if (temp == "not found") { temp = "0.05"; } + convert(temp, experimentwiseAlpha); } } @@ -318,19 +153,14 @@ AnosimCommand::AnosimCommand(string option) { void AnosimCommand::help(){ try { - m->mothurOut("The anosim command can only be executed after a successful read.otu command of a list and group or shared file, or by providing a phylip formatted distance matrix.\n"); + m->mothurOut("Referenced: Clarke, K. R. (1993). Non-parametric multivariate analysis of changes in community structure. _Australian Journal of Ecology_ 18, 117-143.\n"); m->mothurOut("The anosim command outputs a .anosim file. \n"); - m->mothurOut("The anosim command parameters are phylip, iters, groups, label, design, sets and processors. The design parameter is required.\n"); - m->mothurOut("The design parameter allows you to assign your groups to sets when you are running anosim. It is required. \n"); - m->mothurOut("The design file looks like the group file. It is a 2 column tab delimited file, where the first column is the group name and the second column is the set the group belongs to.\n"); - m->mothurOut("The sets parameter allows you to specify which of the sets in your designfile you would like to analyze. The set names are separated by dashes. THe default is all sets in the designfile. To run the pairwise comparisons of all sets use sets=all.\n"); + m->mothurOut("The anosim command parameters are phylip, iters, and alpha. The phylip and design parameters are required.\n"); + m->mothurOut("The design parameter allows you to assign your samples to groups when you are running anosim. It is required. \n"); + m->mothurOut("The design file looks like the group file. It is a 2 column tab delimited file, where the first column is the sample name and the second column is the group the sample belongs to.\n"); m->mothurOut("The iters parameter allows you to set number of randomization for the P value. The default is 1000. \n"); - m->mothurOut("The groups parameter allows you to specify which of the groups you would like included. The group names are separated by dashes. groups=all will find all pairwise comparisons. \n"); - m->mothurOut("The label parameter allows you to select what distance levels you would like, and are also separated by dashes.\n"); - m->mothurOut("The processors parameter allows you to specify how many processors you would like to use. The default is 1. \n"); - m->mothurOut("The anosim command should be in the following format: anosim(design=yourDesignFile).\n"); - m->mothurOut("Example anosim(design=temp.design, groups=A-B-C).\n"); - m->mothurOut("Note: No spaces between parameter labels (i.e. groups), '=' and parameters (i.e.yourGroups).\n\n"); + m->mothurOut("The anosim command should be in the following format: anosim(phylip=file.dist, design=file.design).\n"); + m->mothurOut("Note: No spaces between parameter labels (i.e. iters), '=' and parameters (i.e. 1000).\n\n"); } catch(exception& e) { @@ -351,223 +181,91 @@ int AnosimCommand::execute(){ if (abort == true) { if (calledHelp) { return 0; } return 2; } //read design file - designMap = new GroupMap(designfile); + designMap = new GroupMap(designFileName); designMap->readDesignMap(); - //make sure sets are all in designMap - SharedUtil* util = new SharedUtil(); - util->setGroups(Sets, designMap->namesOfGroups); - delete util; - - //if the user pickedGroups or set sets=all, then we want to figure out how to split the pairwise comparison between processors - int numGroups = Sets.size(); - if (pickedGroups) { - for (int a=0; a groups; groups.push_back(Sets[a]); groups.push_back(Sets[l]); - namesOfGroupCombos.push_back(groups); - } - } - }else { //one giant compare - vector groups; - for (int a=0; ahasPath(phylipFileName); } - //only 1 combo - if (numGroups == 2) { processors = 1; } - else if (numGroups < 2) { m->mothurOut("Not enough sets, I need at least 2 valid sets. Unable to complete command."); m->mothurOutEndLine(); m->control_pressed = true; } + //read in distance matrix and square it + ReadPhylipVector readMatrix(phylipFileName); + vector sampleNames = readMatrix.read(distanceMatrix); -#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) - if(processors != 1){ - int numPairs = namesOfGroupCombos.size(); - int numPairsPerProcessor = numPairs / processors; - - for (int i = 0; i < processors; i++) { - int startPos = i * numPairsPerProcessor; - if(i == processors - 1){ - numPairsPerProcessor = numPairs - i * numPairsPerProcessor; - } - lines.push_back(linePair(startPos, numPairsPerProcessor)); + for(int i=0;ihasPath(sharedfile); } - - //for each calculator - for(int i = 0 ; i < calculators.size(); i++) { - - //create a new filename - ofstream out; - string outputFile = outputDir + m->getRootName(m->getSimpleName(sharedfile)) + calculators[i]->getName() + ".anosim"; - m->openOutputFile(outputFile, out); - outputNames.push_back(outputFile); outputTypes["anosim"].push_back(outputFile); - - //print headers - out << "label\tgroupsCompared\tRValue\tpValue" << endl; - m->mothurOut("\nlabel\tgroupsCompared\tRValue\tpValue"); m->mothurOutEndLine(); - out.close(); - } - - InputData input(sharedfile, "sharedfile"); - vector lookup = input.getSharedRAbundVectors(); - string lastLabel = lookup[0]->getLabel(); - - //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label. - set processedLabels; - set userLabels = labels; - - //sanity check between shared file groups and design file groups - for (int i = 0; i < lookup.size(); i++) { - string group = designMap->getGroup(lookup[i]->getGroup()); - - if (group == "not found") { m->control_pressed = true; m->mothurOut("[ERROR]: " + lookup[i]->getGroup() + " is not in your design file, please correct."); m->mothurOutEndLine(); } - } - - //as long as you are not at the end of the file or done wih the lines you want - while((lookup[0] != NULL) && ((allLines == 1) || (userLabels.size() != 0))) { - - if (m->control_pressed) { for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } globaldata->Groups.clear(); delete designMap; for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; } - - if(allLines == 1 || labels.count(lookup[0]->getLabel()) == 1){ - - m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine(); - process(lookup); - - processedLabels.insert(lookup[0]->getLabel()); - userLabels.erase(lookup[0]->getLabel()); - } - - if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) { - string saveLabel = lookup[0]->getLabel(); - - for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } - lookup = input.getSharedRAbundVectors(lastLabel); - m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine(); - - process(lookup); + //link designMap to rows/columns in distance matrix + map > origGroupSampleMap; + for(int i=0;igetGroup(sampleNames[i])].push_back(i); + } + int numGroups = origGroupSampleMap.size(); + + //create a new filename + ofstream ANOSIMFile; + string ANOSIMFileName = outputDir + m->getRootName(m->getSimpleName(phylipFileName)) + "anosim"; + m->openOutputFile(ANOSIMFileName, ANOSIMFile); + outputNames.push_back(ANOSIMFileName); outputTypes["anosim"].push_back(ANOSIMFileName); + m->mothurOut("\ncomparison\tR-value\tP-value\n"); + ANOSIMFile << "comparison\tR-value\tP-value\n"; + + + double fullANOSIMPValue = runANOSIM(ANOSIMFile, distanceMatrix, origGroupSampleMap, experimentwiseAlpha); + + + if(fullANOSIMPValue <= experimentwiseAlpha && numGroups > 2){ + + int numCombos = numGroups * (numGroups-1) / 2; + double pairwiseAlpha = experimentwiseAlpha / (double) numCombos; + + for(map >::iterator itA=origGroupSampleMap.begin();itA!=origGroupSampleMap.end();itA++){ + map >::iterator itB = itA; + itB++; + for(itB;itB!=origGroupSampleMap.end();itB++){ - processedLabels.insert(lookup[0]->getLabel()); - userLabels.erase(lookup[0]->getLabel()); + map > subGroupSampleMap; - //restore real lastlabel to save below - lookup[0]->setLabel(saveLabel); - } - - lastLabel = lookup[0]->getLabel(); - //prevent memory leak - for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; lookup[i] = NULL; } - - if (m->control_pressed) { globaldata->Groups.clear(); delete designMap; for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; } - - //get next line to process - lookup = input.getSharedRAbundVectors(); - } - - if (m->control_pressed) { globaldata->Groups.clear(); delete designMap; for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; } - - //output error messages about any remaining user labels - set::iterator it; - bool needToRun = false; - for (it = userLabels.begin(); it != userLabels.end(); it++) { - m->mothurOut("Your file does not include the label " + *it); - if (processedLabels.count(lastLabel) != 1) { - m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine(); - needToRun = true; - }else { - m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine(); - } - } - - //run last label if you need to - if (needToRun == true) { - for (int i = 0; i < lookup.size(); i++) { if (lookup[i] != NULL) { delete lookup[i]; } } - lookup = input.getSharedRAbundVectors(lastLabel); - - m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine(); - - process(lookup); - - for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } - } - - //reset groups parameter - globaldata->Groups.clear(); - - - }else { //user provided distance matrix - - if (outputDir == "") { outputDir = m->hasPath(phylipfile); } + subGroupSampleMap[itA->first] = itA->second; string groupA = itA->first; + subGroupSampleMap[itB->first] = itB->second; string groupB = itB->first; - //create a new filename - ofstream out; - string outputFile = outputDir + m->getRootName(m->getSimpleName(phylipfile)) + "anosim"; - m->openOutputFile(outputFile, out); - outputNames.push_back(outputFile); outputTypes["anosim"].push_back(outputFile); - - //print headers - out << "groupsCompared\tRValue\tpValue" << endl; - m->mothurOut("\ngroupsCompared\tRValue\tpValue"); m->mothurOutEndLine(); - out.close(); - - ReadPhylipVector readMatrix(phylipfile); - vector< vector > completeMatrix; - vector names = readMatrix.read(completeMatrix); - -#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) - if(processors == 1){ - driver(0, namesOfGroupCombos.size(), names, "", completeMatrix); - }else{ - int process = 1; - vector processIDS; - - //loop through and create all the processes you want - while (process != processors) { - int pid = fork(); + vector subIndices; + for(map >::iterator it=subGroupSampleMap.begin();it!=subGroupSampleMap.end();it++){ + subIndices.insert(subIndices.end(), it->second.begin(), it->second.end()); + } + int subNumSamples = subIndices.size(); + + sort(subIndices.begin(), subIndices.end()); - if (pid > 0) { - processIDS.push_back(pid); - process++; - }else if (pid == 0){ - driver(lines[process].start, lines[process].num, names, toString(getpid()), completeMatrix); - exit(0); - }else { - m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine(); - for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); } - exit(0); + vector > subDistMatrix(distanceMatrix.size()); + for(int i=0;igetRootName(m->getSimpleName(phylipfile)) + "anosim"; - for (int j = 0; j < processIDS.size(); j++) { - m->appendFiles((outputFile + toString(processIDS[j])), outputFile); - remove((outputFile + toString(processIDS[j])).c_str()); - } - } -#else - driver(0, namesOfGroupCombos.size(), names, "", completeMatrix); -#endif + m->mothurOut("\nExperiment-wise error rate: " + toString(experimentwiseAlpha) + '\n'); + m->mothurOut("Pair-wise error rate (Bonferroni): " + toString(pairwiseAlpha) + '\n'); + } + else{ + m->mothurOut("\nExperiment-wise error rate: " + toString(experimentwiseAlpha) + '\n'); } + m->mothurOut("If you have borderline P-values, you should try increasing the number of iterations\n"); + ANOSIMFile.close(); + delete designMap; - - if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; } - + m->mothurOutEndLine(); m->mothurOut("Output File Names: "); m->mothurOutEndLine(); for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); } @@ -582,366 +280,204 @@ int AnosimCommand::execute(){ } //********************************************************************************************************************** -int AnosimCommand::process(vector thisLookup) { - try{ +double AnosimCommand::runANOSIM(ofstream& ANOSIMFile, vector > dMatrix, map > groupSampleMap, double alpha) { + try { + -#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) - if(processors == 1){ - driver(0, namesOfGroupCombos.size(), thisLookup, ""); - }else{ - int process = 1; - vector processIDS; - - //loop through and create all the processes you want - while (process != processors) { - int pid = fork(); - - if (pid > 0) { - processIDS.push_back(pid); - process++; - }else if (pid == 0){ - driver(lines[process].start, lines[process].num, thisLookup, toString(getpid())); - exit(0); - }else { - m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine(); - for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); } - exit(0); - } - } - - //do my part - driver(lines[0].start, lines[0].num, thisLookup, ""); - - //force parent to wait until all the processes are done - for (int i=0;i<(processors-1);i++) { - int temp = processIDS[i]; - wait(&temp); - } - - //append files - for(int i = 0 ; i < calculators.size(); i++) { - string outputFile = outputDir + m->getRootName(m->getSimpleName(sharedfile)) + calculators[i]->getName() + ".anosim"; - - for (int j = 0; j < processIDS.size(); j++) { - m->appendFiles((outputFile + toString(processIDS[j])), outputFile); - remove((outputFile + toString(processIDS[j])).c_str()); - } - } + vector > rankMatrix = convertToRanks(dMatrix); + double RValue = calcR(rankMatrix, groupSampleMap); + + int pCount = 0; + for(int i=0;i > randGroupSampleMap = getRandomizedGroups(groupSampleMap); + double RValueRand = calcR(rankMatrix, randGroupSampleMap); + if(RValue <= RValueRand){ pCount++; } } -#else - driver(0, namesOfGroupCombos.size(), thisLookUp, ""); -#endif + + double pValue = (double)pCount / (double) iters; + string pString = ""; + if(pValue < 1/(double)iters){ pString = '<' + toString(1/(double)iters); } + else { pString = toString(pValue); } - return 0; + + map >::iterator it=groupSampleMap.begin(); + m->mothurOut(it->first); + ANOSIMFile << it->first; + it++; + for(it;it!=groupSampleMap.end();it++){ + m->mothurOut('-' + it->first); + ANOSIMFile << '-' << it->first; + + } + m->mothurOut('\t' + toString(RValue) + '\t' + pString); + ANOSIMFile << '\t' << RValue << '\t' << pString; + + if(pValue < alpha){ + ANOSIMFile << "*"; + m->mothurOut("*"); + } + ANOSIMFile << endl; + m->mothurOutEndLine(); + + return pValue; } catch(exception& e) { - m->errorOut(e, "AnosimCommand", "process"); + m->errorOut(e, "AnosimCommand", "calcAnisom"); exit(1); } } + //********************************************************************************************************************** -int AnosimCommand::driver(int start, int num, vector thisLookup, string pidValue) { +double AnosimCommand::calcR(vector > rankMatrix, map > groupSampleMap){ try { - vector subset; - EstOutput data; + + int numSamples = 0; + for(map >::iterator it=groupSampleMap.begin();it!=groupSampleMap.end();it++){ + numSamples += it->second.size(); + } - //for each calculator - for(int i = 0 ; i < calculators.size(); i++) { - - //create a new filename - ofstream out; - string outputFile = outputDir + m->getRootName(m->getSimpleName(sharedfile)) + calculators[i]->getName() + ".anosim" + pidValue; - m->openOutputFileAppend(outputFile, out); - out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint); - cout.setf(ios::fixed, ios::floatfield); cout.setf(ios::showpoint); - - //for each combo - for (int c = start; c < (start+num); c++) { - - if (m->control_pressed) { out.close(); return 0; } - - //get set names - vector setNames; - for (int j = 0; j < namesOfGroupCombos[c].size(); j++) { setNames.push_back(namesOfGroupCombos[c][j]); } - - vector thisCombosLookup; - vector thisCombosLookupSets; //what set each sharedRabund is from to be used when calculating SSWithin - for (int k = 0; k < thisLookup.size(); k++) { - string thisGroup = thisLookup[k]->getGroup(); - - //is this group for a set we want to compare?? - if (m->inUsersGroups(designMap->getGroup(thisGroup), setNames)) { - thisCombosLookup.push_back(thisLookup[k]); - thisCombosLookupSets.push_back(designMap->getGroup(thisGroup)); - } - - } - - int numGroups = thisCombosLookup.size(); - - //calc the distance matrix - matrix.clear(); - matrix.resize(numGroups); - for (int k = 0; k < matrix.size(); k++) { for (int j = 0; j < matrix.size(); j++) { matrix[k].push_back(1.0); } } - - if (thisCombosLookup.size() == 0) { - m->mothurOut("[ERROR]: Missing shared info for sets. Skipping comparison."); m->mothurOutEndLine(); - }else{ - - out << thisLookup[0]->getLabel() << '\t'; - m->mothurOut(thisLookup[0]->getLabel() + "\t"); - if (setNames.size() == 2) { out << setNames[0] << "-" << setNames[1] << '\t'; m->mothurOut(setNames[0] + "-" + setNames[1] + "\t"); } - else { out << "all" << '\t'; m->mothurOut("all\t"); } - - for (int k = 0; k < thisCombosLookup.size(); k++) { - for (int l = k; l < thisCombosLookup.size(); l++) { - - if (m->control_pressed) { out.close(); return 0; } - - if (k != l) { //we dont need to similiarity of a groups to itself - //get estimated similarity between 2 groups - subset.clear(); //clear out old pair of sharedrabunds - //add new pair of sharedrabunds - subset.push_back(thisCombosLookup[k]); subset.push_back(thisCombosLookup[l]); - - //if this calc needs all groups to calculate the pair load all groups - if (calculators[i]->getNeedsAll()) { - //load subset with rest of lookup for those calcs that need everyone to calc for a pair - for (int w = 0; w < thisCombosLookup.size(); w++) { - if ((w != k) && (w != l)) { subset.push_back(thisCombosLookup[w]); } - } - } - - data = calculators[i]->getValues(subset); //saves the calculator outputs - - //save values in similarity matrix - matrix[k][l] = 1.0 - data[0]; - matrix[l][k] = 1.0 - data[0]; - } - } - } - - //calc anosim - calcAnosim(out, setNames.size(), thisCombosLookupSets); + + double within = 0.0; + int numWithinComps = 0; + + for(map >::iterator it=groupSampleMap.begin();it!=groupSampleMap.end();it++){ + vector indices = it->second; + for(int i=0;i indices[j]) { within += rankMatrix[indices[i]][indices[j]]; } + else { within += rankMatrix[indices[j]][indices[i]]; } + numWithinComps++; } } - - out.close(); - } + } - return 0; + within /= (float) numWithinComps; - } - catch(exception& e) { - m->errorOut(e, "AnosimCommand", "driver"); - exit(1); - } -} -//********************************************************************************************************************** + double between = 0.0; + int numBetweenComps = 0; -int AnosimCommand::driver(int start, int num, vector names, string pidValue, vector< vector >& completeMatrix) { - try { + map >::iterator itB; - //create a new filename - ofstream out; - string outputFile = outputDir + m->getRootName(m->getSimpleName(phylipfile)) + "anosim" + pidValue; - m->openOutputFileAppend(outputFile, out); - out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint); - cout.setf(ios::fixed, ios::floatfield); cout.setf(ios::showpoint); - - //for each combo - for (int c = start; c < (start+num); c++) { - - if (m->control_pressed) { out.close(); return 0; } - - //get set names - vector setNames; - for (int j = 0; j < namesOfGroupCombos[c].size(); j++) { setNames.push_back(namesOfGroupCombos[c][j]); } - - vector thisCombosSets; //what set each line in the distance matrix is from to be used when calculating SSWithin - set indexes; - for (int k = 0; k < names.size(); k++) { - //is this group for a set we want to compare?? - if (m->inUsersGroups(designMap->getGroup(names[k]), setNames)) { - thisCombosSets.push_back(designMap->getGroup(names[k])); - indexes.insert(k); //save indexes of valid rows in matrix for submatrix - } - } - - int numGroups = thisCombosSets.size(); - - //calc the distance matrix - matrix.clear(); - matrix.resize(numGroups); - - for (int k = 0; k < matrix.size(); k++) { for (int j = 0; j < matrix.size(); j++) { matrix[k].push_back(1.0); } } - - if (thisCombosSets.size() == 0) { - m->mothurOut("[ERROR]: Missing distance info for sets. Skipping comparison."); m->mothurOutEndLine(); - }else{ - - if (setNames.size() == 2) { out << setNames[0] << "-" << setNames[1] << '\t'; m->mothurOut(setNames[0] + "-" + setNames[1] + "\t"); } - else { out << "all" << '\t'; m->mothurOut("all\t"); } - - //fill submatrix - int rowCount = 0; - for (int j = 0; j < completeMatrix.size(); j++) { - - if (indexes.count(j) != 0) { //we want at least part of this row - int columnCount = 0; - for (int k = 0; k < completeMatrix[j].size(); k++) { - - if (indexes.count(k) != 0) { //we want this distance - matrix[rowCount][columnCount] = completeMatrix[j][k]; - matrix[columnCount][rowCount] = completeMatrix[j][k]; - columnCount++; - } - } - rowCount++; - } + for(map >::iterator itA=groupSampleMap.begin();itA!=groupSampleMap.end();itA++){ + + for(int i=0;isecond.size();i++){ + int A = itA->second[i]; + map >::iterator itB = itA; + itB++; + for(itB;itB!=groupSampleMap.end();itB++){ + for(int j=0;jsecond.size();j++){ + int B = itB->second[j]; + if(A>B) { between += rankMatrix[A][B]; } + else { between += rankMatrix[B][A]; } + numBetweenComps++; + } } - - //calc anosim - calcAnosim(out, setNames.size(), thisCombosSets); } } - out.close(); + between /= (float) numBetweenComps; - return 0; - + double Rvalue = (between - within)/(numSamples * (numSamples-1) / 4.0); + + return Rvalue; } catch(exception& e) { - m->errorOut(e, "AnosimCommand", "driver"); + m->errorOut(e, "AnosimCommand", "calcWithinBetween"); exit(1); } } + //********************************************************************************************************************** -int AnosimCommand::calcAnosim(ofstream& out, int numTreatments, vector thisCombosLookupSets) { + +vector > AnosimCommand::convertToRanks(vector > dist) { try { - //rank distances - vector rankMatrix = convertToRanks(); - - double meanBetween, meanWithin, RValue, pValue; + vector cells; + vector > ranks = dist; - meanWithin = calcWithinBetween(rankMatrix, thisCombosLookupSets, meanBetween); - - int N = matrix.size(); - int div = (N * (N-1) / 4); - - //calc RValue - RValue = (meanBetween - meanWithin) / (double) div; - - //calc Pvalue - int count = 0; - for (int i = 0; i < iters; i++) { - if (m->control_pressed) { break; } - - //randomly shuffle names to randomize the matrix - vector copyNames = thisCombosLookupSets; - random_shuffle(copyNames.begin(), copyNames.end()); - - meanWithin = calcWithinBetween(rankMatrix, copyNames, meanBetween); - - //calc RValue - double randomRValue = (meanBetween - meanWithin) / (double) div; - - if (randomRValue >= RValue) { count++; } + for (int i = 0; i < dist.size(); i++) { + for (int j = 0; j < i; j++) { + if(dist[i][j] != -1){ + seqDist member(i, j, dist[i][j]); + cells.push_back(member); + } + } } - pValue = count / (float) iters; - - out << RValue << '\t' << pValue << endl; - cout << RValue << '\t' << pValue << endl; - m->mothurOutJustToLog(toString(RValue) + "\t" + toString(pValue) + "\n"); - return 0; - } - catch(exception& e) { - m->errorOut(e, "AnosimCommand", "calcAnisom"); - exit(1); - } -} -//********************************************************************************************************************** -double AnosimCommand::calcWithinBetween(vector& thisMatrix, vector thisCombosLookupSets, double& between) { - try { - double within = 0.0; - int count = 0; - int count2 = 0; - between = 0.0; - - for (int l = 0; l < thisMatrix.size(); l++) { - //if you are from the same treatment - if (thisCombosLookupSets[thisMatrix[l].seq1] == thisCombosLookupSets[thisMatrix[l].seq2]) { - within += thisMatrix[l].dist; //rank of this distance - count++; - }else { //different treatments - between += thisMatrix[l].dist; //rank of this distance - count2++; + //sort distances + sort(cells.begin(), cells.end(), compareSequenceDistance); + + //find ranks of distances + int index = 0; + int indexSum = 0; + for(int i=0;ierrorOut(e, "AnosimCommand", "calcWithinBetween"); + m->errorOut(e, "AnosimCommand", "convertToRanks"); exit(1); } } + //********************************************************************************************************************** -vector AnosimCommand::convertToRanks() { - try { - vector ranks; - - for (int i = 0; i < matrix.size(); i++) { - for (int j = 0; j < i; j++) { - seqDist member(i, j, matrix[i][j]); - ranks.push_back(member); - } + +map > AnosimCommand::getRandomizedGroups(map > origMapping){ + try{ + vector sampleIndices; + vector samplesPerGroup; + + map >::iterator it; + for(it=origMapping.begin();it!=origMapping.end();it++){ + vector indices = it->second; + samplesPerGroup.push_back(indices.size()); + sampleIndices.insert(sampleIndices.end(), indices.begin(), indices.end()); } - //sort distances - sort(ranks.begin(), ranks.end(), compareSequenceDistance); + random_shuffle(sampleIndices.begin(), sampleIndices.end()); - //find ranks of distances - vector ties; - int rankTotal = 0; - for (int j = 0; j < ranks.size(); j++) { - rankTotal += (j+1); - ties.push_back(&ranks[j]); - - if (j != (ranks.size()-1)) { // you are not the last so you can look ahead - if (ranks[j].dist != ranks[j+1].dist) { // you are done with ties, rank them and continue - - for (int k = 0; k < ties.size(); k++) { - float thisrank = rankTotal / (float) ties.size(); - (*ties[k]).dist = thisrank; - } - ties.clear(); - rankTotal = 0; - } - }else { // you are the last one - - for (int k = 0; k < ties.size(); k++) { - float thisrank = rankTotal / (float) ties.size(); - (*ties[k]).dist = thisrank; - } + int index = 0; + map > randomizedGroups = origMapping; + for(it=randomizedGroups.begin();it!=randomizedGroups.end();it++){ + for(int i=0;isecond.size();i++){ + it->second[i] = sampleIndices[index++]; } } - return ranks; + return randomizedGroups; } - catch(exception& e) { - m->errorOut(e, "AnosimCommand", "convertToRanks"); + catch (exception& e) { + m->errorOut(e, "AnosimCommand", "randomizeGroups"); exit(1); } } + //********************************************************************************************************************** diff --git a/anosimcommand.h b/anosimcommand.h index 0b14563..ab72400 100644 --- a/anosimcommand.h +++ b/anosimcommand.h @@ -12,12 +12,8 @@ #include "command.hpp" -#include "inputdata.h" -#include "sharedrabundvector.h" -#include "validcalculator.h" -#include "readphylipvector.h" -class GlobalData; +class GroupMap; class AnosimCommand : public Command { @@ -33,32 +29,23 @@ public: void help(); private: - struct linePair { - int start; - int num; - linePair(int i, int j) : start(i), num(j) {} - }; - vector lines; - - GlobalData* globaldata; + bool abort; GroupMap* designMap; map > outputTypes; + string outputDir, inputDir, designFileName, phylipFileName; + + vector > convertToRanks(vector >); + double calcR(vector >, map >); + map > getRandomizedGroups(map >); + double runANOSIM(ofstream&, vector >, map >, double); - vector< vector > matrix; - bool abort, allLines, pickedGroups; - set labels; //holds labels to be used - string format, groups, label, outputDir, inputDir, designfile, sets, phylipfile, calc, sharedfile; - vector Groups, outputNames, Sets; + vector< vector > distanceMatrix; + vector outputNames; + int iters; + double experimentwiseAlpha; vector< vector > namesOfGroupCombos; - int iters, processors; - vector calculators; - int driver(int, int, vector, string); - int driver(int, int, vector, string, vector< vector >&); - int process(vector); - int calcAnosim(ofstream&, int, vector); - double calcWithinBetween(vector&, vector, double&); - vector convertToRanks(); + }; #endif diff --git a/homovacommand.cpp b/homovacommand.cpp index 28f560f..52456b4 100644 --- a/homovacommand.cpp +++ b/homovacommand.cpp @@ -8,53 +8,14 @@ */ #include "homovacommand.h" -#include "sharedutilities.h" -#include "sharedsobscollectsummary.h" -#include "sharedchao1.h" -#include "sharedace.h" -#include "sharednseqs.h" -#include "sharedjabund.h" -#include "sharedsorabund.h" -#include "sharedjclass.h" -#include "sharedsorclass.h" -#include "sharedjest.h" -#include "sharedsorest.h" -#include "sharedthetayc.h" -#include "sharedthetan.h" -#include "sharedkstest.h" -#include "whittaker.h" -#include "sharedochiai.h" -#include "sharedanderbergs.h" -#include "sharedkulczynski.h" -#include "sharedkulczynskicody.h" -#include "sharedlennon.h" -#include "sharedmorisitahorn.h" -#include "sharedbraycurtis.h" -#include "sharedjackknife.h" -#include "whittaker.h" -#include "odum.h" -#include "canberra.h" -#include "structeuclidean.h" -#include "structchord.h" -#include "hellinger.h" -#include "manhattan.h" -#include "structpearson.h" -#include "soergel.h" -#include "spearman.h" -#include "structkulczynski.h" -#include "structchi2.h" -#include "speciesprofile.h" -#include "hamming.h" -#include "gower.h" -#include "memchi2.h" -#include "memchord.h" -#include "memeuclidean.h" -#include "mempearson.h" +#include "groupmap.h" +#include "readphylipvector.h" //********************************************************************************************************************** + vector HomovaCommand::getValidParameters(){ try { - string Array[] = {"groups","label","outputdir","iters","phylip","design","sets","processors","inputdir"}; + string Array[] = {"outputdir","iters","phylip","design","alpha", "inputdir"}; vector myArray (Array, Array+(sizeof(Array)/sizeof(string))); return myArray; } @@ -63,7 +24,9 @@ vector HomovaCommand::getValidParameters(){ exit(1); } } + //********************************************************************************************************************** + HomovaCommand::HomovaCommand(){ try { abort = true; calledHelp = true; @@ -75,7 +38,9 @@ HomovaCommand::HomovaCommand(){ exit(1); } } + //********************************************************************************************************************** + vector HomovaCommand::getRequiredParameters(){ try { string Array[] = {"design"}; @@ -87,7 +52,9 @@ vector HomovaCommand::getRequiredParameters(){ exit(1); } } + //********************************************************************************************************************** + vector HomovaCommand::getRequiredFiles(){ try { string Array[] = {}; @@ -99,21 +66,19 @@ vector HomovaCommand::getRequiredFiles(){ exit(1); } } + //********************************************************************************************************************** HomovaCommand::HomovaCommand(string option) { try { - globaldata = GlobalData::getInstance(); abort = false; calledHelp = false; - allLines = 1; - labels.clear(); //allow user to run help if(option == "help") { help(); abort = true; calledHelp = true; } else { //valid paramters for this command - string AlignArray[] = {"groups","label","outputdir","iters","phylip","design","sets","processors","inputdir"}; + string AlignArray[] = {"design","outputdir","iters","phylip","alpha", "inputdir"}; vector myArray (AlignArray, AlignArray+(sizeof(AlignArray)/sizeof(string))); OptionParser parser(option); @@ -156,155 +121,33 @@ HomovaCommand::HomovaCommand(string option) { } } - phylipfile = validParameter.validFile(parameters, "phylip", true); - if (phylipfile == "not open") { phylipfile = ""; abort = true; } - else if (phylipfile == "not found") { phylipfile = ""; } - else { globaldata->newRead(); format = "phylip"; globaldata->setPhylipFile(phylipfile); } + phylipFileName = validParameter.validFile(parameters, "phylip", true); + if (phylipFileName == "not open") { phylipFileName = ""; abort = true; } + else if (phylipFileName == "not found") { phylipFileName = ""; } + else if (designFileName == "not found") { + designFileName = ""; + m->mothurOut("You must provide an phylip file."); + m->mothurOutEndLine(); + abort = true; + } //check for required parameters - designfile = validParameter.validFile(parameters, "design", true); - if (designfile == "not open") { abort = true; } - else if (designfile == "not found") { designfile = ""; m->mothurOut("You must provide an design file."); m->mothurOutEndLine(); abort = true; } - - //make sure the user has already run the read.otu command - if ((globaldata->getSharedFile() == "")) { - if ((phylipfile == "")) { - m->mothurOut("You must read a list and a group, a shared file or provide a distance matrix before you can use the homova command."); m->mothurOutEndLine(); abort = true; - } - }else { sharedfile = globaldata->getSharedFile(); } - - //use distance matrix if one is provided - if ((sharedfile != "") && (phylipfile != "")) { sharedfile = ""; } - - //check for optional parameter and set defaults - // ...at some point should added some additional type checking... - label = validParameter.validFile(parameters, "label", false); - if (label == "not found") { label = ""; } - else { - if(label != "all") { m->splitAtDash(label, labels); allLines = 0; } - else { allLines = 1; } - } - - //if the user has not specified any labels use the ones from read.otu - if (label == "") { - allLines = globaldata->allLines; - labels = globaldata->labels; - } - - groups = validParameter.validFile(parameters, "groups", false); - if (groups == "not found") { groups = ""; } - else { - m->splitAtDash(groups, Groups); - globaldata->Groups = Groups; - } - - sets = validParameter.validFile(parameters, "sets", false); - if (sets == "not found") { sets = ""; pickedGroups = false; } - else { - pickedGroups = true; - m->splitAtDash(sets, Sets); - } - - - string temp = validParameter.validFile(parameters, "iters", false); if (temp == "not found") { temp = "1000"; } + designFileName = validParameter.validFile(parameters, "design", true); + if (designFileName == "not open") { abort = true; } + else if (designFileName == "not found") { + designFileName = ""; + m->mothurOut("You must provide an design file."); + m->mothurOutEndLine(); + abort = true; + } + + string temp = validParameter.validFile(parameters, "iters", false); + if (temp == "not found") { temp = "1000"; } convert(temp, iters); - temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = "1"; } - convert(temp, processors); - - vector Estimators; - calc = validParameter.validFile(parameters, "calc", false); - if (calc == "not found") { calc = "morisitahorn"; } - m->splitAtDash(calc, Estimators); - - if (abort == false) { - - ValidCalculators* validCalculator = new ValidCalculators(); - - for (int i=0; iisValidCalculator("treegroup", Estimators[i]) == true) { - if (Estimators[i] == "sharedsobs") { - calculators.push_back(new SharedSobsCS()); - }else if (Estimators[i] == "sharedchao") { - calculators.push_back(new SharedChao1()); - }else if (Estimators[i] == "sharedace") { - calculators.push_back(new SharedAce()); - }else if (Estimators[i] == "jabund") { - calculators.push_back(new JAbund()); - }else if (Estimators[i] == "sorabund") { - calculators.push_back(new SorAbund()); - }else if (Estimators[i] == "jclass") { - calculators.push_back(new Jclass()); - }else if (Estimators[i] == "sorclass") { - calculators.push_back(new SorClass()); - }else if (Estimators[i] == "jest") { - calculators.push_back(new Jest()); - }else if (Estimators[i] == "sorest") { - calculators.push_back(new SorEst()); - }else if (Estimators[i] == "thetayc") { - calculators.push_back(new ThetaYC()); - }else if (Estimators[i] == "thetan") { - calculators.push_back(new ThetaN()); - }else if (Estimators[i] == "kstest") { - calculators.push_back(new KSTest()); - }else if (Estimators[i] == "sharednseqs") { - calculators.push_back(new SharedNSeqs()); - }else if (Estimators[i] == "ochiai") { - calculators.push_back(new Ochiai()); - }else if (Estimators[i] == "anderberg") { - calculators.push_back(new Anderberg()); - }else if (Estimators[i] == "kulczynski") { - calculators.push_back(new Kulczynski()); - }else if (Estimators[i] == "kulczynskicody") { - calculators.push_back(new KulczynskiCody()); - }else if (Estimators[i] == "lennon") { - calculators.push_back(new Lennon()); - }else if (Estimators[i] == "morisitahorn") { - calculators.push_back(new MorHorn()); - }else if (Estimators[i] == "braycurtis") { - calculators.push_back(new BrayCurtis()); - }else if (Estimators[i] == "whittaker") { - calculators.push_back(new Whittaker()); - }else if (Estimators[i] == "odum") { - calculators.push_back(new Odum()); - }else if (Estimators[i] == "canberra") { - calculators.push_back(new Canberra()); - }else if (Estimators[i] == "structeuclidean") { - calculators.push_back(new StructEuclidean()); - }else if (Estimators[i] == "structchord") { - calculators.push_back(new StructChord()); - }else if (Estimators[i] == "hellinger") { - calculators.push_back(new Hellinger()); - }else if (Estimators[i] == "manhattan") { - calculators.push_back(new Manhattan()); - }else if (Estimators[i] == "structpearson") { - calculators.push_back(new StructPearson()); - }else if (Estimators[i] == "soergel") { - calculators.push_back(new Soergel()); - }else if (Estimators[i] == "spearman") { - calculators.push_back(new Spearman()); - }else if (Estimators[i] == "structkulczynski") { - calculators.push_back(new StructKulczynski()); - }else if (Estimators[i] == "speciesprofile") { - calculators.push_back(new SpeciesProfile()); - }else if (Estimators[i] == "hamming") { - calculators.push_back(new Hamming()); - }else if (Estimators[i] == "structchi2") { - calculators.push_back(new StructChi2()); - }else if (Estimators[i] == "gower") { - calculators.push_back(new Gower()); - }else if (Estimators[i] == "memchi2") { - calculators.push_back(new MemChi2()); - }else if (Estimators[i] == "memchord") { - calculators.push_back(new MemChord()); - }else if (Estimators[i] == "memeuclidean") { - calculators.push_back(new MemEuclidean()); - }else if (Estimators[i] == "mempearson") { - calculators.push_back(new MemPearson()); - } - } - } - } + temp = validParameter.validFile(parameters, "alpha", false); + if (temp == "not found") { temp = "0.05"; } + convert(temp, experimentwiseAlpha); } } @@ -319,20 +162,13 @@ HomovaCommand::HomovaCommand(string option) { void HomovaCommand::help(){ try { m->mothurOut("Referenced: Stewart CN, Excoffier L (1996). Assessing population genetic structure and variability with RAPD data: Application to Vaccinium macrocarpon (American Cranberry). J Evol Biol 9: 153-71.\n"); - m->mothurOut("The homova command can only be executed after a successful read.otu command of a list and group or shared file, or by providing a phylip formatted distance matrix.\n"); m->mothurOut("The homova command outputs a .homova file. \n"); - m->mothurOut("The homova command parameters are phylip, iters, groups, label, design, sets and processors. The design parameter is required.\n"); - m->mothurOut("The design parameter allows you to assign your groups to sets when you are running homova. It is required. \n"); - m->mothurOut("The design file looks like the group file. It is a 2 column tab delimited file, where the first column is the group name and the second column is the set the group belongs to.\n"); - m->mothurOut("The sets parameter allows you to specify which of the sets in your designfile you would like to analyze. The set names are separated by dashes. THe default is all sets in the designfile. To run the pairwise comparisons of all sets use sets=all.\n"); + m->mothurOut("The homova command parameters are phylip, iters, and alpha. The phylip and design parameters are required.\n"); + m->mothurOut("The design parameter allows you to assign your samples to groups when you are running homova. It is required. \n"); + m->mothurOut("The design file looks like the group file. It is a 2 column tab delimited file, where the first column is the sample name and the second column is the group the sample belongs to.\n"); m->mothurOut("The iters parameter allows you to set number of randomization for the P value. The default is 1000. \n"); - m->mothurOut("The groups parameter allows you to specify which of the groups you would like included. The group names are separated by dashes. groups=all will find all pairwise comparisons. \n"); - m->mothurOut("The label parameter allows you to select what distance levels you would like, and are also separated by dashes.\n"); - m->mothurOut("The processors parameter allows you to specify how many processors you would like to use. The default is 1. \n"); - m->mothurOut("The homova command should be in the following format: homova(design=yourDesignFile).\n"); - m->mothurOut("Example amova(design=temp.design, groups=A-B-C).\n"); - m->mothurOut("Note: No spaces between parameter labels (i.e. groups), '=' and parameters (i.e.yourGroups).\n\n"); - + m->mothurOut("The homova command should be in the following format: homova(phylip=file.dist, design=file.design).\n"); + m->mothurOut("Note: No spaces between parameter labels (i.e. iters), '=' and parameters (i.e. 1000).\n\n"); } catch(exception& e) { m->errorOut(e, "HomovaCommand", "help"); @@ -352,220 +188,71 @@ int HomovaCommand::execute(){ if (abort == true) { if (calledHelp) { return 0; } return 2; } //read design file - designMap = new GroupMap(designfile); + designMap = new GroupMap(designFileName); designMap->readDesignMap(); - //make sure sets are all in designMap - SharedUtil* util = new SharedUtil(); - util->setGroups(Sets, designMap->namesOfGroups); - delete util; - - //if the user pickedGroups or set sets=all, then we want to figure out how to split the pairwise comparison between processors - int numGroups = Sets.size(); - if (pickedGroups) { - for (int a=0; a groups; groups.push_back(Sets[a]); groups.push_back(Sets[l]); - namesOfGroupCombos.push_back(groups); - } - } - }else { //one giant compare - vector groups; - for (int a=0; ahasPath(phylipFileName); } - //only 1 combo - if (numGroups == 2) { processors = 1; } - else if (numGroups < 2) { m->mothurOut("Not enough sets, I need at least 2 valid sets. Unable to complete command."); m->mothurOutEndLine(); m->control_pressed = true; } + //read in distance matrix and square it + ReadPhylipVector readMatrix(phylipFileName); + vector sampleNames = readMatrix.read(distanceMatrix); -#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) - if(processors != 1){ - int numPairs = namesOfGroupCombos.size(); - int numPairsPerProcessor = numPairs / processors; - - for (int i = 0; i < processors; i++) { - int startPos = i * numPairsPerProcessor; - if(i == processors - 1){ - numPairsPerProcessor = numPairs - i * numPairsPerProcessor; - } - lines.push_back(linePair(startPos, numPairsPerProcessor)); + for(int i=0;ihasPath(sharedfile); } - - //for each calculator - for(int i = 0 ; i < calculators.size(); i++) { - - //create a new filename - ofstream out; - string outputFile = outputDir + m->getRootName(m->getSimpleName(sharedfile)) + calculators[i]->getName() + ".homova"; - m->openOutputFile(outputFile, out); - outputNames.push_back(outputFile); outputTypes["homova"].push_back(outputFile); - - //print headers - out << "label\tgroupsCompared\tBValue\tpValue" << endl; - out.close(); - } - - InputData input(sharedfile, "sharedfile"); - vector lookup = input.getSharedRAbundVectors(); - string lastLabel = lookup[0]->getLabel(); + //link designMap to rows/columns in distance matrix + map > origGroupSampleMap; + for(int i=0;igetGroup(sampleNames[i])].push_back(i); + } + int numGroups = origGroupSampleMap.size(); + + //create a new filename + ofstream HOMOVAFile; + string HOMOVAFileName = outputDir + m->getRootName(m->getSimpleName(phylipFileName)) + "homova"; + m->openOutputFile(HOMOVAFileName, HOMOVAFile); + outputNames.push_back(HOMOVAFileName); outputTypes["homova"].push_back(HOMOVAFileName); + + HOMOVAFile << "HOMOVA\tBValue\tP-value\tSSwithin_values" << endl; + m->mothurOut("HOMOVA\tBValue\tP-value\tSSwithin_values\n"); + + double fullHOMOVAPValue = runHOMOVA(HOMOVAFile, origGroupSampleMap, experimentwiseAlpha); + + if(fullHOMOVAPValue <= experimentwiseAlpha && numGroups > 2){ - //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label. - set processedLabels; - set userLabels = labels; + int numCombos = numGroups * (numGroups-1) / 2; + double pairwiseAlpha = experimentwiseAlpha / (double) numCombos; - //sanity check between shared file groups and design file groups - for (int i = 0; i < lookup.size(); i++) { - string group = designMap->getGroup(lookup[i]->getGroup()); - - if (group == "not found") { m->control_pressed = true; m->mothurOut("[ERROR]: " + lookup[i]->getGroup() + " is not in your design file, please correct."); m->mothurOutEndLine(); } - } + map >::iterator itA; + map >::iterator itB; - //as long as you are not at the end of the file or done wih the lines you want - while((lookup[0] != NULL) && ((allLines == 1) || (userLabels.size() != 0))) { - - if (m->control_pressed) { for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } globaldata->Groups.clear(); delete designMap; for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; } - - if(allLines == 1 || labels.count(lookup[0]->getLabel()) == 1){ + for(itA=origGroupSampleMap.begin();itA!=origGroupSampleMap.end();itA++){ + itB = itA;itB++; + for(itB;itB!=origGroupSampleMap.end();itB++){ - m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine(); - process(lookup); + map > pairwiseGroupSampleMap; + pairwiseGroupSampleMap[itA->first] = itA->second; + pairwiseGroupSampleMap[itB->first] = itB->second; - processedLabels.insert(lookup[0]->getLabel()); - userLabels.erase(lookup[0]->getLabel()); - } - - if ((m->anyLabelsToProcess(lookup[0]->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) { - string saveLabel = lookup[0]->getLabel(); - - for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } - lookup = input.getSharedRAbundVectors(lastLabel); - m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine(); - - process(lookup); - - processedLabels.insert(lookup[0]->getLabel()); - userLabels.erase(lookup[0]->getLabel()); - - //restore real lastlabel to save below - lookup[0]->setLabel(saveLabel); - } - - lastLabel = lookup[0]->getLabel(); - //prevent memory leak - for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; lookup[i] = NULL; } - - if (m->control_pressed) { globaldata->Groups.clear(); delete designMap; for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; } - - //get next line to process - lookup = input.getSharedRAbundVectors(); - } - - if (m->control_pressed) { globaldata->Groups.clear(); delete designMap; for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; } - - //output error messages about any remaining user labels - set::iterator it; - bool needToRun = false; - for (it = userLabels.begin(); it != userLabels.end(); it++) { - m->mothurOut("Your file does not include the label " + *it); - if (processedLabels.count(lastLabel) != 1) { - m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine(); - needToRun = true; - }else { - m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine(); - } - } - - //run last label if you need to - if (needToRun == true) { - for (int i = 0; i < lookup.size(); i++) { if (lookup[i] != NULL) { delete lookup[i]; } } - lookup = input.getSharedRAbundVectors(lastLabel); - - m->mothurOut(lookup[0]->getLabel()); m->mothurOutEndLine(); - - process(lookup); - - for (int i = 0; i < lookup.size(); i++) { delete lookup[i]; } - } - - //reset groups parameter - globaldata->Groups.clear(); - - - }else { //user provided distance matrix - - if (outputDir == "") { outputDir = m->hasPath(phylipfile); } - - //create a new filename - ofstream out; - string outputFile = outputDir + m->getRootName(m->getSimpleName(phylipfile)) + "homova"; - m->openOutputFile(outputFile, out); - outputNames.push_back(outputFile); outputTypes["homova"].push_back(outputFile); - - //print headers - out << "groupsCompared\tBValue\tpValue" << endl; - out.close(); - - ReadPhylipVector readMatrix(phylipfile); - vector< vector > completeMatrix; - vector names = readMatrix.read(completeMatrix); - -#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) - if(processors == 1){ - driver(0, namesOfGroupCombos.size(), names, "", completeMatrix); - }else{ - int process = 1; - vector processIDS; - - //loop through and create all the processes you want - while (process != processors) { - int pid = fork(); - - if (pid > 0) { - processIDS.push_back(pid); - process++; - }else if (pid == 0){ - driver(lines[process].start, lines[process].num, names, toString(getpid()), completeMatrix); - exit(0); - }else { - m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine(); - for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); } - exit(0); - } - } - - //do my part - driver(lines[0].start, lines[0].num, names, "", completeMatrix); - - //force parent to wait until all the processes are done - for (int i=0;i<(processors-1);i++) { - int temp = processIDS[i]; - wait(&temp); - } - - //append files - string outputFile = outputDir + m->getRootName(m->getSimpleName(phylipfile)) + "homova"; - for (int j = 0; j < processIDS.size(); j++) { - m->appendFiles((outputFile + toString(processIDS[j])), outputFile); - remove((outputFile + toString(processIDS[j])).c_str()); - } - + runHOMOVA(HOMOVAFile, pairwiseGroupSampleMap, pairwiseAlpha); + } } -#else - driver(0, namesOfGroupCombos.size(), names, "", completeMatrix); -#endif + HOMOVAFile << endl; + m->mothurOutEndLine(); + m->mothurOut("Experiment-wise error rate: " + toString(experimentwiseAlpha) + '\n'); + m->mothurOut("Pair-wise error rate (Bonferroni): " + toString(pairwiseAlpha) + '\n'); + } + else{ + m->mothurOut("Experiment-wise error rate: " + toString(experimentwiseAlpha) + '\n'); } - delete designMap; + m->mothurOut("If you have borderline P-values, you should try increasing the number of iterations\n"); - if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; } + delete designMap; m->mothurOutEndLine(); m->mothurOut("Output File Names: "); m->mothurOutEndLine(); @@ -579,384 +266,167 @@ int HomovaCommand::execute(){ exit(1); } } -//********************************************************************************************************************** -int HomovaCommand::process(vector thisLookup) { - try{ - -#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) - if(processors == 1){ - driver(0, namesOfGroupCombos.size(), thisLookup, ""); - }else{ - int process = 1; - vector processIDS; - - //loop through and create all the processes you want - while (process != processors) { - int pid = fork(); - - if (pid > 0) { - processIDS.push_back(pid); - process++; - }else if (pid == 0){ - driver(lines[process].start, lines[process].num, thisLookup, toString(getpid())); - exit(0); - }else { - m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine(); - for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); } - exit(0); - } - } - - //do my part - driver(lines[0].start, lines[0].num, thisLookup, ""); - - //force parent to wait until all the processes are done - for (int i=0;i<(processors-1);i++) { - int temp = processIDS[i]; - wait(&temp); - } - - //append files - for(int i = 0 ; i < calculators.size(); i++) { - string outputFile = outputDir + m->getRootName(m->getSimpleName(sharedfile)) + calculators[i]->getName() + ".homova"; - - for (int j = 0; j < processIDS.size(); j++) { - m->appendFiles((outputFile + toString(processIDS[j])), outputFile); - remove((outputFile + toString(processIDS[j])).c_str()); - } - } - } -#else - driver(0, namesOfGroupCombos.size(), thisLookUp, ""); -#endif - - return 0; - } - catch(exception& e) { - m->errorOut(e, "HomovaCommand", "process"); - exit(1); - } -} //********************************************************************************************************************** -int HomovaCommand::driver(int start, int num, vector thisLookup, string pidValue) { +double HomovaCommand::runHOMOVA(ofstream& HOMOVAFile, map > groupSampleMap, double alpha){ try { - vector subset; - EstOutput data; + map >::iterator it; + int numGroups = groupSampleMap.size(); + + vector ssWithinOrigVector; + double bValueOrig = calcBValue(groupSampleMap, ssWithinOrigVector); + + double counter = 0; + for(int i=0;i ssWithinRandVector; + map > randomizedGroup = getRandomizedGroups(groupSampleMap); + double bValueRand = calcBValue(randomizedGroup, ssWithinRandVector); + if(bValueRand > bValueOrig){ counter++; } + } - //for each calculator - for(int i = 0 ; i < calculators.size(); i++) { - - //create a new filename - ofstream out; - string outputFile = outputDir + m->getRootName(m->getSimpleName(sharedfile)) + calculators[i]->getName() + ".homova" + pidValue; - m->openOutputFileAppend(outputFile, out); - out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint); - - //for each combo - for (int c = start; c < (start+num); c++) { - - if (m->control_pressed) { out.close(); return 0; } - - //get set names - vector setNames; - for (int j = 0; j < namesOfGroupCombos[c].size(); j++) { setNames.push_back(namesOfGroupCombos[c][j]); } - - vector thisCombosLookup; - vector thisCombosLookupSets; //what set each sharedRabund is from to be used when calculating SSWithin - for (int k = 0; k < thisLookup.size(); k++) { - string thisGroup = thisLookup[k]->getGroup(); - - //is this group for a set we want to compare?? - if (m->inUsersGroups(designMap->getGroup(thisGroup), setNames)) { - thisCombosLookup.push_back(thisLookup[k]); - thisCombosLookupSets.push_back(designMap->getGroup(thisGroup)); - } - - } - - int numGroups = thisCombosLookup.size(); - - //calc the distance matrix - matrix.clear(); - matrix.resize(numGroups); - for (int k = 0; k < matrix.size(); k++) { for (int j = 0; j < matrix.size(); j++) { matrix[k].push_back(1.0); } } - - if (thisCombosLookup.size() == 0) { - m->mothurOut("[ERROR]: Missing shared info for sets. Skipping comparison."); m->mothurOutEndLine(); - }else{ - - out << thisLookup[0]->getLabel() << '\t'; - if (setNames.size() == 2) { out << setNames[0] << "-" << setNames[1] << '\t'; } - else { out << "all" << '\t'; } - - for (int k = 0; k < thisCombosLookup.size(); k++) { - for (int l = k; l < thisCombosLookup.size(); l++) { - - if (m->control_pressed) { out.close(); return 0; } - - if (k != l) { //we dont need to similiarity of a groups to itself - //get estimated similarity between 2 groups - subset.clear(); //clear out old pair of sharedrabunds - //add new pair of sharedrabunds - subset.push_back(thisCombosLookup[k]); subset.push_back(thisCombosLookup[l]); - - //if this calc needs all groups to calculate the pair load all groups - if (calculators[i]->getNeedsAll()) { - //load subset with rest of lookup for those calcs that need everyone to calc for a pair - for (int w = 0; w < thisCombosLookup.size(); w++) { - if ((w != k) && (w != l)) { subset.push_back(thisCombosLookup[w]); } - } - } - - data = calculators[i]->getValues(subset); //saves the calculator outputs - - //save values in similarity matrix - matrix[k][l] = 1.0 - data[0]; - matrix[l][k] = 1.0 - data[0]; - } - } - } - - //calc homova - calcHomova(out, setNames.size(), thisCombosLookupSets); - } - } - - out.close(); - } + double pValue = (double) counter / (double) iters; + string pString = ""; + if(pValue < 1/(double)iters){ pString = '<' + toString(1/(double)iters); } + else { pString = toString(pValue); } - return 0; + //print homova table + it = groupSampleMap.begin(); + HOMOVAFile << it->first; + m->mothurOut(it->first); + it++; + for(it;it!=groupSampleMap.end();it++){ + HOMOVAFile << '-' << it->first; + m->mothurOut('-' + it->first); + } + + HOMOVAFile << '\t' << bValueOrig << '\t' << pString; + m->mothurOut('\t' + toString(bValueOrig) + '\t' + pString); + + if(pValue < alpha){ + HOMOVAFile << "*"; + m->mothurOut("*"); + } + + for(int i=0;imothurOut('\t' + toString(ssWithinOrigVector[i])); + } + HOMOVAFile << endl; + m->mothurOutEndLine(); + + return pValue; } catch(exception& e) { - m->errorOut(e, "HomovaCommand", "driver"); + m->errorOut(e, "HomovaCommand", "runHOMOVA"); exit(1); } } + //********************************************************************************************************************** -int HomovaCommand::driver(int start, int num, vector names, string pidValue, vector< vector >& completeMatrix) { +double HomovaCommand::calcSigleSSWithin(vector sampleIndices) { try { + double ssWithin = 0.0; + int numSamplesInGroup = sampleIndices.size(); - //create a new filename - ofstream out; - string outputFile = outputDir + m->getRootName(m->getSimpleName(phylipfile)) + "homova" + pidValue; - m->openOutputFileAppend(outputFile, out); - out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint); - - //for each combo - for (int c = start; c < (start+num); c++) { - - if (m->control_pressed) { out.close(); return 0; } - - //get set names - vector setNames; - for (int j = 0; j < namesOfGroupCombos[c].size(); j++) { setNames.push_back(namesOfGroupCombos[c][j]); } - - vector thisCombosSets; //what set each line in the distance matrix is from to be used when calculating SSWithin - set indexes; - for (int k = 0; k < names.size(); k++) { - //is this group for a set we want to compare?? - if (m->inUsersGroups(designMap->getGroup(names[k]), setNames)) { - thisCombosSets.push_back(designMap->getGroup(names[k])); - indexes.insert(k); //save indexes of valid rows in matrix for submatrix - } - } + for(int i=0;imothurOut("[ERROR]: Missing distance info for sets. Skipping comparison."); m->mothurOutEndLine(); - }else{ + for(int j=0;jerrorOut(e, "HomovaCommand", "driver"); + m->errorOut(e, "HomovaCommand", "calcSigleSSWithin"); exit(1); } } + //********************************************************************************************************************** -int HomovaCommand::calcHomova(ofstream& out, int numTreatments, vector thisCombosLookupSets) { + +double HomovaCommand::calcBValue(map > groupSampleMap, vector& ssWithinVector) { try { + + map >::iterator it; - double SSTotal, BValue, pValue; - map SSWithin; - map::iterator it; - - SSTotal = calcWithin(matrix, numTreatments, thisCombosLookupSets); - int numSamples = matrix.size(); - - //calc BValue - map counts; - SSWithin = calcWithinEach(matrix, numTreatments, thisCombosLookupSets, counts); - - double sum = 0.0; - double sumDenom = 0.0; - for (it = SSWithin.begin(); it != SSWithin.end(); it++) { - double temp2 = (it->second) / (double) (counts[it->first] - 1); - sum += ((counts[it->first] - 1) * log(temp2)); - sumDenom += ((1 / (double) (counts[it->first] - 1)) - ( 1 / (double) (numSamples - numTreatments) )); - } - - double temp = SSTotal / (double) (numSamples - numTreatments); - double numeratorTerm1 = (numSamples - numTreatments) * log(temp); - double numerator = numeratorTerm1 - sum; - double denom = 1 + ((1 / (double) (3 * (numTreatments - 1))) * sumDenom); - - BValue = numerator / denom; + double numGroups = (double)groupSampleMap.size(); + ssWithinVector.resize(numGroups, 0); - //calc Pvalue - int count = 0; - for (int i = 0; i < iters; i++) { - if (m->control_pressed) { break; } - - //randomly shuffle names to randomize the matrix - vector copyNames = thisCombosLookupSets; - random_shuffle(copyNames.begin(), copyNames.end()); - - SSTotal = calcWithin(matrix, numTreatments, copyNames); - - counts.clear(); - map randomSSWithin = calcWithinEach(matrix, numTreatments, copyNames, counts); - - sum = 0.0; - sumDenom = 0.0; - for (it = randomSSWithin.begin(); it != randomSSWithin.end(); it++) { - double temp2 = (it->second) / (double) (counts[it->first] - numTreatments); - sum += ((counts[it->first] - 1) * log(temp2)); - sumDenom += ((1 / (double) (counts[it->first] - 1)) - ( 1 / (double) (numSamples - numTreatments) )); - } - - temp = SSTotal / (double) (numSamples - numTreatments); - numeratorTerm1 = (numSamples - numTreatments) * log(temp); - numerator = numeratorTerm1 - sum; - denom = 1 + ((1 / (double) (3 * (numTreatments - 1))) * sumDenom); - - double randomBValue = numerator / denom; - - if (randomBValue <= BValue) { count++; } - } - - pValue = count / (float) iters; - - out << BValue << '\t' << pValue << endl; - - return 0; - } - catch(exception& e) { - m->errorOut(e, "HomovaCommand", "calcHomova"); - exit(1); - } -} -//********************************************************************************************************************** -map HomovaCommand::calcWithinEach(vector< vector >& thisMatrix, int numTreatments, vector thisCombosLookupSets, map& count) { - try { - map within; //maps treatment to within value - map::iterator it; - map::iterator itCount; + double totalNumSamples = 0; + double ssWithinFull; + double secondTermSum = 0; + double inverseOneMinusSum = 0; + int index = 0; - for (int i = 0; i < thisCombosLookupSets.size(); i++) { - itCount = count.find(thisCombosLookupSets[i]); + ssWithinVector.resize(numGroups, 0); + for(it = groupSampleMap.begin();it!=groupSampleMap.end();it++){ + int numSamplesInGroup = it->second.size(); + totalNumSamples += numSamplesInGroup; - if (itCount == count.end()) { //first time we have seen this treatment - count[thisCombosLookupSets[i]] = 1; - }else { - count[thisCombosLookupSets[i]]++; - } + ssWithinVector[index] = calcSigleSSWithin(it->second); + ssWithinFull += ssWithinVector[index]; - //initialize within - within[thisCombosLookupSets[i]] = 0.0; + secondTermSum += (numSamplesInGroup - 1) * log(ssWithinVector[index] / (double)(numSamplesInGroup - 1)); + inverseOneMinusSum += 1.0 / (double)(numSamplesInGroup - 1); + index++; } - - //traverse lower triangle - for (int k = 0; k < thisMatrix.size(); k++) { - for (int l = k; l < thisMatrix[k].size(); l++) { - - //if you are from the same treatment then eij is 1 so add, else eij = 0 - if (thisCombosLookupSets[k] == thisCombosLookupSets[l]) { - within[thisCombosLookupSets[k]] += (thisMatrix[k][l] * thisMatrix[k][l]); //dij^2 - } - } - } + double B = (totalNumSamples - numGroups) * log(ssWithinFull/(totalNumSamples-numGroups)) - secondTermSum; + double denomintor = 1 + 1.0/(3.0 * (numGroups - 1.0)) * (inverseOneMinusSum - 1.0 / (double) (totalNumSamples - numGroups)); + B /= denomintor; - //1 / (numSamples in this treatment) - for (it = within.begin(); it != within.end(); it++) { - (it->second) *= (1.0 / (float) count[it->first]); - } + return B; - return within; } catch(exception& e) { - m->errorOut(e, "HomovaCommand", "calcWithinEach"); + m->errorOut(e, "HomovaCommand", "calcBValue"); exit(1); } } + //********************************************************************************************************************** -double HomovaCommand::calcWithin(vector< vector >& thisMatrix, int numTreatments, vector thisCombosLookupSets) { - try { - double total = 0.0; + +map > HomovaCommand::getRandomizedGroups(map > origMapping){ + try{ + vector sampleIndices; + vector samplesPerGroup; + + map >::iterator it; + for(it=origMapping.begin();it!=origMapping.end();it++){ + vector indices = it->second; + samplesPerGroup.push_back(indices.size()); + sampleIndices.insert(sampleIndices.end(), indices.begin(), indices.end()); + } - //traverse lower triangle - for (int k = 0; k < thisMatrix.size(); k++) { - for (int l = k; l < thisMatrix[k].size(); l++) { - - //if you are from the same treatment then eij is 1 so add, else eij = 0 - if (thisCombosLookupSets[k] == thisCombosLookupSets[l]) { - total += (thisMatrix[k][l] * thisMatrix[k][l]); //dij^2 - } + random_shuffle(sampleIndices.begin(), sampleIndices.end()); + + int index = 0; + map > randomizedGroups = origMapping; + for(it=randomizedGroups.begin();it!=randomizedGroups.end();it++){ + for(int i=0;isecond.size();i++){ + it->second[i] = sampleIndices[index++]; } } - //1 / (numSamples / numTreatments) - total *= (1.0 / (float) (thisMatrix.size() / (float) numTreatments)); - - return total; + return randomizedGroups; } - catch(exception& e) { - m->errorOut(e, "HomovaCommand", "calcWithin"); + catch (exception& e) { + m->errorOut(e, "AmovaCommand", "randomizeGroups"); exit(1); } } + //********************************************************************************************************************** diff --git a/homovacommand.h b/homovacommand.h index 55e4eac..11f3c07 100644 --- a/homovacommand.h +++ b/homovacommand.h @@ -12,12 +12,9 @@ #include "command.hpp" -#include "inputdata.h" -#include "sharedrabundvector.h" -#include "validcalculator.h" -#include "readphylipvector.h" -class GlobalData; +//class GlobalData; +class GroupMap; class HomovaCommand : public Command { @@ -33,35 +30,20 @@ public: void help(); private: - struct linePair { - int start; - int num; - linePair(int i, int j) : start(i), num(j) {} - }; - vector lines; - - GlobalData* globaldata; - GroupMap* designMap; + double runHOMOVA(ofstream& , map >, double); + double calcSigleSSWithin(vector); + double calcBValue(map >, vector&); + map > getRandomizedGroups(map >); + + bool abort; map > outputTypes; - - vector< vector > matrix; - bool abort, allLines, pickedGroups; - set labels; //holds labels to be used - string format, groups, label, outputDir, inputDir, designfile, sets, phylipfile, calc, sharedfile; - vector Groups, outputNames, Sets; - vector< vector > namesOfGroupCombos; - int iters, processors; - vector calculators; - - int driver(int, int, vector, string); - int driver(int, int, vector, string, vector< vector >&); - int process(vector); - int calcHomova(ofstream&, int, vector); - map calcWithinEach(vector< vector >&, int, vector, map&); - double calcWithin(vector< vector >&, int, vector); - double calcTotal(int); + vector outputNames; + + string outputDir, inputDir, designFileName, phylipFileName; + GroupMap* designMap; + vector< vector > distanceMatrix; + int iters; + double experimentwiseAlpha; }; #endif - - diff --git a/readphylipvector.cpp b/readphylipvector.cpp index f8907cc..c3ffcce 100644 --- a/readphylipvector.cpp +++ b/readphylipvector.cpp @@ -94,6 +94,7 @@ vector ReadPhylipVector::read(vector< vector >& matrix) { } } } + f.close(); return names; } @@ -171,6 +172,7 @@ vector ReadPhylipVector::read(vector& matrix) { } } } + f.close(); return names; } diff --git a/venn.cpp b/venn.cpp index 378e004..cf77a92 100644 --- a/venn.cpp +++ b/venn.cpp @@ -33,7 +33,7 @@ vector Venn::getPic(SAbundVector* sabund, vector vCalcs) { try { vector outputNames; - + for(int i=0;igetRootName(m->getSimpleName(globaldata->inputFileName)) + "." + sabund->getLabel() + "." + vCalcs[i]->getName() + ".svg"; outputNames.push_back(filenamesvg); @@ -78,7 +78,7 @@ vector Venn::getPic(SAbundVector* sabund, vector vCalcs) { //********************************************************************************************************************** vector Venn::getPic(vector lookup, vector vCalcs) { try { - + vector subset; vector outputNames; @@ -168,14 +168,16 @@ vector Venn::getPic(vector lookup, vector data = nseqsCalc->getValues(lookup); - + cout << data[0] << '\t' << data[1] << endl; sharedVal = data[0] + data[1]; numSeqsA = sabundA->getNumSeqs(); numSeqsB = sabundB->getNumSeqs(); + uniqSeqsToA = numSeqsA-data[0]; + uniqSeqsToB = numSeqsB-data[1]; delete nseqsCalc; } @@ -205,21 +207,22 @@ vector Venn::getPic(vector lookup, vector\n"; outsvg << "The number of species in group " + lookup[1]->getGroup() + " is " + toString(numB[0]); if (numB.size() == 3) { outsvg << " the lci is " + toString(numB[1]) + " and the hci is " + toString(numB[2]); } - if (nseqs) { outsvg << ", and the number of squences is " + toString(numSeqsB); } + if (nseqs) { outsvg << ", and the number of squences is " + toString(numSeqsB) + "; " + toString(uniqSeqsToB) + " sequences are not shared"; } outsvg << "\n"; outsvg << "The number of sepecies shared between groups " + lookup[0]->getGroup() + " and " + lookup[1]->getGroup() + " is " + toString(shared[0]); - if (nseqs) { outsvg << ", and the number of squences is " + toString(sharedVal); } + if (nseqs) { outsvg << ", and the number of squences is " + toString(sharedVal) + "; " + toString((sharedVal / (float)(numSeqsA + numSeqsB))*100) + "% of these sequences are shared"; } outsvg << "\n"; outsvg << "Percentage of species that are shared in groups " + lookup[0]->getGroup() + " and " + lookup[1]->getGroup() + " is " + toString((shared[0] / (float)(numA[0] + numB[0] - shared[0]))*100) + "\n"; + outsvg << "The total richness for all groups is " + toString((float)(numA[0] + numB[0] - shared[0]))+ "\n";; @@ -254,7 +257,7 @@ vector Venn::getPic(vector lookup, vectorcontrol_pressed) { outsvg.close(); return outputNames; } - int sharedVal, sharedABVal, sharedACVal, sharedBCVal, numSeqsA, numSeqsB, numSeqsC; + int sharedVal, sharedABVal, sharedACVal, sharedBCVal, numSeqsA, numSeqsB, numSeqsC, uniqSeqsToA, uniqSeqsToB, uniqSeqsToC; if (nseqs) { NSeqs* nseqsCalc = new NSeqs(); @@ -276,6 +279,9 @@ vector Venn::getPic(vector lookup, vectorgetNumSeqs(); numSeqsB = sabundB->getNumSeqs(); numSeqsC = sabundC->getNumSeqs(); + uniqSeqsToA = numSeqsA-sharedData[0]; + uniqSeqsToB = numSeqsC-sharedData[1]; + uniqSeqsToC = numSeqsB-sharedData[1]; delete nseqsCalc; } @@ -409,21 +415,21 @@ vector Venn::getPic(vector lookup, vector\n"; outsvg << "The number of species in group " + lookup[1]->getGroup() + " is " + toString(numB[0]); if (numB.size() == 3) { outsvg << " the lci is " + toString(numB[1]) + " and the hci is " + toString(numB[2]); } - if (nseqs) { outsvg << ", and the number of squences is " + toString(numSeqsB); } + if (nseqs) { outsvg << ", and the number of squences is " + toString(numSeqsB) + "; " + toString(uniqSeqsToB) + " sequences are not shared"; } outsvg << "\n"; outsvg << "The number of species in group " + lookup[2]->getGroup() + " is " + toString(numC[0]); if (numC.size() == 3) { outsvg << " the lci is " + toString(numC[1]) + " and the hci is " + toString(numC[2]); } - if (nseqs) { outsvg << ", and the number of squences is " + toString(numSeqsC); } + if (nseqs) { outsvg << ", and the number of squences is " + toString(numSeqsC) + "; " + toString(uniqSeqsToC) + " sequences are not shared"; } outsvg << "\n"; outsvg << "The total richness of all the groups is " + toString(numA[0] + numB[0] + numC[0] - sharedAB[0] - sharedAC[0] - sharedBC[0] + sharedABC) + "\n"; -- 2.39.2