X-Git-Url: https://git.donarmstrong.com/?p=mothur.git;a=blobdiff_plain;f=seqerrorcommand.cpp;h=41ddf35a7705264445749caf8868d80049ae61ac;hp=4ff1b6f07bbd782c40cfcd69d1032d16ed420956;hb=b206f634aae1b4ce13978d203247fb64757d5482;hpb=4458418562cc9dfc9a29ed4f8f6cfc7bfb927d40 diff --git a/seqerrorcommand.cpp b/seqerrorcommand.cpp index 4ff1b6f..41ddf35 100644 --- a/seqerrorcommand.cpp +++ b/seqerrorcommand.cpp @@ -11,20 +11,23 @@ #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","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 pname("name", "InputTypes", "", "", "namecount", "none", "none","",false,false,true); parameters.push_back(pname); + CommandParameter pcount("count", "InputTypes", "", "", "namecount", "none", "none","",false,false,true); parameters.push_back(pcount); 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 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); @@ -39,7 +42,9 @@ vector SeqErrorCommand::setParameters(){ exit(1); } } + //********************************************************************************************************************** + string SeqErrorCommand::getHelpString(){ try { string helpString = ""; @@ -48,7 +53,8 @@ string SeqErrorCommand::getHelpString(){ helpString += "The reference parameter...\n"; helpString += "The qfile parameter...\n"; helpString += "The report parameter...\n"; - helpString += "The name parameter...\n"; + helpString += "The name parameter allows you to provide a name file associated with the fasta file.\n"; + helpString += "The count parameter allows you to provide a count file associated with the fasta file.\n"; helpString += "The ignorechimeras parameter...\n"; helpString += "The threshold parameter...\n"; helpString += "The processors parameter...\n"; @@ -63,7 +69,9 @@ string SeqErrorCommand::getHelpString(){ exit(1); } } + //********************************************************************************************************************** + string SeqErrorCommand::getOutputPattern(string type) { try { string pattern = ""; @@ -88,7 +96,9 @@ string SeqErrorCommand::getOutputPattern(string type) { exit(1); } } + //********************************************************************************************************************** + SeqErrorCommand::SeqErrorCommand(){ try { abort = true; calledHelp = true; @@ -111,6 +121,7 @@ SeqErrorCommand::SeqErrorCommand(){ exit(1); } } + //*************************************************************************************************************** SeqErrorCommand::SeqErrorCommand(string option) { @@ -181,6 +192,14 @@ SeqErrorCommand::SeqErrorCommand(string option) { //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 names 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; } + } it = parameters.find("qfile"); //user has given a quality score file @@ -218,6 +237,12 @@ SeqErrorCommand::SeqErrorCommand(string option) { if(namesFileName == "not found"){ namesFileName = ""; } else if (namesFileName == "not open") { namesFileName = ""; abort = true; } else { m->setNameFile(namesFileName); } + + //check for optional parameters + countfile = validParameter.validFile(parameters, "count", true); + if(countfile == "not found"){ countfile = ""; } + else if (countfile == "not open") { countfile = ""; abort = true; } + else { m->setCountTableFile(countfile); } qualFileName = validParameter.validFile(parameters, "qfile", true); if(qualFileName == "not found"){ qualFileName = ""; } @@ -228,29 +253,19 @@ 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 = ""; outputDir += m->hasPath(queryFileName); //if user entered a file with a path then preserve it } + if ((countfile != "") && (namesFileName != "")) { m->mothurOut("You must enter ONLY ONE of the following: count or name."); m->mothurOutEndLine(); abort = true; } + //check for optional parameter and set defaults // ...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, "aligned", true); if (temp == "not found"){ temp = "t"; } - aligned = m->isTrue(temp); -// rdb->aligned = aligned; #do we need these lines for aligned? -// if (aligned) { //clear out old references -// rdb->clearMemory(); -// } temp = validParameter.validFile(parameters, "save", false); if (temp == "not found"){ temp = "f"; } save = m->isTrue(temp); @@ -277,6 +292,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); @@ -288,6 +306,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) { @@ -295,6 +329,7 @@ SeqErrorCommand::SeqErrorCommand(string option) { exit(1); } } + //*************************************************************************************************************** int SeqErrorCommand::execute(){ @@ -302,7 +337,7 @@ int SeqErrorCommand::execute(){ if (abort == true) { if (calledHelp) { return 0; } return 2; } int start = time(NULL); - maxLength = 2000; + maxLength = 5000; totalBases = 0; totalMatches = 0; @@ -320,7 +355,12 @@ int SeqErrorCommand::execute(){ getReferences(); //read in reference sequences - make sure there's no ambiguous bases - if(namesFileName != ""){ weights = getWeights(); } + if(namesFileName != "") { weights = getWeights(); } + else if (countfile != "") { + CountTable ct; + ct.readTable(countfile, false, false); + weights = ct.getNameMap(); + } vector fastaFilePos; vector qFilePos; @@ -333,22 +373,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); } @@ -370,23 +411,32 @@ 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 = getOutputFileName("errorref-query",variables); ofstream megAlignmentFile; m->openOutputFile(megAlignmentFileName, megAlignmentFile); 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(); + //set fasta file as new current fastafile + string current = ""; + itTypes = outputTypes.find("errorseq"); + if (itTypes != outputTypes.end()) { + if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setFastaFile(current); } + } + m->mothurOutEndLine(); m->mothurOut("Output File Names: "); m->mothurOutEndLine(); for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); } @@ -399,7 +449,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; @@ -614,7 +666,7 @@ int SeqErrorCommand::createProcesses(string filename, string qFileName, string r int misMatchSize; in >> misMatchSize; m->gobble(in); if (misMatchSize > misMatchCounts.size()) { misMatchCounts.resize(misMatchSize, 0); } - for (int j = 0; j < misMatchCounts.size(); j++) { + for (int j = 0; j < misMatchSize; j++) { in >> tempNum; misMatchCounts[j] += tempNum; } m->gobble(in); @@ -637,7 +689,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 { @@ -649,10 +703,10 @@ int SeqErrorCommand::driver(string filename, string qFileName, string rFileName, int numSeqs = 0; map::iterator it; - qScoreErrorMap['m'].assign(41, 0); - qScoreErrorMap['s'].assign(41, 0); - qScoreErrorMap['i'].assign(41, 0); - qScoreErrorMap['a'].assign(41, 0); + qScoreErrorMap['m'].assign(101, 0); + qScoreErrorMap['s'].assign(101, 0); + qScoreErrorMap['i'].assign(101, 0); + qScoreErrorMap['a'].assign(101, 0); errorForward['m'].assign(maxLength,0); errorForward['s'].assign(maxLength,0); @@ -669,11 +723,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); @@ -687,17 +742,32 @@ int SeqErrorCommand::driver(string filename, string qFileName, string rFileName, qualForwardMap.resize(maxLength); qualReverseMap.resize(maxLength); for(int i=0;iopenInputFile(qFileName, qualFile); + qualFile.seekg(qline.start); + + qualForwardMap.resize(maxLength); + qualReverseMap.resize(maxLength); + for(int i=0;iopenOutputFile(chimeraOutputFileName, outChimeraReport); - RefChimeraTest chimeraTest(referenceSeqs, aligned); - 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); @@ -706,48 +776,58 @@ 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 != ""){ + if((namesFileName != "") || (countfile != "")){ it = weights.find(query.getName()); minCompare.weight = it->second; } else{ minCompare.weight = 1; } + printErrorData(minCompare, numParentSeqs, errorSummaryFile, errorSeqFile); if(!ignoreSeq){ - for(int i=0;i maxMismatch){ @@ -789,15 +881,18 @@ int SeqErrorCommand::driver(string filename, string qFileName, string rFileName, if (queryFile.eof()) { break; } #endif - if(index % 100 == 0){ m->mothurOut(toString(index)); m->mothurOutEndLine(); } + if(index % 100 == 0){ m->mothurOutJustToScreen(toString(index)+"\n"); } } 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)); m->mothurOutEndLine(); } + m->mothurOutJustToScreen(toString(index)+"\n"); return index; } @@ -806,6 +901,7 @@ int SeqErrorCommand::driver(string filename, string qFileName, string rFileName, exit(1); } } + //*************************************************************************************************************** void SeqErrorCommand::getReferences(){ @@ -901,9 +997,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){ @@ -1139,6 +1238,7 @@ void SeqErrorCommand::printSubMatrix(){ exit(1); } } + //*************************************************************************************************************** void SeqErrorCommand::printErrorFRFile(map > errorForward, map > errorReverse){ @@ -1201,18 +1301,17 @@ void SeqErrorCommand::printErrorQuality(map > qScoreErrorMap){ outputNames.push_back(errorQualityFileName); outputTypes["errorquality"].push_back(errorQualityFileName); errorQualityFile << "qscore\tmatches\tsubstitutions\tinsertions\tambiguous" << endl; - for(int i=0;i<41;i++){ + for(int i=0;i<101;i++){ errorQualityFile << i << '\t' << qScoreErrorMap['m'][i] << '\t' << qScoreErrorMap['s'][i] << '\t' << qScoreErrorMap['i'][i] << '\t'<< qScoreErrorMap['a'][i] << endl; } errorQualityFile.close(); } catch(exception& e) { - m->errorOut(e, "SeqErrorCommand", "printErrorFRFile"); + m->errorOut(e, "SeqErrorCommand", "printErrorQuality"); exit(1); } } - //*************************************************************************************************************** void SeqErrorCommand::printQualityFR(vector > qualForwardMap, vector > qualReverseMap){ @@ -1266,11 +1365,12 @@ void SeqErrorCommand::printQualityFR(vector > 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) { @@ -1311,6 +1411,8 @@ int SeqErrorCommand::setLines(string filename, string qfilename, string rfilenam string sname = ""; nameStream >> sname; sname = sname.substr(1); + + m->checkName(sname); map::iterator it = firstSeqNames.find(sname); @@ -1349,65 +1451,67 @@ 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; + + m->checkName(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); + } + + rfileFilePos.push_back(sizeR); } - - rfileFilePos.push_back(sizeR); - return processors; #else @@ -1449,4 +1553,5 @@ int SeqErrorCommand::setLines(string filename, string qfilename, string rfilenam exit(1); } } + //***************************************************************************************************************