#include "homovacommand.h"
#include "groupmap.h"
#include "readphylipvector.h"
+#include "sharedutilities.h"
//**********************************************************************************************************************
vector<string> HomovaCommand::setParameters(){
try {
CommandParameter pdesign("design", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(pdesign);
CommandParameter pphylip("phylip", "InputTypes", "", "", "none", "none", "none",false,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);
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, and alpha. The phylip and design parameters are required, unless valid current files exist.\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";
exit(1);
}
}
-
//**********************************************************************************************************************
-
+string HomovaCommand::getOutputFileNameTag(string type, string inputName=""){
+ try {
+ string outputFileName = "";
+ map<string, vector<string> >::iterator it;
+
+ //is this a type this command creates
+ it = outputTypes.find(type);
+ if (it == outputTypes.end()) { m->mothurOut("[ERROR]: this command doesn't create a " + type + " output file.\n"); }
+ else {
+ if (type == "homova") { outputFileName = "homova"; }
+ else { m->mothurOut("[ERROR]: No definition for type " + type + " output file tag.\n"); m->control_pressed = true; }
+ }
+ return outputFileName;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "HomovaCommand", "getOutputFileNameTag");
+ exit(1);
+ }
+}
+//**********************************************************************************************************************
HomovaCommand::HomovaCommand(){
try {
abort = true; calledHelp = true;
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);
+ }
}
}
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];
//link designMap to rows/columns in distance matrix
map<string, vector<int> > origGroupSampleMap;
for(int i=0;i<sampleNames.size();i++){
- origGroupSampleMap[designMap->getGroup(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";
+ string HOMOVAFileName = outputDir + m->getRootName(m->getSimpleName(phylipFileName)) + getOutputFileNameTag("homova");
m->openOutputFile(HOMOVAFileName, HOMOVAFile);
outputNames.push_back(HOMOVAFileName); outputTypes["homova"].push_back(HOMOVAFileName);
projects / mothur.git / blobdiff