*/
#include "homovacommand.h"
+#include "groupmap.h"
+#include "readphylipvector.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"
//**********************************************************************************************************************
-vector<string> HomovaCommand::getValidParameters(){
+vector<string> HomovaCommand::setParameters(){
try {
- string Array[] = {"groups","label","outputdir","iters","phylip","design","sets","processors","inputdir"};
- vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
+ CommandParameter pdesign("design", "InputTypes", "", "", "none", "none", "none","homova",false,true,true); parameters.push_back(pdesign);
+ CommandParameter pphylip("phylip", "InputTypes", "", "", "none", "none", "none","homova",false,true,true); parameters.push_back(pphylip);
+ CommandParameter psets("sets", "String", "", "", "", "", "","",false,false); parameters.push_back(psets);
+ CommandParameter piters("iters", "Number", "", "1000", "", "", "","",false,false); parameters.push_back(piters);
+ CommandParameter palpha("alpha", "Number", "", "0.05", "", "", "","",false,false); parameters.push_back(palpha);
+ CommandParameter pinputdir("inputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(pinputdir);
+ CommandParameter poutputdir("outputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(poutputdir);
+
+ vector<string> myArray;
+ for (int i = 0; i < parameters.size(); i++) { myArray.push_back(parameters[i].name); }
return myArray;
}
catch(exception& e) {
- m->errorOut(e, "HomovaCommand", "getValidParameters");
+ m->errorOut(e, "HomovaCommand", "setParameters");
exit(1);
}
}
//**********************************************************************************************************************
-HomovaCommand::HomovaCommand(){
+string HomovaCommand::getHelpString(){
try {
- abort = true;
- //initialize outputTypes
- vector<string> tempOutNames;
- outputTypes["homova"] = tempOutNames;
+ string helpString = "";
+ helpString += "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";
+ helpString += "The homova command outputs a .homova file. \n";
+ helpString += "The homova command parameters are phylip, iters, sets and alpha. The phylip and design parameters are required, unless valid current files exist.\n";
+ helpString += "The design parameter allows you to assign your samples to groups when you are running homova. It is required. \n";
+ helpString += "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";
+ helpString += "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.\n";
+ helpString += "The iters parameter allows you to set number of randomization for the P value. The default is 1000. \n";
+ helpString += "The homova command should be in the following format: homova(phylip=file.dist, design=file.design).\n";
+ helpString += "Note: No spaces between parameter labels (i.e. iters), '=' and parameters (i.e. 1000).\n";
+ return helpString;
}
catch(exception& e) {
- m->errorOut(e, "HomovaCommand", "HomovaCommand");
+ m->errorOut(e, "HomovaCommand", "getHelpString");
exit(1);
}
}
//**********************************************************************************************************************
-vector<string> HomovaCommand::getRequiredParameters(){
- try {
- string Array[] = {"design"};
- vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
- return myArray;
- }
- catch(exception& e) {
- m->errorOut(e, "HomovaCommand", "getRequiredParameters");
- exit(1);
- }
+string HomovaCommand::getOutputPattern(string type) {
+ try {
+ string pattern = "";
+
+ if (type == "homova") { pattern = "[filename],homova"; }
+ else { m->mothurOut("[ERROR]: No definition for type " + type + " output pattern.\n"); m->control_pressed = true; }
+
+ return pattern;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "HomovaCommand", "getOutputPattern");
+ exit(1);
+ }
}
//**********************************************************************************************************************
-vector<string> HomovaCommand::getRequiredFiles(){
+HomovaCommand::HomovaCommand(){
try {
- string Array[] = {};
- vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
- return myArray;
+ abort = true; calledHelp = true;
+ setParameters();
+ vector<string> tempOutNames;
+ outputTypes["homova"] = tempOutNames;
}
catch(exception& e) {
- m->errorOut(e, "HomovaCommand", "getRequiredFiles");
+ m->errorOut(e, "HomovaCommand", "HomovaCommand");
exit(1);
}
}
HomovaCommand::HomovaCommand(string option) {
try {
- globaldata = GlobalData::getInstance();
- abort = false;
- allLines = 1;
- labels.clear();
+ abort = false; calledHelp = false;
//allow user to run help
- if(option == "help") { help(); abort = true; }
+ if(option == "help") { help(); abort = true; calledHelp = true; }
+ else if(option == "citation") { citation(); abort = true; calledHelp = true;}
else {
- //valid paramters for this command
- string AlignArray[] = {"groups","label","outputdir","iters","phylip","design","sets","processors","inputdir"};
- vector<string> myArray (AlignArray, AlignArray+(sizeof(AlignArray)/sizeof(string)));
+ vector<string> myArray = setParameters();
OptionParser parser(option);
map<string,string> parameters = parser.getParameters();
}
}
- 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") {
+ //if there is a current phylip file, use it
+ phylipFileName = m->getPhylipFile();
+ if (phylipFileName != "") { m->mothurOut("Using " + phylipFileName + " as input file for the phylip parameter."); m->mothurOutEndLine(); }
+ else { m->mothurOut("You have no current phylip file and the phylip parameter is required."); m->mothurOutEndLine(); abort = true; }
+
+ }else { m->setPhylipFile(phylipFileName); }
//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; }
+ designFileName = validParameter.validFile(parameters, "design", true);
+ if (designFileName == "not open") { abort = true; }
+ else if (designFileName == "not found") {
+ //if there is a current design file, use it
+ designFileName = m->getDesignFile();
+ if (designFileName != "") { m->mothurOut("Using " + designFileName + " as input file for the design parameter."); m->mothurOutEndLine(); }
+ else { m->mothurOut("You have no current design file and the design parameter is required."); m->mothurOutEndLine(); abort = true; }
+ }else { m->setDesignFile(designFileName); }
+
+ string temp = validParameter.validFile(parameters, "iters", false);
+ if (temp == "not found") { temp = "1000"; }
+ m->mothurConvert(temp, iters);
+
+ temp = validParameter.validFile(parameters, "alpha", false);
+ if (temp == "not found") { temp = "0.05"; }
+ m->mothurConvert(temp, experimentwiseAlpha);
+
+ string sets = validParameter.validFile(parameters, "sets", false);
+ if (sets == "not found") { sets = ""; }
else {
- pickedGroups = true;
m->splitAtDash(sets, Sets);
}
-
-
- 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<string> 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; i<Estimators.size(); i++) {
- if (validCalculator->isValidCalculator("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());
- }
- }
- }
- }
}
}
exit(1);
}
}
-
-//**********************************************************************************************************************
-
-void HomovaCommand::help(){
- try {
- 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 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 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 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");
-
- }
- catch(exception& e) {
- m->errorOut(e, "HomovaCommand", "help");
- exit(1);
- }
-}
-
-//**********************************************************************************************************************
-
-HomovaCommand::~HomovaCommand(){}
-
//**********************************************************************************************************************
int HomovaCommand::execute(){
try {
- if (abort == true) { return 0; }
+ 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<numGroups; a++) {
- for (int l = 0; l < a; l++) {
- vector<string> groups; groups.push_back(Sets[a]); groups.push_back(Sets[l]);
- namesOfGroupCombos.push_back(groups);
- }
+ if (outputDir == "") { outputDir = m->hasPath(phylipFileName); }
+
+ //read in distance matrix and square it
+ ReadPhylipVector readMatrix(phylipFileName);
+ vector<string> sampleNames = readMatrix.read(distanceMatrix);
+
+ if (Sets.size() != 0) { //user selected sets, so we want to remove the samples not in those sets
+ SharedUtil util;
+ vector<string> dGroups = designMap->getNamesOfGroups();
+ util.setGroups(Sets, dGroups);
+
+ for(int i=0;i<distanceMatrix.size();i++){
+
+ if (m->control_pressed) { delete designMap; return 0; }
+
+ string group = designMap->getGroup(sampleNames[i]);
+
+ if (group == "not found") {
+ m->mothurOut("[ERROR]: " + sampleNames[i] + " is not in your design file, please correct."); m->mothurOutEndLine(); m->control_pressed = true;
+ }else if (!m->inUsersGroups(group, Sets)){ //not in set we want remove it
+ //remove from all other rows
+ for(int j=0;j<distanceMatrix.size();j++){
+ distanceMatrix[j].erase(distanceMatrix[j].begin()+i);
+ }
+ distanceMatrix.erase(distanceMatrix.begin()+i);
+ sampleNames.erase(sampleNames.begin()+i);
+ i--;
+ }
+ }
+ }
+
+ for(int i=0;i<distanceMatrix.size();i++){
+ for(int j=0;j<i;j++){
+ distanceMatrix[i][j] *= distanceMatrix[i][j];
}
- }else { //one giant compare
- vector<string> groups;
- for (int a=0; a<numGroups; a++) { groups.push_back(Sets[a]); }
- namesOfGroupCombos.push_back(groups);
}
- //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; }
-
-#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
- if(processors != 1){
- int numPairs = namesOfGroupCombos.size();
- int numPairsPerProcessor = numPairs / processors;
+ //link designMap to rows/columns in distance matrix
+ map<string, vector<int> > origGroupSampleMap;
+ for(int i=0;i<sampleNames.size();i++){
+ string group = designMap->getGroup(sampleNames[i]);
- 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));
- }
+ if (group == "not found") {
+ m->mothurOut("[ERROR]: " + sampleNames[i] + " is not in your design file, please correct."); m->mothurOutEndLine(); m->control_pressed = true;
+ }else { origGroupSampleMap[group].push_back(i); }
}
-#endif
+ int numGroups = origGroupSampleMap.size();
- 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() + ".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<SharedRAbundVector*> 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<string> processedLabels;
- set<string> 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<string>::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);
+ if (m->control_pressed) { delete designMap; return 0; }
+
+ //create a new filename
+ ofstream HOMOVAFile;
+ map<string, string> variables;
+ variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(phylipFileName));
+ string HOMOVAFileName = getOutputFileName("homova", variables);
+ m->openOutputFile(HOMOVAFileName, HOMOVAFile);
+ outputNames.push_back(HOMOVAFileName); outputTypes["homova"].push_back(HOMOVAFileName);
+
+ HOMOVAFile << "HOMOVA\tBValue\tP-value\tSSwithin/(Ni-1)_values" << endl;
+ m->mothurOut("HOMOVA\tBValue\tP-value\tSSwithin/(Ni-1)_values\n");
+
+ double fullHOMOVAPValue = runHOMOVA(HOMOVAFile, origGroupSampleMap, experimentwiseAlpha);
+
+ if(fullHOMOVAPValue <= experimentwiseAlpha && numGroups > 2){
- //print headers
- out << "groupsCompared\tBValue\tpValue" << endl;
- out.close();
+ int numCombos = numGroups * (numGroups-1) / 2;
+ double pairwiseAlpha = experimentwiseAlpha / (double) numCombos;
- ReadPhylipVector readMatrix(phylipfile);
- vector< vector<double> > completeMatrix;
- vector<string> names = readMatrix.read(completeMatrix);
+ map<string, vector<int> >::iterator itA;
+ map<string, vector<int> >::iterator itB;
-#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
- if(processors == 1){
- driver(0, namesOfGroupCombos.size(), names, "", completeMatrix);
- }else{
- int process = 1;
- vector<int> processIDS;
-
- //loop through and create all the processes you want
- while (process != processors) {
- int pid = fork();
+ for(itA=origGroupSampleMap.begin();itA!=origGroupSampleMap.end();itA++){
+ itB = itA;itB++;
+ for(;itB!=origGroupSampleMap.end();itB++){
+ map<string, vector<int> > pairwiseGroupSampleMap;
+ pairwiseGroupSampleMap[itA->first] = itA->second;
+ pairwiseGroupSampleMap[itB->first] = itB->second;
- 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();
exit(1);
}
}
-//**********************************************************************************************************************
-int HomovaCommand::process(vector<SharedRAbundVector*> thisLookup) {
- try{
-
-#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
- if(processors == 1){
- driver(0, namesOfGroupCombos.size(), thisLookup, "");
- }else{
- int process = 1;
- vector<int> 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<SharedRAbundVector*> thisLookup, string pidValue) {
+double HomovaCommand::runHOMOVA(ofstream& HOMOVAFile, map<string, vector<int> > groupSampleMap, double alpha){
try {
- vector<SharedRAbundVector*> subset;
- EstOutput data;
+ map<string, vector<int> >::iterator it;
+ int numGroups = groupSampleMap.size();
+
+ vector<double> ssWithinOrigVector;
+ double bValueOrig = calcBValue(groupSampleMap, ssWithinOrigVector);
+
+ double counter = 0;
+ for(int i=0;i<iters;i++){
+ vector<double> ssWithinRandVector;
+ map<string, vector<int> > 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<string> setNames;
- for (int j = 0; j < namesOfGroupCombos[c].size(); j++) { setNames.push_back(namesOfGroupCombos[c][j]); }
-
- vector<SharedRAbundVector*> thisCombosLookup;
- vector<string> 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!=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;i<numGroups;i++){
+ HOMOVAFile << '\t' << ssWithinOrigVector[i];
+ m->mothurOut('\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<string> names, string pidValue, vector< vector<double> >& completeMatrix) {
+double HomovaCommand::calcSigleSSWithin(vector<int> 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<string> setNames;
- for (int j = 0; j < namesOfGroupCombos[c].size(); j++) { setNames.push_back(namesOfGroupCombos[c][j]); }
+ for(int i=0;i<numSamplesInGroup;i++){
+ int row = sampleIndices[i];
- vector<string> thisCombosSets; //what set each line in the distance matrix is from to be used when calculating SSWithin
- set<int> 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'; }
- else { out << "all" << '\t'; }
+ for(int j=0;j<numSamplesInGroup;j++){
+ int col = sampleIndices[j];
- //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++;
- }
+ if(col < row){
+ ssWithin += distanceMatrix[row][col];
}
- //calc homova
- calcHomova(out, setNames.size(), thisCombosSets);
}
}
- out.close();
-
-
- return 0;
-
+ ssWithin /= numSamplesInGroup;
+ return ssWithin;
}
catch(exception& e) {
- m->errorOut(e, "HomovaCommand", "driver");
+ m->errorOut(e, "HomovaCommand", "calcSigleSSWithin");
exit(1);
}
}
+
//**********************************************************************************************************************
-int HomovaCommand::calcHomova(ofstream& out, int numTreatments, vector<string> thisCombosLookupSets) {
+
+double HomovaCommand::calcBValue(map<string, vector<int> > groupSampleMap, vector<double>& ssWithinVector) {
try {
+
+ map<string, vector<int> >::iterator it;
- double SSTotal, BValue, pValue;
- map<string, double> SSWithin;
- map<string, double>::iterator it;
-
- //SSTotal = 0.0;
- SSTotal = calcTotal(numTreatments);
- cout << "sstotal = " << SSTotal << " (n-P) = " << (matrix.size() - numTreatments) << endl;
- int numSamples = matrix.size();
-
- //calc BValue
- map<string, int> counts;
- SSWithin = calcWithin(matrix, numTreatments, thisCombosLookupSets, counts);
-
- double sum = 0.0;
- double sumDenom = 0.0;
- for (it = SSWithin.begin(); it != SSWithin.end(); it++) {
- cout << it->first << '\t' << it->second << '\t' << (counts[it->first] - numTreatments) << endl;
- double temp2 = (it->second) / (double) (counts[it->first] - 1);
- cout << "sumTerm = " << temp2 << '\t' << "ln sumTerm = " << log(temp2) << endl;
- 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);
+ double numGroups = (double)groupSampleMap.size();
+ ssWithinVector.resize(numGroups, 0);
- BValue = numerator / denom;
- cout << "numeratorTerm1 = " << numeratorTerm1 << " sum = " << sum << " sumDenom = " << sumDenom << " numerator = " << numerator << " denom = " << denom << " B = " << BValue << endl;
+ double totalNumSamples = 0;
+ double ssWithinFull;
+ double secondTermSum = 0;
+ double inverseOneMinusSum = 0;
+ int index = 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<string> copyNames = thisCombosLookupSets;
- random_shuffle(copyNames.begin(), copyNames.end());
+ ssWithinVector.resize(numGroups, 0);
+ for(it = groupSampleMap.begin();it!=groupSampleMap.end();it++){
+ int numSamplesInGroup = it->second.size();
+ totalNumSamples += numSamplesInGroup;
- counts.clear();
- map<string, double> randomSSWithin = calcWithin(matrix, numTreatments, copyNames, counts);
+ ssWithinVector[index] = calcSigleSSWithin(it->second);
+ ssWithinFull += ssWithinVector[index];
- 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) ));
- }
-
- numerator = numeratorTerm1 - sum;
- denom = 1 + ((1 / (double) (3 * (numTreatments - 1))) * sumDenom);
+ secondTermSum += (numSamplesInGroup - 1) * log(ssWithinVector[index] / (double)(numSamplesInGroup - 1));
+ inverseOneMinusSum += 1.0 / (double)(numSamplesInGroup - 1);
- double randomBValue = numerator / denom;
-
- if (randomBValue >= BValue) { count++; }
+ ssWithinVector[index] /= (double)(numSamplesInGroup - 1); //this line is only for output purposes to scale SSw by the number of samples in the group
+ index++;
}
-
- pValue = count / (float) iters;
- out << BValue << '\t' << pValue << endl;
-
- return 0;
- }
- catch(exception& e) {
- m->errorOut(e, "HomovaCommand", "calcHomova");
- exit(1);
- }
-}
-//**********************************************************************************************************************
-map<string, double> HomovaCommand::calcWithin(vector< vector<double> >& thisMatrix, int numTreatments, vector<string> thisCombosLookupSets, map<string, int>& count) {
- try {
- map<string, double> within; //maps treatment to within value
- map<string, double>::iterator it;
- map<string, int>::iterator itCount;
+ 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;
- for (int i = 0; i < thisCombosLookupSets.size(); i++) {
- itCount = count.find(thisCombosLookupSets[i]);
-
- if (itCount == count.end()) { //first time we have seen this treatment
- count[thisCombosLookupSets[i]] = 1;
- }else {
- count[thisCombosLookupSets[i]]++;
- }
-
- //initialize within
- within[thisCombosLookupSets[i]] = 0.0;
- }
-
+ return B;
- //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
- }
- }
- }
-
- //1 / (numSamples in this treatment)
- for (it = within.begin(); it != within.end(); it++) {
- (it->second) *= (1.0 / (float) count[it->first]);
- }
-
- return within;
}
catch(exception& e) {
- m->errorOut(e, "HomovaCommand", "calcWithin");
+ m->errorOut(e, "HomovaCommand", "calcBValue");
exit(1);
}
}
+
//**********************************************************************************************************************
-double HomovaCommand::calcTotal(int numTreatments) {
- try {
- double total = 0.0;
+
+map<string, vector<int> > HomovaCommand::getRandomizedGroups(map<string, vector<int> > origMapping){
+ try{
+ vector<int> sampleIndices;
+ vector<int> samplesPerGroup;
+
+ map<string, vector<int> >::iterator it;
+ for(it=origMapping.begin();it!=origMapping.end();it++){
+ vector<int> indices = it->second;
+ samplesPerGroup.push_back(indices.size());
+ sampleIndices.insert(sampleIndices.end(), indices.begin(), indices.end());
+ }
+
+ random_shuffle(sampleIndices.begin(), sampleIndices.end());
- //traverse lower triangle
- for (int k = 0; k < matrix.size(); k++) {
- for (int l = k; l < matrix[k].size(); l++) {
- total += (matrix[k][l] * matrix[k][l]); //dij^2
+ int index = 0;
+ map<string, vector<int> > randomizedGroups = origMapping;
+ for(it=randomizedGroups.begin();it!=randomizedGroups.end();it++){
+ for(int i=0;i<it->second.size();i++){
+ it->second[i] = sampleIndices[index++];
}
}
- //1 / numSamples
- total *= (1.0 / (float) matrix.size());
-
- return total;
+ return randomizedGroups;
}
- catch(exception& e) {
- m->errorOut(e, "HomovaCommand", "calcTotal");
+ catch (exception& e) {
+ m->errorOut(e, "AmovaCommand", "randomizeGroups");
exit(1);
}
}
+
//**********************************************************************************************************************
projects / mothur.git / blobdiff