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"
14 //**********************************************************************************************************************
15 vector<string> SeqErrorCommand::getValidParameters(){
17 string Array[] = {"query", "reference", "name", "qfile", "report", "threshold", "inputdir", "outputdir"};
18 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
22 m->errorOut(e, "SeqErrorCommand", "getValidParameters");
26 //**********************************************************************************************************************
27 SeqErrorCommand::SeqErrorCommand(){
30 //initialize outputTypes
31 vector<string> tempOutNames;
32 outputTypes["error"] = tempOutNames;
33 outputTypes["count"] = tempOutNames;
36 m->errorOut(e, "SeqErrorCommand", "SeqErrorCommand");
40 //**********************************************************************************************************************
41 vector<string> SeqErrorCommand::getRequiredParameters(){
43 string Array[] = {"query","reference"};
44 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
48 m->errorOut(e, "SeqErrorCommand", "getRequiredParameters");
52 //**********************************************************************************************************************
53 vector<string> SeqErrorCommand::getRequiredFiles(){
55 vector<string> myArray;
59 m->errorOut(e, "SeqErrorCommand", "getRequiredFiles");
63 //***************************************************************************************************************
65 SeqErrorCommand::SeqErrorCommand(string option) {
70 //allow user to run help
71 if(option == "help") { help(); abort = true; }
76 //valid paramters for this command
77 string AlignArray[] = {"query", "reference", "name", "qfile", "report", "threshold", "inputdir", "outputdir"};
79 vector<string> myArray (AlignArray, AlignArray+(sizeof(AlignArray)/sizeof(string)));
81 OptionParser parser(option);
82 map<string,string> parameters = parser.getParameters();
84 ValidParameters validParameter;
85 map<string,string>::iterator it;
87 //check to make sure all parameters are valid for command
88 for (it = parameters.begin(); it != parameters.end(); it++) {
89 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
92 //initialize outputTypes
93 vector<string> tempOutNames;
94 outputTypes["error"] = tempOutNames;
95 outputTypes["count"] = tempOutNames;
97 //if the user changes the input directory command factory will send this info to us in the output parameter
98 string inputDir = validParameter.validFile(parameters, "inputdir", false);
99 if (inputDir == "not found"){ inputDir = ""; }
102 it = parameters.find("query");
103 //user has given a template file
104 if(it != parameters.end()){
105 path = m->hasPath(it->second);
106 //if the user has not given a path then, add inputdir. else leave path alone.
107 if (path == "") { parameters["query"] = inputDir + it->second; }
110 it = parameters.find("reference");
111 //user has given a template file
112 if(it != parameters.end()){
113 path = m->hasPath(it->second);
114 //if the user has not given a path then, add inputdir. else leave path alone.
115 if (path == "") { parameters["reference"] = inputDir + it->second; }
118 it = parameters.find("name");
119 //user has given a names file
120 if(it != parameters.end()){
121 path = m->hasPath(it->second);
122 //if the user has not given a path then, add inputdir. else leave path alone.
123 if (path == "") { parameters["name"] = inputDir + it->second; }
126 it = parameters.find("qfile");
127 //user has given a quality score file
128 if(it != parameters.end()){
129 path = m->hasPath(it->second);
130 //if the user has not given a path then, add inputdir. else leave path alone.
131 if (path == "") { parameters["qfile"] = inputDir + it->second; }
134 it = parameters.find("report");
135 //user has given a alignment report file
136 if(it != parameters.end()){
137 path = m->hasPath(it->second);
138 //if the user has not given a path then, add inputdir. else leave path alone.
139 if (path == "") { parameters["report"] = inputDir + it->second; }
143 //check for required parameters
144 queryFileName = validParameter.validFile(parameters, "query", true);
145 if (queryFileName == "not found") { m->mothurOut("query is a required parameter for the seq.error command."); m->mothurOutEndLine(); abort = true; }
146 else if (queryFileName == "not open") { abort = true; }
148 referenceFileName = validParameter.validFile(parameters, "reference", true);
149 if (referenceFileName == "not found") { m->mothurOut("reference is a required parameter for the seq.error command."); m->mothurOutEndLine(); abort = true; }
150 else if (referenceFileName == "not open") { abort = true; }
153 //check for optional parameters
154 namesFileName = validParameter.validFile(parameters, "name", true);
155 if(namesFileName == "not found"){ namesFileName = ""; }
157 qualFileName = validParameter.validFile(parameters, "qfile", true);
158 if(qualFileName == "not found"){ qualFileName = ""; }
160 reportFileName = validParameter.validFile(parameters, "report", true);
161 if(reportFileName == "not found"){ reportFileName = ""; }
163 if((reportFileName != "" && qualFileName == "") || (reportFileName == "" && qualFileName != "")){
164 m->mothurOut("if you use either a qual file or a report file, you have to have both.");
165 m->mothurOutEndLine();
169 outputDir = validParameter.validFile(parameters, "outputdir", false);
170 if (outputDir == "not found"){
172 outputDir += m->hasPath(queryFileName); //if user entered a file with a path then preserve it
175 //check for optional parameter and set defaults
176 // ...at some point should added some additional type checking...
177 temp = validParameter.validFile(parameters, "threshold", false); if (temp == "not found") { temp = "1.00"; }
178 convert(temp, threshold);
180 errorSummaryFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.summary";
181 m->openOutputFile(errorSummaryFileName, errorSummaryFile);
182 outputNames.push_back(errorSummaryFileName); outputTypes["error.summary"].push_back(errorSummaryFileName);
185 errorSeqFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.seq";
186 m->openOutputFile(errorSeqFileName, errorSeqFile);
187 outputNames.push_back(errorSeqFileName); outputTypes["error.seq"].push_back(errorSeqFileName);
189 substitutionMatrix.resize(6);
190 for(int i=0;i<6;i++){ substitutionMatrix[i].assign(6,0); }
193 catch(exception& e) {
194 m->errorOut(e, "SeqErrorCommand", "SeqErrorCommand");
199 //**********************************************************************************************************************
201 void SeqErrorCommand::help(){
203 m->mothurOut("The seq.error command reads a query alignment file and a reference alignment file and creates .....\n");
207 m->mothurOut("Example seq.error(...).\n");
208 m->mothurOut("Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFasta).\n");
209 m->mothurOut("For more details please check out the wiki http://www.mothur.org/wiki/seq.error .\n\n");
211 catch(exception& e) {
212 m->errorOut(e, "SeqErrorCommand", "help");
217 //***************************************************************************************************************
219 SeqErrorCommand::~SeqErrorCommand(){
220 errorSummaryFile.close();
221 errorSeqFile.close();
224 //***************************************************************************************************************
226 int SeqErrorCommand::execute(){
228 if (abort == true) { return 0; }
230 getReferences(); //read in reference sequences - make sure there's no ambiguous bases
232 map<string, int> weights;
233 if(namesFileName != ""){ weights = getWeights(); }
236 m->openInputFile(queryFileName, queryFile);
242 QualityScores quality;
243 vector<vector<int> > qualForwardMap;
244 vector<vector<int> > qualReverseMap;
246 if(qualFileName != "" && reportFileName != ""){
247 m->openInputFile(qualFileName, qualFile);
248 report = ReportFile(reportFile, reportFileName);
250 qualForwardMap.resize(1000);
251 qualReverseMap.resize(1000);
252 for(int i=0;i<1000;i++){
253 qualForwardMap[i].assign(100,0);
254 qualReverseMap[i].assign(100,0);
259 int totalMatches = 0;
261 vector<int> misMatchCounts(11, 0);
265 map<string, int>::iterator it;
266 map<char, vector<int> > qScoreErrorMap;
267 qScoreErrorMap['m'].assign(41, 0);
268 qScoreErrorMap['s'].assign(41, 0);
269 qScoreErrorMap['i'].assign(41, 0);
270 qScoreErrorMap['a'].assign(41, 0);
274 map<char, vector<int> > errorForward;
275 errorForward['m'].assign(1000,0);
276 errorForward['s'].assign(1000,0);
277 errorForward['i'].assign(1000,0);
278 errorForward['d'].assign(1000,0);
279 errorForward['a'].assign(1000,0);
281 map<char, vector<int> > errorReverse;
282 errorReverse['m'].assign(1000,0);
283 errorReverse['s'].assign(1000,0);
284 errorReverse['i'].assign(1000,0);
285 errorReverse['d'].assign(1000,0);
286 errorReverse['a'].assign(1000,0);
292 Sequence query(queryFile);
294 for(int i=0;i<numRefs;i++){
295 Compare currCompare = getErrors(query, referenceSeqs[i]);
297 if(currCompare.errorRate < minCompare.errorRate){
298 minCompare = currCompare;
302 if(namesFileName != ""){
303 it = weights.find(query.getName());
304 minCompare.weight = it->second;
306 else { minCompare.weight = 1; }
308 printErrorData(minCompare);
310 for(int i=0;i<minCompare.total;i++){
311 char letter = minCompare.sequence[i];
312 errorForward[letter][i] += minCompare.weight;
313 errorReverse[letter][minCompare.total-i-1] += minCompare.weight;
316 if(qualFileName != "" && reportFileName != ""){
317 report = ReportFile(reportFile);
319 int origLength = report.getQueryLength();
320 int startBase = report.getQueryStart();
321 int endBase = report.getQueryEnd();
323 quality = QualityScores(qualFile, origLength);
324 quality.updateQScoreErrorMap(qScoreErrorMap, minCompare.sequence, startBase, endBase, minCompare.weight);
325 quality.updateForwardMap(qualForwardMap, startBase, endBase, minCompare.weight);
326 quality.updateReverseMap(qualReverseMap, startBase, endBase, minCompare.weight);
329 if(minCompare.errorRate < threshold){
330 totalBases += (minCompare.total * minCompare.weight);
331 totalMatches += minCompare.matches * minCompare.weight;
332 if(minCompare.mismatches > maxMismatch){
333 maxMismatch = minCompare.mismatches;
334 misMatchCounts.resize(maxMismatch + 1, 0);
336 misMatchCounts[minCompare.mismatches] += minCompare.weight;
346 if(qualFileName != "" && reportFileName != ""){
347 string errorQualityFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.quality";
348 ofstream errorQualityFile;
349 m->openOutputFile(errorQualityFileName, errorQualityFile);
350 outputNames.push_back(errorQualityFileName); outputTypes["error.quality"].push_back(errorQualityFileName);
352 errorQualityFile << "qscore\tmatches\tsubstitutions\tinsertions\tambiguous" << endl;
353 for(int i=0;i<41;i++){
354 errorQualityFile << i << '\t' << qScoreErrorMap['m'][i] << '\t' << qScoreErrorMap['s'][i] << '\t' << qScoreErrorMap['i'][i] << '\t'<< qScoreErrorMap['a'][i] << endl;
356 errorQualityFile.close();
363 for(int i=0;i<qualForwardMap.size();i++){
364 for(int j=0;j<qualForwardMap[i].size();j++){
365 if(qualForwardMap[i][j] != 0){
366 if(lastRow < i) { lastRow = i+2; }
367 if(lastColumn < j) { lastColumn = j+2; }
374 string qualityForwardFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.qual.forward";
375 ofstream qualityForwardFile;
376 m->openOutputFile(qualityForwardFileName, qualityForwardFile);
377 outputNames.push_back(errorQualityFileName); outputTypes["error.qual.forward"].push_back(qualityForwardFileName);
379 for(int i=0;i<lastColumn;i++){ qualityForwardFile << '\t' << i; } qualityForwardFile << endl;
380 for(int i=0;i<lastRow;i++){
381 qualityForwardFile << i+1;
382 for(int j=0;j<lastColumn;j++){
383 qualityForwardFile << '\t' << qualForwardMap[i][j];
385 qualityForwardFile << endl;
387 qualityForwardFile.close();
390 string qualityReverseFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.qual.reverse";
391 ofstream qualityReverseFile;
392 m->openOutputFile(qualityReverseFileName, qualityReverseFile);
393 outputNames.push_back(errorQualityFileName); outputTypes["error.qual.reverse"].push_back(qualityReverseFileName);
395 for(int i=0;i<lastColumn;i++){ qualityReverseFile << '\t' << i; } qualityReverseFile << endl;
396 for(int i=0;i<lastRow;i++){
397 qualityReverseFile << i+1;
398 for(int j=0;j<lastColumn;j++){
399 qualityReverseFile << '\t' << qualReverseMap[i][j];
401 qualityReverseFile << endl;
407 string errorForwardFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.seq.forward";
408 ofstream errorForwardFile;
409 m->openOutputFile(errorForwardFileName, errorForwardFile);
410 outputNames.push_back(errorForwardFileName); outputTypes["error.forward"].push_back(errorForwardFileName);
412 errorForwardFile << "position\ttotalseqs\tmatch\tsubstitution\tinsertion\tdeletion\tambiguous" << endl;
413 for(int i=0;i<1000;i++){
414 float match = (float)errorForward['m'][i];
415 float subst = (float)errorForward['s'][i];
416 float insert = (float)errorForward['i'][i];
417 float del = (float)errorForward['d'][i];
418 float amb = (float)errorForward['a'][i];
419 float total = match + subst + insert + del + amb;
420 if(total == 0){ break; }
421 errorForwardFile << i+1 << '\t' << total << '\t' << match/total << '\t' << subst/total << '\t' << insert/total << '\t' << del/total << '\t' << amb/total << endl;
423 errorForwardFile.close();
426 string errorReverseFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.seq.reverse";
427 ofstream errorReverseFile;
428 m->openOutputFile(errorReverseFileName, errorReverseFile);
429 outputNames.push_back(errorReverseFileName); outputTypes["error.reverse"].push_back(errorReverseFileName);
431 errorReverseFile << "position\ttotalseqs\tmatch\tsubstitution\tinsertion\tdeletion\tambiguous" << endl;
432 for(int i=0;i<1000;i++){
433 float match = (float)errorReverse['m'][i];
434 float subst = (float)errorReverse['s'][i];
435 float insert = (float)errorReverse['i'][i];
436 float del = (float)errorReverse['d'][i];
437 float amb = (float)errorReverse['a'][i];
438 float total = match + subst + insert + del + amb;
439 if(total == 0){ break; }
440 errorReverseFile << i+1 << '\t' << total << '\t' << match/total << '\t' << subst/total << '\t' << insert/total << '\t' << del/total << '\t' << amb/total << endl;
442 errorReverseFile.close();
446 string errorCountFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.count";
447 ofstream errorCountFile;
448 m->openOutputFile(errorCountFileName, errorCountFile);
449 outputNames.push_back(errorCountFileName); outputTypes["error.count"].push_back(errorCountFileName);
450 m->mothurOut("Overall error rate:\t" + toString((double)(totalBases - totalMatches) / (double)totalBases) + "\n\n");
451 m->mothurOut("Errors\tSequences\n");
452 errorCountFile << "Errors\tSequences\n";
453 for(int i=0;i<misMatchCounts.size();i++){
454 m->mothurOut(toString(i) + '\t' + toString(misMatchCounts[i]) + '\n');
455 errorCountFile << i << '\t' << misMatchCounts[i] << endl;
457 errorCountFile.close();
463 string subMatrixFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.matrix";
464 ofstream subMatrixFile;
465 m->openOutputFile(subMatrixFileName, subMatrixFile);
466 outputNames.push_back(subMatrixFileName); outputTypes["error.matrix"].push_back(subMatrixFileName);
467 vector<string> bases(6);
474 vector<int> refSums(5,1);
476 for(int i=0;i<5;i++){
477 subMatrixFile << "\tr" << bases[i];
479 for(int j=0;j<6;j++){
480 refSums[i] += substitutionMatrix[i][j];
484 subMatrixFile << endl;
486 for(int i=0;i<6;i++){
487 subMatrixFile << 'q' << bases[i];
488 for(int j=0;j<5;j++){
489 subMatrixFile << '\t' << substitutionMatrix[j][i];
491 subMatrixFile << endl;
493 subMatrixFile << "total";
494 for(int i=0;i<5;i++){
495 subMatrixFile << '\t' << refSums[i];
497 subMatrixFile << endl;
498 subMatrixFile.close();
504 catch(exception& e) {
505 m->errorOut(e, "SeqErrorCommand", "execute");
510 //***************************************************************************************************************
512 void SeqErrorCommand::getReferences(){
515 ifstream referenceFile;
516 m->openInputFile(referenceFileName, referenceFile);
518 while(referenceFile){
519 Sequence currentSeq(referenceFile);
520 int numAmbigs = currentSeq.getAmbigBases();
523 m->mothurOut("Warning: " + toString(currentSeq.getName()) + " has " + toString(numAmbigs) + " ambiguous bases, these bases will be removed\n");
524 currentSeq.removeAmbigBases();
526 referenceSeqs.push_back(currentSeq);
527 m->gobble(referenceFile);
529 numRefs = referenceSeqs.size();
531 referenceFile.close();
533 catch(exception& e) {
534 m->errorOut(e, "SeqErrorCommand", "getReferences");
539 //***************************************************************************************************************
541 Compare SeqErrorCommand::getErrors(Sequence query, Sequence reference){
543 if(query.getAlignLength() != reference.getAlignLength()){
544 m->mothurOut("Warning: " + toString(query.getName()) + " and " + toString(reference.getName()) + " are different lengths\n");
546 int alignLength = query.getAlignLength();
548 string q = query.getAligned();
549 string r = reference.getAligned();
554 for(int i=0;i<alignLength;i++){
555 if(q[i] != '.' && r[i] != '.' && (q[i] != '-' || r[i] != '-')){ // no missing data and no double gaps
559 if(r[i] == 'A'){ errors.AA++; errors.matches++; errors.sequence += 'm'; }
560 if(r[i] == 'T'){ errors.AT++; errors.sequence += 's'; }
561 if(r[i] == 'G'){ errors.AG++; errors.sequence += 's'; }
562 if(r[i] == 'C'){ errors.AC++; errors.sequence += 's'; }
563 if(r[i] == '-'){ errors.Ai++; errors.sequence += 'i'; }
565 else if(q[i] == 'T'){
566 if(r[i] == 'A'){ errors.TA++; errors.sequence += 's'; }
567 if(r[i] == 'T'){ errors.TT++; errors.matches++; errors.sequence += 'm'; }
568 if(r[i] == 'G'){ errors.TG++; errors.sequence += 's'; }
569 if(r[i] == 'C'){ errors.TC++; errors.sequence += 's'; }
570 if(r[i] == '-'){ errors.Ti++; errors.sequence += 'i'; }
572 else if(q[i] == 'G'){
573 if(r[i] == 'A'){ errors.GA++; errors.sequence += 's'; }
574 if(r[i] == 'T'){ errors.GT++; errors.sequence += 's'; }
575 if(r[i] == 'G'){ errors.GG++; errors.matches++; errors.sequence += 'm'; }
576 if(r[i] == 'C'){ errors.GC++; errors.sequence += 's'; }
577 if(r[i] == '-'){ errors.Gi++; errors.sequence += 'i'; }
579 else if(q[i] == 'C'){
580 if(r[i] == 'A'){ errors.CA++; errors.sequence += 's'; }
581 if(r[i] == 'T'){ errors.CT++; errors.sequence += 's'; }
582 if(r[i] == 'G'){ errors.CG++; errors.sequence += 's'; }
583 if(r[i] == 'C'){ errors.CC++; errors.matches++; errors.sequence += 'm'; }
584 if(r[i] == '-'){ errors.Ci++; errors.sequence += 'i'; }
586 else if(q[i] == 'N'){
587 if(r[i] == 'A'){ errors.NA++; errors.sequence += 'a'; }
588 if(r[i] == 'T'){ errors.NT++; errors.sequence += 'a'; }
589 if(r[i] == 'G'){ errors.NG++; errors.sequence += 'a'; }
590 if(r[i] == 'C'){ errors.NC++; errors.sequence += 'a'; }
591 if(r[i] == '-'){ errors.Ni++; errors.sequence += 'a'; }
593 else if(q[i] == '-' && r[i] != '-'){
594 if(r[i] == 'A'){ errors.dA++; errors.sequence += 'd'; }
595 if(r[i] == 'T'){ errors.dT++; errors.sequence += 'd'; }
596 if(r[i] == 'G'){ errors.dG++; errors.sequence += 'd'; }
597 if(r[i] == 'C'){ errors.dC++; errors.sequence += 'd'; }
602 else if(q[i] == '.' && r[i] != '.'){ // reference extends beyond query
603 if(started == 1){ break; }
605 else if(q[i] != '.' && r[i] == '.'){ // query extends beyond reference
606 m->mothurOut("Warning: " + toString(query.getName()) + " extend beyond " + toString(reference.getName()) + ". Ignoring the extra bases in the query\n");
607 if(started == 1){ break; }
609 else if(q[i] == '.' && r[i] == '.'){ // both are missing data
610 if(started == 1){ break; }
614 errors.mismatches = errors.total-errors.matches;
615 errors.errorRate = (double)(errors.total-errors.matches) / (double)errors.total;
616 errors.queryName = query.getName();
617 errors.refName = reference.getName();
621 catch(exception& e) {
622 m->errorOut(e, "SeqErrorCommand", "getErrors");
627 //***************************************************************************************************************
629 map<string, int> SeqErrorCommand::getWeights(){
631 m->openInputFile(namesFileName, nameFile);
634 string redundantSeqs;
635 map<string, int> nameCountMap;
638 nameFile >> seqName >> redundantSeqs;
639 nameCountMap[seqName] = m->getNumNames(redundantSeqs);
646 //***************************************************************************************************************
648 void SeqErrorCommand::printErrorHeader(){
650 errorSummaryFile << "query\treference\tweight\t";
651 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";
652 errorSummaryFile << "insertions\tdeletions\tsubstitutions\tambig\tmatches\tmismatches\ttotal\terror\n";
654 errorSummaryFile << setprecision(6);
655 errorSummaryFile.setf(ios::fixed);
657 catch(exception& e) {
658 m->errorOut(e, "SeqErrorCommand", "printErrorHeader");
663 //***************************************************************************************************************
665 void SeqErrorCommand::printErrorData(Compare error){
667 errorSummaryFile << error.queryName << '\t' << error.refName << '\t' << error.weight << '\t';
668 errorSummaryFile << error.AA << '\t' << error.AT << '\t' << error.AG << '\t' << error.AC << '\t';
669 errorSummaryFile << error.TA << '\t' << error.TT << '\t' << error.TG << '\t' << error.TC << '\t';
670 errorSummaryFile << error.GA << '\t' << error.GT << '\t' << error.GG << '\t' << error.GC << '\t';
671 errorSummaryFile << error.CA << '\t' << error.CT << '\t' << error.CG << '\t' << error.CC << '\t';
672 errorSummaryFile << error.NA << '\t' << error.NT << '\t' << error.NG << '\t' << error.NC << '\t';
673 errorSummaryFile << error.Ai << '\t' << error.Ti << '\t' << error.Gi << '\t' << error.Ci << '\t' << error.Ni << '\t' ;
674 errorSummaryFile << error.dA << '\t' << error.dT << '\t' << error.dG << '\t' << error.dC << '\t';
676 errorSummaryFile << error.Ai + error.Ti + error.Gi + error.Ci << '\t'; //insertions
677 errorSummaryFile << error.dA + error.dT + error.dG + error.dC << '\t'; //deletions
678 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
679 errorSummaryFile << error.NA + error.NT + error.NG + error.NC + error.Ni << '\t'; //ambiguities
680 errorSummaryFile << error.matches << '\t' << error.mismatches << '\t' << error.total << '\t' << error.errorRate << endl;
682 errorSeqFile << '>' << error.queryName << "\tref:" << error.refName << '\n' << error.sequence << endl;
685 int a=0; int t=1; int g=2; int c=3;
688 substitutionMatrix[a][a] += error.weight * error.AA;
689 substitutionMatrix[a][t] += error.weight * error.TA;
690 substitutionMatrix[a][g] += error.weight * error.GA;
691 substitutionMatrix[a][c] += error.weight * error.CA;
692 substitutionMatrix[a][gap] += error.weight * error.dA;
693 substitutionMatrix[a][n] += error.weight * error.NA;
695 substitutionMatrix[t][a] += error.weight * error.AT;
696 substitutionMatrix[t][t] += error.weight * error.TT;
697 substitutionMatrix[t][g] += error.weight * error.GT;
698 substitutionMatrix[t][c] += error.weight * error.CT;
699 substitutionMatrix[t][gap] += error.weight * error.dT;
700 substitutionMatrix[t][n] += error.weight * error.NT;
702 substitutionMatrix[g][a] += error.weight * error.AG;
703 substitutionMatrix[g][t] += error.weight * error.TG;
704 substitutionMatrix[g][g] += error.weight * error.GG;
705 substitutionMatrix[g][c] += error.weight * error.CG;
706 substitutionMatrix[g][gap] += error.weight * error.dG;
707 substitutionMatrix[g][n] += error.weight * error.NG;
709 substitutionMatrix[c][a] += error.weight * error.AC;
710 substitutionMatrix[c][t] += error.weight * error.TC;
711 substitutionMatrix[c][g] += error.weight * error.GC;
712 substitutionMatrix[c][c] += error.weight * error.CC;
713 substitutionMatrix[c][gap] += error.weight * error.dC;
714 substitutionMatrix[c][n] += error.weight * error.NC;
716 substitutionMatrix[gap][a] += error.weight * error.Ai;
717 substitutionMatrix[gap][t] += error.weight * error.Ti;
718 substitutionMatrix[gap][g] += error.weight * error.Gi;
719 substitutionMatrix[gap][c] += error.weight * error.Ci;
720 substitutionMatrix[gap][n] += error.weight * error.Ni;
723 catch(exception& e) {
724 m->errorOut(e, "SeqErrorCommand", "printErrorData");
729 //***************************************************************************************************************