X-Git-Url: https://git.donarmstrong.com/?p=mothur.git;a=blobdiff_plain;f=homovacommand.cpp;h=488b68e0b6b69a8eb8e351d39be9ff09193fb36f;hp=52456b4ddedcd7b1db62a4fa8d03c5082eedf932;hb=cf9987b67aa49777a4c91c2d21f96e58bf17aa82;hpb=7b287636ea67fb2254b119c21b4057a177d3ce19 diff --git a/homovacommand.cpp b/homovacommand.cpp index 52456b4..488b68e 100644 --- a/homovacommand.cpp +++ b/homovacommand.cpp @@ -10,63 +10,76 @@ #include "homovacommand.h" #include "groupmap.h" #include "readphylipvector.h" +#include "sharedutilities.h" //********************************************************************************************************************** - -vector HomovaCommand::getValidParameters(){ +vector HomovaCommand::setParameters(){ try { - string Array[] = {"outputdir","iters","phylip","design","alpha", "inputdir"}; - vector 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 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; calledHelp = true; - vector 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 HomovaCommand::getRequiredParameters(){ - try { - string Array[] = {"design"}; - vector 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 HomovaCommand::getRequiredFiles(){ +HomovaCommand::HomovaCommand(){ try { - string Array[] = {}; - vector myArray (Array, Array+(sizeof(Array)/sizeof(string))); - return myArray; + abort = true; calledHelp = true; + setParameters(); + vector tempOutNames; + outputTypes["homova"] = tempOutNames; } catch(exception& e) { - m->errorOut(e, "HomovaCommand", "getRequiredFiles"); + m->errorOut(e, "HomovaCommand", "HomovaCommand"); exit(1); } } - //********************************************************************************************************************** HomovaCommand::HomovaCommand(string option) { @@ -75,11 +88,10 @@ HomovaCommand::HomovaCommand(string option) { //allow user to run help if(option == "help") { help(); abort = true; calledHelp = true; } + else if(option == "citation") { citation(); abort = true; calledHelp = true;} else { - //valid paramters for this command - string AlignArray[] = {"design","outputdir","iters","phylip","alpha", "inputdir"}; - vector myArray (AlignArray, AlignArray+(sizeof(AlignArray)/sizeof(string))); + vector myArray = setParameters(); OptionParser parser(option); map parameters = parser.getParameters(); @@ -123,31 +135,37 @@ HomovaCommand::HomovaCommand(string option) { 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; - } + 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 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; - } + //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"; } - convert(temp, iters); + m->mothurConvert(temp, iters); temp = validParameter.validFile(parameters, "alpha", false); if (temp == "not found") { temp = "0.05"; } - convert(temp, experimentwiseAlpha); + m->mothurConvert(temp, experimentwiseAlpha); + + string sets = validParameter.validFile(parameters, "sets", false); + if (sets == "not found") { sets = ""; } + else { + m->splitAtDash(sets, Sets); + } } } @@ -156,30 +174,6 @@ HomovaCommand::HomovaCommand(string option) { exit(1); } } - -//********************************************************************************************************************** - -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 outputs a .homova file. \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 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"); - exit(1); - } -} - -//********************************************************************************************************************** - -HomovaCommand::~HomovaCommand(){} - //********************************************************************************************************************** int HomovaCommand::execute(){ @@ -197,6 +191,31 @@ int HomovaCommand::execute(){ ReadPhylipVector readMatrix(phylipFileName); vector 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 dGroups = designMap->getNamesOfGroups(); + util.setGroups(Sets, dGroups); + + for(int i=0;icontrol_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 > origGroupSampleMap; for(int i=0;igetGroup(sampleNames[i])].push_back(i); + 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 { origGroupSampleMap[group].push_back(i); } } int numGroups = origGroupSampleMap.size(); + if (m->control_pressed) { delete designMap; return 0; } + //create a new filename ofstream HOMOVAFile; - string HOMOVAFileName = outputDir + m->getRootName(m->getSimpleName(phylipFileName)) + "homova"; + map 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_values" << endl; - m->mothurOut("HOMOVA\tBValue\tP-value\tSSwithin_values\n"); + 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); @@ -231,8 +258,7 @@ int HomovaCommand::execute(){ for(itA=origGroupSampleMap.begin();itA!=origGroupSampleMap.end();itA++){ itB = itA;itB++; - for(itB;itB!=origGroupSampleMap.end();itB++){ - + for(;itB!=origGroupSampleMap.end();itB++){ map > pairwiseGroupSampleMap; pairwiseGroupSampleMap[itA->first] = itA->second; pairwiseGroupSampleMap[itB->first] = itB->second; @@ -282,7 +308,7 @@ double HomovaCommand::runHOMOVA(ofstream& HOMOVAFile, map > vector ssWithinRandVector; map > randomizedGroup = getRandomizedGroups(groupSampleMap); double bValueRand = calcBValue(randomizedGroup, ssWithinRandVector); - if(bValueRand > bValueOrig){ counter++; } + if(bValueRand >= bValueOrig){ counter++; } } double pValue = (double) counter / (double) iters; @@ -296,7 +322,7 @@ double HomovaCommand::runHOMOVA(ofstream& HOMOVAFile, map > HOMOVAFile << it->first; m->mothurOut(it->first); it++; - for(it;it!=groupSampleMap.end();it++){ + for(;it!=groupSampleMap.end();it++){ HOMOVAFile << '-' << it->first; m->mothurOut('-' + it->first); } @@ -379,6 +405,8 @@ double HomovaCommand::calcBValue(map > groupSampleMap, vecto secondTermSum += (numSamplesInGroup - 1) * log(ssWithinVector[index] / (double)(numSamplesInGroup - 1)); inverseOneMinusSum += 1.0 / (double)(numSamplesInGroup - 1); + + ssWithinVector[index] /= (double)(numSamplesInGroup - 1); //this line is only for output purposes to scale SSw by the number of samples in the group index++; }