X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=seqerrorcommand.cpp;h=0a6eae93df72e8556554bef32c5a52074a4d0a54;hb=01d5a60fc5f396339acf529151fa8186ee7c1a41;hp=0dc15f59ee98d0ee2d93659c4c3ad1f033f8a193;hpb=a33a385cc5b7481488f92f794425f01fbf40a543;p=mothur.git diff --git a/seqerrorcommand.cpp b/seqerrorcommand.cpp index 0dc15f5..0a6eae9 100644 --- a/seqerrorcommand.cpp +++ b/seqerrorcommand.cpp @@ -11,23 +11,26 @@ #include "reportfile.h" #include "qualityscores.h" #include "refchimeratest.h" +#include "myPerseus.h" #include "filterseqscommand.h" //********************************************************************************************************************** + vector SeqErrorCommand::setParameters(){ try { - CommandParameter pquery("fasta", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(pquery); - CommandParameter preference("reference", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(preference); - CommandParameter pqfile("qfile", "InputTypes", "", "", "none", "none", "QualReport",false,false); parameters.push_back(pqfile); - CommandParameter preport("report", "InputTypes", "", "", "none", "none", "QualReport",false,false); parameters.push_back(preport); - CommandParameter pname("name", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(pname); - 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 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); + CommandParameter pquery("fasta", "InputTypes", "", "", "none", "none", "none","errorType",false,true,true); parameters.push_back(pquery); + CommandParameter preference("reference", "InputTypes", "", "", "none", "none", "none","",false,true,true); parameters.push_back(preference); + CommandParameter pqfile("qfile", "InputTypes", "", "", "none", "none", "QualReport","",false,false); parameters.push_back(pqfile); + CommandParameter preport("report", "InputTypes", "", "", "none", "none", "QualReport","",false,false); parameters.push_back(preport); + CommandParameter pname("name", "InputTypes", "", "", "none", "none", "none","",false,false); parameters.push_back(pname); + CommandParameter pignorechimeras("ignorechimeras", "Boolean", "", "T", "", "", "","",false,false); parameters.push_back(pignorechimeras); + CommandParameter pthreshold("threshold", "Number", "", "1.0", "", "", "","",false,false); parameters.push_back(pthreshold); + CommandParameter paligned("aligned", "Boolean", "", "T", "", "", "","",false,false); parameters.push_back(paligned); + CommandParameter pprocessors("processors", "Number", "", "1", "", "", "","",false,false,true); parameters.push_back(pprocessors); + 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); vector myArray; for (int i = 0; i < parameters.size(); i++) { myArray.push_back(parameters[i].name); } @@ -38,7 +41,9 @@ vector SeqErrorCommand::setParameters(){ exit(1); } } + //********************************************************************************************************************** + string SeqErrorCommand::getHelpString(){ try { string helpString = ""; @@ -62,27 +67,59 @@ string SeqErrorCommand::getHelpString(){ exit(1); } } + //********************************************************************************************************************** + +string SeqErrorCommand::getOutputPattern(string type) { + try { + string pattern = ""; + + if (type == "errorsummary") { pattern = "[filename],error.summary"; } + else if (type == "errorseq") { pattern = "[filename],error.seq"; } + else if (type == "errorquality") { pattern = "[filename],error.quality"; } + else if (type == "errorqualforward") { pattern = "[filename],error.qual.forward"; } + else if (type == "errorqualreverse") { pattern = "[filename],error.qual.reverse"; } + else if (type == "errorforward") { pattern = "[filename],error.seq.forward"; } + else if (type == "errorreverse") { pattern = "[filename],error.seq.reverse"; } + else if (type == "errorcount") { pattern = "[filename],error.count"; } + else if (type == "errormatrix") { pattern = "[filename],error.matrix"; } + else if (type == "errorchimera") { pattern = "[filename],error.chimera"; } + else if (type == "errorref-query") { pattern = "[filename],error.ref-query"; } + else { m->mothurOut("[ERROR]: No definition for type " + type + " output pattern.\n"); m->control_pressed = true; } + + return pattern; + } + catch(exception& e) { + m->errorOut(e, "SeqErrorCommand", "getOutputPattern"); + exit(1); + } +} + +//********************************************************************************************************************** + SeqErrorCommand::SeqErrorCommand(){ try { abort = true; calledHelp = true; setParameters(); vector tempOutNames; - outputTypes["error.summary"] = tempOutNames; - outputTypes["error.seq"] = tempOutNames; - outputTypes["error.quality"] = tempOutNames; - outputTypes["error.qual.forward"] = tempOutNames; - outputTypes["error.qual.reverse"] = tempOutNames; - outputTypes["error.forward"] = tempOutNames; - outputTypes["error.reverse"] = tempOutNames; - outputTypes["error.count"] = tempOutNames; - outputTypes["error.matrix"] = tempOutNames; + outputTypes["errorsummary"] = tempOutNames; + outputTypes["errorseq"] = tempOutNames; + outputTypes["errorquality"] = tempOutNames; + outputTypes["errorqualforward"] = tempOutNames; + outputTypes["errorqualreverse"] = tempOutNames; + outputTypes["errorforward"] = tempOutNames; + outputTypes["errorreverse"] = tempOutNames; + outputTypes["errorcount"] = tempOutNames; + outputTypes["errormatrix"] = tempOutNames; + outputTypes["errorchimera"] = tempOutNames; + outputTypes["errorref-query"] = tempOutNames; } catch(exception& e) { m->errorOut(e, "SeqErrorCommand", "SeqErrorCommand"); exit(1); } } + //*************************************************************************************************************** SeqErrorCommand::SeqErrorCommand(string option) { @@ -112,15 +149,17 @@ SeqErrorCommand::SeqErrorCommand(string option) { //initialize outputTypes vector tempOutNames; - outputTypes["error.summary"] = tempOutNames; - outputTypes["error.seq"] = tempOutNames; - outputTypes["error.quality"] = tempOutNames; - outputTypes["error.qual.forward"] = tempOutNames; - outputTypes["error.qual.reverse"] = tempOutNames; - outputTypes["error.forward"] = tempOutNames; - outputTypes["error.reverse"] = tempOutNames; - outputTypes["error.count"] = tempOutNames; - outputTypes["error.matrix"] = tempOutNames; + outputTypes["errorsummary"] = tempOutNames; + outputTypes["errorseq"] = tempOutNames; + outputTypes["errorquality"] = tempOutNames; + outputTypes["errorqualforward"] = tempOutNames; + outputTypes["errorqualreverse"] = tempOutNames; + outputTypes["errorforward"] = tempOutNames; + outputTypes["errorreverse"] = tempOutNames; + outputTypes["errorcount"] = tempOutNames; + outputTypes["errormatrix"] = tempOutNames; + outputTypes["errorchimera"] = tempOutNames; + outputTypes["errorref-query"] = tempOutNames; //if the user changes the input directory command factory will send this info to us in the output parameter @@ -198,12 +237,6 @@ SeqErrorCommand::SeqErrorCommand(string option) { if(reportFileName == "not found"){ reportFileName = ""; } else if (reportFileName == "not open") { reportFileName = ""; abort = true; } - if((reportFileName != "" && qualFileName == "") || (reportFileName == "" && qualFileName != "")){ - m->mothurOut("if you use either a qual file or a report file, you have to have both."); - m->mothurOutEndLine(); - abort = true; - } - outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){ outputDir = ""; @@ -214,6 +247,7 @@ SeqErrorCommand::SeqErrorCommand(string option) { // ...at some point should added some additional type checking... temp = validParameter.validFile(parameters, "threshold", false); if (temp == "not found") { temp = "1.00"; } m->mothurConvert(temp, threshold); + temp = validParameter.validFile(parameters, "save", false); if (temp == "not found"){ temp = "f"; } save = m->isTrue(temp); @@ -240,6 +274,9 @@ SeqErrorCommand::SeqErrorCommand(string option) { temp = validParameter.validFile(parameters, "ignorechimeras", false); if (temp == "not found") { temp = "T"; } ignoreChimeras = m->isTrue(temp); + temp = validParameter.validFile(parameters, "aligned", false); if (temp == "not found"){ temp = "t"; } + aligned = m->isTrue(temp); + temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); } m->setProcessors(temp); m->mothurConvert(temp, processors); @@ -251,6 +288,22 @@ SeqErrorCommand::SeqErrorCommand(string option) { vector files; files.push_back(queryFileName); parser.getNameFile(files); } + + if(aligned == true){ + if((reportFileName != "" && qualFileName == "") || (reportFileName == "" && qualFileName != "")){ + m->mothurOut("if you use either a qual file or a report file, you have to have both."); + m->mothurOutEndLine(); + abort = true; + } + } + else{ + if(reportFileName != ""){ + m->mothurOut("we are ignoring the report file if your sequences are not aligned. we will check that the sequences in your fasta and and qual file are the same length."); + m->mothurOutEndLine(); + } + } + + } } catch(exception& e) { @@ -258,6 +311,7 @@ SeqErrorCommand::SeqErrorCommand(string option) { exit(1); } } + //*************************************************************************************************************** int SeqErrorCommand::execute(){ @@ -265,18 +319,21 @@ int SeqErrorCommand::execute(){ if (abort == true) { if (calledHelp) { return 0; } return 2; } int start = time(NULL); - maxLength = 2000; + maxLength = 5000; totalBases = 0; totalMatches = 0; - string errorSummaryFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.summary"; - outputNames.push_back(errorSummaryFileName); outputTypes["error.summary"].push_back(errorSummaryFileName); + string fileNameRoot = outputDir + m->getRootName(m->getSimpleName(queryFileName)); + map variables; + variables["[filename]"] = fileNameRoot; + string errorSummaryFileName = getOutputFileName("errorsummary",variables); + outputNames.push_back(errorSummaryFileName); outputTypes["errorsummary"].push_back(errorSummaryFileName); - string errorSeqFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.seq"; - outputNames.push_back(errorSeqFileName); outputTypes["error.seq"].push_back(errorSeqFileName); + string errorSeqFileName = getOutputFileName("errorseq",variables); + outputNames.push_back(errorSeqFileName); outputTypes["errorseq"].push_back(errorSeqFileName); - string errorChimeraFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.chimera"; - outputNames.push_back(errorChimeraFileName); outputTypes["error.chimera"].push_back(errorChimeraFileName); + string errorChimeraFileName = getOutputFileName("errorchimera",variables); + outputNames.push_back(errorChimeraFileName); outputTypes["errorchimera"].push_back(errorChimeraFileName); getReferences(); //read in reference sequences - make sure there's no ambiguous bases @@ -293,22 +350,23 @@ int SeqErrorCommand::execute(){ for (int i = 0; i < (fastaFilePos.size()-1); i++) { lines.push_back(linePair(fastaFilePos[i], fastaFilePos[(i+1)])); if (qualFileName != "") { qLines.push_back(linePair(qFilePos[i], qFilePos[(i+1)])); } - if (reportFileName != "") { rLines.push_back(linePair(reportFilePos[i], reportFilePos[(i+1)])); } + if (reportFileName != "" && aligned == true) { rLines.push_back(linePair(reportFilePos[i], reportFilePos[(i+1)])); } } if(qualFileName == "") { qLines = lines; rLines = lines; } //fills with duds - + if(aligned == false){ rLines = lines; } int numSeqs = 0; #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix) if(processors == 1){ numSeqs = driver(queryFileName, qualFileName, reportFileName, errorSummaryFileName, errorSeqFileName, errorChimeraFileName, lines[0], qLines[0], rLines[0]); + }else{ numSeqs = createProcesses(queryFileName, qualFileName, reportFileName, errorSummaryFileName, errorSeqFileName, errorChimeraFileName); } #else numSeqs = driver(queryFileName, qualFileName, reportFileName, errorSummaryFileName, errorSeqFileName, errorChimeraFileName, lines[0], qLines[0], rLines[0]); #endif - - if(qualFileName != "" && reportFileName != ""){ + + if(qualFileName != ""){ printErrorQuality(qScoreErrorMap); printQualityFR(qualForwardMap, qualReverseMap); } @@ -317,10 +375,10 @@ int SeqErrorCommand::execute(){ 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"; + string errorCountFileName = getOutputFileName("errorcount",variables); ofstream errorCountFile; m->openOutputFile(errorCountFileName, errorCountFile); - outputNames.push_back(errorCountFileName); outputTypes["error.count"].push_back(errorCountFileName); + outputNames.push_back(errorCountFileName); outputTypes["errorcount"].push_back(errorCountFileName); m->mothurOut("Overall error rate:\t" + toString((double)(totalBases - totalMatches) / (double)totalBases) + "\n"); m->mothurOut("Errors\tSequences\n"); errorCountFile << "Errors\tSequences\n"; @@ -330,19 +388,21 @@ int SeqErrorCommand::execute(){ } errorCountFile.close(); - if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; } +// if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; } printSubMatrix(); - - string megAlignmentFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.ref-query"; + + string megAlignmentFileName = getOutputFileName("errorref-query",variables); ofstream megAlignmentFile; m->openOutputFile(megAlignmentFileName, megAlignmentFile); - outputNames.push_back(megAlignmentFileName); outputTypes["error.ref-query"].push_back(megAlignmentFileName); - + outputNames.push_back(megAlignmentFileName); outputTypes["errorref-query"].push_back(megAlignmentFileName); + for(int i=0;imothurOut("It took " + toString(time(NULL) - start) + " secs to check " + toString(numSeqs) + " sequences."); m->mothurOutEndLine(); @@ -359,7 +419,9 @@ int SeqErrorCommand::execute(){ exit(1); } } + //********************************************************************************************************************** + int SeqErrorCommand::createProcesses(string filename, string qFileName, string rFileName, string summaryFileName, string errorOutputFileName, string chimeraOutputFileName) { try { int process = 1; @@ -597,7 +659,9 @@ int SeqErrorCommand::createProcesses(string filename, string qFileName, string r exit(1); } } + //********************************************************************************************************************** + int SeqErrorCommand::driver(string filename, string qFileName, string rFileName, string summaryFileName, string errorOutputFileName, string chimeraOutputFileName, linePair line, linePair qline, linePair rline) { try { @@ -629,11 +693,12 @@ int SeqErrorCommand::driver(string filename, string qFileName, string rFileName, //open inputfiles and go to beginning place for this processor ifstream queryFile; m->openInputFile(filename, queryFile); + queryFile.seekg(line.start); ifstream reportFile; ifstream qualFile; - if(qFileName != "" && rFileName != ""){ + if((qFileName != "" && rFileName != "" && aligned)){ m->openInputFile(qFileName, qualFile); qualFile.seekg(qline.start); @@ -651,11 +716,28 @@ int SeqErrorCommand::driver(string filename, string qFileName, string rFileName, qualReverseMap[i].assign(41,0); } } - + else if(qFileName != "" && !aligned){ + + m->openInputFile(qFileName, qualFile); + qualFile.seekg(qline.start); + + qualForwardMap.resize(maxLength); + qualReverseMap.resize(maxLength); + for(int i=0;iopenOutputFile(chimeraOutputFileName, outChimeraReport); - RefChimeraTest chimeraTest(referenceSeqs); - if (line.start == 0) { chimeraTest.printHeader(outChimeraReport); } + + + RefChimeraTest chimeraTest; + + chimeraTest = RefChimeraTest(referenceSeqs, aligned); + if (line.start == 0) { chimeraTest.printHeader(outChimeraReport); } + ofstream errorSummaryFile; m->openOutputFile(summaryFileName, errorSummaryFile); @@ -664,26 +746,37 @@ int SeqErrorCommand::driver(string filename, string qFileName, string rFileName, ofstream errorSeqFile; m->openOutputFile(errorOutputFileName, errorSeqFile); - megaAlignVector.resize(numRefs, ""); + megaAlignVector.assign(numRefs, ""); int index = 0; bool ignoreSeq = 0; bool moreSeqs = 1; while (moreSeqs) { - - if (m->control_pressed) { queryFile.close(); if(qFileName != "" && rFileName != ""){ reportFile.close(); qualFile.close(); } outChimeraReport.close(); errorSummaryFile.close();errorSeqFile.close(); return 0; } - + Sequence query(queryFile); - - int numParentSeqs = chimeraTest.analyzeQuery(query.getName(), query.getAligned(), outChimeraReport); - int closestRefIndex = chimeraTest.getClosestRefIndex(); - + Sequence reference; + int numParentSeqs = -1; + int closestRefIndex = -1; + + string querySeq = query.getAligned(); + if (!aligned) { querySeq = query.getUnaligned(); } + + numParentSeqs = chimeraTest.analyzeQuery(query.getName(), querySeq, outChimeraReport); + + closestRefIndex = chimeraTest.getClosestRefIndex(); + + reference = referenceSeqs[closestRefIndex]; + + reference.setAligned(chimeraTest.getClosestRefAlignment()); + query.setAligned(chimeraTest.getQueryAlignment()); + if(numParentSeqs > 1 && ignoreChimeras == 1) { ignoreSeq = 1; } - else { ignoreSeq = 0; } - + else { ignoreSeq = 0; } + Compare minCompare; - getErrors(query, referenceSeqs[closestRefIndex], minCompare); + + getErrors(query, reference, minCompare); if(namesFileName != ""){ it = weights.find(query.getName()); @@ -691,21 +784,20 @@ int SeqErrorCommand::driver(string filename, string qFileName, string rFileName, } else{ minCompare.weight = 1; } + printErrorData(minCompare, numParentSeqs, errorSummaryFile, errorSeqFile); if(!ignoreSeq){ - for(int i=0;i maxMismatch){ @@ -747,15 +851,18 @@ int SeqErrorCommand::driver(string filename, string qFileName, string rFileName, if (queryFile.eof()) { break; } #endif - if(index % 100 == 0){ m->mothurOut(toString(index) + '\n'); } + if(index % 100 == 0){ m->mothurOut(toString(index)); m->mothurOutEndLine(); } } queryFile.close(); - if(qFileName != "" && rFileName != ""){ reportFile.close(); qualFile.close(); } - errorSummaryFile.close(); + outChimeraReport.close(); + errorSummaryFile.close(); errorSeqFile.close(); - + + if(qFileName != "" && rFileName != "") { reportFile.close(); qualFile.close(); } + else if(qFileName != "" && aligned == false){ qualFile.close(); } + //report progress - if(index % 100 != 0){ m->mothurOut(toString(index) + '\n'); } + m->mothurOut(toString(index)); m->mothurOutEndLine(); return index; } @@ -764,6 +871,7 @@ int SeqErrorCommand::driver(string filename, string qFileName, string rFileName, exit(1); } } + //*************************************************************************************************************** void SeqErrorCommand::getReferences(){ @@ -786,12 +894,16 @@ void SeqErrorCommand::getReferences(){ // // int endPos = rdb->referenceSeqs[i].getEndPos(); // if(endPos < minEndPos) { minEndPos = endPos; } + if (rdb->referenceSeqs[i].getNumBases() == 0) { + m->mothurOut("[WARNING]: " + rdb->referenceSeqs[i].getName() + " is blank, ignoring.");m->mothurOutEndLine(); + }else { + referenceSeqs.push_back(rdb->referenceSeqs[i]); + } - 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(); + m->mothurOut("It took " + toString(time(NULL) - start) + " to load " + toString(referenceSeqs.size()) + " sequences.");m->mothurOutEndLine(); }else { int start = time(NULL); @@ -809,9 +921,12 @@ void SeqErrorCommand::getReferences(){ // // int endPos = currentSeq.getEndPos(); // if(endPos < minEndPos) { minEndPos = endPos; } - referenceSeqs.push_back(currentSeq); - - if (rdb->save) { rdb->referenceSeqs.push_back(currentSeq); } + if (currentSeq.getNumBases() == 0) { + m->mothurOut("[WARNING]: " + currentSeq.getName() + " is blank, ignoring.");m->mothurOutEndLine(); + }else { + referenceSeqs.push_back(currentSeq); + if (rdb->save) { rdb->referenceSeqs.push_back(currentSeq); } + } m->gobble(referenceFile); } @@ -825,7 +940,7 @@ void SeqErrorCommand::getReferences(){ for(int i=0;imothurOut("Warning: " + toString(numAmbigSeqs) + " reference sequences have ambiguous bases, these bases will be ignored\n"); @@ -852,9 +967,10 @@ int SeqErrorCommand::getErrors(Sequence query, Sequence reference, Compare& erro int started = 0; //Compare errors; - + + errors.sequence = ""; for(int i=0;i SeqErrorCommand::getWeights(){ nameCountMap[seqName] = m->getNumNames(redundantSeqs); m->gobble(nameFile); } + + nameFile.close(); + return nameCountMap; } - //*************************************************************************************************************** void SeqErrorCommand::printErrorHeader(ofstream& errorSummaryFile){ @@ -1046,10 +1163,13 @@ void SeqErrorCommand::printErrorData(Compare error, int numParentSeqs, ofstream& void SeqErrorCommand::printSubMatrix(){ try { - string subMatrixFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.matrix"; + string fileNameRoot = outputDir + m->getRootName(m->getSimpleName(queryFileName)); + map variables; + variables["[filename]"] = fileNameRoot; + string subMatrixFileName = getOutputFileName("errormatrix",variables); ofstream subMatrixFile; m->openOutputFile(subMatrixFileName, subMatrixFile); - outputNames.push_back(subMatrixFileName); outputTypes["error.matrix"].push_back(subMatrixFileName); + outputNames.push_back(subMatrixFileName); outputTypes["errormatrix"].push_back(subMatrixFileName); vector bases(6); bases[0] = "A"; bases[1] = "T"; @@ -1088,14 +1208,18 @@ void SeqErrorCommand::printSubMatrix(){ exit(1); } } + //*************************************************************************************************************** void SeqErrorCommand::printErrorFRFile(map > errorForward, map > errorReverse){ try{ - string errorForwardFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.seq.forward"; + string fileNameRoot = outputDir + m->getRootName(m->getSimpleName(queryFileName)); + map variables; + variables["[filename]"] = fileNameRoot; + string errorForwardFileName = getOutputFileName("errorforward",variables); ofstream errorForwardFile; m->openOutputFile(errorForwardFileName, errorForwardFile); - outputNames.push_back(errorForwardFileName); outputTypes["error.forward"].push_back(errorForwardFileName); + outputNames.push_back(errorForwardFileName); outputTypes["errorforward"].push_back(errorForwardFileName); errorForwardFile << "position\ttotalseqs\tmatch\tsubstitution\tinsertion\tdeletion\tambiguous" << endl; for(int i=0;i > errorForward, map } errorForwardFile.close(); - string errorReverseFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.seq.reverse"; + string errorReverseFileName = getOutputFileName("errorreverse",variables); ofstream errorReverseFile; m->openOutputFile(errorReverseFileName, errorReverseFile); - outputNames.push_back(errorReverseFileName); outputTypes["error.reverse"].push_back(errorReverseFileName); + outputNames.push_back(errorReverseFileName); outputTypes["errorreverse"].push_back(errorReverseFileName); errorReverseFile << "position\ttotalseqs\tmatch\tsubstitution\tinsertion\tdeletion\tambiguous" << endl; for(int i=0;i > errorForward, map void SeqErrorCommand::printErrorQuality(map > qScoreErrorMap){ try{ - - string errorQualityFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.quality"; + string fileNameRoot = outputDir + m->getRootName(m->getSimpleName(queryFileName)); + map variables; + variables["[filename]"] = fileNameRoot; + string errorQualityFileName = getOutputFileName("errorquality",variables); ofstream errorQualityFile; m->openOutputFile(errorQualityFileName, errorQualityFile); - outputNames.push_back(errorQualityFileName); outputTypes["error.quality"].push_back(errorQualityFileName); + outputNames.push_back(errorQualityFileName); outputTypes["errorquality"].push_back(errorQualityFileName); errorQualityFile << "qscore\tmatches\tsubstitutions\tinsertions\tambiguous" << endl; for(int i=0;i<41;i++){ @@ -1151,12 +1277,11 @@ void SeqErrorCommand::printErrorQuality(map > qScoreErrorMap){ errorQualityFile.close(); } catch(exception& e) { - m->errorOut(e, "SeqErrorCommand", "printErrorFRFile"); + m->errorOut(e, "SeqErrorCommand", "printErrorQuality"); exit(1); } } - //*************************************************************************************************************** void SeqErrorCommand::printQualityFR(vector > qualForwardMap, vector > qualReverseMap){ @@ -1172,11 +1297,13 @@ void SeqErrorCommand::printQualityFR(vector > qualForwardMap, vector } } } - - string qualityForwardFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.qual.forward"; + string fileNameRoot = outputDir + m->getRootName(m->getSimpleName(queryFileName)); + map variables; + variables["[filename]"] = fileNameRoot; + string qualityForwardFileName = getOutputFileName("errorqualforward",variables); ofstream qualityForwardFile; m->openOutputFile(qualityForwardFileName, qualityForwardFile); - outputNames.push_back(qualityForwardFileName); outputTypes["error.qual.forward"].push_back(qualityForwardFileName); + outputNames.push_back(qualityForwardFileName); outputTypes["errorqualforward"].push_back(qualityForwardFileName); for(int i=0;i > qualForwardMap, vector qualityForwardFile.close(); - string qualityReverseFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.qual.reverse"; + string qualityReverseFileName = getOutputFileName("errorqualreverse",variables); ofstream qualityReverseFile; m->openOutputFile(qualityReverseFileName, qualityReverseFile); - outputNames.push_back(qualityReverseFileName); outputTypes["error.qual.reverse"].push_back(qualityReverseFileName); + outputNames.push_back(qualityReverseFileName); outputTypes["errorqualreverse"].push_back(qualityReverseFileName); for(int i=0;i > qualForwardMap, vector qualityReverseFile.close(); } catch(exception& e) { - m->errorOut(e, "SeqErrorCommand", "printErrorFRFile"); + m->errorOut(e, "SeqErrorCommand", "printQualityFR"); exit(1); } } + /**************************************************************************************************/ int SeqErrorCommand::setLines(string filename, string qfilename, string rfilename, vector& fastaFilePos, vector& qfileFilePos, vector& rfileFilePos) { @@ -1253,6 +1381,10 @@ int SeqErrorCommand::setLines(string filename, string qfilename, string rfilenam string sname = ""; nameStream >> sname; sname = sname.substr(1); + + for (int i = 0; i < sname.length(); i++) { + if (sname[i] == ':') { sname[i] = '_'; m->changedSeqNames = true; } + } map::iterator it = firstSeqNames.find(sname); @@ -1291,65 +1423,69 @@ int SeqErrorCommand::setLines(string filename, string qfilename, string rfilenam qfileFilePos.push_back(size); - //seach for filePos of each first name in the rfile and save in rfileFilePos - string junk; - ifstream inR; - m->openInputFile(rfilename, inR); - - //read column headers - for (int i = 0; i < 16; i++) { - if (!inR.eof()) { inR >> junk; } - else { break; } - } - - while(!inR.eof()){ - - if (m->control_pressed) { inR.close(); return processors; } - - input = m->getline(inR); - - if (input.length() != 0) { - - istringstream nameStream(input); - string sname = ""; nameStream >> sname; - - map::iterator it = firstSeqNamesReport.find(sname); - - if(it != firstSeqNamesReport.end()) { //this is the start of a new chunk - unsigned long long pos = inR.tellg(); - rfileFilePos.push_back(pos - input.length() - 1); - firstSeqNamesReport.erase(it); - } - } - - if (firstSeqNamesReport.size() == 0) { break; } - m->gobble(inR); - } - inR.close(); - - if (firstSeqNamesReport.size() != 0) { - for (map::iterator it = firstSeqNamesReport.begin(); it != firstSeqNamesReport.end(); it++) { - m->mothurOut(it->first + " is in your fasta file and not in your report file, aborting."); m->mothurOutEndLine(); - } - m->control_pressed = true; - return processors; - } - - //get last file position of qfile - FILE * rFile; - unsigned long long sizeR; - - //get num bytes in file - rFile = fopen (rfilename.c_str(),"rb"); - if (rFile==NULL) perror ("Error opening file"); - else{ - fseek (rFile, 0, SEEK_END); - sizeR=ftell (rFile); - fclose (rFile); + if(aligned){ + //seach for filePos of each first name in the rfile and save in rfileFilePos + string junk; + ifstream inR; + + m->openInputFile(rfilename, inR); + + //read column headers + for (int i = 0; i < 16; i++) { + if (!inR.eof()) { inR >> junk; } + else { break; } + } + + while(!inR.eof()){ + + input = m->getline(inR); + + if (input.length() != 0) { + + istringstream nameStream(input); + string sname = ""; nameStream >> sname; + + for (int i = 0; i < sname.length(); i++) { + if (sname[i] == ':') { sname[i] = '_'; m->changedSeqNames = true; } + } + + map::iterator it = firstSeqNamesReport.find(sname); + + if(it != firstSeqNamesReport.end()) { //this is the start of a new chunk + unsigned long long pos = inR.tellg(); + rfileFilePos.push_back(pos - input.length() - 1); + firstSeqNamesReport.erase(it); + } + } + + if (firstSeqNamesReport.size() == 0) { break; } + m->gobble(inR); + } + inR.close(); + + if (firstSeqNamesReport.size() != 0) { + for (map::iterator it = firstSeqNamesReport.begin(); it != firstSeqNamesReport.end(); it++) { + m->mothurOut(it->first + " is in your fasta file and not in your report file, aborting."); m->mothurOutEndLine(); + } + m->control_pressed = true; + return processors; + } + + //get last file position of qfile + FILE * rFile; + unsigned long long sizeR; + + //get num bytes in file + rFile = fopen (rfilename.c_str(),"rb"); + if (rFile==NULL) perror ("Error opening file"); + else{ + fseek (rFile, 0, SEEK_END); + sizeR=ftell (rFile); + fclose (rFile); + } + + rfileFilePos.push_back(sizeR); } - - rfileFilePos.push_back(sizeR); - return processors; #else @@ -1391,4 +1527,5 @@ int SeqErrorCommand::setLines(string filename, string qfilename, string rfilenam exit(1); } } + //***************************************************************************************************************