//**********************************************************************************************************************
vector<string> MGClusterCommand::setParameters(){
try {
- CommandParameter pblast("blast", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(pblast);
- CommandParameter pname("name", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(pname);
- CommandParameter plength("length", "Number", "", "5", "", "", "",false,false); parameters.push_back(plength);
- CommandParameter ppenalty("penalty", "Number", "", "0.10", "", "", "",false,false); parameters.push_back(ppenalty);
- CommandParameter pcutoff("cutoff", "Number", "", "0.70", "", "", "",false,false); parameters.push_back(pcutoff);
- CommandParameter pprecision("precision", "Number", "", "100", "", "", "",false,false); parameters.push_back(pprecision);
- CommandParameter pmethod("method", "Multiple", "furthest-nearest-average", "average", "", "", "",false,false); parameters.push_back(pmethod);
- CommandParameter phard("hard", "Boolean", "", "T", "", "", "",false,false); parameters.push_back(phard);
- CommandParameter pmin("min", "Boolean", "", "T", "", "", "",false,false); parameters.push_back(pmin);
- CommandParameter pmerge("merge", "Boolean", "", "T", "", "", "",false,false); parameters.push_back(pmerge);
- CommandParameter phcluster("hcluster", "Boolean", "", "F", "", "", "",false,false); parameters.push_back(phcluster);
- CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
- CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
+ CommandParameter pblast("blast", "InputTypes", "", "", "none", "none", "none","list",false,true,true); parameters.push_back(pblast);
+ CommandParameter pname("name", "InputTypes", "", "", "NameCount", "none", "ColumnName","rabund-sabund",false,false,true); parameters.push_back(pname);
+ CommandParameter pcount("count", "InputTypes", "", "", "NameCount", "none", "none","",false,false,true); parameters.push_back(pcount);
+ CommandParameter plength("length", "Number", "", "5", "", "", "","",false,false); parameters.push_back(plength);
+ CommandParameter ppenalty("penalty", "Number", "", "0.10", "", "", "","",false,false); parameters.push_back(ppenalty);
+ CommandParameter pcutoff("cutoff", "Number", "", "0.70", "", "", "","",false,false,true); parameters.push_back(pcutoff);
+ CommandParameter pprecision("precision", "Number", "", "100", "", "", "","",false,false); parameters.push_back(pprecision);
+ CommandParameter pmethod("method", "Multiple", "furthest-nearest-average", "average", "", "", "","",false,false); parameters.push_back(pmethod);
+ CommandParameter phard("hard", "Boolean", "", "T", "", "", "","",false,false); parameters.push_back(phard);
+ CommandParameter pmin("min", "Boolean", "", "T", "", "", "","",false,false); parameters.push_back(pmin);
+ CommandParameter pmerge("merge", "Boolean", "", "T", "", "", "","",false,false); parameters.push_back(pmerge);
+ CommandParameter padjust("adjust", "String", "", "F", "", "", "","",false,false); parameters.push_back(padjust);
+ CommandParameter phcluster("hcluster", "Boolean", "", "F", "", "", "","",false,false); parameters.push_back(phcluster);
+ 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); }
string MGClusterCommand::getHelpString(){
try {
string helpString = "";
- helpString += "The mgcluster command parameter options are blast, name, cutoff, precision, hard, method, merge, min, length, penalty and hcluster. The blast parameter is required.\n";
+ helpString += "The mgcluster command parameter options are blast, name, cutoff, precision, hard, method, merge, min, length, penalty, adjust and hcluster. The blast parameter is required.\n";
helpString += "The mgcluster command reads a blast and name file and clusters the sequences into OPF units similiar to the OTUs.\n";
helpString += "This command outputs a .list, .rabund and .sabund file that can be used with mothur other commands to estimate richness.\n";
helpString += "The cutoff parameter is used to specify the maximum distance you would like to cluster to. The default is 0.70.\n";
helpString += "The acceptable mgcluster methods are furthest, nearest and average. If no method is provided then average is assumed.\n";
helpString += "The min parameter allows you to specify is you want the minimum or maximum blast score ratio used in calculating the distance. The default is true, meaning you want the minimum.\n";
helpString += "The length parameter is used to specify the minimum overlap required. The default is 5.\n";
+ helpString += "The adjust parameter is used to handle missing distances. If you set a cutoff, adjust=f by default. If not, adjust=t by default. Adjust=f, means ignore missing distances and adjust cutoff as needed with the average neighbor method. Adjust=t, will treat missing distances as 1.0. You can also set the value the missing distances should be set to, adjust=0.5 would give missing distances a value of 0.5.\n";
helpString += "The penalty parameter is used to adjust the error rate. The default is 0.10.\n";
helpString += "The merge parameter allows you to shut off merging based on overlaps and just cluster. By default merge is true, meaning you want to merge.\n";
helpString += "The hcluster parameter allows you to use the hcluster algorithm when clustering. This may be neccessary if your file is too large to fit into RAM. The default is false.\n";
}
}
//**********************************************************************************************************************
-string MGClusterCommand::getOutputFileNameTag(string type, string inputName=""){
- try {
- string outputFileName = "";
- map<string, vector<string> >::iterator it;
+string MGClusterCommand::getOutputPattern(string type) {
+ try {
+ string pattern = "";
- //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 == "list") { outputFileName = "list"; }
- else if (type == "rabund") { outputFileName = "rabund"; }
- else if (type == "sabund") { outputFileName = "sabund"; }
- 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, "MGClusterCommand", "getOutputFileNameTag");
- exit(1);
- }
+ if (type == "list") { pattern = "[filename],[clustertag],list-[filename],[clustertag],[tag2],list"; }
+ else if (type == "rabund") { pattern = "[filename],[clustertag],rabund"; }
+ else if (type == "sabund") { pattern = "[filename],[clustertag],sabund"; }
+ else { m->mothurOut("[ERROR]: No definition for type " + type + " output pattern.\n"); m->control_pressed = true; }
+
+ return pattern;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "MGClusterCommand", "getOutputPattern");
+ exit(1);
+ }
}
-//**********************************************************************************************************************
+//*******************************************************************************************************************
MGClusterCommand::MGClusterCommand(){
try {
abort = true; calledHelp = true;
//if the user has not given a path then, add inputdir. else leave path alone.
if (path == "") { parameters["name"] = inputDir + it->second; }
}
+
+ it = parameters.find("count");
+ //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["count"] = inputDir + it->second; }
+ }
}
if (namefile == "not open") { abort = true; }
else if (namefile == "not found") { namefile = ""; }
else { m->setNameFile(namefile); }
+
+ countfile = validParameter.validFile(parameters, "count", true);
+ if (countfile == "not open") { abort = true; }
+ else if (countfile == "not found") { countfile = ""; }
+ else { m->setCountTableFile(countfile); }
+
+ if (countfile != "" && namefile != "") { m->mothurOut("[ERROR]: Cannot have both a name file and count file. Please use one or the other."); m->mothurOutEndLine(); abort = true; }
if ((blastfile == "")) { m->mothurOut("When executing a mgcluster command you must provide a blastfile."); m->mothurOutEndLine(); abort = true; }
//check for optional parameter and set defaults
string temp;
- temp = validParameter.validFile(parameters, "precision", false); if (temp == "not found") { temp = "100"; }
+ temp = validParameter.validFile(parameters, "precision", false); if (temp == "not found") { temp = "100"; }
precisionLength = temp.length();
m->mothurConvert(temp, precision);
- temp = validParameter.validFile(parameters, "cutoff", false); if (temp == "not found") { temp = "0.70"; }
+ cutoffSet = false;
+ temp = validParameter.validFile(parameters, "cutoff", false);
+ if (temp == "not found") { temp = "0.70"; }
+ else { cutoffSet = true; }
m->mothurConvert(temp, cutoff);
cutoff += (5 / (precision * 10.0));
temp = validParameter.validFile(parameters, "hard", false); if (temp == "not found") { temp = "T"; }
hard = m->isTrue(temp);
+
+ temp = validParameter.validFile(parameters, "adjust", false); if (temp == "not found") { if (cutoffSet) { temp = "F"; }else { temp="T"; } }
+ if (m->isNumeric1(temp)) { m->mothurConvert(temp, adjust); }
+ else if (m->isTrue(temp)) { adjust = 1.0; }
+ else { adjust = -1.0; }
}
}
//**********************************************************************************************************************
int MGClusterCommand::execute(){
try {
-
if (abort == true) { if (calledHelp) { return 0; } return 2; }
//read names file
if (namefile != "") {
nameMap = new NameAssignment(namefile);
nameMap->readMap();
- }else{ nameMap= new NameAssignment(); }
+ }else if (countfile != "") {
+ ct = new CountTable();
+ ct->readTable(countfile, false);
+ nameMap= new NameAssignment();
+ vector<string> tempNames = ct->getNamesOfSeqs();
+ for (int i = 0; i < tempNames.size(); i++) { nameMap->push_back(tempNames[i]); }
+ }else{ nameMap= new NameAssignment(); }
string fileroot = outputDir + m->getRootName(m->getSimpleName(blastfile));
string tag = "";
time_t start;
float previousDist = 0.00000;
float rndPreviousDist = 0.00000;
-
+
//read blastfile - creates sparsematrices for the distances and overlaps as well as a listvector
//must remember to delete those objects here since readBlast does not
read = new ReadBlast(blastfile, cutoff, penalty, length, minWanted, hclusterWanted);
read->read(nameMap);
-
- list = new ListVector(nameMap->getListVector());
- RAbundVector* rabund = new RAbundVector(list->getRAbundVector());
+
+ list = new ListVector(nameMap->getListVector());
+ RAbundVector* rabund = NULL;
+
+ if(countfile != "") {
+ rabund = new RAbundVector();
+ createRabund(ct, list, rabund);
+ }else {
+ rabund = new RAbundVector(list->getRAbundVector());
+ }
+
+
+ //list = new ListVector(nameMap->getListVector());
+ //rabund = new RAbundVector(list->getRAbundVector());
if (m->control_pressed) { outputTypes.clear(); delete nameMap; delete read; delete list; delete rabund; return 0; }
else if (method == "nearest") { tag = "nn"; }
else { tag = "an"; }
- string sabundFileName = fileroot+ tag + "." + getOutputFileNameTag("sabund");
- string rabundFileName = fileroot+ tag + "." + getOutputFileNameTag("rabund");
- string listFileName = fileroot+ tag + "." + getOutputFileNameTag("list");
+ map<string, string> variables;
+ variables["[filename]"] = fileroot;
+ variables["[clustertag]"] = tag;
+ string sabundFileName = getOutputFileName("sabund", variables);
+ string rabundFileName = getOutputFileName("rabund", variables);
+ if (countfile != "") { variables["[tag2]"] = "unique_list"; }
+ string listFileName = getOutputFileName("list", variables);
- m->openOutputFile(sabundFileName, sabundFile);
- m->openOutputFile(rabundFileName, rabundFile);
+ if (countfile == "") {
+ m->openOutputFile(sabundFileName, sabundFile);
+ m->openOutputFile(rabundFileName, rabundFile);
+ }
m->openOutputFile(listFileName, listFile);
if (m->control_pressed) {
delete nameMap; delete read; delete list; delete rabund;
- listFile.close(); rabundFile.close(); sabundFile.close(); m->mothurRemove((fileroot+ tag + ".list")); m->mothurRemove((fileroot+ tag + ".rabund")); m->mothurRemove((fileroot+ tag + ".sabund"));
+ listFile.close(); if (countfile == "") { rabundFile.close(); sabundFile.close(); m->mothurRemove((fileroot+ tag + ".rabund")); m->mothurRemove((fileroot+ tag + ".sabund")); } m->mothurRemove((fileroot+ tag + ".list"));
outputTypes.clear();
return 0;
}
if (!hclusterWanted) {
//get distmatrix and overlap
- SparseMatrix* distMatrix = read->getDistMatrix();
+ SparseDistanceMatrix* distMatrix = read->getDistMatrix();
overlapMatrix = read->getOverlapMatrix(); //already sorted by read
delete read;
//create cluster
- if (method == "furthest") { cluster = new CompleteLinkage(rabund, list, distMatrix, cutoff, method); }
- else if(method == "nearest"){ cluster = new SingleLinkage(rabund, list, distMatrix, cutoff, method); }
- else if(method == "average"){ cluster = new AverageLinkage(rabund, list, distMatrix, cutoff, method); }
+ if (method == "furthest") { cluster = new CompleteLinkage(rabund, list, distMatrix, cutoff, method, adjust); }
+ else if(method == "nearest"){ cluster = new SingleLinkage(rabund, list, distMatrix, cutoff, method, adjust); }
+ else if(method == "average"){ cluster = new AverageLinkage(rabund, list, distMatrix, cutoff, method, adjust); }
cluster->setMapWanted(true);
Seq2Bin = cluster->getSeqtoBin();
oldSeq2Bin = Seq2Bin;
if (m->control_pressed) {
delete nameMap; delete distMatrix; delete list; delete rabund; delete cluster;
- listFile.close(); rabundFile.close(); sabundFile.close(); m->mothurRemove((fileroot+ tag + ".list")); m->mothurRemove((fileroot+ tag + ".rabund")); m->mothurRemove((fileroot+ tag + ".sabund"));
+ listFile.close(); if (countfile == "") { rabundFile.close(); sabundFile.close(); m->mothurRemove((fileroot+ tag + ".rabund")); m->mothurRemove((fileroot+ tag + ".sabund")); } m->mothurRemove((fileroot+ tag + ".list"));
outputTypes.clear();
return 0;
}
-
+
+
//cluster using cluster classes
while (distMatrix->getSmallDist() < cutoff && distMatrix->getNNodes() > 0){
+ if (m->debug) { cout << "numNodes=" << distMatrix->getNNodes() << " smallDist = " << distMatrix->getSmallDist() << endl; }
+
cluster->update(cutoff);
if (m->control_pressed) {
delete nameMap; delete distMatrix; delete list; delete rabund; delete cluster;
- listFile.close(); rabundFile.close(); sabundFile.close(); m->mothurRemove((fileroot+ tag + ".list")); m->mothurRemove((fileroot+ tag + ".rabund")); m->mothurRemove((fileroot+ tag + ".sabund"));
+ listFile.close(); if (countfile == "") { rabundFile.close(); sabundFile.close(); m->mothurRemove((fileroot+ tag + ".rabund")); m->mothurRemove((fileroot+ tag + ".sabund")); } m->mothurRemove((fileroot+ tag + ".list"));
outputTypes.clear();
return 0;
}
if (m->control_pressed) {
delete nameMap; delete distMatrix; delete list; delete rabund; delete cluster; delete temp;
- listFile.close(); rabundFile.close(); sabundFile.close(); m->mothurRemove((fileroot+ tag + ".list")); m->mothurRemove((fileroot+ tag + ".rabund")); m->mothurRemove((fileroot+ tag + ".sabund"));
+ listFile.close(); if (countfile == "") { rabundFile.close(); sabundFile.close(); m->mothurRemove((fileroot+ tag + ".rabund")); m->mothurRemove((fileroot+ tag + ".sabund")); } m->mothurRemove((fileroot+ tag + ".list"));
outputTypes.clear();
return 0;
}
if (m->control_pressed) {
delete nameMap; delete distMatrix; delete list; delete rabund; delete cluster; delete temp;
- listFile.close(); rabundFile.close(); sabundFile.close(); m->mothurRemove((fileroot+ tag + ".list")); m->mothurRemove((fileroot+ tag + ".rabund")); m->mothurRemove((fileroot+ tag + ".sabund"));
+ listFile.close(); if (countfile == "") { rabundFile.close(); sabundFile.close(); m->mothurRemove((fileroot+ tag + ".rabund")); m->mothurRemove((fileroot+ tag + ".sabund")); } m->mothurRemove((fileroot+ tag + ".list"));
outputTypes.clear();
return 0;
}
if (m->control_pressed) {
delete nameMap; delete list; delete rabund;
- listFile.close(); rabundFile.close(); sabundFile.close(); m->mothurRemove((fileroot+ tag + ".list")); m->mothurRemove((fileroot+ tag + ".rabund")); m->mothurRemove((fileroot+ tag + ".sabund"));
+ listFile.close(); if (countfile == "") { rabundFile.close(); sabundFile.close(); m->mothurRemove((fileroot+ tag + ".rabund")); m->mothurRemove((fileroot+ tag + ".sabund")); } m->mothurRemove((fileroot+ tag + ".list"));
outputTypes.clear();
return 0;
}
if (m->control_pressed) {
delete nameMap; delete list; delete rabund; delete hcluster;
- listFile.close(); rabundFile.close(); sabundFile.close(); m->mothurRemove((fileroot+ tag + ".list")); m->mothurRemove((fileroot+ tag + ".rabund")); m->mothurRemove((fileroot+ tag + ".sabund"));
+ listFile.close(); if (countfile == "") { rabundFile.close(); sabundFile.close(); m->mothurRemove((fileroot+ tag + ".rabund")); m->mothurRemove((fileroot+ tag + ".sabund")); } m->mothurRemove((fileroot+ tag + ".list"));
outputTypes.clear();
return 0;
}
if (m->control_pressed) {
delete nameMap; delete list; delete rabund; delete hcluster;
- listFile.close(); rabundFile.close(); sabundFile.close(); m->mothurRemove((fileroot+ tag + ".list")); m->mothurRemove((fileroot+ tag + ".rabund")); m->mothurRemove((fileroot+ tag + ".sabund"));
+ listFile.close(); if (countfile == "") { rabundFile.close(); sabundFile.close(); m->mothurRemove((fileroot+ tag + ".rabund")); m->mothurRemove((fileroot+ tag + ".sabund")); } m->mothurRemove((fileroot+ tag + ".list"));
m->mothurRemove(distFile);
m->mothurRemove(overlapFile);
outputTypes.clear();
if (m->control_pressed) {
delete nameMap; delete list; delete rabund; delete hcluster;
- listFile.close(); rabundFile.close(); sabundFile.close(); m->mothurRemove((fileroot+ tag + ".list")); m->mothurRemove((fileroot+ tag + ".rabund")); m->mothurRemove((fileroot+ tag + ".sabund"));
+ listFile.close(); if (countfile == "") { rabundFile.close(); sabundFile.close(); m->mothurRemove((fileroot+ tag + ".rabund")); m->mothurRemove((fileroot+ tag + ".sabund")); } m->mothurRemove((fileroot+ tag + ".list"));
m->mothurRemove(distFile);
m->mothurRemove(overlapFile);
outputTypes.clear();
if (m->control_pressed) {
delete nameMap; delete list; delete rabund; delete hcluster; delete temp;
- listFile.close(); rabundFile.close(); sabundFile.close(); m->mothurRemove((fileroot+ tag + ".list")); m->mothurRemove((fileroot+ tag + ".rabund")); m->mothurRemove((fileroot+ tag + ".sabund"));
+ listFile.close(); if (countfile == "") { rabundFile.close(); sabundFile.close(); m->mothurRemove((fileroot+ tag + ".rabund")); m->mothurRemove((fileroot+ tag + ".sabund")); } m->mothurRemove((fileroot+ tag + ".list"));
m->mothurRemove(distFile);
m->mothurRemove(overlapFile);
outputTypes.clear();
if (m->control_pressed) {
delete nameMap; delete list; delete rabund; delete hcluster; delete temp;
- listFile.close(); rabundFile.close(); sabundFile.close(); m->mothurRemove((fileroot+ tag + ".list")); m->mothurRemove((fileroot+ tag + ".rabund")); m->mothurRemove((fileroot+ tag + ".sabund"));
+ listFile.close(); if (countfile == "") { rabundFile.close(); sabundFile.close(); m->mothurRemove((fileroot+ tag + ".rabund")); m->mothurRemove((fileroot+ tag + ".sabund")); } m->mothurRemove((fileroot+ tag + ".list"));
m->mothurRemove(distFile);
m->mothurRemove(overlapFile);
outputTypes.clear();
m->mothurRemove(overlapFile);
}
- delete list;
+ delete list;
delete rabund;
listFile.close();
- sabundFile.close();
- rabundFile.close();
-
+ if (countfile == "") {
+ sabundFile.close();
+ rabundFile.close();
+ }
if (m->control_pressed) {
delete nameMap;
- listFile.close(); rabundFile.close(); sabundFile.close(); m->mothurRemove((fileroot+ tag + ".list")); m->mothurRemove((fileroot+ tag + ".rabund")); m->mothurRemove((fileroot+ tag + ".sabund"));
+ listFile.close(); if (countfile == "") { rabundFile.close(); sabundFile.close(); m->mothurRemove((fileroot+ tag + ".rabund")); m->mothurRemove((fileroot+ tag + ".sabund")); } m->mothurRemove((fileroot+ tag + ".list"));
outputTypes.clear();
return 0;
}
m->mothurOutEndLine();
m->mothurOut("Output File Names: "); m->mothurOutEndLine();
m->mothurOut(listFileName); m->mothurOutEndLine(); outputNames.push_back(listFileName); outputTypes["list"].push_back(listFileName);
- m->mothurOut(rabundFileName); m->mothurOutEndLine(); outputNames.push_back(rabundFileName); outputTypes["rabund"].push_back(rabundFileName);
- m->mothurOut(sabundFileName); m->mothurOutEndLine(); outputNames.push_back(sabundFileName); outputTypes["sabund"].push_back(sabundFileName);
+ if (countfile == "") {
+ m->mothurOut(rabundFileName); m->mothurOutEndLine(); outputNames.push_back(rabundFileName); outputTypes["rabund"].push_back(rabundFileName);
+ m->mothurOut(sabundFileName); m->mothurOutEndLine(); outputNames.push_back(sabundFileName); outputTypes["sabund"].push_back(sabundFileName);
+ }
m->mothurOutEndLine();
if (saveCutoff != cutoff) {
void MGClusterCommand::printData(ListVector* mergedList){
try {
mergedList->print(listFile);
- mergedList->getRAbundVector().print(rabundFile);
-
- SAbundVector sabund = mergedList->getSAbundVector();
+ SAbundVector sabund = mergedList->getSAbundVector();
+
+ if (countfile == "") {
+ mergedList->getRAbundVector().print(rabundFile);
+ sabund.print(sabundFile);
+ }
sabund.print(cout);
- sabund.print(sabundFile);
}
catch(exception& e) {
m->errorOut(e, "MGClusterCommand", "printData");
//**********************************************************************************************************************
+void MGClusterCommand::createRabund(CountTable*& ct, ListVector*& list, RAbundVector*& rabund){
+ try {
+ //vector<string> names = ct.getNamesOfSeqs();
+ //for ( int i; i < ct.getNumGroups(); i++ ) { rav.push_back( ct.getNumSeqs(names[i]) ); }
+ //return rav;
+
+ for(int i = 0; i < list->getNumBins(); i++) {
+ vector<string> binNames;
+ string bin = list->get(i);
+ m->splitAtComma(bin, binNames);
+ int total = 0;
+ for (int j = 0; j < binNames.size(); j++) {
+ total += ct->getNumSeqs(binNames[j]);
+ }
+ rabund->push_back(total);
+ }
+
+
+ }
+ catch(exception& e) {
+ m->errorOut(e, "MGClusterCommand", "createRabund");
+ exit(1);
+ }
+
+}
+//**********************************************************************************************************************
projects / mothur.git / blobdiff