From ea34218a94af7cb96d4a1b44e444a87c7ab7456c Mon Sep 17 00:00:00 2001 From: westcott Date: Tue, 8 Feb 2011 12:03:48 +0000 Subject: [PATCH] added amova command --- Mothur.xcodeproj/project.pbxproj | 8 +- amovacommand.cpp | 912 +++++++++++++++++++++++++++++++ amovacommand.h | 64 +++ chimeraslayercommand.cpp | 2 +- commandfactory.cpp | 5 + makefile | 4 +- summarysharedcommand.cpp | 2 +- 7 files changed, 992 insertions(+), 5 deletions(-) create mode 100644 amovacommand.cpp create mode 100644 amovacommand.h diff --git a/Mothur.xcodeproj/project.pbxproj b/Mothur.xcodeproj/project.pbxproj index 849d202..c3a3d00 100644 --- a/Mothur.xcodeproj/project.pbxproj +++ b/Mothur.xcodeproj/project.pbxproj @@ -14,6 +14,7 @@ A713EBED12DC7C5E000092AC /* nmdscommand.cpp in Sources */ = {isa = PBXBuildFile; fileRef = A713EBEC12DC7C5E000092AC /* nmdscommand.cpp */; }; A71FE12C12EDF72400963CA7 /* mergegroupscommand.cpp in Sources */ = {isa = PBXBuildFile; fileRef = A71FE12B12EDF72400963CA7 /* mergegroupscommand.cpp */; }; A727864412E9E28C00F86ABA /* removerarecommand.cpp in Sources */ = {isa = PBXBuildFile; fileRef = A727864312E9E28C00F86ABA /* removerarecommand.cpp */; }; + A7A61F2D130062E000E05B6B /* amovacommand.cpp in Sources */ = {isa = PBXBuildFile; fileRef = A7A61F2C130062E000E05B6B /* amovacommand.cpp */; }; A7E9B88112D37EC400DA6239 /* ace.cpp in Sources */ = {isa = PBXBuildFile; fileRef = A7E9B64F12D37EC300DA6239 /* ace.cpp */; }; A7E9B88212D37EC400DA6239 /* aligncommand.cpp in Sources */ = {isa = PBXBuildFile; fileRef = A7E9B65112D37EC300DA6239 /* aligncommand.cpp */; }; A7E9B88312D37EC400DA6239 /* alignment.cpp in Sources */ = {isa = PBXBuildFile; fileRef = A7E9B65312D37EC300DA6239 /* alignment.cpp */; }; @@ -316,6 +317,8 @@ A727864212E9E28C00F86ABA /* removerarecommand.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = removerarecommand.h; sourceTree = ""; }; A727864312E9E28C00F86ABA /* removerarecommand.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = removerarecommand.cpp; sourceTree = ""; }; A7A61F1A130035C800E05B6B /* LICENSE */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = text; path = LICENSE; sourceTree = ""; }; + A7A61F2B130062E000E05B6B /* amovacommand.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = amovacommand.h; sourceTree = ""; }; + A7A61F2C130062E000E05B6B /* amovacommand.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = amovacommand.cpp; sourceTree = ""; }; A7E9B64F12D37EC300DA6239 /* ace.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = ace.cpp; sourceTree = ""; }; A7E9B65012D37EC300DA6239 /* ace.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = ace.h; sourceTree = ""; }; A7E9B65112D37EC300DA6239 /* aligncommand.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = aligncommand.cpp; sourceTree = ""; }; @@ -1021,6 +1024,8 @@ A7E9B6AE12D37EC400DA6239 /* command.hpp */, A7E9B65112D37EC300DA6239 /* aligncommand.cpp */, A7E9B65212D37EC300DA6239 /* aligncommand.h */, + A7A61F2B130062E000E05B6B /* amovacommand.h */, + A7A61F2C130062E000E05B6B /* amovacommand.cpp */, A7E9B66012D37EC300DA6239 /* binsequencecommand.cpp */, A7E9B66112D37EC300DA6239 /* binsequencecommand.h */, A7E9B66A12D37EC400DA6239 /* bootstrapsharedcommand.cpp */, @@ -1903,6 +1908,7 @@ A727864412E9E28C00F86ABA /* removerarecommand.cpp in Sources */, A71FE12C12EDF72400963CA7 /* mergegroupscommand.cpp in Sources */, 7E6BE10A12F710D8007ADDBE /* refchimeratest.cpp in Sources */, + A7A61F2D130062E000E05B6B /* amovacommand.cpp in Sources */, ); runOnlyForDeploymentPostprocessing = 0; }; @@ -1946,7 +1952,7 @@ GCC_PREPROCESSOR_DEFINITIONS = ( "MOTHUR_FILES=\"\\\"../release\\\"\"", "VERSION=\"\\\"1.16.0\\\"\"", - "RELEASE_DATE=\"\\\"1/26/2010\\\"\"", + "RELEASE_DATE=\"\\\"2/7/2011\\\"\"", ); GCC_WARN_ABOUT_MISSING_NEWLINE = YES; GCC_WARN_ABOUT_RETURN_TYPE = YES; diff --git a/amovacommand.cpp b/amovacommand.cpp new file mode 100644 index 0000000..05c462d --- /dev/null +++ b/amovacommand.cpp @@ -0,0 +1,912 @@ +/* + * amovacommand.cpp + * mothur + * + * Created by westcott on 2/7/11. + * Copyright 2011 Schloss Lab. All rights reserved. + * + */ + +#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" + +//********************************************************************************************************************** +vector AmovaCommand::getValidParameters(){ + try { + string Array[] = {"groups","label","outputdir","iters","phylip","design","sets","processors","inputdir"}; + vector myArray (Array, Array+(sizeof(Array)/sizeof(string))); + return myArray; + } + catch(exception& e) { + m->errorOut(e, "AmovaCommand", "getValidParameters"); + exit(1); + } +} +//********************************************************************************************************************** +AmovaCommand::AmovaCommand(){ + try { + abort = true; + //initialize outputTypes + vector tempOutNames; + outputTypes["amova"] = tempOutNames; + } + catch(exception& e) { + m->errorOut(e, "AmovaCommand", "AmovaCommand"); + exit(1); + } +} +//********************************************************************************************************************** +vector AmovaCommand::getRequiredParameters(){ + try { + string Array[] = {"design"}; + vector myArray (Array, Array+(sizeof(Array)/sizeof(string))); + return myArray; + } + catch(exception& e) { + m->errorOut(e, "AmovaCommand", "getRequiredParameters"); + exit(1); + } +} +//********************************************************************************************************************** +vector AmovaCommand::getRequiredFiles(){ + try { + string Array[] = {}; + vector myArray (Array, Array+(sizeof(Array)/sizeof(string))); + return myArray; + } + catch(exception& e) { + m->errorOut(e, "AmovaCommand", "getRequiredFiles"); + exit(1); + } +} +//********************************************************************************************************************** + +AmovaCommand::AmovaCommand(string option) { + try { + globaldata = GlobalData::getInstance(); + abort = false; + allLines = 1; + labels.clear(); + + //allow user to run help + if(option == "help") { help(); abort = true; } + + else { + //valid paramters for this command + string AlignArray[] = {"groups","label","outputdir","iters","phylip","design","sets","processors","inputdir"}; + vector myArray (AlignArray, AlignArray+(sizeof(AlignArray)/sizeof(string))); + + OptionParser parser(option); + map parameters = parser.getParameters(); + + ValidParameters validParameter; + + //check to make sure all parameters are valid for command + map::iterator it; + for (it = parameters.begin(); it != parameters.end(); it++) { + if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; } + } + + //initialize outputTypes + vector tempOutNames; + outputTypes["amova"] = tempOutNames; + + //if the user changes the output directory command factory will send this info to us in the output parameter + outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){ outputDir = ""; } + + //if the user changes the input directory command factory will send this info to us in the output parameter + string inputDir = validParameter.validFile(parameters, "inputdir", false); + if (inputDir == "not found"){ inputDir = ""; } + else { + string path; + it = parameters.find("design"); + //user has given a template file + if(it != parameters.end()){ + path = m->hasPath(it->second); + //if the user has not given a path then, add inputdir. else leave path alone. + if (path == "") { parameters["design"] = inputDir + it->second; } + } + + it = parameters.find("phylip"); + //user has given a template file + if(it != parameters.end()){ + path = m->hasPath(it->second); + //if the user has not given a path then, add inputdir. else leave path alone. + if (path == "") { parameters["phylip"] = inputDir + it->second; } + } + } + + 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); } + + //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 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); + } + + + 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(); + + if (sharedfile != "") { + + 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"); + exit(1); + } +} + +//********************************************************************************************************************** + +void AmovaCommand::help(){ + try { + 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, 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.\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(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, "AmovaCommand", "help"); + exit(1); + } +} + +//********************************************************************************************************************** + +AmovaCommand::~AmovaCommand(){} + +//********************************************************************************************************************** + +int AmovaCommand::execute(){ + try { + + if (abort == true) { return 0; } + + //read design file + designMap = new GroupMap(designfile); + 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; 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(); + } + + 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 + + //for each calculator + for(int i = 0 ; i < calculators.size(); i++) { + //create a new filename + ofstream out; + string outputFile = outputDir + m->getRootName(m->getSimpleName(phylipfile)) + calculators[i]->getName() + ".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(); + } + + 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 + for(int i = 0 ; i < calculators.size(); i++) { + string outputFile = outputDir + m->getRootName(m->getSimpleName(phylipfile)) + 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(), names, "", completeMatrix); + #endif + + } + + 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(); } + m->mothurOutEndLine(); + + return 0; + } + catch(exception& e) { + m->errorOut(e, "AmovaCommand", "execute"); + 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{ + + 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 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 { + + //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 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'; } + else { out << "all" << '\t'; } + + //fill submatrix + int rowCount = 0; + int columnCount = 0; + for (int j = 0; j < completeMatrix.size(); j++) { + + if (indexes.count(j) != 0) { //we want at least part of this row + 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); + } +} +//********************************************************************************************************************** +int AmovaCommand::calcAmova(ofstream& out, int numTreatments, vector thisCombosLookupSets) { + try { + + double SSWithin, SSTotal, SSAmoung, MSWithin, MSAmoung, FStatistic, pValue; + + SSWithin = calcWithin(matrix, numTreatments, thisCombosLookupSets); + SSTotal = calcTotal(numTreatments); + + 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()); + + double randomSSWithin = calcWithin(matrix, numTreatments, copyNames); + + if (randomSSWithin <= SSWithin) { count++; } + } + + SSAmoung = SSTotal - SSWithin; + MSWithin = SSWithin / (double) (matrix.size() - numTreatments); + MSAmoung = SSAmoung / (double) (numTreatments - 1); + FStatistic = MSAmoung / MSWithin; + + pValue = count / (float) iters; + + out << MSWithin << '\t' << MSAmoung << '\t' << FStatistic << '\t' << pValue << endl; + + return 0; + } + catch(exception& e) { + m->errorOut(e, "AmovaCommand", "calcAmova"); + exit(1); + } +} +//********************************************************************************************************************** +double AmovaCommand::calcWithin(vector< vector >& thisMatrix, int numTreatments, vector thisCombosLookupSets) { + try { + double within = 0.0; + + //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 += (thisMatrix[k][l] * thisMatrix[k][l]); //dij^2 + } + } + } + + //1 / (numSamples / numTreatments) + within *= (1.0 / (float) (thisMatrix.size() / (float) numTreatments)); + + return within; + } + catch(exception& e) { + m->errorOut(e, "AmovaCommand", "calcWithin"); + exit(1); + } +} +//********************************************************************************************************************** +double AmovaCommand::calcTotal(int numTreatments) { + try { + double total = 0.0; + + //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 + } + } + + //1 / numSamples + total *= (1.0 / (float) matrix.size()); + + return total; + } + catch(exception& e) { + m->errorOut(e, "AmovaCommand", "calcTotal"); + exit(1); + } +} +//********************************************************************************************************************** + diff --git a/amovacommand.h b/amovacommand.h new file mode 100644 index 0000000..3b90b96 --- /dev/null +++ b/amovacommand.h @@ -0,0 +1,64 @@ +#ifndef AMOVACOMMAND_H +#define AMOVACOMMAND_H + +/* + * amovacommand.h + * mothur + * + * Created by westcott on 2/7/11. + * Copyright 2011 Schloss Lab. All rights reserved. + * + */ + +#include "command.hpp" +#include "inputdata.h" +#include "sharedrabundvector.h" +#include "validcalculator.h" +#include "readphylipvector.h" + +class GlobalData; + +class AmovaCommand : public Command { + +public: + AmovaCommand(string); + AmovaCommand(); + ~AmovaCommand(); + vector getRequiredParameters(); + vector getValidParameters(); + vector getRequiredFiles(); + map > getOutputFiles() { return outputTypes; } + int execute(); + void help(); + +private: + struct linePair { + int start; + int num; + linePair(int i, int j) : start(i), num(j) {} + }; + vector lines; + + GlobalData* globaldata; + GroupMap* designMap; + 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 calcAmova(ofstream&, int, vector); + double calcWithin(vector< vector >&, int, vector); + double calcTotal(int); +}; + +#endif + diff --git a/chimeraslayercommand.cpp b/chimeraslayercommand.cpp index bd9ff4d..a7ed94f 100644 --- a/chimeraslayercommand.cpp +++ b/chimeraslayercommand.cpp @@ -307,7 +307,7 @@ void ChimeraSlayerCommand::help(){ #ifdef USE_MPI m->mothurOut("When using MPI, the processors parameter is set to the number of MPI processes running. \n"); #endif - m->mothurOut("The trim parameter allows you to output a new fasta file containing your sequences with the chimeric ones trimmed to include only their longest peice, default=F. \n"); + m->mothurOut("The trim parameter allows you to output a new fasta file containing your sequences with the chimeric ones trimmed to include only their longest piece, default=F. \n"); m->mothurOut("The window parameter allows you to specify the window size for searching for chimeras, default=50. \n"); m->mothurOut("The increment parameter allows you to specify how far you move each window while finding chimeric sequences, default=5.\n"); m->mothurOut("The numwanted parameter allows you to specify how many sequences you would each query sequence compared with, default=15.\n"); diff --git a/commandfactory.cpp b/commandfactory.cpp index 250e4f6..42567e4 100644 --- a/commandfactory.cpp +++ b/commandfactory.cpp @@ -109,6 +109,7 @@ #include "nmdscommand.h" #include "removerarecommand.h" #include "mergegroupscommand.h" +#include "amovacommand.h" /*******************************************************/ @@ -220,6 +221,7 @@ CommandFactory::CommandFactory(){ commands["pca"] = "pca"; commands["nmds"] = "nmds"; commands["remove.rare"] = "remove.rare"; + commands["amova"] = "amova"; commands["merge.groups"] = "merge.groups"; commands["pairwise.seqs"] = "MPIEnabled"; commands["pipeline.pds"] = "MPIEnabled"; @@ -382,6 +384,7 @@ Command* CommandFactory::getCommand(string commandName, string optionString){ else if(commandName == "corr.axes") { command = new CorrAxesCommand(optionString); } else if(commandName == "remove.rare") { command = new RemoveRareCommand(optionString); } else if(commandName == "merge.groups") { command = new MergeGroupsCommand(optionString); } + else if(commandName == "amova") { command = new AmovaCommand(optionString); } else { command = new NoCommand(optionString); } return command; @@ -509,6 +512,7 @@ Command* CommandFactory::getCommand(string commandName, string optionString, str else if(commandName == "corr.axes") { pipecommand = new CorrAxesCommand(optionString); } else if(commandName == "remove.rare") { pipecommand = new RemoveRareCommand(optionString); } else if(commandName == "merge.groups") { pipecommand = new MergeGroupsCommand(optionString); } + else if(commandName == "amova") { pipecommand = new AmovaCommand(optionString); } else { pipecommand = new NoCommand(optionString); } return pipecommand; @@ -624,6 +628,7 @@ Command* CommandFactory::getCommand(string commandName){ else if(commandName == "corr.axes") { shellcommand = new CorrAxesCommand(); } else if(commandName == "remove.rare") { shellcommand = new RemoveRareCommand(); } else if(commandName == "merge.groups") { shellcommand = new MergeGroupsCommand(); } + else if(commandName == "amova") { shellcommand = new AmovaCommand(); } else { shellcommand = new NoCommand(); } return shellcommand; diff --git a/makefile b/makefile index 266c0c3..caa7357 100644 --- a/makefile +++ b/makefile @@ -30,8 +30,8 @@ endif MOTHUR_FILES = "\"../release\"" -RELEASE_DATE = "\"12/10/2010\"" -VERSION = "\"1.15.0\"" +RELEASE_DATE = "\"2/7/2011\"" +VERSION = "\"1.16.0\"" CXXFLAGS += -DRELEASE_DATE=${RELEASE_DATE} -DVERSION=${VERSION} diff --git a/summarysharedcommand.cpp b/summarysharedcommand.cpp index 029ca26..1211ea9 100644 --- a/summarysharedcommand.cpp +++ b/summarysharedcommand.cpp @@ -744,7 +744,7 @@ int SummarySharedCommand::driver(vector thisLookup, int sta outputFileHandle << '\t'; sumCalculators[i]->print(outputFileHandle); - seqDist temp(l, k, tempdata[0]); + seqDist temp(l, k, (1.0 - tempdata[0])); calcDists[i].push_back(temp); } outputFileHandle << endl; -- 2.39.2