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 "filterseqscommand.h"
16 //**********************************************************************************************************************
17 vector<string> SeqErrorCommand::setParameters(){
19 CommandParameter pquery("fasta", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(pquery);
20 CommandParameter preference("reference", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(preference);
21 CommandParameter pqfile("qfile", "InputTypes", "", "", "none", "none", "QualReport",false,false); parameters.push_back(pqfile);
22 CommandParameter preport("report", "InputTypes", "", "", "none", "none", "QualReport",false,false); parameters.push_back(preport);
23 CommandParameter pname("name", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(pname);
24 CommandParameter pignorechimeras("ignorechimeras", "Boolean", "", "T", "", "", "",false,false); parameters.push_back(pignorechimeras);
25 CommandParameter pthreshold("threshold", "Number", "", "1.0", "", "", "",false,false); parameters.push_back(pthreshold);
26 CommandParameter pprocessors("processors", "Number", "", "1", "", "", "",false,false); parameters.push_back(pprocessors);
27 CommandParameter pfilter("filter", "Boolean", "", "T", "", "", "",false,false); parameters.push_back(pfilter);
28 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
29 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
31 vector<string> myArray;
32 for (int i = 0; i < parameters.size(); i++) { myArray.push_back(parameters[i].name); }
36 m->errorOut(e, "SeqErrorCommand", "setParameters");
40 //**********************************************************************************************************************
41 string SeqErrorCommand::getHelpString(){
43 string helpString = "";
44 helpString += "The seq.error command reads a query alignment file and a reference alignment file and creates .....\n";
45 helpString += "Example seq.error(...).\n";
46 helpString += "Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFasta).\n";
47 helpString += "For more details please check out the wiki http://www.mothur.org/wiki/seq.error .\n";
51 m->errorOut(e, "SeqErrorCommand", "getHelpString");
55 //**********************************************************************************************************************
56 SeqErrorCommand::SeqErrorCommand(){
58 abort = true; calledHelp = true;
59 vector<string> tempOutNames;
60 outputTypes["error.summary"] = tempOutNames;
61 outputTypes["error.seq"] = tempOutNames;
62 outputTypes["error.quality"] = tempOutNames;
63 outputTypes["error.qual.forward"] = tempOutNames;
64 outputTypes["error.qual.reverse"] = tempOutNames;
65 outputTypes["error.forward"] = tempOutNames;
66 outputTypes["error.reverse"] = tempOutNames;
67 outputTypes["error.count"] = tempOutNames;
68 outputTypes["error.matrix"] = tempOutNames;
71 m->errorOut(e, "SeqErrorCommand", "SeqErrorCommand");
75 //***************************************************************************************************************
77 SeqErrorCommand::SeqErrorCommand(string option) {
80 abort = false; calledHelp = false;
82 //allow user to run help
83 if(option == "help") { help(); abort = true; calledHelp = true; }
84 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
88 vector<string> myArray = setParameters();
90 OptionParser parser(option);
91 map<string,string> parameters = parser.getParameters();
93 ValidParameters validParameter;
94 map<string,string>::iterator it;
96 //check to make sure all parameters are valid for command
97 for (it = parameters.begin(); it != parameters.end(); it++) {
98 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
101 //initialize outputTypes
102 vector<string> tempOutNames;
103 outputTypes["error.summary"] = tempOutNames;
104 outputTypes["error.seq"] = tempOutNames;
105 outputTypes["error.quality"] = tempOutNames;
106 outputTypes["error.qual.forward"] = tempOutNames;
107 outputTypes["error.qual.reverse"] = tempOutNames;
108 outputTypes["error.forward"] = tempOutNames;
109 outputTypes["error.reverse"] = tempOutNames;
110 outputTypes["error.count"] = tempOutNames;
111 outputTypes["error.matrix"] = tempOutNames;
114 //if the user changes the input directory command factory will send this info to us in the output parameter
115 string inputDir = validParameter.validFile(parameters, "inputdir", false);
116 if (inputDir == "not found"){ inputDir = ""; }
119 it = parameters.find("fasta");
120 //user has given a template file
121 if(it != parameters.end()){
122 path = m->hasPath(it->second);
123 //if the user has not given a path then, add inputdir. else leave path alone.
124 if (path == "") { parameters["fasta"] = inputDir + it->second; }
127 it = parameters.find("reference");
128 //user has given a template file
129 if(it != parameters.end()){
130 path = m->hasPath(it->second);
131 //if the user has not given a path then, add inputdir. else leave path alone.
132 if (path == "") { parameters["reference"] = inputDir + it->second; }
135 it = parameters.find("name");
136 //user has given a names file
137 if(it != parameters.end()){
138 path = m->hasPath(it->second);
139 //if the user has not given a path then, add inputdir. else leave path alone.
140 if (path == "") { parameters["name"] = inputDir + it->second; }
143 it = parameters.find("qfile");
144 //user has given a quality score file
145 if(it != parameters.end()){
146 path = m->hasPath(it->second);
147 //if the user has not given a path then, add inputdir. else leave path alone.
148 if (path == "") { parameters["qfile"] = inputDir + it->second; }
151 it = parameters.find("report");
152 //user has given a alignment report file
153 if(it != parameters.end()){
154 path = m->hasPath(it->second);
155 //if the user has not given a path then, add inputdir. else leave path alone.
156 if (path == "") { parameters["report"] = inputDir + it->second; }
160 //check for required parameters
161 queryFileName = validParameter.validFile(parameters, "fasta", true);
162 if (queryFileName == "not found") {
163 queryFileName = m->getFastaFile();
164 if (queryFileName != "") { m->mothurOut("Using " + queryFileName + " as input file for the fasta parameter."); m->mothurOutEndLine(); }
165 else { m->mothurOut("You have no current fasta file and the fasta parameter is required."); m->mothurOutEndLine(); abort = true; }
167 else if (queryFileName == "not open") { abort = true; }
169 referenceFileName = validParameter.validFile(parameters, "reference", true);
170 if (referenceFileName == "not found") { m->mothurOut("reference is a required parameter for the seq.error command."); m->mothurOutEndLine(); abort = true; }
171 else if (referenceFileName == "not open") { abort = true; }
174 //check for optional parameters
175 namesFileName = validParameter.validFile(parameters, "name", true);
176 if(namesFileName == "not found"){ namesFileName = ""; }
177 else if (namesFileName == "not open") { namesFileName = ""; abort = true; }
179 qualFileName = validParameter.validFile(parameters, "qfile", true);
180 if(qualFileName == "not found"){ qualFileName = ""; }
181 else if (qualFileName == "not open") { qualFileName = ""; abort = true; }
183 reportFileName = validParameter.validFile(parameters, "report", true);
184 if(reportFileName == "not found"){ reportFileName = ""; }
185 else if (reportFileName == "not open") { reportFileName = ""; abort = true; }
187 if((reportFileName != "" && qualFileName == "") || (reportFileName == "" && qualFileName != "")){
188 m->mothurOut("if you use either a qual file or a report file, you have to have both.");
189 m->mothurOutEndLine();
193 outputDir = validParameter.validFile(parameters, "outputdir", false);
194 if (outputDir == "not found"){
196 outputDir += m->hasPath(queryFileName); //if user entered a file with a path then preserve it
199 //check for optional parameter and set defaults
200 // ...at some point should added some additional type checking...
201 temp = validParameter.validFile(parameters, "threshold", false); if (temp == "not found") { temp = "1.00"; }
202 convert(temp, threshold);
204 temp = validParameter.validFile(parameters, "ignorechimeras", false); if (temp == "not found") { temp = "T"; }
205 ignoreChimeras = m->isTrue(temp);
207 temp = validParameter.validFile(parameters, "filter", false); if (temp == "not found") { temp = "T"; }
208 filter = m->isTrue(temp);
210 temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); }
211 m->setProcessors(temp);
212 convert(temp, processors);
214 substitutionMatrix.resize(6);
215 for(int i=0;i<6;i++){ substitutionMatrix[i].resize(6,0); }
218 catch(exception& e) {
219 m->errorOut(e, "SeqErrorCommand", "SeqErrorCommand");
223 //***************************************************************************************************************
225 int SeqErrorCommand::execute(){
227 if (abort == true) { if (calledHelp) { return 0; } return 2; }
229 int start = time(NULL);
234 //run vertical filter on query and reference files.
236 string inputString = "fasta=" + queryFileName + "-" + referenceFileName;
237 m->mothurOut("/******************************************/"); m->mothurOutEndLine();
238 m->mothurOut("Running command: filter.seqs(" + inputString + ") to improve processing time."); m->mothurOutEndLine();
240 Command* filterCommand = new FilterSeqsCommand(inputString);
241 filterCommand->execute();
243 map<string, vector<string> > filenames = filterCommand->getOutputFiles();
245 delete filterCommand;
247 m->mothurOut("/******************************************/"); m->mothurOutEndLine();
249 queryFileName = filenames["fasta"][0];
250 referenceFileName = filenames["fasta"][1];
253 string errorSummaryFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.summary";
254 outputNames.push_back(errorSummaryFileName); outputTypes["error.summary"].push_back(errorSummaryFileName);
256 string errorSeqFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.seq";
257 outputNames.push_back(errorSeqFileName); outputTypes["error.seq"].push_back(errorSeqFileName);
259 string errorChimeraFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.chimera";
260 outputNames.push_back(errorChimeraFileName); outputTypes["error.chimera"].push_back(errorChimeraFileName);
262 getReferences(); //read in reference sequences - make sure there's no ambiguous bases
264 if(namesFileName != ""){ weights = getWeights(); }
266 vector<unsigned long int> fastaFilePos;
267 vector<unsigned long int> qFilePos;
268 vector<unsigned long int> reportFilePos;
270 setLines(queryFileName, qualFileName, reportFileName, fastaFilePos, qFilePos, reportFilePos);
272 if (m->control_pressed) { return 0; }
274 for (int i = 0; i < (fastaFilePos.size()-1); i++) {
275 lines.push_back(linePair(fastaFilePos[i], fastaFilePos[(i+1)]));
276 if (qualFileName != "") { qLines.push_back(linePair(qFilePos[i], qFilePos[(i+1)])); }
277 if (reportFileName != "") { rLines.push_back(linePair(reportFilePos[i], reportFilePos[(i+1)])); }
279 if(qualFileName == "") { qLines = lines; rLines = lines; } //fills with duds
282 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
284 numSeqs = driver(queryFileName, qualFileName, reportFileName, errorSummaryFileName, errorSeqFileName, errorChimeraFileName, lines[0], qLines[0], rLines[0]);
286 numSeqs = createProcesses(queryFileName, qualFileName, reportFileName, errorSummaryFileName, errorSeqFileName, errorChimeraFileName);
289 numSeqs = driver(queryFileName, qualFileName, reportFileName, errorSummaryFileName, errorSeqFileName, errorChimeraFileName, lines[0], qLines[0], rLines[0]);
292 if(qualFileName != "" && reportFileName != ""){
293 printErrorQuality(qScoreErrorMap);
294 printQualityFR(qualForwardMap, qualReverseMap);
297 printErrorFRFile(errorForward, errorReverse);
299 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
301 string errorCountFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.count";
302 ofstream errorCountFile;
303 m->openOutputFile(errorCountFileName, errorCountFile);
304 outputNames.push_back(errorCountFileName); outputTypes["error.count"].push_back(errorCountFileName);
305 m->mothurOut("Overall error rate:\t" + toString((double)(totalBases - totalMatches) / (double)totalBases) + "\n");
306 m->mothurOut("Errors\tSequences\n");
307 errorCountFile << "Errors\tSequences\n";
308 for(int i=0;i<misMatchCounts.size();i++){
309 m->mothurOut(toString(i) + '\t' + toString(misMatchCounts[i]) + '\n');
310 errorCountFile << i << '\t' << misMatchCounts[i] << endl;
312 errorCountFile.close();
314 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
318 string megAlignmentFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.ref-query";
319 ofstream megAlignmentFile;
320 m->openOutputFile(megAlignmentFileName, megAlignmentFile);
321 outputNames.push_back(megAlignmentFileName); outputTypes["error.ref-query"].push_back(megAlignmentFileName);
323 for(int i=0;i<numRefs;i++){
324 megAlignmentFile << referenceSeqs[i].getInlineSeq() << endl;
325 megAlignmentFile << megaAlignVector[i] << endl;
328 m->mothurOut("It took " + toString(time(NULL) - start) + " secs to check " + toString(numSeqs) + " sequences.");
329 m->mothurOutEndLine();
331 m->mothurOutEndLine();
332 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
333 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
334 m->mothurOutEndLine();
338 catch(exception& e) {
339 m->errorOut(e, "SeqErrorCommand", "execute");
343 //**********************************************************************************************************************
344 int SeqErrorCommand::createProcesses(string filename, string qFileName, string rFileName, string summaryFileName, string errorOutputFileName, string chimeraOutputFileName) {
348 map<char, vector<int> >::iterator it;
350 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
352 //loop through and create all the processes you want
353 while (process != processors) {
357 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
361 num = driver(filename, qFileName, rFileName, summaryFileName + toString(getpid()) + ".temp", errorOutputFileName+ toString(getpid()) + ".temp", chimeraOutputFileName + toString(getpid()) + ".temp", lines[process], qLines[process], rLines[process]);
363 //pass groupCounts to parent
365 string tempFile = filename + toString(getpid()) + ".info.temp";
366 m->openOutputFile(tempFile, out);
368 //output totalBases and totalMatches
369 out << num << '\t' << totalBases << '\t' << totalMatches << endl << endl;
371 //output substitutionMatrix
372 for(int i = 0; i < substitutionMatrix.size(); i++) {
373 for (int j = 0; j < substitutionMatrix[i].size(); j++) {
374 out << substitutionMatrix[i][j] << '\t';
380 //output qScoreErrorMap
381 for (it = qScoreErrorMap.begin(); it != qScoreErrorMap.end(); it++) {
382 vector<int> thisScoreErrorMap = it->second;
383 out << it->first << '\t';
384 for (int i = 0; i < thisScoreErrorMap.size(); i++) {
385 out << thisScoreErrorMap[i] << '\t';
391 //output qualForwardMap
392 for(int i = 0; i < qualForwardMap.size(); i++) {
393 for (int j = 0; j < qualForwardMap[i].size(); j++) {
394 out << qualForwardMap[i][j] << '\t';
400 //output qualReverseMap
401 for(int i = 0; i < qualReverseMap.size(); i++) {
402 for (int j = 0; j < qualReverseMap[i].size(); j++) {
403 out << qualReverseMap[i][j] << '\t';
410 //output errorForward
411 for (it = errorForward.begin(); it != errorForward.end(); it++) {
412 vector<int> thisErrorForward = it->second;
413 out << it->first << '\t';
414 for (int i = 0; i < thisErrorForward.size(); i++) {
415 out << thisErrorForward[i] << '\t';
421 //output errorReverse
422 for (it = errorReverse.begin(); it != errorReverse.end(); it++) {
423 vector<int> thisErrorReverse = it->second;
424 out << it->first << '\t';
425 for (int i = 0; i < thisErrorReverse.size(); i++) {
426 out << thisErrorReverse[i] << '\t';
432 //output misMatchCounts
433 out << misMatchCounts.size() << endl;
434 for (int j = 0; j < misMatchCounts.size(); j++) {
435 out << misMatchCounts[j] << '\t';
440 //output megaAlignVector
441 for (int j = 0; j < megaAlignVector.size(); j++) {
442 out << megaAlignVector[j] << endl;
450 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
451 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
457 num = driver(filename, qFileName, rFileName, summaryFileName, errorOutputFileName, chimeraOutputFileName, lines[0], qLines[0], rLines[0]);
459 //force parent to wait until all the processes are done
460 for (int i=0;i<processIDS.size();i++) {
461 int temp = processIDS[i];
466 for(int i=0;i<processIDS.size();i++){
468 m->mothurOut("Appending files from process " + toString(processIDS[i])); m->mothurOutEndLine();
470 m->appendFiles((summaryFileName + toString(processIDS[i]) + ".temp"), summaryFileName);
471 remove((summaryFileName + toString(processIDS[i]) + ".temp").c_str());
472 m->appendFiles((errorOutputFileName + toString(processIDS[i]) + ".temp"), errorOutputFileName);
473 remove((errorOutputFileName + toString(processIDS[i]) + ".temp").c_str());
474 m->appendFiles((chimeraOutputFileName + toString(processIDS[i]) + ".temp"), chimeraOutputFileName);
475 remove((chimeraOutputFileName + toString(processIDS[i]) + ".temp").c_str());
478 string tempFile = filename + toString(processIDS[i]) + ".info.temp";
479 m->openInputFile(tempFile, in);
481 //input totalBases and totalMatches
482 int tempBases, tempMatches, tempNumSeqs;
483 in >> tempNumSeqs >> tempBases >> tempMatches; m->gobble(in);
484 totalBases += tempBases; totalMatches += tempMatches; num += tempNumSeqs;
486 //input substitutionMatrix
488 for(int i = 0; i < substitutionMatrix.size(); i++) {
489 for (int j = 0; j < substitutionMatrix[i].size(); j++) {
490 in >> tempNum; substitutionMatrix[i][j] += tempNum;
496 //input qScoreErrorMap
498 for (int i = 0; i < qScoreErrorMap.size(); i++) {
500 vector<int> thisScoreErrorMap = qScoreErrorMap[first];
502 for (int i = 0; i < thisScoreErrorMap.size(); i++) {
503 in >> tempNum; thisScoreErrorMap[i] += tempNum;
505 qScoreErrorMap[first] = thisScoreErrorMap;
510 //input qualForwardMap
511 for(int i = 0; i < qualForwardMap.size(); i++) {
512 for (int j = 0; j < qualForwardMap[i].size(); j++) {
513 in >> tempNum; qualForwardMap[i][j] += tempNum;
519 //input qualReverseMap
520 for(int i = 0; i < qualReverseMap.size(); i++) {
521 for (int j = 0; j < qualReverseMap[i].size(); j++) {
522 in >> tempNum; qualReverseMap[i][j] += tempNum;
529 for (int i = 0; i < errorForward.size(); i++) {
531 vector<int> thisErrorForward = errorForward[first];
533 for (int i = 0; i < thisErrorForward.size(); i++) {
534 in >> tempNum; thisErrorForward[i] += tempNum;
536 errorForward[first] = thisErrorForward;
542 for (int i = 0; i < errorReverse.size(); i++) {
544 vector<int> thisErrorReverse = errorReverse[first];
546 for (int i = 0; i < thisErrorReverse.size(); i++) {
547 in >> tempNum; thisErrorReverse[i] += tempNum;
549 errorReverse[first] = thisErrorReverse;
554 //input misMatchCounts
556 in >> misMatchSize; m->gobble(in);
557 if (misMatchSize > misMatchCounts.size()) { misMatchCounts.resize(misMatchSize, 0); }
558 for (int j = 0; j < misMatchCounts.size(); j++) {
559 in >> tempNum; misMatchCounts[j] += tempNum;
563 //input megaAlignVector
565 for (int j = 0; j < megaAlignVector.size(); j++) {
566 thisLine = m->getline(in); m->gobble(in); megaAlignVector[j] += thisLine + '\n';
570 in.close(); remove(tempFile.c_str());
576 catch(exception& e) {
577 m->errorOut(e, "SeqErrorCommand", "createProcesses");
581 //**********************************************************************************************************************
582 int SeqErrorCommand::driver(string filename, string qFileName, string rFileName, string summaryFileName, string errorOutputFileName, string chimeraOutputFileName, linePair line, linePair qline, linePair rline) {
586 QualityScores quality;
588 misMatchCounts.resize(11, 0);
592 map<string, int>::iterator it;
593 qScoreErrorMap['m'].assign(41, 0);
594 qScoreErrorMap['s'].assign(41, 0);
595 qScoreErrorMap['i'].assign(41, 0);
596 qScoreErrorMap['a'].assign(41, 0);
598 errorForward['m'].assign(maxLength,0);
599 errorForward['s'].assign(maxLength,0);
600 errorForward['i'].assign(maxLength,0);
601 errorForward['d'].assign(maxLength,0);
602 errorForward['a'].assign(maxLength,0);
604 errorReverse['m'].assign(maxLength,0);
605 errorReverse['s'].assign(maxLength,0);
606 errorReverse['i'].assign(maxLength,0);
607 errorReverse['d'].assign(maxLength,0);
608 errorReverse['a'].assign(maxLength,0);
610 //open inputfiles and go to beginning place for this processor
612 m->openInputFile(filename, queryFile);
613 queryFile.seekg(line.start);
617 if(qFileName != "" && rFileName != ""){
618 m->openInputFile(qFileName, qualFile);
619 qualFile.seekg(qline.start);
622 if (rline.start == 0) { report = ReportFile(reportFile, rFileName); }
624 m->openInputFile(rFileName, reportFile);
625 reportFile.seekg(rline.start);
628 qualForwardMap.resize(maxLength);
629 qualReverseMap.resize(maxLength);
630 for(int i=0;i<maxLength;i++){
631 qualForwardMap[i].assign(41,0);
632 qualReverseMap[i].assign(41,0);
636 ofstream outChimeraReport;
637 m->openOutputFile(chimeraOutputFileName, outChimeraReport);
638 RefChimeraTest chimeraTest(referenceSeqs);
639 if (line.start == 0) { chimeraTest.printHeader(outChimeraReport); }
641 ofstream errorSummaryFile;
642 m->openOutputFile(summaryFileName, errorSummaryFile);
643 if (line.start == 0) { printErrorHeader(errorSummaryFile); }
645 ofstream errorSeqFile;
646 m->openOutputFile(errorOutputFileName, errorSeqFile);
648 megaAlignVector.resize(numRefs, "");
656 if (m->control_pressed) { queryFile.close(); if(qFileName != "" && rFileName != ""){ reportFile.close(); qualFile.close(); } outChimeraReport.close(); errorSummaryFile.close();errorSeqFile.close(); return 0; }
658 Sequence query(queryFile);
660 int numParentSeqs = chimeraTest.analyzeQuery(query.getName(), query.getAligned(), outChimeraReport);
661 int closestRefIndex = chimeraTest.getClosestRefIndex();
663 if(numParentSeqs > 1 && ignoreChimeras == 1) { ignoreSeq = 1; }
664 else { ignoreSeq = 0; }
666 Compare minCompare = getErrors(query, referenceSeqs[closestRefIndex]);
668 if(namesFileName != ""){
669 it = weights.find(query.getName());
670 minCompare.weight = it->second;
672 else{ minCompare.weight = 1; }
674 printErrorData(minCompare, numParentSeqs, errorSummaryFile, errorSeqFile);
678 for(int i=0;i<minCompare.sequence.length();i++){
679 char letter = minCompare.sequence[i];
681 errorForward[letter][i] += minCompare.weight;
682 errorReverse[letter][minCompare.total-i-1] += minCompare.weight;
686 if(qualFileName != "" && reportFileName != ""){
687 report = ReportFile(reportFile);
689 // int origLength = report.getQueryLength();
690 int startBase = report.getQueryStart();
691 int endBase = report.getQueryEnd();
693 quality = QualityScores(qualFile);
696 quality.updateQScoreErrorMap(qScoreErrorMap, minCompare.sequence, startBase, endBase, minCompare.weight);
697 quality.updateForwardMap(qualForwardMap, startBase, endBase, minCompare.weight);
698 quality.updateReverseMap(qualReverseMap, startBase, endBase, minCompare.weight);
702 if(minCompare.errorRate < threshold && !ignoreSeq){
703 totalBases += (minCompare.total * minCompare.weight);
704 totalMatches += minCompare.matches * minCompare.weight;
705 if(minCompare.mismatches > maxMismatch){
706 maxMismatch = minCompare.mismatches;
707 misMatchCounts.resize(maxMismatch + 1, 0);
709 misMatchCounts[minCompare.mismatches] += minCompare.weight;
712 megaAlignVector[closestRefIndex] += query.getInlineSeq() + '\n';
717 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
718 unsigned long int pos = queryFile.tellg();
719 if ((pos == -1) || (pos >= line.end)) { break; }
721 if (queryFile.eof()) { break; }
724 if(index % 100 == 0){ m->mothurOut(toString(index) + '\n'); }
727 if(qFileName != "" && rFileName != ""){ reportFile.close(); qualFile.close(); }
728 errorSummaryFile.close();
729 errorSeqFile.close();
732 if(index % 100 != 0){ m->mothurOut(toString(index) + '\n'); }
736 catch(exception& e) {
737 m->errorOut(e, "SeqErrorCommand", "driver");
741 //***************************************************************************************************************
743 void SeqErrorCommand::getReferences(){
746 ifstream referenceFile;
747 m->openInputFile(referenceFileName, referenceFile);
749 int numAmbigSeqs = 0;
752 int minEndPos = 100000;
754 while(referenceFile){
755 Sequence currentSeq(referenceFile);
756 int numAmbigs = currentSeq.getAmbigBases();
757 if(numAmbigs > 0){ numAmbigSeqs++; }
759 // int startPos = currentSeq.getStartPos();
760 // if(startPos > maxStartPos) { maxStartPos = startPos; }
762 // int endPos = currentSeq.getEndPos();
763 // if(endPos < minEndPos) { minEndPos = endPos; }
764 referenceSeqs.push_back(currentSeq);
766 m->gobble(referenceFile);
768 referenceFile.close();
769 numRefs = referenceSeqs.size();
772 for(int i=0;i<numRefs;i++){
773 referenceSeqs[i].padToPos(maxStartPos);
774 referenceSeqs[i].padFromPos(minEndPos);
777 if(numAmbigSeqs != 0){
778 m->mothurOut("Warning: " + toString(numAmbigSeqs) + " reference sequences have ambiguous bases, these bases will be ignored\n");
782 catch(exception& e) {
783 m->errorOut(e, "SeqErrorCommand", "getReferences");
788 //***************************************************************************************************************
790 Compare SeqErrorCommand::getErrors(Sequence query, Sequence reference){
792 if(query.getAlignLength() != reference.getAlignLength()){
793 m->mothurOut("Warning: " + toString(query.getName()) + " and " + toString(reference.getName()) + " are different lengths\n");
795 int alignLength = query.getAlignLength();
797 string q = query.getAligned();
798 string r = reference.getAligned();
803 for(int i=0;i<alignLength;i++){
804 if(r[i] != 'N' && q[i] != '.' && r[i] != '.' && (q[i] != '-' || r[i] != '-')){ // no missing data and no double gaps
808 if(r[i] == 'A'){ errors.AA++; errors.matches++; errors.sequence += 'm'; }
809 if(r[i] == 'T'){ errors.AT++; errors.sequence += 's'; }
810 if(r[i] == 'G'){ errors.AG++; errors.sequence += 's'; }
811 if(r[i] == 'C'){ errors.AC++; errors.sequence += 's'; }
812 if(r[i] == '-'){ errors.Ai++; errors.sequence += 'i'; }
814 else if(q[i] == 'T'){
815 if(r[i] == 'A'){ errors.TA++; errors.sequence += 's'; }
816 if(r[i] == 'T'){ errors.TT++; errors.matches++; errors.sequence += 'm'; }
817 if(r[i] == 'G'){ errors.TG++; errors.sequence += 's'; }
818 if(r[i] == 'C'){ errors.TC++; errors.sequence += 's'; }
819 if(r[i] == '-'){ errors.Ti++; errors.sequence += 'i'; }
821 else if(q[i] == 'G'){
822 if(r[i] == 'A'){ errors.GA++; errors.sequence += 's'; }
823 if(r[i] == 'T'){ errors.GT++; errors.sequence += 's'; }
824 if(r[i] == 'G'){ errors.GG++; errors.matches++; errors.sequence += 'm'; }
825 if(r[i] == 'C'){ errors.GC++; errors.sequence += 's'; }
826 if(r[i] == '-'){ errors.Gi++; errors.sequence += 'i'; }
828 else if(q[i] == 'C'){
829 if(r[i] == 'A'){ errors.CA++; errors.sequence += 's'; }
830 if(r[i] == 'T'){ errors.CT++; errors.sequence += 's'; }
831 if(r[i] == 'G'){ errors.CG++; errors.sequence += 's'; }
832 if(r[i] == 'C'){ errors.CC++; errors.matches++; errors.sequence += 'm'; }
833 if(r[i] == '-'){ errors.Ci++; errors.sequence += 'i'; }
835 else if(q[i] == 'N'){
836 if(r[i] == 'A'){ errors.NA++; errors.sequence += 'a'; }
837 if(r[i] == 'T'){ errors.NT++; errors.sequence += 'a'; }
838 if(r[i] == 'G'){ errors.NG++; errors.sequence += 'a'; }
839 if(r[i] == 'C'){ errors.NC++; errors.sequence += 'a'; }
840 if(r[i] == '-'){ errors.Ni++; errors.sequence += 'a'; }
842 else if(q[i] == '-' && r[i] != '-'){
843 if(r[i] == 'A'){ errors.dA++; errors.sequence += 'd'; }
844 if(r[i] == 'T'){ errors.dT++; errors.sequence += 'd'; }
845 if(r[i] == 'G'){ errors.dG++; errors.sequence += 'd'; }
846 if(r[i] == 'C'){ errors.dC++; errors.sequence += 'd'; }
851 else if(q[i] == '.' && r[i] != '.'){ // reference extends beyond query
852 if(started == 1){ break; }
854 else if(q[i] != '.' && r[i] == '.'){ // query extends beyond reference
855 if(started == 1){ break; }
857 else if(q[i] == '.' && r[i] == '.'){ // both are missing data
858 if(started == 1){ break; }
862 errors.mismatches = errors.total-errors.matches;
863 errors.errorRate = (double)(errors.total-errors.matches) / (double)errors.total;
864 errors.queryName = query.getName();
865 errors.refName = reference.getName();
869 catch(exception& e) {
870 m->errorOut(e, "SeqErrorCommand", "getErrors");
875 //***************************************************************************************************************
877 map<string, int> SeqErrorCommand::getWeights(){
879 m->openInputFile(namesFileName, nameFile);
882 string redundantSeqs;
883 map<string, int> nameCountMap;
886 nameFile >> seqName >> redundantSeqs;
887 nameCountMap[seqName] = m->getNumNames(redundantSeqs);
894 //***************************************************************************************************************
896 void SeqErrorCommand::printErrorHeader(ofstream& errorSummaryFile){
898 errorSummaryFile << "query\treference\tweight\t";
899 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";
900 errorSummaryFile << "insertions\tdeletions\tsubstitutions\tambig\tmatches\tmismatches\ttotal\terror\tnumparents\n";
902 errorSummaryFile << setprecision(6);
903 errorSummaryFile.setf(ios::fixed);
905 catch(exception& e) {
906 m->errorOut(e, "SeqErrorCommand", "printErrorHeader");
911 //***************************************************************************************************************
913 void SeqErrorCommand::printErrorData(Compare error, int numParentSeqs, ofstream& errorSummaryFile, ofstream& errorSeqFile){
916 errorSummaryFile << error.queryName << '\t' << error.refName << '\t' << error.weight << '\t';
917 errorSummaryFile << error.AA << '\t' << error.AT << '\t' << error.AG << '\t' << error.AC << '\t';
918 errorSummaryFile << error.TA << '\t' << error.TT << '\t' << error.TG << '\t' << error.TC << '\t';
919 errorSummaryFile << error.GA << '\t' << error.GT << '\t' << error.GG << '\t' << error.GC << '\t';
920 errorSummaryFile << error.CA << '\t' << error.CT << '\t' << error.CG << '\t' << error.CC << '\t';
921 errorSummaryFile << error.NA << '\t' << error.NT << '\t' << error.NG << '\t' << error.NC << '\t';
922 errorSummaryFile << error.Ai << '\t' << error.Ti << '\t' << error.Gi << '\t' << error.Ci << '\t' << error.Ni << '\t';
923 errorSummaryFile << error.dA << '\t' << error.dT << '\t' << error.dG << '\t' << error.dC << '\t';
925 errorSummaryFile << error.Ai + error.Ti + error.Gi + error.Ci << '\t'; //insertions
926 errorSummaryFile << error.dA + error.dT + error.dG + error.dC << '\t'; //deletions
927 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
928 errorSummaryFile << error.NA + error.NT + error.NG + error.NC + error.Ni << '\t'; //ambiguities
929 errorSummaryFile << error.matches << '\t' << error.mismatches << '\t' << error.total << '\t' << error.errorRate << '\t' << numParentSeqs << endl;
931 errorSeqFile << '>' << error.queryName << "\tref:" << error.refName << '\n' << error.sequence << endl;
933 int a=0; int t=1; int g=2; int c=3;
936 if(numParentSeqs == 1 || ignoreChimeras == 0){
937 substitutionMatrix[a][a] += error.weight * error.AA;
938 substitutionMatrix[a][t] += error.weight * error.TA;
939 substitutionMatrix[a][g] += error.weight * error.GA;
940 substitutionMatrix[a][c] += error.weight * error.CA;
941 substitutionMatrix[a][gap] += error.weight * error.dA;
942 substitutionMatrix[a][n] += error.weight * error.NA;
944 substitutionMatrix[t][a] += error.weight * error.AT;
945 substitutionMatrix[t][t] += error.weight * error.TT;
946 substitutionMatrix[t][g] += error.weight * error.GT;
947 substitutionMatrix[t][c] += error.weight * error.CT;
948 substitutionMatrix[t][gap] += error.weight * error.dT;
949 substitutionMatrix[t][n] += error.weight * error.NT;
951 substitutionMatrix[g][a] += error.weight * error.AG;
952 substitutionMatrix[g][t] += error.weight * error.TG;
953 substitutionMatrix[g][g] += error.weight * error.GG;
954 substitutionMatrix[g][c] += error.weight * error.CG;
955 substitutionMatrix[g][gap] += error.weight * error.dG;
956 substitutionMatrix[g][n] += error.weight * error.NG;
958 substitutionMatrix[c][a] += error.weight * error.AC;
959 substitutionMatrix[c][t] += error.weight * error.TC;
960 substitutionMatrix[c][g] += error.weight * error.GC;
961 substitutionMatrix[c][c] += error.weight * error.CC;
962 substitutionMatrix[c][gap] += error.weight * error.dC;
963 substitutionMatrix[c][n] += error.weight * error.NC;
965 substitutionMatrix[gap][a] += error.weight * error.Ai;
966 substitutionMatrix[gap][t] += error.weight * error.Ti;
967 substitutionMatrix[gap][g] += error.weight * error.Gi;
968 substitutionMatrix[gap][c] += error.weight * error.Ci;
969 substitutionMatrix[gap][n] += error.weight * error.Ni;
972 catch(exception& e) {
973 m->errorOut(e, "SeqErrorCommand", "printErrorData");
978 //***************************************************************************************************************
980 void SeqErrorCommand::printSubMatrix(){
982 string subMatrixFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.matrix";
983 ofstream subMatrixFile;
984 m->openOutputFile(subMatrixFileName, subMatrixFile);
985 outputNames.push_back(subMatrixFileName); outputTypes["error.matrix"].push_back(subMatrixFileName);
986 vector<string> bases(6);
993 vector<int> refSums(5,1);
995 for(int i=0;i<5;i++){
996 subMatrixFile << "\tr" << bases[i];
998 for(int j=0;j<6;j++){
999 refSums[i] += substitutionMatrix[i][j];
1002 subMatrixFile << endl;
1004 for(int i=0;i<6;i++){
1005 subMatrixFile << 'q' << bases[i];
1006 for(int j=0;j<5;j++){
1007 subMatrixFile << '\t' << substitutionMatrix[j][i];
1009 subMatrixFile << endl;
1012 subMatrixFile << "total";
1013 for(int i=0;i<5;i++){
1014 subMatrixFile << '\t' << refSums[i];
1016 subMatrixFile << endl;
1017 subMatrixFile.close();
1019 catch(exception& e) {
1020 m->errorOut(e, "SeqErrorCommand", "printSubMatrix");
1024 //***************************************************************************************************************
1026 void SeqErrorCommand::printErrorFRFile(map<char, vector<int> > errorForward, map<char, vector<int> > errorReverse){
1028 string errorForwardFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.seq.forward";
1029 ofstream errorForwardFile;
1030 m->openOutputFile(errorForwardFileName, errorForwardFile);
1031 outputNames.push_back(errorForwardFileName); outputTypes["error.forward"].push_back(errorForwardFileName);
1033 errorForwardFile << "position\ttotalseqs\tmatch\tsubstitution\tinsertion\tdeletion\tambiguous" << endl;
1034 for(int i=0;i<maxLength;i++){
1035 float match = (float)errorForward['m'][i];
1036 float subst = (float)errorForward['s'][i];
1037 float insert = (float)errorForward['i'][i];
1038 float del = (float)errorForward['d'][i];
1039 float amb = (float)errorForward['a'][i];
1040 float total = match + subst + insert + del + amb;
1041 if(total == 0){ break; }
1042 errorForwardFile << i+1 << '\t' << total << '\t' << match/total << '\t' << subst/total << '\t' << insert/total << '\t' << del/total << '\t' << amb/total << endl;
1044 errorForwardFile.close();
1046 string errorReverseFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.seq.reverse";
1047 ofstream errorReverseFile;
1048 m->openOutputFile(errorReverseFileName, errorReverseFile);
1049 outputNames.push_back(errorReverseFileName); outputTypes["error.reverse"].push_back(errorReverseFileName);
1051 errorReverseFile << "position\ttotalseqs\tmatch\tsubstitution\tinsertion\tdeletion\tambiguous" << endl;
1052 for(int i=0;i<maxLength;i++){
1053 float match = (float)errorReverse['m'][i];
1054 float subst = (float)errorReverse['s'][i];
1055 float insert = (float)errorReverse['i'][i];
1056 float del = (float)errorReverse['d'][i];
1057 float amb = (float)errorReverse['a'][i];
1058 float total = match + subst + insert + del + amb;
1059 if(total == 0){ break; }
1060 errorReverseFile << i+1 << '\t' << total << '\t' << match/total << '\t' << subst/total << '\t' << insert/total << '\t' << del/total << '\t' << amb/total << endl;
1062 errorReverseFile.close();
1064 catch(exception& e) {
1065 m->errorOut(e, "SeqErrorCommand", "printErrorFRFile");
1070 //***************************************************************************************************************
1072 void SeqErrorCommand::printErrorQuality(map<char, vector<int> > qScoreErrorMap){
1075 string errorQualityFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.quality";
1076 ofstream errorQualityFile;
1077 m->openOutputFile(errorQualityFileName, errorQualityFile);
1078 outputNames.push_back(errorQualityFileName); outputTypes["error.quality"].push_back(errorQualityFileName);
1080 errorQualityFile << "qscore\tmatches\tsubstitutions\tinsertions\tambiguous" << endl;
1081 for(int i=0;i<41;i++){
1082 errorQualityFile << i << '\t' << qScoreErrorMap['m'][i] << '\t' << qScoreErrorMap['s'][i] << '\t' << qScoreErrorMap['i'][i] << '\t'<< qScoreErrorMap['a'][i] << endl;
1084 errorQualityFile.close();
1086 catch(exception& e) {
1087 m->errorOut(e, "SeqErrorCommand", "printErrorFRFile");
1093 //***************************************************************************************************************
1095 void SeqErrorCommand::printQualityFR(vector<vector<int> > qualForwardMap, vector<vector<int> > qualReverseMap){
1099 int numColumns = qualForwardMap[0].size();
1101 for(int i=0;i<qualForwardMap.size();i++){
1102 for(int j=0;j<numColumns;j++){
1103 if(qualForwardMap[i][j] != 0){
1104 if(numRows < i) { numRows = i+20; }
1109 string qualityForwardFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.qual.forward";
1110 ofstream qualityForwardFile;
1111 m->openOutputFile(qualityForwardFileName, qualityForwardFile);
1112 outputNames.push_back(qualityForwardFileName); outputTypes["error.qual.forward"].push_back(qualityForwardFileName);
1114 for(int i=0;i<numColumns;i++){ qualityForwardFile << '\t' << i; } qualityForwardFile << endl;
1116 for(int i=0;i<numRows;i++){
1117 qualityForwardFile << i+1;
1118 for(int j=0;j<numColumns;j++){
1119 qualityForwardFile << '\t' << qualForwardMap[i][j];
1122 qualityForwardFile << endl;
1124 qualityForwardFile.close();
1127 string qualityReverseFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.qual.reverse";
1128 ofstream qualityReverseFile;
1129 m->openOutputFile(qualityReverseFileName, qualityReverseFile);
1130 outputNames.push_back(qualityReverseFileName); outputTypes["error.qual.reverse"].push_back(qualityReverseFileName);
1132 for(int i=0;i<numColumns;i++){ qualityReverseFile << '\t' << i; } qualityReverseFile << endl;
1133 for(int i=0;i<numRows;i++){
1135 qualityReverseFile << i+1;
1136 for(int j=0;j<numColumns;j++){
1137 qualityReverseFile << '\t' << qualReverseMap[i][j];
1139 qualityReverseFile << endl;
1141 qualityReverseFile.close();
1143 catch(exception& e) {
1144 m->errorOut(e, "SeqErrorCommand", "printErrorFRFile");
1149 /**************************************************************************************************/
1151 int SeqErrorCommand::setLines(string filename, string qfilename, string rfilename, vector<unsigned long int>& fastaFilePos, vector<unsigned long int>& qfileFilePos, vector<unsigned long int>& rfileFilePos) {
1153 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
1154 //set file positions for fasta file
1155 fastaFilePos = m->divideFile(filename, processors);
1157 if (qfilename == "") { return processors; }
1159 //get name of first sequence in each chunk
1160 map<string, int> firstSeqNames;
1161 for (int i = 0; i < (fastaFilePos.size()-1); i++) {
1163 m->openInputFile(filename, in);
1164 in.seekg(fastaFilePos[i]);
1167 firstSeqNames[temp.getName()] = i;
1172 //make copy to use below
1173 map<string, int> firstSeqNamesReport = firstSeqNames;
1175 //seach for filePos of each first name in the qfile and save in qfileFilePos
1177 m->openInputFile(qfilename, inQual);
1180 while(!inQual.eof()){
1181 input = m->getline(inQual);
1183 if (input.length() != 0) {
1184 if(input[0] == '>'){ //this is a sequence name line
1185 istringstream nameStream(input);
1187 string sname = ""; nameStream >> sname;
1188 sname = sname.substr(1);
1190 map<string, int>::iterator it = firstSeqNames.find(sname);
1192 if(it != firstSeqNames.end()) { //this is the start of a new chunk
1193 unsigned long int pos = inQual.tellg();
1194 qfileFilePos.push_back(pos - input.length() - 1);
1195 firstSeqNames.erase(it);
1200 if (firstSeqNames.size() == 0) { break; }
1204 if (firstSeqNames.size() != 0) {
1205 for (map<string, int>::iterator it = firstSeqNames.begin(); it != firstSeqNames.end(); it++) {
1206 m->mothurOut(it->first + " is in your fasta file and not in your quality file, aborting."); m->mothurOutEndLine();
1208 m->control_pressed = true;
1212 //get last file position of qfile
1214 unsigned long int size;
1216 //get num bytes in file
1217 pFile = fopen (qfilename.c_str(),"rb");
1218 if (pFile==NULL) perror ("Error opening file");
1220 fseek (pFile, 0, SEEK_END);
1225 qfileFilePos.push_back(size);
1227 //seach for filePos of each first name in the rfile and save in rfileFilePos
1230 m->openInputFile(rfilename, inR);
1232 //read column headers
1233 for (int i = 0; i < 16; i++) {
1234 if (!inR.eof()) { inR >> junk; }
1240 if (m->control_pressed) { inR.close(); return processors; }
1242 input = m->getline(inR);
1244 if (input.length() != 0) {
1246 istringstream nameStream(input);
1247 string sname = ""; nameStream >> sname;
1249 map<string, int>::iterator it = firstSeqNamesReport.find(sname);
1251 if(it != firstSeqNamesReport.end()) { //this is the start of a new chunk
1252 unsigned long int pos = inR.tellg();
1253 rfileFilePos.push_back(pos - input.length() - 1);
1254 firstSeqNamesReport.erase(it);
1258 if (firstSeqNamesReport.size() == 0) { break; }
1263 if (firstSeqNamesReport.size() != 0) {
1264 for (map<string, int>::iterator it = firstSeqNamesReport.begin(); it != firstSeqNamesReport.end(); it++) {
1265 m->mothurOut(it->first + " is in your fasta file and not in your report file, aborting."); m->mothurOutEndLine();
1267 m->control_pressed = true;
1271 //get last file position of qfile
1273 unsigned long int sizeR;
1275 //get num bytes in file
1276 rFile = fopen (rfilename.c_str(),"rb");
1277 if (rFile==NULL) perror ("Error opening file");
1279 fseek (rFile, 0, SEEK_END);
1280 sizeR=ftell (rFile);
1284 rfileFilePos.push_back(sizeR);
1290 fastaFilePos.push_back(0); qfileFilePos.push_back(0);
1291 //get last file position of fastafile
1293 unsigned long int size;
1295 //get num bytes in file
1296 pFile = fopen (filename.c_str(),"rb");
1297 if (pFile==NULL) perror ("Error opening file");
1299 fseek (pFile, 0, SEEK_END);
1303 fastaFilePos.push_back(size);
1305 //get last file position of fastafile
1308 //get num bytes in file
1309 qFile = fopen (qfilename.c_str(),"rb");
1310 if (qFile==NULL) perror ("Error opening file");
1312 fseek (qFile, 0, SEEK_END);
1316 qfileFilePos.push_back(size);
1322 catch(exception& e) {
1323 m->errorOut(e, "SeqErrorCommand", "setLines");
1327 //***************************************************************************************************************