5 * Created by Pat Schloss on 7/15/10.
6 * Copyright 2010 Schloss Lab. All rights reserved.
10 #include "seqerrorcommand.h"
11 #include "reportfile.h"
12 #include "qualityscores.h"
13 #include "refchimeratest.h"
14 #include "myPerseus.h"
15 #include "filterseqscommand.h"
18 //**********************************************************************************************************************
20 vector<string> SeqErrorCommand::setParameters(){
22 CommandParameter pquery("fasta", "InputTypes", "", "", "none", "none", "none","errorType",false,true,true); parameters.push_back(pquery);
23 CommandParameter preference("reference", "InputTypes", "", "", "none", "none", "none","",false,true,true); parameters.push_back(preference);
24 CommandParameter pqfile("qfile", "InputTypes", "", "", "none", "none", "QualReport","",false,false); parameters.push_back(pqfile);
25 CommandParameter preport("report", "InputTypes", "", "", "none", "none", "QualReport","",false,false); parameters.push_back(preport);
26 CommandParameter pname("name", "InputTypes", "", "", "namecount", "none", "none","",false,false,true); parameters.push_back(pname);
27 CommandParameter pcount("count", "InputTypes", "", "", "namecount", "none", "none","",false,false,true); parameters.push_back(pcount);
28 CommandParameter pignorechimeras("ignorechimeras", "Boolean", "", "T", "", "", "","",false,false); parameters.push_back(pignorechimeras);
29 CommandParameter pthreshold("threshold", "Number", "", "1.0", "", "", "","",false,false); parameters.push_back(pthreshold);
30 CommandParameter paligned("aligned", "Boolean", "", "T", "", "", "","",false,false); parameters.push_back(paligned);
31 CommandParameter pprocessors("processors", "Number", "", "1", "", "", "","",false,false,true); parameters.push_back(pprocessors);
32 CommandParameter psave("save", "Boolean", "", "F", "", "", "","",false,false); parameters.push_back(psave);
33 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(pinputdir);
34 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(poutputdir);
36 vector<string> myArray;
37 for (int i = 0; i < parameters.size(); i++) { myArray.push_back(parameters[i].name); }
41 m->errorOut(e, "SeqErrorCommand", "setParameters");
46 //**********************************************************************************************************************
48 string SeqErrorCommand::getHelpString(){
50 string helpString = "";
51 helpString += "The seq.error command reads a query alignment file and a reference alignment file and creates .....\n";
52 helpString += "The fasta parameter...\n";
53 helpString += "The reference parameter...\n";
54 helpString += "The qfile parameter...\n";
55 helpString += "The report parameter...\n";
56 helpString += "The name parameter allows you to provide a name file associated with the fasta file.\n";
57 helpString += "The count parameter allows you to provide a count file associated with the fasta file.\n";
58 helpString += "The ignorechimeras parameter...\n";
59 helpString += "The threshold parameter...\n";
60 helpString += "The processors parameter...\n";
61 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.";
62 helpString += "Example seq.error(...).\n";
63 helpString += "Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFasta).\n";
64 helpString += "For more details please check out the wiki http://www.mothur.org/wiki/seq.error .\n";
68 m->errorOut(e, "SeqErrorCommand", "getHelpString");
73 //**********************************************************************************************************************
75 string SeqErrorCommand::getOutputPattern(string type) {
79 if (type == "errorsummary") { pattern = "[filename],error.summary"; }
80 else if (type == "errorseq") { pattern = "[filename],error.seq"; }
81 else if (type == "errorquality") { pattern = "[filename],error.quality"; }
82 else if (type == "errorqualforward") { pattern = "[filename],error.qual.forward"; }
83 else if (type == "errorqualreverse") { pattern = "[filename],error.qual.reverse"; }
84 else if (type == "errorforward") { pattern = "[filename],error.seq.forward"; }
85 else if (type == "errorreverse") { pattern = "[filename],error.seq.reverse"; }
86 else if (type == "errorcount") { pattern = "[filename],error.count"; }
87 else if (type == "errormatrix") { pattern = "[filename],error.matrix"; }
88 else if (type == "errorchimera") { pattern = "[filename],error.chimera"; }
89 else if (type == "errorref-query") { pattern = "[filename],error.ref-query"; }
90 else { m->mothurOut("[ERROR]: No definition for type " + type + " output pattern.\n"); m->control_pressed = true; }
95 m->errorOut(e, "SeqErrorCommand", "getOutputPattern");
100 //**********************************************************************************************************************
102 SeqErrorCommand::SeqErrorCommand(){
104 abort = true; calledHelp = true;
106 vector<string> tempOutNames;
107 outputTypes["errorsummary"] = tempOutNames;
108 outputTypes["errorseq"] = tempOutNames;
109 outputTypes["errorquality"] = tempOutNames;
110 outputTypes["errorqualforward"] = tempOutNames;
111 outputTypes["errorqualreverse"] = tempOutNames;
112 outputTypes["errorforward"] = tempOutNames;
113 outputTypes["errorreverse"] = tempOutNames;
114 outputTypes["errorcount"] = tempOutNames;
115 outputTypes["errormatrix"] = tempOutNames;
116 outputTypes["errorchimera"] = tempOutNames;
117 outputTypes["errorref-query"] = tempOutNames;
119 catch(exception& e) {
120 m->errorOut(e, "SeqErrorCommand", "SeqErrorCommand");
125 //***************************************************************************************************************
127 SeqErrorCommand::SeqErrorCommand(string option) {
130 abort = false; calledHelp = false;
131 rdb = ReferenceDB::getInstance();
133 //allow user to run help
134 if(option == "help") { help(); abort = true; calledHelp = true; }
135 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
139 vector<string> myArray = setParameters();
141 OptionParser parser(option);
142 map<string,string> parameters = parser.getParameters();
144 ValidParameters validParameter;
145 map<string,string>::iterator it;
147 //check to make sure all parameters are valid for command
148 for (it = parameters.begin(); it != parameters.end(); it++) {
149 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
152 //initialize outputTypes
153 vector<string> tempOutNames;
154 outputTypes["errorsummary"] = tempOutNames;
155 outputTypes["errorseq"] = tempOutNames;
156 outputTypes["errorquality"] = tempOutNames;
157 outputTypes["errorqualforward"] = tempOutNames;
158 outputTypes["errorqualreverse"] = tempOutNames;
159 outputTypes["errorforward"] = tempOutNames;
160 outputTypes["errorreverse"] = tempOutNames;
161 outputTypes["errorcount"] = tempOutNames;
162 outputTypes["errormatrix"] = tempOutNames;
163 outputTypes["errorchimera"] = tempOutNames;
164 outputTypes["errorref-query"] = tempOutNames;
167 //if the user changes the input directory command factory will send this info to us in the output parameter
168 string inputDir = validParameter.validFile(parameters, "inputdir", false);
169 if (inputDir == "not found"){ inputDir = ""; }
172 it = parameters.find("fasta");
173 //user has given a template file
174 if(it != parameters.end()){
175 path = m->hasPath(it->second);
176 //if the user has not given a path then, add inputdir. else leave path alone.
177 if (path == "") { parameters["fasta"] = inputDir + it->second; }
180 it = parameters.find("reference");
181 //user has given a template file
182 if(it != parameters.end()){
183 path = m->hasPath(it->second);
184 //if the user has not given a path then, add inputdir. else leave path alone.
185 if (path == "") { parameters["reference"] = inputDir + it->second; }
188 it = parameters.find("name");
189 //user has given a names file
190 if(it != parameters.end()){
191 path = m->hasPath(it->second);
192 //if the user has not given a path then, add inputdir. else leave path alone.
193 if (path == "") { parameters["name"] = inputDir + it->second; }
196 it = parameters.find("count");
197 //user has given a names file
198 if(it != parameters.end()){
199 path = m->hasPath(it->second);
200 //if the user has not given a path then, add inputdir. else leave path alone.
201 if (path == "") { parameters["count"] = inputDir + it->second; }
204 it = parameters.find("qfile");
205 //user has given a quality score file
206 if(it != parameters.end()){
207 path = m->hasPath(it->second);
208 //if the user has not given a path then, add inputdir. else leave path alone.
209 if (path == "") { parameters["qfile"] = inputDir + it->second; }
212 it = parameters.find("report");
213 //user has given a alignment report file
214 if(it != parameters.end()){
215 path = m->hasPath(it->second);
216 //if the user has not given a path then, add inputdir. else leave path alone.
217 if (path == "") { parameters["report"] = inputDir + it->second; }
221 //check for required parameters
222 queryFileName = validParameter.validFile(parameters, "fasta", true);
223 if (queryFileName == "not found") {
224 queryFileName = m->getFastaFile();
225 if (queryFileName != "") { m->mothurOut("Using " + queryFileName + " as input file for the fasta parameter."); m->mothurOutEndLine(); }
226 else { m->mothurOut("You have no current fasta file and the fasta parameter is required."); m->mothurOutEndLine(); abort = true; }
228 else if (queryFileName == "not open") { queryFileName = ""; abort = true; }
229 else { m->setFastaFile(queryFileName); }
231 referenceFileName = validParameter.validFile(parameters, "reference", true);
232 if (referenceFileName == "not found") { m->mothurOut("reference is a required parameter for the seq.error command."); m->mothurOutEndLine(); abort = true; }
233 else if (referenceFileName == "not open") { abort = true; }
235 //check for optional parameters
236 namesFileName = validParameter.validFile(parameters, "name", true);
237 if(namesFileName == "not found"){ namesFileName = ""; }
238 else if (namesFileName == "not open") { namesFileName = ""; abort = true; }
239 else { m->setNameFile(namesFileName); }
241 //check for optional parameters
242 countfile = validParameter.validFile(parameters, "count", true);
243 if(countfile == "not found"){ countfile = ""; }
244 else if (countfile == "not open") { countfile = ""; abort = true; }
245 else { m->setCountTableFile(countfile); }
247 qualFileName = validParameter.validFile(parameters, "qfile", true);
248 if(qualFileName == "not found"){ qualFileName = ""; }
249 else if (qualFileName == "not open") { qualFileName = ""; abort = true; }
250 else { m->setQualFile(qualFileName); }
252 reportFileName = validParameter.validFile(parameters, "report", true);
253 if(reportFileName == "not found"){ reportFileName = ""; }
254 else if (reportFileName == "not open") { reportFileName = ""; abort = true; }
256 outputDir = validParameter.validFile(parameters, "outputdir", false);
257 if (outputDir == "not found"){
259 outputDir += m->hasPath(queryFileName); //if user entered a file with a path then preserve it
262 if ((countfile != "") && (namesFileName != "")) { m->mothurOut("You must enter ONLY ONE of the following: count or name."); m->mothurOutEndLine(); abort = true; }
264 //check for optional parameter and set defaults
265 // ...at some point should added some additional type checking...
266 temp = validParameter.validFile(parameters, "threshold", false); if (temp == "not found") { temp = "1.00"; }
267 m->mothurConvert(temp, threshold);
270 temp = validParameter.validFile(parameters, "save", false); if (temp == "not found"){ temp = "f"; }
271 save = m->isTrue(temp);
273 if (save) { //clear out old references
277 //this has to go after save so that if the user sets save=t and provides no reference we abort
278 referenceFileName = validParameter.validFile(parameters, "reference", true);
279 if (referenceFileName == "not found") {
280 //check for saved reference sequences
281 if (rdb->referenceSeqs.size() != 0) {
282 referenceFileName = "saved";
284 m->mothurOut("[ERROR]: You don't have any saved reference sequences and the reference parameter is a required.");
285 m->mothurOutEndLine();
288 }else if (referenceFileName == "not open") { abort = true; }
289 else { if (save) { rdb->setSavedReference(referenceFileName); } }
292 temp = validParameter.validFile(parameters, "ignorechimeras", false); if (temp == "not found") { temp = "T"; }
293 ignoreChimeras = m->isTrue(temp);
295 temp = validParameter.validFile(parameters, "aligned", false); if (temp == "not found"){ temp = "t"; }
296 aligned = m->isTrue(temp);
298 temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); }
299 m->setProcessors(temp);
300 m->mothurConvert(temp, processors);
302 substitutionMatrix.resize(6);
303 for(int i=0;i<6;i++){ substitutionMatrix[i].resize(6,0); }
305 if ((namesFileName == "") && (queryFileName != "")){
306 vector<string> files; files.push_back(queryFileName);
307 parser.getNameFile(files);
311 if((reportFileName != "" && qualFileName == "") || (reportFileName == "" && qualFileName != "")){
312 m->mothurOut("if you use either a qual file or a report file, you have to have both.");
313 m->mothurOutEndLine();
318 if(reportFileName != ""){
319 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.");
320 m->mothurOutEndLine();
327 catch(exception& e) {
328 m->errorOut(e, "SeqErrorCommand", "SeqErrorCommand");
333 //***************************************************************************************************************
335 int SeqErrorCommand::execute(){
337 if (abort == true) { if (calledHelp) { return 0; } return 2; }
339 int start = time(NULL);
344 string fileNameRoot = outputDir + m->getRootName(m->getSimpleName(queryFileName));
345 map<string, string> variables;
346 variables["[filename]"] = fileNameRoot;
347 string errorSummaryFileName = getOutputFileName("errorsummary",variables);
348 outputNames.push_back(errorSummaryFileName); outputTypes["errorsummary"].push_back(errorSummaryFileName);
350 string errorSeqFileName = getOutputFileName("errorseq",variables);
351 outputNames.push_back(errorSeqFileName); outputTypes["errorseq"].push_back(errorSeqFileName);
353 string errorChimeraFileName = getOutputFileName("errorchimera",variables);
354 outputNames.push_back(errorChimeraFileName); outputTypes["errorchimera"].push_back(errorChimeraFileName);
356 getReferences(); //read in reference sequences - make sure there's no ambiguous bases
358 if(namesFileName != "") { weights = getWeights(); }
359 else if (countfile != "") {
361 ct.readTable(countfile, false);
362 weights = ct.getNameMap();
365 vector<unsigned long long> fastaFilePos;
366 vector<unsigned long long> qFilePos;
367 vector<unsigned long long> reportFilePos;
369 setLines(queryFileName, qualFileName, reportFileName, fastaFilePos, qFilePos, reportFilePos);
371 if (m->control_pressed) { return 0; }
373 for (int i = 0; i < (fastaFilePos.size()-1); i++) {
374 lines.push_back(linePair(fastaFilePos[i], fastaFilePos[(i+1)]));
375 if (qualFileName != "") { qLines.push_back(linePair(qFilePos[i], qFilePos[(i+1)])); }
376 if (reportFileName != "" && aligned == true) { rLines.push_back(linePair(reportFilePos[i], reportFilePos[(i+1)])); }
378 if(qualFileName == "") { qLines = lines; rLines = lines; } //fills with duds
379 if(aligned == false){ rLines = lines; }
381 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
383 numSeqs = driver(queryFileName, qualFileName, reportFileName, errorSummaryFileName, errorSeqFileName, errorChimeraFileName, lines[0], qLines[0], rLines[0]);
386 numSeqs = createProcesses(queryFileName, qualFileName, reportFileName, errorSummaryFileName, errorSeqFileName, errorChimeraFileName);
389 numSeqs = driver(queryFileName, qualFileName, reportFileName, errorSummaryFileName, errorSeqFileName, errorChimeraFileName, lines[0], qLines[0], rLines[0]);
392 if(qualFileName != ""){
393 printErrorQuality(qScoreErrorMap);
394 printQualityFR(qualForwardMap, qualReverseMap);
397 printErrorFRFile(errorForward, errorReverse);
399 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; }
401 string errorCountFileName = getOutputFileName("errorcount",variables);
402 ofstream errorCountFile;
403 m->openOutputFile(errorCountFileName, errorCountFile);
404 outputNames.push_back(errorCountFileName); outputTypes["errorcount"].push_back(errorCountFileName);
405 m->mothurOut("Overall error rate:\t" + toString((double)(totalBases - totalMatches) / (double)totalBases) + "\n");
406 m->mothurOut("Errors\tSequences\n");
407 errorCountFile << "Errors\tSequences\n";
408 for(int i=0;i<misMatchCounts.size();i++){
409 m->mothurOut(toString(i) + '\t' + toString(misMatchCounts[i]) + '\n');
410 errorCountFile << i << '\t' << misMatchCounts[i] << endl;
412 errorCountFile.close();
414 // if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; }
418 string megAlignmentFileName = getOutputFileName("errorref-query",variables);
419 ofstream megAlignmentFile;
420 m->openOutputFile(megAlignmentFileName, megAlignmentFile);
421 outputNames.push_back(megAlignmentFileName); outputTypes["errorref-query"].push_back(megAlignmentFileName);
423 for(int i=0;i<numRefs;i++){
424 megAlignmentFile << referenceSeqs[i].getInlineSeq() << endl;
425 megAlignmentFile << megaAlignVector[i] << endl;
427 megAlignmentFile.close();
430 m->mothurOut("It took " + toString(time(NULL) - start) + " secs to check " + toString(numSeqs) + " sequences.");
431 m->mothurOutEndLine();
433 //set fasta file as new current fastafile
435 itTypes = outputTypes.find("errorseq");
436 if (itTypes != outputTypes.end()) {
437 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setFastaFile(current); }
440 m->mothurOutEndLine();
441 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
442 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
443 m->mothurOutEndLine();
447 catch(exception& e) {
448 m->errorOut(e, "SeqErrorCommand", "execute");
453 //**********************************************************************************************************************
455 int SeqErrorCommand::createProcesses(string filename, string qFileName, string rFileName, string summaryFileName, string errorOutputFileName, string chimeraOutputFileName) {
459 map<char, vector<int> >::iterator it;
461 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
463 //loop through and create all the processes you want
464 while (process != processors) {
468 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
472 num = driver(filename, qFileName, rFileName, summaryFileName + toString(getpid()) + ".temp", errorOutputFileName+ toString(getpid()) + ".temp", chimeraOutputFileName + toString(getpid()) + ".temp", lines[process], qLines[process], rLines[process]);
474 //pass groupCounts to parent
476 string tempFile = filename + toString(getpid()) + ".info.temp";
477 m->openOutputFile(tempFile, out);
479 //output totalBases and totalMatches
480 out << num << '\t' << totalBases << '\t' << totalMatches << endl << endl;
482 //output substitutionMatrix
483 for(int i = 0; i < substitutionMatrix.size(); i++) {
484 for (int j = 0; j < substitutionMatrix[i].size(); j++) {
485 out << substitutionMatrix[i][j] << '\t';
491 //output qScoreErrorMap
492 for (it = qScoreErrorMap.begin(); it != qScoreErrorMap.end(); it++) {
493 vector<int> thisScoreErrorMap = it->second;
494 out << it->first << '\t';
495 for (int i = 0; i < thisScoreErrorMap.size(); i++) {
496 out << thisScoreErrorMap[i] << '\t';
502 //output qualForwardMap
503 for(int i = 0; i < qualForwardMap.size(); i++) {
504 for (int j = 0; j < qualForwardMap[i].size(); j++) {
505 out << qualForwardMap[i][j] << '\t';
511 //output qualReverseMap
512 for(int i = 0; i < qualReverseMap.size(); i++) {
513 for (int j = 0; j < qualReverseMap[i].size(); j++) {
514 out << qualReverseMap[i][j] << '\t';
521 //output errorForward
522 for (it = errorForward.begin(); it != errorForward.end(); it++) {
523 vector<int> thisErrorForward = it->second;
524 out << it->first << '\t';
525 for (int i = 0; i < thisErrorForward.size(); i++) {
526 out << thisErrorForward[i] << '\t';
532 //output errorReverse
533 for (it = errorReverse.begin(); it != errorReverse.end(); it++) {
534 vector<int> thisErrorReverse = it->second;
535 out << it->first << '\t';
536 for (int i = 0; i < thisErrorReverse.size(); i++) {
537 out << thisErrorReverse[i] << '\t';
543 //output misMatchCounts
544 out << misMatchCounts.size() << endl;
545 for (int j = 0; j < misMatchCounts.size(); j++) {
546 out << misMatchCounts[j] << '\t';
551 //output megaAlignVector
552 for (int j = 0; j < megaAlignVector.size(); j++) {
553 out << megaAlignVector[j] << endl;
561 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
562 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
568 num = driver(filename, qFileName, rFileName, summaryFileName, errorOutputFileName, chimeraOutputFileName, lines[0], qLines[0], rLines[0]);
570 //force parent to wait until all the processes are done
571 for (int i=0;i<processIDS.size();i++) {
572 int temp = processIDS[i];
577 for(int i=0;i<processIDS.size();i++){
579 m->mothurOut("Appending files from process " + toString(processIDS[i])); m->mothurOutEndLine();
581 m->appendFiles((summaryFileName + toString(processIDS[i]) + ".temp"), summaryFileName);
582 m->mothurRemove((summaryFileName + toString(processIDS[i]) + ".temp"));
583 m->appendFiles((errorOutputFileName + toString(processIDS[i]) + ".temp"), errorOutputFileName);
584 m->mothurRemove((errorOutputFileName + toString(processIDS[i]) + ".temp"));
585 m->appendFiles((chimeraOutputFileName + toString(processIDS[i]) + ".temp"), chimeraOutputFileName);
586 m->mothurRemove((chimeraOutputFileName + toString(processIDS[i]) + ".temp"));
589 string tempFile = filename + toString(processIDS[i]) + ".info.temp";
590 m->openInputFile(tempFile, in);
592 //input totalBases and totalMatches
593 int tempBases, tempMatches, tempNumSeqs;
594 in >> tempNumSeqs >> tempBases >> tempMatches; m->gobble(in);
595 totalBases += tempBases; totalMatches += tempMatches; num += tempNumSeqs;
597 //input substitutionMatrix
599 for(int i = 0; i < substitutionMatrix.size(); i++) {
600 for (int j = 0; j < substitutionMatrix[i].size(); j++) {
601 in >> tempNum; substitutionMatrix[i][j] += tempNum;
607 //input qScoreErrorMap
609 for (int i = 0; i < qScoreErrorMap.size(); i++) {
611 vector<int> thisScoreErrorMap = qScoreErrorMap[first];
613 for (int i = 0; i < thisScoreErrorMap.size(); i++) {
614 in >> tempNum; thisScoreErrorMap[i] += tempNum;
616 qScoreErrorMap[first] = thisScoreErrorMap;
621 //input qualForwardMap
622 for(int i = 0; i < qualForwardMap.size(); i++) {
623 for (int j = 0; j < qualForwardMap[i].size(); j++) {
624 in >> tempNum; qualForwardMap[i][j] += tempNum;
630 //input qualReverseMap
631 for(int i = 0; i < qualReverseMap.size(); i++) {
632 for (int j = 0; j < qualReverseMap[i].size(); j++) {
633 in >> tempNum; qualReverseMap[i][j] += tempNum;
640 for (int i = 0; i < errorForward.size(); i++) {
642 vector<int> thisErrorForward = errorForward[first];
644 for (int i = 0; i < thisErrorForward.size(); i++) {
645 in >> tempNum; thisErrorForward[i] += tempNum;
647 errorForward[first] = thisErrorForward;
653 for (int i = 0; i < errorReverse.size(); i++) {
655 vector<int> thisErrorReverse = errorReverse[first];
657 for (int i = 0; i < thisErrorReverse.size(); i++) {
658 in >> tempNum; thisErrorReverse[i] += tempNum;
660 errorReverse[first] = thisErrorReverse;
665 //input misMatchCounts
667 in >> misMatchSize; m->gobble(in);
668 if (misMatchSize > misMatchCounts.size()) { misMatchCounts.resize(misMatchSize, 0); }
669 for (int j = 0; j < misMatchSize; j++) {
670 in >> tempNum; misMatchCounts[j] += tempNum;
674 //input megaAlignVector
676 for (int j = 0; j < megaAlignVector.size(); j++) {
677 thisLine = m->getline(in); m->gobble(in); megaAlignVector[j] += thisLine + '\n';
681 in.close(); m->mothurRemove(tempFile);
687 catch(exception& e) {
688 m->errorOut(e, "SeqErrorCommand", "createProcesses");
693 //**********************************************************************************************************************
695 int SeqErrorCommand::driver(string filename, string qFileName, string rFileName, string summaryFileName, string errorOutputFileName, string chimeraOutputFileName, linePair line, linePair qline, linePair rline) {
699 QualityScores quality;
701 misMatchCounts.resize(11, 0);
705 map<string, int>::iterator it;
706 qScoreErrorMap['m'].assign(101, 0);
707 qScoreErrorMap['s'].assign(101, 0);
708 qScoreErrorMap['i'].assign(101, 0);
709 qScoreErrorMap['a'].assign(101, 0);
711 errorForward['m'].assign(maxLength,0);
712 errorForward['s'].assign(maxLength,0);
713 errorForward['i'].assign(maxLength,0);
714 errorForward['d'].assign(maxLength,0);
715 errorForward['a'].assign(maxLength,0);
717 errorReverse['m'].assign(maxLength,0);
718 errorReverse['s'].assign(maxLength,0);
719 errorReverse['i'].assign(maxLength,0);
720 errorReverse['d'].assign(maxLength,0);
721 errorReverse['a'].assign(maxLength,0);
723 //open inputfiles and go to beginning place for this processor
725 m->openInputFile(filename, queryFile);
727 queryFile.seekg(line.start);
731 if((qFileName != "" && rFileName != "" && aligned)){
732 m->openInputFile(qFileName, qualFile);
733 qualFile.seekg(qline.start);
736 if (rline.start == 0) { report = ReportFile(reportFile, rFileName); }
738 m->openInputFile(rFileName, reportFile);
739 reportFile.seekg(rline.start);
742 qualForwardMap.resize(maxLength);
743 qualReverseMap.resize(maxLength);
744 for(int i=0;i<maxLength;i++){
745 qualForwardMap[i].assign(101,0);
746 qualReverseMap[i].assign(101,0);
749 else if(qFileName != "" && !aligned){
751 m->openInputFile(qFileName, qualFile);
752 qualFile.seekg(qline.start);
754 qualForwardMap.resize(maxLength);
755 qualReverseMap.resize(maxLength);
756 for(int i=0;i<maxLength;i++){
757 qualForwardMap[i].assign(101,0);
758 qualReverseMap[i].assign(101,0);
762 ofstream outChimeraReport;
763 m->openOutputFile(chimeraOutputFileName, outChimeraReport);
766 RefChimeraTest chimeraTest;
768 chimeraTest = RefChimeraTest(referenceSeqs, aligned);
769 if (line.start == 0) { chimeraTest.printHeader(outChimeraReport); }
772 ofstream errorSummaryFile;
773 m->openOutputFile(summaryFileName, errorSummaryFile);
774 if (line.start == 0) { printErrorHeader(errorSummaryFile); }
776 ofstream errorSeqFile;
777 m->openOutputFile(errorOutputFileName, errorSeqFile);
779 megaAlignVector.assign(numRefs, "");
787 Sequence query(queryFile);
789 int numParentSeqs = -1;
790 int closestRefIndex = -1;
792 string querySeq = query.getAligned();
793 if (!aligned) { querySeq = query.getUnaligned(); }
795 numParentSeqs = chimeraTest.analyzeQuery(query.getName(), querySeq, outChimeraReport);
797 closestRefIndex = chimeraTest.getClosestRefIndex();
799 reference = referenceSeqs[closestRefIndex];
801 reference.setAligned(chimeraTest.getClosestRefAlignment());
802 query.setAligned(chimeraTest.getQueryAlignment());
804 if(numParentSeqs > 1 && ignoreChimeras == 1) { ignoreSeq = 1; }
805 else { ignoreSeq = 0; }
809 getErrors(query, reference, minCompare);
811 if((namesFileName != "") || (countfile != "")){
812 it = weights.find(query.getName());
813 minCompare.weight = it->second;
815 else{ minCompare.weight = 1; }
818 printErrorData(minCompare, numParentSeqs, errorSummaryFile, errorSeqFile);
821 for(int i=0;i<minCompare.sequence.length();i++){
822 char letter = minCompare.sequence[i];
824 errorForward[letter][i] += minCompare.weight;
825 errorReverse[letter][minCompare.total-i-1] += minCompare.weight;
830 if(aligned && qualFileName != "" && reportFileName != ""){
831 report = ReportFile(reportFile);
833 // int origLength = report.getQueryLength();
834 int startBase = report.getQueryStart();
835 int endBase = report.getQueryEnd();
837 quality = QualityScores(qualFile);
840 quality.updateQScoreErrorMap(qScoreErrorMap, minCompare.sequence, startBase, endBase, minCompare.weight);
841 quality.updateForwardMap(qualForwardMap, startBase, endBase, minCompare.weight);
842 quality.updateReverseMap(qualReverseMap, startBase, endBase, minCompare.weight);
845 else if(aligned == false && qualFileName != ""){
847 quality = QualityScores(qualFile);
848 int qualityLength = quality.getLength();
850 if(qualityLength != query.getNumBases()){ cout << "warning - quality and fasta sequence files do not match at " << query.getName() << '\t' << qualityLength <<'\t' << query.getNumBases() << endl; }
853 int endBase = qualityLength;
856 quality.updateQScoreErrorMap(qScoreErrorMap, minCompare.sequence, startBase, endBase, minCompare.weight);
857 quality.updateForwardMap(qualForwardMap, startBase, endBase, minCompare.weight);
858 quality.updateReverseMap(qualReverseMap, startBase, endBase, minCompare.weight);
862 if(minCompare.errorRate <= threshold && !ignoreSeq){
863 totalBases += (minCompare.total * minCompare.weight);
864 totalMatches += minCompare.matches * minCompare.weight;
865 if(minCompare.mismatches > maxMismatch){
866 maxMismatch = minCompare.mismatches;
867 misMatchCounts.resize(maxMismatch + 1, 0);
869 misMatchCounts[minCompare.mismatches] += minCompare.weight;
872 megaAlignVector[closestRefIndex] += query.getInlineSeq() + '\n';
877 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
878 unsigned long long pos = queryFile.tellg();
879 if ((pos == -1) || (pos >= line.end)) { break; }
881 if (queryFile.eof()) { break; }
884 if(index % 100 == 0){ m->mothurOutJustToScreen(toString(index)+"\n"); }
887 outChimeraReport.close();
888 errorSummaryFile.close();
889 errorSeqFile.close();
891 if(qFileName != "" && rFileName != "") { reportFile.close(); qualFile.close(); }
892 else if(qFileName != "" && aligned == false){ qualFile.close(); }
895 m->mothurOutJustToScreen(toString(index)+"\n");
899 catch(exception& e) {
900 m->errorOut(e, "SeqErrorCommand", "driver");
905 //***************************************************************************************************************
907 void SeqErrorCommand::getReferences(){
909 int numAmbigSeqs = 0;
912 int minEndPos = 100000;
914 if (referenceFileName == "saved") {
915 int start = time(NULL);
916 m->mothurOutEndLine(); m->mothurOut("Using sequences from " + rdb->getSavedReference() + " that are saved in memory."); m->mothurOutEndLine();
918 for (int i = 0; i < rdb->referenceSeqs.size(); i++) {
919 int numAmbigs = rdb->referenceSeqs[i].getAmbigBases();
920 if(numAmbigs > 0){ numAmbigSeqs++; }
922 // int startPos = rdb->referenceSeqs[i].getStartPos();
923 // if(startPos > maxStartPos) { maxStartPos = startPos; }
925 // int endPos = rdb->referenceSeqs[i].getEndPos();
926 // if(endPos < minEndPos) { minEndPos = endPos; }
927 if (rdb->referenceSeqs[i].getNumBases() == 0) {
928 m->mothurOut("[WARNING]: " + rdb->referenceSeqs[i].getName() + " is blank, ignoring.");m->mothurOutEndLine();
930 referenceSeqs.push_back(rdb->referenceSeqs[i]);
934 referenceFileName = rdb->getSavedReference();
936 m->mothurOut("It took " + toString(time(NULL) - start) + " to load " + toString(referenceSeqs.size()) + " sequences.");m->mothurOutEndLine();
939 int start = time(NULL);
941 ifstream referenceFile;
942 m->openInputFile(referenceFileName, referenceFile);
944 while(referenceFile){
945 Sequence currentSeq(referenceFile);
946 int numAmbigs = currentSeq.getAmbigBases();
947 if(numAmbigs > 0){ numAmbigSeqs++; }
949 // int startPos = currentSeq.getStartPos();
950 // if(startPos > maxStartPos) { maxStartPos = startPos; }
952 // int endPos = currentSeq.getEndPos();
953 // if(endPos < minEndPos) { minEndPos = endPos; }
954 if (currentSeq.getNumBases() == 0) {
955 m->mothurOut("[WARNING]: " + currentSeq.getName() + " is blank, ignoring.");m->mothurOutEndLine();
957 referenceSeqs.push_back(currentSeq);
958 if (rdb->save) { rdb->referenceSeqs.push_back(currentSeq); }
961 m->gobble(referenceFile);
963 referenceFile.close();
965 m->mothurOut("It took " + toString(time(NULL) - start) + " to read " + toString(referenceSeqs.size()) + " sequences.");m->mothurOutEndLine();
968 numRefs = referenceSeqs.size();
970 for(int i=0;i<numRefs;i++){
971 referenceSeqs[i].padToPos(maxStartPos);
972 referenceSeqs[i].padFromPos(minEndPos);
975 if(numAmbigSeqs != 0){
976 m->mothurOut("Warning: " + toString(numAmbigSeqs) + " reference sequences have ambiguous bases, these bases will be ignored\n");
980 catch(exception& e) {
981 m->errorOut(e, "SeqErrorCommand", "getReferences");
986 //***************************************************************************************************************
988 int SeqErrorCommand::getErrors(Sequence query, Sequence reference, Compare& errors){
990 if(query.getAlignLength() != reference.getAlignLength()){
991 m->mothurOut("Warning: " + toString(query.getName()) + " and " + toString(reference.getName()) + " are different lengths\n");
993 int alignLength = query.getAlignLength();
995 string q = query.getAligned();
996 string r = reference.getAligned();
1001 errors.sequence = "";
1002 for(int i=0;i<alignLength;i++){
1004 if(q[i] != '.' && r[i] != '.' && (q[i] != '-' || r[i] != '-')){ // no missing data and no double gaps
1009 if(r[i] == 'A'){ errors.AA++; errors.matches++; errors.sequence += 'm'; }
1010 if(r[i] == 'T'){ errors.AT++; errors.sequence += 's'; }
1011 if(r[i] == 'G'){ errors.AG++; errors.sequence += 's'; }
1012 if(r[i] == 'C'){ errors.AC++; errors.sequence += 's'; }
1013 if(r[i] == '-'){ errors.Ai++; errors.sequence += 'i'; }
1015 else if(q[i] == 'T'){
1016 if(r[i] == 'A'){ errors.TA++; errors.sequence += 's'; }
1017 if(r[i] == 'T'){ errors.TT++; errors.matches++; errors.sequence += 'm'; }
1018 if(r[i] == 'G'){ errors.TG++; errors.sequence += 's'; }
1019 if(r[i] == 'C'){ errors.TC++; errors.sequence += 's'; }
1020 if(r[i] == '-'){ errors.Ti++; errors.sequence += 'i'; }
1022 else if(q[i] == 'G'){
1023 if(r[i] == 'A'){ errors.GA++; errors.sequence += 's'; }
1024 if(r[i] == 'T'){ errors.GT++; errors.sequence += 's'; }
1025 if(r[i] == 'G'){ errors.GG++; errors.matches++; errors.sequence += 'm'; }
1026 if(r[i] == 'C'){ errors.GC++; errors.sequence += 's'; }
1027 if(r[i] == '-'){ errors.Gi++; errors.sequence += 'i'; }
1029 else if(q[i] == 'C'){
1030 if(r[i] == 'A'){ errors.CA++; errors.sequence += 's'; }
1031 if(r[i] == 'T'){ errors.CT++; errors.sequence += 's'; }
1032 if(r[i] == 'G'){ errors.CG++; errors.sequence += 's'; }
1033 if(r[i] == 'C'){ errors.CC++; errors.matches++; errors.sequence += 'm'; }
1034 if(r[i] == '-'){ errors.Ci++; errors.sequence += 'i'; }
1036 else if(q[i] == 'N'){
1037 if(r[i] == 'A'){ errors.NA++; errors.sequence += 'a'; }
1038 if(r[i] == 'T'){ errors.NT++; errors.sequence += 'a'; }
1039 if(r[i] == 'G'){ errors.NG++; errors.sequence += 'a'; }
1040 if(r[i] == 'C'){ errors.NC++; errors.sequence += 'a'; }
1041 if(r[i] == '-'){ errors.Ni++; errors.sequence += 'a'; }
1043 else if(q[i] == '-' && r[i] != '-'){
1044 if(r[i] == 'A'){ errors.dA++; errors.sequence += 'd'; }
1045 if(r[i] == 'T'){ errors.dT++; errors.sequence += 'd'; }
1046 if(r[i] == 'G'){ errors.dG++; errors.sequence += 'd'; }
1047 if(r[i] == 'C'){ errors.dC++; errors.sequence += 'd'; }
1054 errors.sequence += 'd'; errors.total++;
1057 errors.sequence += 'r';
1061 else if(q[i] == '.' && r[i] != '.'){ // reference extends beyond query
1062 if(started == 1){ break; }
1064 else if(q[i] != '.' && r[i] == '.'){ // query extends beyond reference
1065 if(started == 1){ break; }
1067 else if(q[i] == '.' && r[i] == '.'){ // both are missing data
1068 if(started == 1){ break; }
1072 errors.mismatches = errors.total-errors.matches;
1073 if(errors.total != 0){ errors.errorRate = (double)(errors.total-errors.matches) / (double)errors.total; }
1074 else{ errors.errorRate = 0; }
1076 errors.queryName = query.getName();
1077 errors.refName = reference.getName();
1081 catch(exception& e) {
1082 m->errorOut(e, "SeqErrorCommand", "getErrors");
1087 //***************************************************************************************************************
1089 map<string, int> SeqErrorCommand::getWeights(){
1091 m->openInputFile(namesFileName, nameFile);
1094 string redundantSeqs;
1095 map<string, int> nameCountMap;
1098 nameFile >> seqName >> redundantSeqs;
1099 nameCountMap[seqName] = m->getNumNames(redundantSeqs);
1100 m->gobble(nameFile);
1105 return nameCountMap;
1108 //***************************************************************************************************************
1110 void SeqErrorCommand::printErrorHeader(ofstream& errorSummaryFile){
1112 errorSummaryFile << "query\treference\tweight\t";
1113 errorSummaryFile << "AA\tAT\tAG\tAC\tTA\tTT\tTG\tTC\tGA\tGT\tGG\tGC\tCA\tCT\tCG\tCC\tNA\tNT\tNG\tNC\tAi\tTi\tGi\tCi\tNi\tdA\tdT\tdG\tdC\t";
1114 errorSummaryFile << "insertions\tdeletions\tsubstitutions\tambig\tmatches\tmismatches\ttotal\terror\tnumparents\n";
1116 errorSummaryFile << setprecision(6);
1117 errorSummaryFile.setf(ios::fixed);
1119 catch(exception& e) {
1120 m->errorOut(e, "SeqErrorCommand", "printErrorHeader");
1125 //***************************************************************************************************************
1127 void SeqErrorCommand::printErrorData(Compare error, int numParentSeqs, ofstream& errorSummaryFile, ofstream& errorSeqFile){
1130 errorSummaryFile << error.queryName << '\t' << error.refName << '\t' << error.weight << '\t';
1131 errorSummaryFile << error.AA << '\t' << error.AT << '\t' << error.AG << '\t' << error.AC << '\t';
1132 errorSummaryFile << error.TA << '\t' << error.TT << '\t' << error.TG << '\t' << error.TC << '\t';
1133 errorSummaryFile << error.GA << '\t' << error.GT << '\t' << error.GG << '\t' << error.GC << '\t';
1134 errorSummaryFile << error.CA << '\t' << error.CT << '\t' << error.CG << '\t' << error.CC << '\t';
1135 errorSummaryFile << error.NA << '\t' << error.NT << '\t' << error.NG << '\t' << error.NC << '\t';
1136 errorSummaryFile << error.Ai << '\t' << error.Ti << '\t' << error.Gi << '\t' << error.Ci << '\t' << error.Ni << '\t';
1137 errorSummaryFile << error.dA << '\t' << error.dT << '\t' << error.dG << '\t' << error.dC << '\t';
1139 errorSummaryFile << error.Ai + error.Ti + error.Gi + error.Ci << '\t'; //insertions
1140 errorSummaryFile << error.dA + error.dT + error.dG + error.dC << '\t'; //deletions
1141 errorSummaryFile << error.mismatches - (error.Ai + error.Ti + error.Gi + error.Ci) - (error.dA + error.dT + error.dG + error.dC) - (error.NA + error.NT + error.NG + error.NC + error.Ni) << '\t'; //substitutions
1142 errorSummaryFile << error.NA + error.NT + error.NG + error.NC + error.Ni << '\t'; //ambiguities
1143 errorSummaryFile << error.matches << '\t' << error.mismatches << '\t' << error.total << '\t' << error.errorRate << '\t' << numParentSeqs << endl;
1145 errorSeqFile << '>' << error.queryName << "\tref:" << error.refName << '\n' << error.sequence << endl;
1147 int a=0; int t=1; int g=2; int c=3;
1150 if(numParentSeqs == 1 || ignoreChimeras == 0){
1151 substitutionMatrix[a][a] += error.weight * error.AA;
1152 substitutionMatrix[a][t] += error.weight * error.TA;
1153 substitutionMatrix[a][g] += error.weight * error.GA;
1154 substitutionMatrix[a][c] += error.weight * error.CA;
1155 substitutionMatrix[a][gap] += error.weight * error.dA;
1156 substitutionMatrix[a][n] += error.weight * error.NA;
1158 substitutionMatrix[t][a] += error.weight * error.AT;
1159 substitutionMatrix[t][t] += error.weight * error.TT;
1160 substitutionMatrix[t][g] += error.weight * error.GT;
1161 substitutionMatrix[t][c] += error.weight * error.CT;
1162 substitutionMatrix[t][gap] += error.weight * error.dT;
1163 substitutionMatrix[t][n] += error.weight * error.NT;
1165 substitutionMatrix[g][a] += error.weight * error.AG;
1166 substitutionMatrix[g][t] += error.weight * error.TG;
1167 substitutionMatrix[g][g] += error.weight * error.GG;
1168 substitutionMatrix[g][c] += error.weight * error.CG;
1169 substitutionMatrix[g][gap] += error.weight * error.dG;
1170 substitutionMatrix[g][n] += error.weight * error.NG;
1172 substitutionMatrix[c][a] += error.weight * error.AC;
1173 substitutionMatrix[c][t] += error.weight * error.TC;
1174 substitutionMatrix[c][g] += error.weight * error.GC;
1175 substitutionMatrix[c][c] += error.weight * error.CC;
1176 substitutionMatrix[c][gap] += error.weight * error.dC;
1177 substitutionMatrix[c][n] += error.weight * error.NC;
1179 substitutionMatrix[gap][a] += error.weight * error.Ai;
1180 substitutionMatrix[gap][t] += error.weight * error.Ti;
1181 substitutionMatrix[gap][g] += error.weight * error.Gi;
1182 substitutionMatrix[gap][c] += error.weight * error.Ci;
1183 substitutionMatrix[gap][n] += error.weight * error.Ni;
1186 catch(exception& e) {
1187 m->errorOut(e, "SeqErrorCommand", "printErrorData");
1192 //***************************************************************************************************************
1194 void SeqErrorCommand::printSubMatrix(){
1196 string fileNameRoot = outputDir + m->getRootName(m->getSimpleName(queryFileName));
1197 map<string, string> variables;
1198 variables["[filename]"] = fileNameRoot;
1199 string subMatrixFileName = getOutputFileName("errormatrix",variables);
1200 ofstream subMatrixFile;
1201 m->openOutputFile(subMatrixFileName, subMatrixFile);
1202 outputNames.push_back(subMatrixFileName); outputTypes["errormatrix"].push_back(subMatrixFileName);
1203 vector<string> bases(6);
1210 vector<int> refSums(5,1);
1212 for(int i=0;i<5;i++){
1213 subMatrixFile << "\tr" << bases[i];
1215 for(int j=0;j<6;j++){
1216 refSums[i] += substitutionMatrix[i][j];
1219 subMatrixFile << endl;
1221 for(int i=0;i<6;i++){
1222 subMatrixFile << 'q' << bases[i];
1223 for(int j=0;j<5;j++){
1224 subMatrixFile << '\t' << substitutionMatrix[j][i];
1226 subMatrixFile << endl;
1229 subMatrixFile << "total";
1230 for(int i=0;i<5;i++){
1231 subMatrixFile << '\t' << refSums[i];
1233 subMatrixFile << endl;
1234 subMatrixFile.close();
1236 catch(exception& e) {
1237 m->errorOut(e, "SeqErrorCommand", "printSubMatrix");
1242 //***************************************************************************************************************
1244 void SeqErrorCommand::printErrorFRFile(map<char, vector<int> > errorForward, map<char, vector<int> > errorReverse){
1246 string fileNameRoot = outputDir + m->getRootName(m->getSimpleName(queryFileName));
1247 map<string, string> variables;
1248 variables["[filename]"] = fileNameRoot;
1249 string errorForwardFileName = getOutputFileName("errorforward",variables);
1250 ofstream errorForwardFile;
1251 m->openOutputFile(errorForwardFileName, errorForwardFile);
1252 outputNames.push_back(errorForwardFileName); outputTypes["errorforward"].push_back(errorForwardFileName);
1254 errorForwardFile << "position\ttotalseqs\tmatch\tsubstitution\tinsertion\tdeletion\tambiguous" << endl;
1255 for(int i=0;i<maxLength;i++){
1256 float match = (float)errorForward['m'][i];
1257 float subst = (float)errorForward['s'][i];
1258 float insert = (float)errorForward['i'][i];
1259 float del = (float)errorForward['d'][i];
1260 float amb = (float)errorForward['a'][i];
1261 float total = match + subst + insert + del + amb;
1262 if(total == 0){ break; }
1263 errorForwardFile << i+1 << '\t' << total << '\t' << match/total << '\t' << subst/total << '\t' << insert/total << '\t' << del/total << '\t' << amb/total << endl;
1265 errorForwardFile.close();
1267 string errorReverseFileName = getOutputFileName("errorreverse",variables);
1268 ofstream errorReverseFile;
1269 m->openOutputFile(errorReverseFileName, errorReverseFile);
1270 outputNames.push_back(errorReverseFileName); outputTypes["errorreverse"].push_back(errorReverseFileName);
1272 errorReverseFile << "position\ttotalseqs\tmatch\tsubstitution\tinsertion\tdeletion\tambiguous" << endl;
1273 for(int i=0;i<maxLength;i++){
1274 float match = (float)errorReverse['m'][i];
1275 float subst = (float)errorReverse['s'][i];
1276 float insert = (float)errorReverse['i'][i];
1277 float del = (float)errorReverse['d'][i];
1278 float amb = (float)errorReverse['a'][i];
1279 float total = match + subst + insert + del + amb;
1280 if(total == 0){ break; }
1281 errorReverseFile << i+1 << '\t' << total << '\t' << match/total << '\t' << subst/total << '\t' << insert/total << '\t' << del/total << '\t' << amb/total << endl;
1283 errorReverseFile.close();
1285 catch(exception& e) {
1286 m->errorOut(e, "SeqErrorCommand", "printErrorFRFile");
1291 //***************************************************************************************************************
1293 void SeqErrorCommand::printErrorQuality(map<char, vector<int> > qScoreErrorMap){
1295 string fileNameRoot = outputDir + m->getRootName(m->getSimpleName(queryFileName));
1296 map<string, string> variables;
1297 variables["[filename]"] = fileNameRoot;
1298 string errorQualityFileName = getOutputFileName("errorquality",variables);
1299 ofstream errorQualityFile;
1300 m->openOutputFile(errorQualityFileName, errorQualityFile);
1301 outputNames.push_back(errorQualityFileName); outputTypes["errorquality"].push_back(errorQualityFileName);
1303 errorQualityFile << "qscore\tmatches\tsubstitutions\tinsertions\tambiguous" << endl;
1304 for(int i=0;i<101;i++){
1305 errorQualityFile << i << '\t' << qScoreErrorMap['m'][i] << '\t' << qScoreErrorMap['s'][i] << '\t' << qScoreErrorMap['i'][i] << '\t'<< qScoreErrorMap['a'][i] << endl;
1307 errorQualityFile.close();
1309 catch(exception& e) {
1310 m->errorOut(e, "SeqErrorCommand", "printErrorQuality");
1315 //***************************************************************************************************************
1317 void SeqErrorCommand::printQualityFR(vector<vector<int> > qualForwardMap, vector<vector<int> > qualReverseMap){
1321 int numColumns = qualForwardMap[0].size();
1323 for(int i=0;i<qualForwardMap.size();i++){
1324 for(int j=0;j<numColumns;j++){
1325 if(qualForwardMap[i][j] != 0){
1326 if(numRows < i) { numRows = i+20; }
1330 string fileNameRoot = outputDir + m->getRootName(m->getSimpleName(queryFileName));
1331 map<string, string> variables;
1332 variables["[filename]"] = fileNameRoot;
1333 string qualityForwardFileName = getOutputFileName("errorqualforward",variables);
1334 ofstream qualityForwardFile;
1335 m->openOutputFile(qualityForwardFileName, qualityForwardFile);
1336 outputNames.push_back(qualityForwardFileName); outputTypes["errorqualforward"].push_back(qualityForwardFileName);
1338 for(int i=0;i<numColumns;i++){ qualityForwardFile << '\t' << i; } qualityForwardFile << endl;
1340 for(int i=0;i<numRows;i++){
1341 qualityForwardFile << i+1;
1342 for(int j=0;j<numColumns;j++){
1343 qualityForwardFile << '\t' << qualForwardMap[i][j];
1346 qualityForwardFile << endl;
1348 qualityForwardFile.close();
1351 string qualityReverseFileName = getOutputFileName("errorqualreverse",variables);
1352 ofstream qualityReverseFile;
1353 m->openOutputFile(qualityReverseFileName, qualityReverseFile);
1354 outputNames.push_back(qualityReverseFileName); outputTypes["errorqualreverse"].push_back(qualityReverseFileName);
1356 for(int i=0;i<numColumns;i++){ qualityReverseFile << '\t' << i; } qualityReverseFile << endl;
1357 for(int i=0;i<numRows;i++){
1359 qualityReverseFile << i+1;
1360 for(int j=0;j<numColumns;j++){
1361 qualityReverseFile << '\t' << qualReverseMap[i][j];
1363 qualityReverseFile << endl;
1365 qualityReverseFile.close();
1367 catch(exception& e) {
1368 m->errorOut(e, "SeqErrorCommand", "printQualityFR");
1374 /**************************************************************************************************/
1376 int SeqErrorCommand::setLines(string filename, string qfilename, string rfilename, vector<unsigned long long>& fastaFilePos, vector<unsigned long long>& qfileFilePos, vector<unsigned long long>& rfileFilePos) {
1378 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
1379 //set file positions for fasta file
1380 fastaFilePos = m->divideFile(filename, processors);
1382 if (qfilename == "") { return processors; }
1384 //get name of first sequence in each chunk
1385 map<string, int> firstSeqNames;
1386 for (int i = 0; i < (fastaFilePos.size()-1); i++) {
1388 m->openInputFile(filename, in);
1389 in.seekg(fastaFilePos[i]);
1392 firstSeqNames[temp.getName()] = i;
1397 //make copy to use below
1398 map<string, int> firstSeqNamesReport = firstSeqNames;
1400 //seach for filePos of each first name in the qfile and save in qfileFilePos
1402 m->openInputFile(qfilename, inQual);
1405 while(!inQual.eof()){
1406 input = m->getline(inQual);
1408 if (input.length() != 0) {
1409 if(input[0] == '>'){ //this is a sequence name line
1410 istringstream nameStream(input);
1412 string sname = ""; nameStream >> sname;
1413 sname = sname.substr(1);
1415 m->checkName(sname);
1417 map<string, int>::iterator it = firstSeqNames.find(sname);
1419 if(it != firstSeqNames.end()) { //this is the start of a new chunk
1420 unsigned long long pos = inQual.tellg();
1421 qfileFilePos.push_back(pos - input.length() - 1);
1422 firstSeqNames.erase(it);
1427 if (firstSeqNames.size() == 0) { break; }
1431 if (firstSeqNames.size() != 0) {
1432 for (map<string, int>::iterator it = firstSeqNames.begin(); it != firstSeqNames.end(); it++) {
1433 m->mothurOut(it->first + " is in your fasta file and not in your quality file, aborting."); m->mothurOutEndLine();
1435 m->control_pressed = true;
1439 //get last file position of qfile
1441 unsigned long long size;
1443 //get num bytes in file
1444 pFile = fopen (qfilename.c_str(),"rb");
1445 if (pFile==NULL) perror ("Error opening file");
1447 fseek (pFile, 0, SEEK_END);
1452 qfileFilePos.push_back(size);
1455 //seach for filePos of each first name in the rfile and save in rfileFilePos
1459 m->openInputFile(rfilename, inR);
1461 //read column headers
1462 for (int i = 0; i < 16; i++) {
1463 if (!inR.eof()) { inR >> junk; }
1469 input = m->getline(inR);
1471 if (input.length() != 0) {
1473 istringstream nameStream(input);
1474 string sname = ""; nameStream >> sname;
1476 m->checkName(sname);
1478 map<string, int>::iterator it = firstSeqNamesReport.find(sname);
1480 if(it != firstSeqNamesReport.end()) { //this is the start of a new chunk
1481 unsigned long long pos = inR.tellg();
1482 rfileFilePos.push_back(pos - input.length() - 1);
1483 firstSeqNamesReport.erase(it);
1487 if (firstSeqNamesReport.size() == 0) { break; }
1492 if (firstSeqNamesReport.size() != 0) {
1493 for (map<string, int>::iterator it = firstSeqNamesReport.begin(); it != firstSeqNamesReport.end(); it++) {
1494 m->mothurOut(it->first + " is in your fasta file and not in your report file, aborting."); m->mothurOutEndLine();
1496 m->control_pressed = true;
1500 //get last file position of qfile
1502 unsigned long long sizeR;
1504 //get num bytes in file
1505 rFile = fopen (rfilename.c_str(),"rb");
1506 if (rFile==NULL) perror ("Error opening file");
1508 fseek (rFile, 0, SEEK_END);
1509 sizeR=ftell (rFile);
1513 rfileFilePos.push_back(sizeR);
1519 fastaFilePos.push_back(0); qfileFilePos.push_back(0);
1520 //get last file position of fastafile
1522 unsigned long long size;
1524 //get num bytes in file
1525 pFile = fopen (filename.c_str(),"rb");
1526 if (pFile==NULL) perror ("Error opening file");
1528 fseek (pFile, 0, SEEK_END);
1532 fastaFilePos.push_back(size);
1534 //get last file position of fastafile
1537 //get num bytes in file
1538 qFile = fopen (qfilename.c_str(),"rb");
1539 if (qFile==NULL) perror ("Error opening file");
1541 fseek (qFile, 0, SEEK_END);
1545 qfileFilePos.push_back(size);
1551 catch(exception& e) {
1552 m->errorOut(e, "SeqErrorCommand", "setLines");
1557 //***************************************************************************************************************