#include "refchimeratest.h"
#include "filterseqscommand.h"
+
//**********************************************************************************************************************
vector<string> SeqErrorCommand::setParameters(){
try {
CommandParameter pignorechimeras("ignorechimeras", "Boolean", "", "T", "", "", "",false,false); parameters.push_back(pignorechimeras);
CommandParameter pthreshold("threshold", "Number", "", "1.0", "", "", "",false,false); parameters.push_back(pthreshold);
CommandParameter pprocessors("processors", "Number", "", "1", "", "", "",false,false); parameters.push_back(pprocessors);
- CommandParameter pfilter("filter", "Boolean", "", "T", "", "", "",false,false); parameters.push_back(pfilter);
+ CommandParameter psave("save", "Boolean", "", "F", "", "", "",false,false); parameters.push_back(psave);
CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
try {
string helpString = "";
helpString += "The seq.error command reads a query alignment file and a reference alignment file and creates .....\n";
+ helpString += "The fasta parameter...\n";
+ helpString += "The reference parameter...\n";
+ helpString += "The qfile parameter...\n";
+ helpString += "The report parameter...\n";
+ helpString += "The name parameter...\n";
+ helpString += "The ignorechimeras parameter...\n";
+ helpString += "The threshold parameter...\n";
+ helpString += "The processors parameter...\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 += "Example seq.error(...).\n";
helpString += "Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFasta).\n";
helpString += "For more details please check out the wiki http://www.mothur.org/wiki/seq.error .\n";
SeqErrorCommand::SeqErrorCommand(){
try {
abort = true; calledHelp = true;
+ setParameters();
vector<string> tempOutNames;
outputTypes["error.summary"] = tempOutNames;
outputTypes["error.seq"] = tempOutNames;
try {
abort = false; calledHelp = false;
+ rdb = ReferenceDB::getInstance();
//allow user to run help
if(option == "help") { help(); abort = true; calledHelp = true; }
else { m->mothurOut("You have no current fasta file and the fasta parameter is required."); m->mothurOutEndLine(); abort = true; }
}
else if (queryFileName == "not open") { abort = true; }
-
- referenceFileName = validParameter.validFile(parameters, "reference", true);
- if (referenceFileName == "not found") { m->mothurOut("reference is a required parameter for the seq.error command."); m->mothurOutEndLine(); abort = true; }
- else if (referenceFileName == "not open") { abort = true; }
-
+ else { m->setFastaFile(queryFileName); }
//check for optional parameters
namesFileName = validParameter.validFile(parameters, "name", true);
if(namesFileName == "not found"){ namesFileName = ""; }
else if (namesFileName == "not open") { namesFileName = ""; abort = true; }
+ else { m->setNameFile(namesFileName); }
qualFileName = validParameter.validFile(parameters, "qfile", true);
if(qualFileName == "not found"){ qualFileName = ""; }
else if (qualFileName == "not open") { qualFileName = ""; abort = true; }
-
+ else { m->setQualFile(qualFileName); }
+
reportFileName = validParameter.validFile(parameters, "report", true);
if(reportFileName == "not found"){ reportFileName = ""; }
else if (reportFileName == "not open") { reportFileName = ""; abort = true; }
temp = validParameter.validFile(parameters, "threshold", false); if (temp == "not found") { temp = "1.00"; }
convert(temp, threshold);
+ 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();
+ }
+
+ //this has to go after save so that if the user sets save=t and provides no reference we abort
+ referenceFileName = validParameter.validFile(parameters, "reference", true);
+ if (referenceFileName == "not found") {
+ //check for saved reference sequences
+ if (rdb->referenceSeqs.size() != 0) {
+ referenceFileName = "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 (referenceFileName == "not open") { abort = true; }
+ else { if (save) { rdb->setSavedReference(referenceFileName); } }
+
+
temp = validParameter.validFile(parameters, "ignorechimeras", false); if (temp == "not found") { temp = "T"; }
ignoreChimeras = m->isTrue(temp);
- temp = validParameter.validFile(parameters, "filter", false); if (temp == "not found") { temp = "T"; }
- filter = m->isTrue(temp);
-
temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); }
m->setProcessors(temp);
convert(temp, processors);
maxLength = 2000;
totalBases = 0;
totalMatches = 0;
-
- //run vertical filter on query and reference files.
- if (filter) {
- string inputString = "fasta=" + queryFileName + "-" + referenceFileName;
- m->mothurOut("/******************************************/"); m->mothurOutEndLine();
- m->mothurOut("Running command: filter.seqs(" + inputString + ") to improve processing time."); m->mothurOutEndLine();
-
- Command* filterCommand = new FilterSeqsCommand(inputString);
- filterCommand->execute();
-
- map<string, vector<string> > filenames = filterCommand->getOutputFiles();
-
- delete filterCommand;
-
- m->mothurOut("/******************************************/"); m->mothurOutEndLine();
-
- queryFileName = filenames["fasta"][0];
- referenceFileName = filenames["fasta"][1];
- }
string errorSummaryFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.summary";
outputNames.push_back(errorSummaryFileName); outputTypes["error.summary"].push_back(errorSummaryFileName);
printErrorFRFile(errorForward, errorReverse);
- if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
+ if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; }
string errorCountFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.count";
ofstream errorCountFile;
}
errorCountFile.close();
- if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
+ if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; }
printSubMatrix();
m->mothurOut("Appending files from process " + toString(processIDS[i])); m->mothurOutEndLine();
m->appendFiles((summaryFileName + toString(processIDS[i]) + ".temp"), summaryFileName);
- remove((summaryFileName + toString(processIDS[i]) + ".temp").c_str());
+ m->mothurRemove((summaryFileName + toString(processIDS[i]) + ".temp"));
m->appendFiles((errorOutputFileName + toString(processIDS[i]) + ".temp"), errorOutputFileName);
- remove((errorOutputFileName + toString(processIDS[i]) + ".temp").c_str());
+ m->mothurRemove((errorOutputFileName + toString(processIDS[i]) + ".temp"));
m->appendFiles((chimeraOutputFileName + toString(processIDS[i]) + ".temp"), chimeraOutputFileName);
- remove((chimeraOutputFileName + toString(processIDS[i]) + ".temp").c_str());
+ m->mothurRemove((chimeraOutputFileName + toString(processIDS[i]) + ".temp"));
ifstream in;
string tempFile = filename + toString(processIDS[i]) + ".info.temp";
}
m->gobble(in);
- in.close(); remove(tempFile.c_str());
+ in.close(); m->mothurRemove(tempFile);
}
#endif
for(int i=0;i<minCompare.sequence.length();i++){
char letter = minCompare.sequence[i];
- errorForward[letter][i] += minCompare.weight;
- errorReverse[letter][minCompare.total-i-1] += minCompare.weight;
+ if(letter != 'r'){
+ errorForward[letter][i] += minCompare.weight;
+ errorReverse[letter][minCompare.total-i-1] += minCompare.weight;
+ }
}
}
quality = QualityScores(qualFile);
if(!ignoreSeq){
+// cout << query.getName() << '\t';
+
quality.updateQScoreErrorMap(qScoreErrorMap, minCompare.sequence, startBase, endBase, minCompare.weight);
quality.updateForwardMap(qualForwardMap, startBase, endBase, minCompare.weight);
quality.updateReverseMap(qualReverseMap, startBase, endBase, minCompare.weight);
+
+// cout << endl;
}
}
void SeqErrorCommand::getReferences(){
try {
-
- ifstream referenceFile;
- m->openInputFile(referenceFileName, referenceFile);
-
int numAmbigSeqs = 0;
int maxStartPos = 0;
int minEndPos = 100000;
- while(referenceFile){
- Sequence currentSeq(referenceFile);
- int numAmbigs = currentSeq.getAmbigBases();
- if(numAmbigs > 0){ numAmbigSeqs++; }
+ if (referenceFileName == "saved") {
+ int start = time(NULL);
+ m->mothurOutEndLine(); m->mothurOut("Using sequences from " + rdb->getSavedReference() + " that are saved in memory."); m->mothurOutEndLine();
+
+ for (int i = 0; i < rdb->referenceSeqs.size(); i++) {
+ int numAmbigs = rdb->referenceSeqs[i].getAmbigBases();
+ if(numAmbigs > 0){ numAmbigSeqs++; }
+
+ // int startPos = rdb->referenceSeqs[i].getStartPos();
+ // if(startPos > maxStartPos) { maxStartPos = startPos; }
+ //
+ // int endPos = rdb->referenceSeqs[i].getEndPos();
+ // if(endPos < minEndPos) { minEndPos = endPos; }
+
+ referenceSeqs.push_back(rdb->referenceSeqs[i]);
+ }
+ referenceFileName = rdb->getSavedReference();
+
+ m->mothurOut("It took " + toString(time(NULL) - start) + " to load " + toString(rdb->referenceSeqs.size()) + " sequences.");m->mothurOutEndLine();
+
+ }else {
+ int start = time(NULL);
+
+ ifstream referenceFile;
+ m->openInputFile(referenceFileName, referenceFile);
-// int startPos = currentSeq.getStartPos();
-// if(startPos > maxStartPos) { maxStartPos = startPos; }
-//
-// int endPos = currentSeq.getEndPos();
-// if(endPos < minEndPos) { minEndPos = endPos; }
- referenceSeqs.push_back(currentSeq);
+ while(referenceFile){
+ Sequence currentSeq(referenceFile);
+ int numAmbigs = currentSeq.getAmbigBases();
+ if(numAmbigs > 0){ numAmbigSeqs++; }
- m->gobble(referenceFile);
+ // int startPos = currentSeq.getStartPos();
+ // if(startPos > maxStartPos) { maxStartPos = startPos; }
+ //
+ // int endPos = currentSeq.getEndPos();
+ // if(endPos < minEndPos) { minEndPos = endPos; }
+ referenceSeqs.push_back(currentSeq);
+
+ if (rdb->save) { rdb->referenceSeqs.push_back(currentSeq); }
+
+ m->gobble(referenceFile);
+ }
+ referenceFile.close();
+
+ m->mothurOut("It took " + toString(time(NULL) - start) + " to read " + toString(referenceSeqs.size()) + " sequences.");m->mothurOutEndLine();
}
- referenceFile.close();
+
numRefs = referenceSeqs.size();
-
for(int i=0;i<numRefs;i++){
referenceSeqs[i].padToPos(maxStartPos);
if(numAmbigSeqs != 0){
m->mothurOut("Warning: " + toString(numAmbigSeqs) + " reference sequences have ambiguous bases, these bases will be ignored\n");
- }
+ }
}
catch(exception& e) {
Compare errors;
for(int i=0;i<alignLength;i++){
- if(r[i] != 'N' && q[i] != '.' && r[i] != '.' && (q[i] != '-' || r[i] != '-')){ // no missing data and no double gaps
- started = 1;
-
- if(q[i] == 'A'){
- if(r[i] == 'A'){ errors.AA++; errors.matches++; errors.sequence += 'm'; }
- if(r[i] == 'T'){ errors.AT++; errors.sequence += 's'; }
- if(r[i] == 'G'){ errors.AG++; errors.sequence += 's'; }
- if(r[i] == 'C'){ errors.AC++; errors.sequence += 's'; }
- if(r[i] == '-'){ errors.Ai++; errors.sequence += 'i'; }
- }
- else if(q[i] == 'T'){
- if(r[i] == 'A'){ errors.TA++; errors.sequence += 's'; }
- if(r[i] == 'T'){ errors.TT++; errors.matches++; errors.sequence += 'm'; }
- if(r[i] == 'G'){ errors.TG++; errors.sequence += 's'; }
- if(r[i] == 'C'){ errors.TC++; errors.sequence += 's'; }
- if(r[i] == '-'){ errors.Ti++; errors.sequence += 'i'; }
- }
- else if(q[i] == 'G'){
- if(r[i] == 'A'){ errors.GA++; errors.sequence += 's'; }
- if(r[i] == 'T'){ errors.GT++; errors.sequence += 's'; }
- if(r[i] == 'G'){ errors.GG++; errors.matches++; errors.sequence += 'm'; }
- if(r[i] == 'C'){ errors.GC++; errors.sequence += 's'; }
- if(r[i] == '-'){ errors.Gi++; errors.sequence += 'i'; }
- }
- else if(q[i] == 'C'){
- if(r[i] == 'A'){ errors.CA++; errors.sequence += 's'; }
- if(r[i] == 'T'){ errors.CT++; errors.sequence += 's'; }
- if(r[i] == 'G'){ errors.CG++; errors.sequence += 's'; }
- if(r[i] == 'C'){ errors.CC++; errors.matches++; errors.sequence += 'm'; }
- if(r[i] == '-'){ errors.Ci++; errors.sequence += 'i'; }
- }
- else if(q[i] == 'N'){
- if(r[i] == 'A'){ errors.NA++; errors.sequence += 'a'; }
- if(r[i] == 'T'){ errors.NT++; errors.sequence += 'a'; }
- if(r[i] == 'G'){ errors.NG++; errors.sequence += 'a'; }
- if(r[i] == 'C'){ errors.NC++; errors.sequence += 'a'; }
- if(r[i] == '-'){ errors.Ni++; errors.sequence += 'a'; }
+// cout << r[i] << '\t' << q[i] << '\t';
+ if(q[i] != '.' && r[i] != '.' && (q[i] != '-' || r[i] != '-')){ // no missing data and no double gaps
+ if(r[i] != 'N'){
+ started = 1;
+
+ if(q[i] == 'A'){
+ if(r[i] == 'A'){ errors.AA++; errors.matches++; errors.sequence += 'm'; }
+ if(r[i] == 'T'){ errors.AT++; errors.sequence += 's'; }
+ if(r[i] == 'G'){ errors.AG++; errors.sequence += 's'; }
+ if(r[i] == 'C'){ errors.AC++; errors.sequence += 's'; }
+ if(r[i] == '-'){ errors.Ai++; errors.sequence += 'i'; }
+ }
+ else if(q[i] == 'T'){
+ if(r[i] == 'A'){ errors.TA++; errors.sequence += 's'; }
+ if(r[i] == 'T'){ errors.TT++; errors.matches++; errors.sequence += 'm'; }
+ if(r[i] == 'G'){ errors.TG++; errors.sequence += 's'; }
+ if(r[i] == 'C'){ errors.TC++; errors.sequence += 's'; }
+ if(r[i] == '-'){ errors.Ti++; errors.sequence += 'i'; }
+ }
+ else if(q[i] == 'G'){
+ if(r[i] == 'A'){ errors.GA++; errors.sequence += 's'; }
+ if(r[i] == 'T'){ errors.GT++; errors.sequence += 's'; }
+ if(r[i] == 'G'){ errors.GG++; errors.matches++; errors.sequence += 'm'; }
+ if(r[i] == 'C'){ errors.GC++; errors.sequence += 's'; }
+ if(r[i] == '-'){ errors.Gi++; errors.sequence += 'i'; }
+ }
+ else if(q[i] == 'C'){
+ if(r[i] == 'A'){ errors.CA++; errors.sequence += 's'; }
+ if(r[i] == 'T'){ errors.CT++; errors.sequence += 's'; }
+ if(r[i] == 'G'){ errors.CG++; errors.sequence += 's'; }
+ if(r[i] == 'C'){ errors.CC++; errors.matches++; errors.sequence += 'm'; }
+ if(r[i] == '-'){ errors.Ci++; errors.sequence += 'i'; }
+ }
+ else if(q[i] == 'N'){
+ if(r[i] == 'A'){ errors.NA++; errors.sequence += 'a'; }
+ if(r[i] == 'T'){ errors.NT++; errors.sequence += 'a'; }
+ if(r[i] == 'G'){ errors.NG++; errors.sequence += 'a'; }
+ if(r[i] == 'C'){ errors.NC++; errors.sequence += 'a'; }
+ if(r[i] == '-'){ errors.Ni++; errors.sequence += 'a'; }
+ }
+ else if(q[i] == '-' && r[i] != '-'){
+ if(r[i] == 'A'){ errors.dA++; errors.sequence += 'd'; }
+ if(r[i] == 'T'){ errors.dT++; errors.sequence += 'd'; }
+ if(r[i] == 'G'){ errors.dG++; errors.sequence += 'd'; }
+ if(r[i] == 'C'){ errors.dC++; errors.sequence += 'd'; }
+ }
+ errors.total++;
}
- else if(q[i] == '-' && r[i] != '-'){
- if(r[i] == 'A'){ errors.dA++; errors.sequence += 'd'; }
- if(r[i] == 'T'){ errors.dT++; errors.sequence += 'd'; }
- if(r[i] == 'G'){ errors.dG++; errors.sequence += 'd'; }
- if(r[i] == 'C'){ errors.dC++; errors.sequence += 'd'; }
+ else{
+
+ if(q[i] == '-'){
+ errors.sequence += 'd'; errors.total++;
+ }
+ else{
+ errors.sequence += 'r';
+ }
}
- errors.total++;
-
}
+
else if(q[i] == '.' && r[i] != '.'){ // reference extends beyond query
if(started == 1){ break; }
}
else if(q[i] == '.' && r[i] == '.'){ // both are missing data
if(started == 1){ break; }
}
-
+// cout << errors.sequence[errors.sequence.length()-1] << endl;
}
+// cout << errors.sequence << endl;
errors.mismatches = errors.total-errors.matches;
errors.errorRate = (double)(errors.total-errors.matches) / (double)errors.total;
errors.queryName = query.getName();
projects / mothur.git / blobdiff