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"
15 //**********************************************************************************************************************
16 vector<string> SeqErrorCommand::getValidParameters(){
18 string Array[] = {"query", "reference", "name", "qfile", "report", "threshold", "inputdir", "outputdir"};
19 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
23 m->errorOut(e, "SeqErrorCommand", "getValidParameters");
27 //**********************************************************************************************************************
28 SeqErrorCommand::SeqErrorCommand(){
30 abort = true; calledHelp = true;
31 vector<string> tempOutNames;
32 outputTypes["error.summary"] = tempOutNames;
33 outputTypes["error.seq"] = tempOutNames;
34 outputTypes["error.quality"] = tempOutNames;
35 outputTypes["error.qual.forward"] = tempOutNames;
36 outputTypes["error.qual.reverse"] = tempOutNames;
37 outputTypes["error.forward"] = tempOutNames;
38 outputTypes["error.reverse"] = tempOutNames;
39 outputTypes["error.count"] = tempOutNames;
40 outputTypes["error.matrix"] = tempOutNames;
43 m->errorOut(e, "SeqErrorCommand", "SeqErrorCommand");
47 //**********************************************************************************************************************
48 vector<string> SeqErrorCommand::getRequiredParameters(){
50 string Array[] = {"query","reference"};
51 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
55 m->errorOut(e, "SeqErrorCommand", "getRequiredParameters");
59 //**********************************************************************************************************************
60 vector<string> SeqErrorCommand::getRequiredFiles(){
62 vector<string> myArray;
66 m->errorOut(e, "SeqErrorCommand", "getRequiredFiles");
70 //***************************************************************************************************************
72 SeqErrorCommand::SeqErrorCommand(string option) {
75 abort = false; calledHelp = false;
77 //allow user to run help
78 if(option == "help") { help(); abort = true; calledHelp = true; }
83 //valid paramters for this command
84 string AlignArray[] = {"query", "reference", "name", "qfile", "report", "threshold", "inputdir", "ignorechimeras", "outputdir"};
86 vector<string> myArray (AlignArray, AlignArray+(sizeof(AlignArray)/sizeof(string)));
88 OptionParser parser(option);
89 map<string,string> parameters = parser.getParameters();
91 ValidParameters validParameter;
92 map<string,string>::iterator it;
94 //check to make sure all parameters are valid for command
95 for (it = parameters.begin(); it != parameters.end(); it++) {
96 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
99 //initialize outputTypes
100 vector<string> tempOutNames;
101 outputTypes["error.summary"] = tempOutNames;
102 outputTypes["error.seq"] = tempOutNames;
103 outputTypes["error.quality"] = tempOutNames;
104 outputTypes["error.qual.forward"] = tempOutNames;
105 outputTypes["error.qual.reverse"] = tempOutNames;
106 outputTypes["error.forward"] = tempOutNames;
107 outputTypes["error.reverse"] = tempOutNames;
108 outputTypes["error.count"] = tempOutNames;
109 outputTypes["error.matrix"] = tempOutNames;
112 //if the user changes the input directory command factory will send this info to us in the output parameter
113 string inputDir = validParameter.validFile(parameters, "inputdir", false);
114 if (inputDir == "not found"){ inputDir = ""; }
117 it = parameters.find("query");
118 //user has given a template file
119 if(it != parameters.end()){
120 path = m->hasPath(it->second);
121 //if the user has not given a path then, add inputdir. else leave path alone.
122 if (path == "") { parameters["query"] = inputDir + it->second; }
125 it = parameters.find("reference");
126 //user has given a template file
127 if(it != parameters.end()){
128 path = m->hasPath(it->second);
129 //if the user has not given a path then, add inputdir. else leave path alone.
130 if (path == "") { parameters["reference"] = inputDir + it->second; }
133 it = parameters.find("name");
134 //user has given a names file
135 if(it != parameters.end()){
136 path = m->hasPath(it->second);
137 //if the user has not given a path then, add inputdir. else leave path alone.
138 if (path == "") { parameters["name"] = inputDir + it->second; }
141 it = parameters.find("qfile");
142 //user has given a quality score file
143 if(it != parameters.end()){
144 path = m->hasPath(it->second);
145 //if the user has not given a path then, add inputdir. else leave path alone.
146 if (path == "") { parameters["qfile"] = inputDir + it->second; }
149 it = parameters.find("report");
150 //user has given a alignment report file
151 if(it != parameters.end()){
152 path = m->hasPath(it->second);
153 //if the user has not given a path then, add inputdir. else leave path alone.
154 if (path == "") { parameters["report"] = inputDir + it->second; }
158 //check for required parameters
159 queryFileName = validParameter.validFile(parameters, "query", true);
160 if (queryFileName == "not found") { m->mothurOut("query is a required parameter for the seq.error command."); m->mothurOutEndLine(); abort = true; }
161 else if (queryFileName == "not open") { abort = true; }
163 referenceFileName = validParameter.validFile(parameters, "reference", true);
164 if (referenceFileName == "not found") { m->mothurOut("reference is a required parameter for the seq.error command."); m->mothurOutEndLine(); abort = true; }
165 else if (referenceFileName == "not open") { abort = true; }
168 //check for optional parameters
169 namesFileName = validParameter.validFile(parameters, "name", true);
170 if(namesFileName == "not found"){ namesFileName = ""; }
172 qualFileName = validParameter.validFile(parameters, "qfile", true);
173 if(qualFileName == "not found"){ qualFileName = ""; }
175 reportFileName = validParameter.validFile(parameters, "report", true);
176 if(reportFileName == "not found"){ reportFileName = ""; }
178 if((reportFileName != "" && qualFileName == "") || (reportFileName == "" && qualFileName != "")){
179 m->mothurOut("if you use either a qual file or a report file, you have to have both.");
180 m->mothurOutEndLine();
184 outputDir = validParameter.validFile(parameters, "outputdir", false);
185 if (outputDir == "not found"){
187 outputDir += m->hasPath(queryFileName); //if user entered a file with a path then preserve it
190 //check for optional parameter and set defaults
191 // ...at some point should added some additional type checking...
192 temp = validParameter.validFile(parameters, "threshold", false); if (temp == "not found") { temp = "1.00"; }
193 convert(temp, threshold);
195 temp = validParameter.validFile(parameters, "ignorechimeras", false); if (temp == "not found") { temp = "1"; }
196 convert(temp, ignoreChimeras);
198 substitutionMatrix.resize(6);
199 for(int i=0;i<6;i++){ substitutionMatrix[i].assign(6,0); }
202 catch(exception& e) {
203 m->errorOut(e, "SeqErrorCommand", "SeqErrorCommand");
208 //**********************************************************************************************************************
210 void SeqErrorCommand::help(){
212 m->mothurOut("The seq.error command reads a query alignment file and a reference alignment file and creates .....\n");
213 m->mothurOut("Example seq.error(...).\n");
214 m->mothurOut("Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFasta).\n");
215 m->mothurOut("For more details please check out the wiki http://www.mothur.org/wiki/seq.error .\n\n");
217 catch(exception& e) {
218 m->errorOut(e, "SeqErrorCommand", "help");
223 //***************************************************************************************************************
225 SeqErrorCommand::~SeqErrorCommand(){
229 //***************************************************************************************************************
231 int SeqErrorCommand::execute(){
233 if (abort == true) { if (calledHelp) { return 0; } return 2; }
235 string errorSummaryFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.summary";
236 m->openOutputFile(errorSummaryFileName, errorSummaryFile);
237 outputNames.push_back(errorSummaryFileName); outputTypes["error.summary"].push_back(errorSummaryFileName);
240 string errorSeqFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.seq";
241 m->openOutputFile(errorSeqFileName, errorSeqFile);
242 outputNames.push_back(errorSeqFileName); outputTypes["error.seq"].push_back(errorSeqFileName);
244 getReferences(); //read in reference sequences - make sure there's no ambiguous bases
246 map<string, int> weights;
247 if(namesFileName != ""){ weights = getWeights(); }
250 m->openInputFile(queryFileName, queryFile);
256 QualityScores quality;
257 vector<vector<int> > qualForwardMap;
258 vector<vector<int> > qualReverseMap;
260 if(qualFileName != "" && reportFileName != ""){
261 m->openInputFile(qualFileName, qualFile);
262 report = ReportFile(reportFile, reportFileName);
264 qualForwardMap.resize(1000);
265 qualReverseMap.resize(1000);
266 for(int i=0;i<1000;i++){
267 qualForwardMap[i].assign(100,0);
268 qualReverseMap[i].assign(100,0);
273 int totalMatches = 0;
275 vector<int> misMatchCounts(11, 0);
279 map<string, int>::iterator it;
280 map<char, vector<int> > qScoreErrorMap;
281 qScoreErrorMap['m'].assign(41, 0);
282 qScoreErrorMap['s'].assign(41, 0);
283 qScoreErrorMap['i'].assign(41, 0);
284 qScoreErrorMap['a'].assign(41, 0);
286 map<char, vector<int> > errorForward;
287 errorForward['m'].assign(1000,0);
288 errorForward['s'].assign(1000,0);
289 errorForward['i'].assign(1000,0);
290 errorForward['d'].assign(1000,0);
291 errorForward['a'].assign(1000,0);
293 map<char, vector<int> > errorReverse;
294 errorReverse['m'].assign(1000,0);
295 errorReverse['s'].assign(1000,0);
296 errorReverse['i'].assign(1000,0);
297 errorReverse['d'].assign(1000,0);
298 errorReverse['a'].assign(1000,0);
301 string errorChimeraFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.chimera";
302 RefChimeraTest chimeraTest(referenceSeqs, errorChimeraFileName);
303 outputNames.push_back(errorChimeraFileName); outputTypes["error.chimera"].push_back(errorChimeraFileName);
310 if (m->control_pressed) { errorSummaryFile.close(); errorSeqFile.close(); for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
312 Sequence query(queryFile);
314 int numParentSeqs = chimeraTest.analyzeQuery(query.getName(), query.getAligned());
315 int closestRefIndex = chimeraTest.getClosestRefIndex();
317 if(numParentSeqs > 1 && ignoreChimeras == 1) { ignoreSeq = 1; }
318 else { ignoreSeq = 0; }
321 Compare minCompare = getErrors(query, referenceSeqs[closestRefIndex]);
323 if(namesFileName != ""){
324 it = weights.find(query.getName());
325 minCompare.weight = it->second;
327 else { minCompare.weight = 1; }
329 printErrorData(minCompare, numParentSeqs);
332 for(int i=0;i<minCompare.total;i++){
333 char letter = minCompare.sequence[i];
334 errorForward[letter][i] += minCompare.weight;
335 errorReverse[letter][minCompare.total-i-1] += minCompare.weight;
339 if(qualFileName != "" && reportFileName != ""){
340 report = ReportFile(reportFile);
342 int origLength = report.getQueryLength();
343 int startBase = report.getQueryStart();
344 int endBase = report.getQueryEnd();
346 quality = QualityScores(qualFile, origLength);
349 quality.updateQScoreErrorMap(qScoreErrorMap, minCompare.sequence, startBase, endBase, minCompare.weight);
350 quality.updateForwardMap(qualForwardMap, startBase, endBase, minCompare.weight);
351 quality.updateReverseMap(qualReverseMap, startBase, endBase, minCompare.weight);
355 if(minCompare.errorRate < threshold && !ignoreSeq){
356 totalBases += (minCompare.total * minCompare.weight);
357 totalMatches += minCompare.matches * minCompare.weight;
358 if(minCompare.mismatches > maxMismatch){
359 maxMismatch = minCompare.mismatches;
360 misMatchCounts.resize(maxMismatch + 1, 0);
362 misMatchCounts[minCompare.mismatches] += minCompare.weight;
367 if(index % 1000 == 0){ cout << index << endl; }
370 errorSummaryFile.close();
371 errorSeqFile.close();
373 if(qualFileName != "" && reportFileName != ""){
374 printErrorQuality(qScoreErrorMap);
375 printQualityFR(qualForwardMap, qualReverseMap);
378 printErrorFRFile(errorForward, errorReverse);
380 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
382 string errorCountFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.count";
383 ofstream errorCountFile;
384 m->openOutputFile(errorCountFileName, errorCountFile);
385 outputNames.push_back(errorCountFileName); outputTypes["error.count"].push_back(errorCountFileName);
386 m->mothurOut("Overall error rate:\t" + toString((double)(totalBases - totalMatches) / (double)totalBases) + "\n\n");
387 m->mothurOut("Errors\tSequences\n");
388 errorCountFile << "Errors\tSequences\n";
389 for(int i=0;i<misMatchCounts.size();i++){
390 m->mothurOut(toString(i) + '\t' + toString(misMatchCounts[i]) + '\n');
391 errorCountFile << i << '\t' << misMatchCounts[i] << endl;
393 errorCountFile.close();
395 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
399 m->mothurOutEndLine();
400 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
401 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
402 m->mothurOutEndLine();
406 catch(exception& e) {
407 m->errorOut(e, "SeqErrorCommand", "execute");
412 //***************************************************************************************************************
414 void SeqErrorCommand::getReferences(){
417 ifstream referenceFile;
418 m->openInputFile(referenceFileName, referenceFile);
420 while(referenceFile){
421 Sequence currentSeq(referenceFile);
422 int numAmbigs = currentSeq.getAmbigBases();
425 m->mothurOut("Warning: " + toString(currentSeq.getName()) + " has " + toString(numAmbigs) + " ambiguous bases, these bases will be removed\n");
426 currentSeq.removeAmbigBases();
428 referenceSeqs.push_back(currentSeq);
429 m->gobble(referenceFile);
431 numRefs = referenceSeqs.size();
433 referenceFile.close();
435 catch(exception& e) {
436 m->errorOut(e, "SeqErrorCommand", "getReferences");
441 //***************************************************************************************************************
443 Compare SeqErrorCommand::getErrors(Sequence query, Sequence reference){
445 if(query.getAlignLength() != reference.getAlignLength()){
446 m->mothurOut("Warning: " + toString(query.getName()) + " and " + toString(reference.getName()) + " are different lengths\n");
448 int alignLength = query.getAlignLength();
450 string q = query.getAligned();
451 string r = reference.getAligned();
456 for(int i=0;i<alignLength;i++){
457 if(q[i] != '.' && r[i] != '.' && (q[i] != '-' || r[i] != '-')){ // no missing data and no double gaps
461 if(r[i] == 'A'){ errors.AA++; errors.matches++; errors.sequence += 'm'; }
462 if(r[i] == 'T'){ errors.AT++; errors.sequence += 's'; }
463 if(r[i] == 'G'){ errors.AG++; errors.sequence += 's'; }
464 if(r[i] == 'C'){ errors.AC++; errors.sequence += 's'; }
465 if(r[i] == '-'){ errors.Ai++; errors.sequence += 'i'; }
467 else if(q[i] == 'T'){
468 if(r[i] == 'A'){ errors.TA++; errors.sequence += 's'; }
469 if(r[i] == 'T'){ errors.TT++; errors.matches++; errors.sequence += 'm'; }
470 if(r[i] == 'G'){ errors.TG++; errors.sequence += 's'; }
471 if(r[i] == 'C'){ errors.TC++; errors.sequence += 's'; }
472 if(r[i] == '-'){ errors.Ti++; errors.sequence += 'i'; }
474 else if(q[i] == 'G'){
475 if(r[i] == 'A'){ errors.GA++; errors.sequence += 's'; }
476 if(r[i] == 'T'){ errors.GT++; errors.sequence += 's'; }
477 if(r[i] == 'G'){ errors.GG++; errors.matches++; errors.sequence += 'm'; }
478 if(r[i] == 'C'){ errors.GC++; errors.sequence += 's'; }
479 if(r[i] == '-'){ errors.Gi++; errors.sequence += 'i'; }
481 else if(q[i] == 'C'){
482 if(r[i] == 'A'){ errors.CA++; errors.sequence += 's'; }
483 if(r[i] == 'T'){ errors.CT++; errors.sequence += 's'; }
484 if(r[i] == 'G'){ errors.CG++; errors.sequence += 's'; }
485 if(r[i] == 'C'){ errors.CC++; errors.matches++; errors.sequence += 'm'; }
486 if(r[i] == '-'){ errors.Ci++; errors.sequence += 'i'; }
488 else if(q[i] == 'N'){
489 if(r[i] == 'A'){ errors.NA++; errors.sequence += 'a'; }
490 if(r[i] == 'T'){ errors.NT++; errors.sequence += 'a'; }
491 if(r[i] == 'G'){ errors.NG++; errors.sequence += 'a'; }
492 if(r[i] == 'C'){ errors.NC++; errors.sequence += 'a'; }
493 if(r[i] == '-'){ errors.Ni++; errors.sequence += 'a'; }
495 else if(q[i] == '-' && r[i] != '-'){
496 if(r[i] == 'A'){ errors.dA++; errors.sequence += 'd'; }
497 if(r[i] == 'T'){ errors.dT++; errors.sequence += 'd'; }
498 if(r[i] == 'G'){ errors.dG++; errors.sequence += 'd'; }
499 if(r[i] == 'C'){ errors.dC++; errors.sequence += 'd'; }
504 else if(q[i] == '.' && r[i] != '.'){ // reference extends beyond query
505 if(started == 1){ break; }
507 else if(q[i] != '.' && r[i] == '.'){ // query extends beyond reference
508 m->mothurOut("Warning: " + toString(query.getName()) + " extend beyond " + toString(reference.getName()) + ". Ignoring the extra bases in the query\n");
509 if(started == 1){ break; }
511 else if(q[i] == '.' && r[i] == '.'){ // both are missing data
512 if(started == 1){ break; }
516 errors.mismatches = errors.total-errors.matches;
517 errors.errorRate = (double)(errors.total-errors.matches) / (double)errors.total;
518 errors.queryName = query.getName();
519 errors.refName = reference.getName();
523 catch(exception& e) {
524 m->errorOut(e, "SeqErrorCommand", "getErrors");
529 //***************************************************************************************************************
531 map<string, int> SeqErrorCommand::getWeights(){
533 m->openInputFile(namesFileName, nameFile);
536 string redundantSeqs;
537 map<string, int> nameCountMap;
540 nameFile >> seqName >> redundantSeqs;
541 nameCountMap[seqName] = m->getNumNames(redundantSeqs);
548 //***************************************************************************************************************
550 void SeqErrorCommand::printErrorHeader(){
552 errorSummaryFile << "query\treference\tweight\t";
553 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";
554 errorSummaryFile << "insertions\tdeletions\tsubstitutions\tambig\tmatches\tmismatches\ttotal\terror\tnumparents\n";
556 errorSummaryFile << setprecision(6);
557 errorSummaryFile.setf(ios::fixed);
559 catch(exception& e) {
560 m->errorOut(e, "SeqErrorCommand", "printErrorHeader");
565 //***************************************************************************************************************
567 void SeqErrorCommand::printErrorData(Compare error, int numParentSeqs){
569 errorSummaryFile << error.queryName << '\t' << error.refName << '\t' << error.weight << '\t';
570 errorSummaryFile << error.AA << '\t' << error.AT << '\t' << error.AG << '\t' << error.AC << '\t';
571 errorSummaryFile << error.TA << '\t' << error.TT << '\t' << error.TG << '\t' << error.TC << '\t';
572 errorSummaryFile << error.GA << '\t' << error.GT << '\t' << error.GG << '\t' << error.GC << '\t';
573 errorSummaryFile << error.CA << '\t' << error.CT << '\t' << error.CG << '\t' << error.CC << '\t';
574 errorSummaryFile << error.NA << '\t' << error.NT << '\t' << error.NG << '\t' << error.NC << '\t';
575 errorSummaryFile << error.Ai << '\t' << error.Ti << '\t' << error.Gi << '\t' << error.Ci << '\t' << error.Ni << '\t' ;
576 errorSummaryFile << error.dA << '\t' << error.dT << '\t' << error.dG << '\t' << error.dC << '\t';
578 errorSummaryFile << error.Ai + error.Ti + error.Gi + error.Ci << '\t'; //insertions
579 errorSummaryFile << error.dA + error.dT + error.dG + error.dC << '\t'; //deletions
580 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
581 errorSummaryFile << error.NA + error.NT + error.NG + error.NC + error.Ni << '\t'; //ambiguities
582 errorSummaryFile << error.matches << '\t' << error.mismatches << '\t' << error.total << '\t' << error.errorRate << '\t' << numParentSeqs << endl;
584 errorSeqFile << '>' << error.queryName << "\tref:" << error.refName << '\n' << error.sequence << endl;
587 int a=0; int t=1; int g=2; int c=3;
589 if(numParentSeqs == 1 || ignoreChimeras == 0){
590 substitutionMatrix[a][a] += error.weight * error.AA;
591 substitutionMatrix[a][t] += error.weight * error.TA;
592 substitutionMatrix[a][g] += error.weight * error.GA;
593 substitutionMatrix[a][c] += error.weight * error.CA;
594 substitutionMatrix[a][gap] += error.weight * error.dA;
595 substitutionMatrix[a][n] += error.weight * error.NA;
597 substitutionMatrix[t][a] += error.weight * error.AT;
598 substitutionMatrix[t][t] += error.weight * error.TT;
599 substitutionMatrix[t][g] += error.weight * error.GT;
600 substitutionMatrix[t][c] += error.weight * error.CT;
601 substitutionMatrix[t][gap] += error.weight * error.dT;
602 substitutionMatrix[t][n] += error.weight * error.NT;
604 substitutionMatrix[g][a] += error.weight * error.AG;
605 substitutionMatrix[g][t] += error.weight * error.TG;
606 substitutionMatrix[g][g] += error.weight * error.GG;
607 substitutionMatrix[g][c] += error.weight * error.CG;
608 substitutionMatrix[g][gap] += error.weight * error.dG;
609 substitutionMatrix[g][n] += error.weight * error.NG;
611 substitutionMatrix[c][a] += error.weight * error.AC;
612 substitutionMatrix[c][t] += error.weight * error.TC;
613 substitutionMatrix[c][g] += error.weight * error.GC;
614 substitutionMatrix[c][c] += error.weight * error.CC;
615 substitutionMatrix[c][gap] += error.weight * error.dC;
616 substitutionMatrix[c][n] += error.weight * error.NC;
618 substitutionMatrix[gap][a] += error.weight * error.Ai;
619 substitutionMatrix[gap][t] += error.weight * error.Ti;
620 substitutionMatrix[gap][g] += error.weight * error.Gi;
621 substitutionMatrix[gap][c] += error.weight * error.Ci;
622 substitutionMatrix[gap][n] += error.weight * error.Ni;
625 catch(exception& e) {
626 m->errorOut(e, "SeqErrorCommand", "printErrorData");
631 //***************************************************************************************************************
633 void SeqErrorCommand::printSubMatrix(){
635 string subMatrixFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.matrix";
636 ofstream subMatrixFile;
637 m->openOutputFile(subMatrixFileName, subMatrixFile);
638 outputNames.push_back(subMatrixFileName); outputTypes["error.matrix"].push_back(subMatrixFileName);
639 vector<string> bases(6);
646 vector<int> refSums(5,1);
648 for(int i=0;i<5;i++){
649 subMatrixFile << "\tr" << bases[i];
651 for(int j=0;j<6;j++){
652 refSums[i] += substitutionMatrix[i][j];
655 subMatrixFile << endl;
657 for(int i=0;i<6;i++){
658 subMatrixFile << 'q' << bases[i];
659 for(int j=0;j<5;j++){
660 subMatrixFile << '\t' << substitutionMatrix[j][i];
662 subMatrixFile << endl;
665 subMatrixFile << "total";
666 for(int i=0;i<5;i++){
667 subMatrixFile << '\t' << refSums[i];
669 subMatrixFile << endl;
670 subMatrixFile.close();
672 catch(exception& e) {
673 m->errorOut(e, "SeqErrorCommand", "printSubMatrix");
677 //***************************************************************************************************************
679 void SeqErrorCommand::printErrorFRFile(map<char, vector<int> > errorForward, map<char, vector<int> > errorReverse){
681 string errorForwardFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.seq.forward";
682 ofstream errorForwardFile;
683 m->openOutputFile(errorForwardFileName, errorForwardFile);
684 outputNames.push_back(errorForwardFileName); outputTypes["error.forward"].push_back(errorForwardFileName);
686 errorForwardFile << "position\ttotalseqs\tmatch\tsubstitution\tinsertion\tdeletion\tambiguous" << endl;
687 for(int i=0;i<1000;i++){
688 float match = (float)errorForward['m'][i];
689 float subst = (float)errorForward['s'][i];
690 float insert = (float)errorForward['i'][i];
691 float del = (float)errorForward['d'][i];
692 float amb = (float)errorForward['a'][i];
693 float total = match + subst + insert + del + amb;
694 if(total == 0){ break; }
695 errorForwardFile << i+1 << '\t' << total << '\t' << match/total << '\t' << subst/total << '\t' << insert/total << '\t' << del/total << '\t' << amb/total << endl;
697 errorForwardFile.close();
699 string errorReverseFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.seq.reverse";
700 ofstream errorReverseFile;
701 m->openOutputFile(errorReverseFileName, errorReverseFile);
702 outputNames.push_back(errorReverseFileName); outputTypes["error.reverse"].push_back(errorReverseFileName);
704 errorReverseFile << "position\ttotalseqs\tmatch\tsubstitution\tinsertion\tdeletion\tambiguous" << endl;
705 for(int i=0;i<1000;i++){
706 float match = (float)errorReverse['m'][i];
707 float subst = (float)errorReverse['s'][i];
708 float insert = (float)errorReverse['i'][i];
709 float del = (float)errorReverse['d'][i];
710 float amb = (float)errorReverse['a'][i];
711 float total = match + subst + insert + del + amb;
712 if(total == 0){ break; }
713 errorReverseFile << i+1 << '\t' << total << '\t' << match/total << '\t' << subst/total << '\t' << insert/total << '\t' << del/total << '\t' << amb/total << endl;
715 errorReverseFile.close();
717 catch(exception& e) {
718 m->errorOut(e, "SeqErrorCommand", "printErrorFRFile");
723 //***************************************************************************************************************
725 void SeqErrorCommand::printErrorQuality(map<char, vector<int> > qScoreErrorMap){
728 string errorQualityFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.quality";
729 ofstream errorQualityFile;
730 m->openOutputFile(errorQualityFileName, errorQualityFile);
731 outputNames.push_back(errorQualityFileName); outputTypes["error.quality"].push_back(errorQualityFileName);
733 errorQualityFile << "qscore\tmatches\tsubstitutions\tinsertions\tambiguous" << endl;
734 for(int i=0;i<41;i++){
735 errorQualityFile << i << '\t' << qScoreErrorMap['m'][i] << '\t' << qScoreErrorMap['s'][i] << '\t' << qScoreErrorMap['i'][i] << '\t'<< qScoreErrorMap['a'][i] << endl;
737 errorQualityFile.close();
739 catch(exception& e) {
740 m->errorOut(e, "SeqErrorCommand", "printErrorFRFile");
746 //***************************************************************************************************************
748 void SeqErrorCommand::printQualityFR(vector<vector<int> > qualForwardMap, vector<vector<int> > qualReverseMap){
755 for(int i=0;i<qualForwardMap.size();i++){
756 for(int j=0;j<qualForwardMap[i].size();j++){
757 if(qualForwardMap[i][j] != 0){
758 if(lastRow < i) { lastRow = i+2; }
759 if(lastColumn < j) { lastColumn = j+2; }
764 string qualityForwardFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.qual.forward";
765 ofstream qualityForwardFile;
766 m->openOutputFile(qualityForwardFileName, qualityForwardFile);
767 outputNames.push_back(qualityForwardFileName); outputTypes["error.qual.forward"].push_back(qualityForwardFileName);
769 for(int i=0;i<lastColumn;i++){ qualityForwardFile << '\t' << i; } qualityForwardFile << endl;
771 for(int i=0;i<lastRow;i++){
772 qualityForwardFile << i+1;
773 for(int j=0;j<lastColumn;j++){
774 qualityForwardFile << '\t' << qualForwardMap[i][j];
777 qualityForwardFile << endl;
779 qualityForwardFile.close();
782 string qualityReverseFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.qual.reverse";
783 ofstream qualityReverseFile;
784 m->openOutputFile(qualityReverseFileName, qualityReverseFile);
785 outputNames.push_back(qualityReverseFileName); outputTypes["error.qual.reverse"].push_back(qualityReverseFileName);
787 for(int i=0;i<lastColumn;i++){ qualityReverseFile << '\t' << i; } qualityReverseFile << endl;
788 for(int i=0;i<lastRow;i++){
790 qualityReverseFile << i+1;
791 for(int j=0;j<lastColumn;j++){
792 qualityReverseFile << '\t' << qualReverseMap[i][j];
794 qualityReverseFile << endl;
796 qualityReverseFile.close();
798 catch(exception& e) {
799 m->errorOut(e, "SeqErrorCommand", "printErrorFRFile");
805 //***************************************************************************************************************