#include "chimeraslayercommand.h"
#include "chimeraslayer.h"
#include "deconvolutecommand.h"
+#include "referencedb.h"
//**********************************************************************************************************************
vector<string> ChimeraSlayerCommand::setParameters(){
CommandParameter pmincov("mincov", "Number", "", "70", "", "", "",false,false); parameters.push_back(pmincov);
CommandParameter pminsnp("minsnp", "Number", "", "10", "", "", "",false,false); parameters.push_back(pminsnp);
CommandParameter pminbs("minbs", "Number", "", "90", "", "", "",false,false); parameters.push_back(pminbs);
- CommandParameter psearch("search", "Multiple", "kmer-blast-distance", "blast", "", "", "",false,false); parameters.push_back(psearch);
+ CommandParameter psearch("search", "Multiple", "kmer-blast", "blast", "", "", "",false,false); parameters.push_back(psearch);
CommandParameter pprocessors("processors", "Number", "", "1", "", "", "",false,false); parameters.push_back(pprocessors);
CommandParameter prealign("realign", "Boolean", "", "T", "", "", "",false,false); parameters.push_back(prealign);
CommandParameter ptrim("trim", "Boolean", "", "F", "", "", "",false,false); parameters.push_back(ptrim);
CommandParameter pincrement("increment", "Number", "", "5", "", "", "",false,false); parameters.push_back(pincrement);
CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
-
+ CommandParameter psave("save", "Boolean", "", "F", "", "", "",false,false); parameters.push_back(psave);
+
vector<string> myArray;
for (int i = 0; i < parameters.size(); i++) { myArray.push_back(parameters[i].name); }
return myArray;
helpString += "This command was modeled after the chimeraSlayer written by the Broad Institute.\n";
helpString += "The chimera.slayer command parameters are fasta, name, template, processors, trim, ksize, window, match, mismatch, divergence. minsim, mincov, minbs, minsnp, parents, search, iters, increment, numwanted, and realign.\n";
helpString += "The fasta parameter allows you to enter the fasta file containing your potentially chimeric sequences, and is required, unless you have a valid current fasta file. \n";
- helpString += "The name parameter allows you to provide a name file, if you are using template=self. \n";
+ helpString += "The name parameter allows you to provide a name file, if you are using reference=self. \n";
helpString += "You may enter multiple fasta files by separating their names with dashes. ie. fasta=abrecovery.fasta-amazon.fasta \n";
helpString += "The reference parameter allows you to enter a reference file containing known non-chimeric sequences, and is required. You may also set template=self, in this case the abundant sequences will be used as potential parents. \n";
helpString += "The processors parameter allows you to specify how many processors you would like to use. The default is 1. \n";
helpString += "The mincov parameter allows you to specify minimum coverage by closest matches found in template. Default is 70, meaning 70%. \n";
helpString += "The minbs parameter allows you to specify minimum bootstrap support for calling a sequence chimeric. Default is 90, meaning 90%. \n";
helpString += "The minsnp parameter allows you to specify percent of SNPs to sample on each side of breakpoint for computing bootstrap support (default: 10) \n";
- helpString += "The search parameter allows you to specify search method for finding the closest parent. Choices are distance, blast, and kmer, default distance. \n";
+ helpString += "The search parameter allows you to specify search method for finding the closest parent. Choices are blast, and kmer, default blast. \n";
helpString += "The realign parameter allows you to realign the query to the potential parents. Choices are true or false, default true. \n";
+ helpString += "If the save parameter is set to true the reference sequences will be saved in memory, to clear them later you can use the clear.memory command. Default=f.";
helpString += "The chimera.slayer command should be in the following format: \n";
helpString += "chimera.slayer(fasta=yourFastaFile, reference=yourTemplate, search=yourSearch) \n";
helpString += "Example: chimera.slayer(fasta=AD.align, reference=core_set_aligned.imputed.fasta, search=kmer) \n";
ChimeraSlayerCommand::ChimeraSlayerCommand(string option) {
try {
abort = false; calledHelp = false;
+ ReferenceDB* rdb = ReferenceDB::getInstance();
//allow user to run help
if(option == "help") { help(); abort = true; calledHelp = true; }
+ else if(option == "citation") { citation(); abort = true; calledHelp = true;}
else {
vector<string> myArray = setParameters();
//erase from file list
fastaFileNames.erase(fastaFileNames.begin()+i);
i--;
+ }else {
+ m->setFastaFile(fastaFileNames[i]);
}
}
}
//erase from file list
nameFileNames.erase(nameFileNames.begin()+i);
i--;
+ }else {
+ m->setNameFile(nameFileNames[i]);
}
}
}
//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 = ""; }
+ string temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); }
+ m->setProcessors(temp);
+ convert(temp, processors);
+
+ temp = validParameter.validFile(parameters, "save", false); if (temp == "not found"){ temp = "f"; }
+ save = m->isTrue(temp);
+ rdb->save = save;
+ if (save) { //clear out old references
+ rdb->clearMemory();
+ }
string path;
it = parameters.find("reference");
//user has given a template file
if(it != parameters.end()){
- if (it->second == "self") { templatefile = "self"; }
+ if (it->second == "self") {
+ templatefile = "self";
+ if (save) {
+ m->mothurOut("[WARNING]: You can't save reference=self, ignoring save.");
+ m->mothurOutEndLine();
+ save = false;
+ }
+ }
else {
path = m->hasPath(it->second);
//if the user has not given a path then, add inputdir. else leave path alone.
templatefile = validParameter.validFile(parameters, "reference", true);
if (templatefile == "not open") { abort = true; }
- else if (templatefile == "not found") { templatefile = ""; m->mothurOut("reference is a required parameter for the chimera.slayer command."); m->mothurOutEndLine(); abort = true; }
+ else if (templatefile == "not found") { //check for saved reference sequences
+ if (rdb->referenceSeqs.size() != 0) {
+ templatefile = "saved";
+ }else {
+ m->mothurOut("[ERROR]: You don't have any saved reference sequences and the reference parameter is a required.");
+ m->mothurOutEndLine();
+ abort = true;
+ }
+ }else { if (save) { rdb->setSavedReference(templatefile); } }
+ }
+ }else if (hasName) { templatefile = "self";
+ if (save) {
+ m->mothurOut("[WARNING]: You can't save reference=self, ignoring save.");
+ m->mothurOutEndLine();
+ save = false;
}
}
+ else {
+ if (rdb->referenceSeqs.size() != 0) {
+ templatefile = "saved";
+ }else {
+ m->mothurOut("[ERROR]: You don't have any saved reference sequences and the reference parameter is a required.");
+ m->mothurOutEndLine();
+ templatefile = ""; abort = true;
+ }
+ }
+
- string temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); }
- m->setProcessors(temp);
- convert(temp, processors);
temp = validParameter.validFile(parameters, "ksize", false); if (temp == "not found") { temp = "7"; }
convert(temp, ksize);
temp = validParameter.validFile(parameters, "split", false); if (temp == "not found") { temp = "f"; }
trimera = m->isTrue(temp);
- search = validParameter.validFile(parameters, "search", false); if (search == "not found") { search = "distance"; }
+ search = validParameter.validFile(parameters, "search", false); if (search == "not found") { search = "blast"; }
temp = validParameter.validFile(parameters, "iters", false); if (temp == "not found") { temp = "1000"; }
convert(temp, iters);
temp = validParameter.validFile(parameters, "numwanted", false); if (temp == "not found") { temp = "15"; }
convert(temp, numwanted);
- if ((search != "distance") && (search != "blast") && (search != "kmer")) { m->mothurOut(search + " is not a valid search."); m->mothurOutEndLine(); abort = true; }
+ if ((search != "blast") && (search != "kmer")) { m->mothurOut(search + " is not a valid search."); m->mothurOutEndLine(); abort = true; }
+
+ if (hasName && (templatefile != "self")) { m->mothurOut("You have provided a namefile and the reference parameter is not set to self. I am not sure what reference you are trying to use, aborting."); m->mothurOutEndLine(); abort=true; }
}
}
catch(exception& e) {
int ChimeraSlayerCommand::execute(){
try{
if (abort == true) { if (calledHelp) { return 0; } return 2; }
-
+
for (int s = 0; s < fastaFileNames.size(); s++) {
m->mothurOut("Checking sequences from " + fastaFileNames[s] + " ..." ); m->mothurOutEndLine();
//do your part
driverMPI(startIndex, numSeqsPerProcessor, inMPI, outMPI, outMPIAccnos, outMPIFasta, MPIPos);
+ int numNoParents = chimera->getNumNoParents();
+ int temp;
+ for(int i = 1; i < processors; i++) {
+ MPI_Recv(&temp, 1, MPI_INT, 1, tag, MPI_COMM_WORLD, &status);
+ numNoParents += temp;
+ }
+
+
+ if (numSeqs == numNoParents) { m->mothurOut("[WARNING]: megablast returned 0 potential parents for all your sequences. This could be due to formatdb.exe not being setup properly, please check formatdb.log for errors."); m->mothurOutEndLine(); }
+
if (m->control_pressed) { outputTypes.clear(); MPI_File_close(&inMPI); MPI_File_close(&outMPI); if (trim) { MPI_File_close(&outMPIFasta); } MPI_File_close(&outMPIAccnos); for (int j = 0; j < outputNames.size(); j++) { remove(outputNames[j].c_str()); } remove(outputFileName.c_str()); remove(accnosFileName.c_str()); delete chimera; return 0; }
}else{ //you are a child process
//do your part
driverMPI(startIndex, numSeqsPerProcessor, inMPI, outMPI, outMPIAccnos, outMPIFasta, MPIPos);
+
+ int numNoParents = chimera->getNumNoParents();
+ MPI_Send(&numNoParents, 1, MPI_INT, 0, tag, MPI_COMM_WORLD);
if (m->control_pressed) { outputTypes.clear(); MPI_File_close(&inMPI); MPI_File_close(&outMPI); if (trim) { MPI_File_close(&outMPIFasta); } MPI_File_close(&outMPIAccnos); for (int j = 0; j < outputNames.size(); j++) { remove(outputNames[j].c_str()); } delete chimera; return 0; }
if(processors == 1){
numSeqs = driver(lines[0], outputFileName, fastaFileNames[s], accnosFileName, trimFastaFileName);
+ int numNoParents = chimera->getNumNoParents();
+ if (numNoParents == numSeqs) { m->mothurOut("[WARNING]: megablast returned 0 potential parents for all your sequences. This could be due to formatdb.exe not being setup properly, please check formatdb.log for errors."); m->mothurOutEndLine(); }
+
if (m->control_pressed) { outputTypes.clear(); if (trim) { remove(trimFastaFileName.c_str()); } remove(outputFileName.c_str()); remove(tempHeader.c_str()); remove(accnosFileName.c_str()); for (int j = 0; j < outputNames.size(); j++) { remove(outputNames[j].c_str()); } for (int i = 0; i < lines.size(); i++) { delete lines[i]; } lines.clear(); delete chimera; return 0; }
}else{
#else
numSeqs = driver(lines[0], outputFileName, fastaFileNames[s], accnosFileName, trimFastaFileName);
+ int numNoParents = chimera->getNumNoParents();
+ if (numNoParents == numSeqs) { m->mothurOut("[WARNING]: megablast returned 0 potential parents for all your sequences. This could be due to formatdb.exe not being setup properly, please check formatdb.log for errors."); m->mothurOutEndLine(); }
+
+
if (m->control_pressed) { outputTypes.clear(); if (trim) { remove(trimFastaFileName.c_str()); } remove(outputFileName.c_str()); remove(tempHeader.c_str()); remove(accnosFileName.c_str()); for (int j = 0; j < outputNames.size(); j++) { remove(outputNames[j].c_str()); } for (int i = 0; i < lines.size(); i++) { delete lines[i]; } lines.clear(); delete chimera; return 0; }
#endif
string candidateAligned = candidateSeq->getAligned();
if (candidateSeq->getName() != "") { //incase there is a commented sequence at the end of a file
-
if (candidateSeq->getAligned().length() != templateSeqsLength) {
m->mothurOut(candidateSeq->getName() + " is not the same length as the template sequences. Skipping."); m->mothurOutEndLine();
}else{
data_results rightResults = chimera->getResults();
//if either piece is chimeric then report
- Sequence* trimmed = chimera->print(out, out2, leftResults, rightResults);
- if (trim) { trimmed->printSequence(out3); delete trimmed; }
+ Sequence trimmed = chimera->print(out, out2, leftResults, rightResults);
+ if (trim) { trimmed.printSequence(out3); }
delete left; delete right;
}else { //already chimeric
//print results
- Sequence* trimmed = chimera->print(out, out2);
- if (trim) { trimmed->printSequence(out3); delete trimmed; }
+ Sequence trimmed = chimera->print(out, out2);
+ if (trim) { trimmed.printSequence(out3); }
}
#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
unsigned long int pos = inFASTA.tellg();
- //cout << candidateSeq->getName() << '\t' << pos << endl;
if ((pos == -1) || (pos >= filePos->end)) { break; }
#else
if (inFASTA.eof()) { break; }
data_results rightResults = chimera->getResults();
//if either piece is chimeric then report
- Sequence* trimmed = chimera->print(outMPI, outAccMPI, leftResults, rightResults);
+ Sequence trimmed = chimera->print(outMPI, outAccMPI, leftResults, rightResults);
if (trim) {
- string outputString = ">" + trimmed->getName() + "\n" + trimmed->getAligned() + "\n";
- delete trimmed;
+ string outputString = ">" + trimmed.getName() + "\n" + trimmed.getAligned() + "\n";
//write to accnos file
int length = outputString.length();
}else {
//print results
- Sequence* trimmed = chimera->print(outMPI, outAccMPI);
+ Sequence trimmed = chimera->print(outMPI, outAccMPI);
if (trim) {
- string outputString = ">" + trimmed->getName() + "\n" + trimmed->getAligned() + "\n";
- delete trimmed;
+ string outputString = ">" + trimmed.getName() + "\n" + trimmed.getAligned() + "\n";
//write to accnos file
int length = outputString.length();
ofstream out;
string tempFile = outputFileName + toString(getpid()) + ".num.temp";
m->openOutputFile(tempFile, out);
- out << num << endl;
+ out << num << '\t' << chimera->getNumNoParents() << endl;
out.close();
-
exit(0);
}else {
m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
wait(&temp);
}
+ int numNoParents = 0;
for (int i = 0; i < processIDS.size(); i++) {
ifstream in;
string tempFile = outputFileName + toString(processIDS[i]) + ".num.temp";
m->openInputFile(tempFile, in);
- if (!in.eof()) { int tempNum = 0; in >> tempNum; num += tempNum; }
+ if (!in.eof()) { int tempNum = 0; int tempNumParents = 0; in >> tempNum >> tempNumParents; num += tempNum; numNoParents += tempNumParents; }
in.close(); remove(tempFile.c_str());
}
+ if (num == numNoParents) { m->mothurOut("[WARNING]: megablast returned 0 potential parents for all your sequences. This could be due to formatdb.exe not being setup properly, please check formatdb.log for errors."); m->mothurOutEndLine(); }
+
return num;
#endif
}
projects / mothur.git / blobdiff