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"
17 //**********************************************************************************************************************
18 vector<string> SeqErrorCommand::setParameters(){
20 CommandParameter pquery("fasta", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(pquery);
21 CommandParameter preference("reference", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(preference);
22 CommandParameter pqfile("qfile", "InputTypes", "", "", "none", "none", "QualReport",false,false); parameters.push_back(pqfile);
23 CommandParameter preport("report", "InputTypes", "", "", "none", "none", "QualReport",false,false); parameters.push_back(preport);
24 CommandParameter pname("name", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(pname);
25 CommandParameter pignorechimeras("ignorechimeras", "Boolean", "", "T", "", "", "",false,false); parameters.push_back(pignorechimeras);
26 CommandParameter pthreshold("threshold", "Number", "", "1.0", "", "", "",false,false); parameters.push_back(pthreshold);
27 CommandParameter pprocessors("processors", "Number", "", "1", "", "", "",false,false); parameters.push_back(pprocessors);
28 CommandParameter psave("save", "Boolean", "", "F", "", "", "",false,false); parameters.push_back(psave);
29 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
30 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
32 vector<string> myArray;
33 for (int i = 0; i < parameters.size(); i++) { myArray.push_back(parameters[i].name); }
37 m->errorOut(e, "SeqErrorCommand", "setParameters");
41 //**********************************************************************************************************************
42 string SeqErrorCommand::getHelpString(){
44 string helpString = "";
45 helpString += "The seq.error command reads a query alignment file and a reference alignment file and creates .....\n";
46 helpString += "The fasta parameter...\n";
47 helpString += "The reference parameter...\n";
48 helpString += "The qfile parameter...\n";
49 helpString += "The report parameter...\n";
50 helpString += "The name parameter...\n";
51 helpString += "The ignorechimeras parameter...\n";
52 helpString += "The threshold parameter...\n";
53 helpString += "The processors parameter...\n";
54 helpString += "If the save parameter is set to true the reference sequences will be saved in memory, to clear them later you can use the clear.memory command. Default=f.";
55 helpString += "Example seq.error(...).\n";
56 helpString += "Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFasta).\n";
57 helpString += "For more details please check out the wiki http://www.mothur.org/wiki/seq.error .\n";
61 m->errorOut(e, "SeqErrorCommand", "getHelpString");
65 //**********************************************************************************************************************
66 SeqErrorCommand::SeqErrorCommand(){
68 abort = true; calledHelp = true;
70 vector<string> tempOutNames;
71 outputTypes["error.summary"] = tempOutNames;
72 outputTypes["error.seq"] = tempOutNames;
73 outputTypes["error.quality"] = tempOutNames;
74 outputTypes["error.qual.forward"] = tempOutNames;
75 outputTypes["error.qual.reverse"] = tempOutNames;
76 outputTypes["error.forward"] = tempOutNames;
77 outputTypes["error.reverse"] = tempOutNames;
78 outputTypes["error.count"] = tempOutNames;
79 outputTypes["error.matrix"] = tempOutNames;
82 m->errorOut(e, "SeqErrorCommand", "SeqErrorCommand");
86 //***************************************************************************************************************
88 SeqErrorCommand::SeqErrorCommand(string option) {
91 abort = false; calledHelp = false;
92 rdb = ReferenceDB::getInstance();
94 //allow user to run help
95 if(option == "help") { help(); abort = true; calledHelp = true; }
96 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
100 vector<string> myArray = setParameters();
102 OptionParser parser(option);
103 map<string,string> parameters = parser.getParameters();
105 ValidParameters validParameter;
106 map<string,string>::iterator it;
108 //check to make sure all parameters are valid for command
109 for (it = parameters.begin(); it != parameters.end(); it++) {
110 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
113 //initialize outputTypes
114 vector<string> tempOutNames;
115 outputTypes["error.summary"] = tempOutNames;
116 outputTypes["error.seq"] = tempOutNames;
117 outputTypes["error.quality"] = tempOutNames;
118 outputTypes["error.qual.forward"] = tempOutNames;
119 outputTypes["error.qual.reverse"] = tempOutNames;
120 outputTypes["error.forward"] = tempOutNames;
121 outputTypes["error.reverse"] = tempOutNames;
122 outputTypes["error.count"] = tempOutNames;
123 outputTypes["error.matrix"] = tempOutNames;
126 //if the user changes the input directory command factory will send this info to us in the output parameter
127 string inputDir = validParameter.validFile(parameters, "inputdir", false);
128 if (inputDir == "not found"){ inputDir = ""; }
131 it = parameters.find("fasta");
132 //user has given a template file
133 if(it != parameters.end()){
134 path = m->hasPath(it->second);
135 //if the user has not given a path then, add inputdir. else leave path alone.
136 if (path == "") { parameters["fasta"] = inputDir + it->second; }
139 it = parameters.find("reference");
140 //user has given a template file
141 if(it != parameters.end()){
142 path = m->hasPath(it->second);
143 //if the user has not given a path then, add inputdir. else leave path alone.
144 if (path == "") { parameters["reference"] = inputDir + it->second; }
147 it = parameters.find("name");
148 //user has given a names file
149 if(it != parameters.end()){
150 path = m->hasPath(it->second);
151 //if the user has not given a path then, add inputdir. else leave path alone.
152 if (path == "") { parameters["name"] = inputDir + it->second; }
155 it = parameters.find("qfile");
156 //user has given a quality score file
157 if(it != parameters.end()){
158 path = m->hasPath(it->second);
159 //if the user has not given a path then, add inputdir. else leave path alone.
160 if (path == "") { parameters["qfile"] = inputDir + it->second; }
163 it = parameters.find("report");
164 //user has given a alignment report file
165 if(it != parameters.end()){
166 path = m->hasPath(it->second);
167 //if the user has not given a path then, add inputdir. else leave path alone.
168 if (path == "") { parameters["report"] = inputDir + it->second; }
172 //check for required parameters
173 queryFileName = validParameter.validFile(parameters, "fasta", true);
174 if (queryFileName == "not found") {
175 queryFileName = m->getFastaFile();
176 if (queryFileName != "") { m->mothurOut("Using " + queryFileName + " as input file for the fasta parameter."); m->mothurOutEndLine(); }
177 else { m->mothurOut("You have no current fasta file and the fasta parameter is required."); m->mothurOutEndLine(); abort = true; }
179 else if (queryFileName == "not open") { queryFileName = ""; abort = true; }
180 else { m->setFastaFile(queryFileName); }
182 referenceFileName = validParameter.validFile(parameters, "reference", true);
183 if (referenceFileName == "not found") { m->mothurOut("reference is a required parameter for the seq.error command."); m->mothurOutEndLine(); abort = true; }
184 else if (referenceFileName == "not open") { abort = true; }
186 //check for optional parameters
187 namesFileName = validParameter.validFile(parameters, "name", true);
188 if(namesFileName == "not found"){ namesFileName = ""; }
189 else if (namesFileName == "not open") { namesFileName = ""; abort = true; }
190 else { m->setNameFile(namesFileName); }
192 qualFileName = validParameter.validFile(parameters, "qfile", true);
193 if(qualFileName == "not found"){ qualFileName = ""; }
194 else if (qualFileName == "not open") { qualFileName = ""; abort = true; }
195 else { m->setQualFile(qualFileName); }
197 reportFileName = validParameter.validFile(parameters, "report", true);
198 if(reportFileName == "not found"){ reportFileName = ""; }
199 else if (reportFileName == "not open") { reportFileName = ""; abort = true; }
201 if((reportFileName != "" && qualFileName == "") || (reportFileName == "" && qualFileName != "")){
202 m->mothurOut("if you use either a qual file or a report file, you have to have both.");
203 m->mothurOutEndLine();
207 outputDir = validParameter.validFile(parameters, "outputdir", false);
208 if (outputDir == "not found"){
210 outputDir += m->hasPath(queryFileName); //if user entered a file with a path then preserve it
213 //check for optional parameter and set defaults
214 // ...at some point should added some additional type checking...
215 temp = validParameter.validFile(parameters, "threshold", false); if (temp == "not found") { temp = "1.00"; }
216 m->mothurConvert(temp, threshold);
218 temp = validParameter.validFile(parameters, "save", false); if (temp == "not found"){ temp = "f"; }
219 save = m->isTrue(temp);
221 if (save) { //clear out old references
225 //this has to go after save so that if the user sets save=t and provides no reference we abort
226 referenceFileName = validParameter.validFile(parameters, "reference", true);
227 if (referenceFileName == "not found") {
228 //check for saved reference sequences
229 if (rdb->referenceSeqs.size() != 0) {
230 referenceFileName = "saved";
232 m->mothurOut("[ERROR]: You don't have any saved reference sequences and the reference parameter is a required.");
233 m->mothurOutEndLine();
236 }else if (referenceFileName == "not open") { abort = true; }
237 else { if (save) { rdb->setSavedReference(referenceFileName); } }
240 temp = validParameter.validFile(parameters, "ignorechimeras", false); if (temp == "not found") { temp = "T"; }
241 ignoreChimeras = m->isTrue(temp);
243 temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); }
244 m->setProcessors(temp);
245 m->mothurConvert(temp, processors);
247 substitutionMatrix.resize(6);
248 for(int i=0;i<6;i++){ substitutionMatrix[i].resize(6,0); }
250 if ((namesFileName == "") && (queryFileName != "")){
251 vector<string> files; files.push_back(queryFileName);
252 parser.getNameFile(files);
256 catch(exception& e) {
257 m->errorOut(e, "SeqErrorCommand", "SeqErrorCommand");
261 //***************************************************************************************************************
263 int SeqErrorCommand::execute(){
265 if (abort == true) { if (calledHelp) { return 0; } return 2; }
267 int start = time(NULL);
272 string errorSummaryFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.summary";
273 outputNames.push_back(errorSummaryFileName); outputTypes["error.summary"].push_back(errorSummaryFileName);
275 string errorSeqFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.seq";
276 outputNames.push_back(errorSeqFileName); outputTypes["error.seq"].push_back(errorSeqFileName);
278 string errorChimeraFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.chimera";
279 outputNames.push_back(errorChimeraFileName); outputTypes["error.chimera"].push_back(errorChimeraFileName);
281 getReferences(); //read in reference sequences - make sure there's no ambiguous bases
283 if(namesFileName != ""){ weights = getWeights(); }
285 vector<unsigned long long> fastaFilePos;
286 vector<unsigned long long> qFilePos;
287 vector<unsigned long long> reportFilePos;
289 setLines(queryFileName, qualFileName, reportFileName, fastaFilePos, qFilePos, reportFilePos);
291 if (m->control_pressed) { return 0; }
293 for (int i = 0; i < (fastaFilePos.size()-1); i++) {
294 lines.push_back(linePair(fastaFilePos[i], fastaFilePos[(i+1)]));
295 if (qualFileName != "") { qLines.push_back(linePair(qFilePos[i], qFilePos[(i+1)])); }
296 if (reportFileName != "") { rLines.push_back(linePair(reportFilePos[i], reportFilePos[(i+1)])); }
298 if(qualFileName == "") { qLines = lines; rLines = lines; } //fills with duds
301 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
303 numSeqs = driver(queryFileName, qualFileName, reportFileName, errorSummaryFileName, errorSeqFileName, errorChimeraFileName, lines[0], qLines[0], rLines[0]);
305 numSeqs = createProcesses(queryFileName, qualFileName, reportFileName, errorSummaryFileName, errorSeqFileName, errorChimeraFileName);
308 numSeqs = driver(queryFileName, qualFileName, reportFileName, errorSummaryFileName, errorSeqFileName, errorChimeraFileName, lines[0], qLines[0], rLines[0]);
311 if(qualFileName != "" && reportFileName != ""){
312 printErrorQuality(qScoreErrorMap);
313 printQualityFR(qualForwardMap, qualReverseMap);
316 printErrorFRFile(errorForward, errorReverse);
318 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; }
320 string errorCountFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.count";
321 ofstream errorCountFile;
322 m->openOutputFile(errorCountFileName, errorCountFile);
323 outputNames.push_back(errorCountFileName); outputTypes["error.count"].push_back(errorCountFileName);
324 m->mothurOut("Overall error rate:\t" + toString((double)(totalBases - totalMatches) / (double)totalBases) + "\n");
325 m->mothurOut("Errors\tSequences\n");
326 errorCountFile << "Errors\tSequences\n";
327 for(int i=0;i<misMatchCounts.size();i++){
328 m->mothurOut(toString(i) + '\t' + toString(misMatchCounts[i]) + '\n');
329 errorCountFile << i << '\t' << misMatchCounts[i] << endl;
331 errorCountFile.close();
333 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; }
337 string megAlignmentFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.ref-query";
338 ofstream megAlignmentFile;
339 m->openOutputFile(megAlignmentFileName, megAlignmentFile);
340 outputNames.push_back(megAlignmentFileName); outputTypes["error.ref-query"].push_back(megAlignmentFileName);
342 for(int i=0;i<numRefs;i++){
343 megAlignmentFile << referenceSeqs[i].getInlineSeq() << endl;
344 megAlignmentFile << megaAlignVector[i] << endl;
347 m->mothurOut("It took " + toString(time(NULL) - start) + " secs to check " + toString(numSeqs) + " sequences.");
348 m->mothurOutEndLine();
350 m->mothurOutEndLine();
351 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
352 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
353 m->mothurOutEndLine();
357 catch(exception& e) {
358 m->errorOut(e, "SeqErrorCommand", "execute");
362 //**********************************************************************************************************************
363 int SeqErrorCommand::createProcesses(string filename, string qFileName, string rFileName, string summaryFileName, string errorOutputFileName, string chimeraOutputFileName) {
367 map<char, vector<int> >::iterator it;
369 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
371 //loop through and create all the processes you want
372 while (process != processors) {
376 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
380 num = driver(filename, qFileName, rFileName, summaryFileName + toString(getpid()) + ".temp", errorOutputFileName+ toString(getpid()) + ".temp", chimeraOutputFileName + toString(getpid()) + ".temp", lines[process], qLines[process], rLines[process]);
382 //pass groupCounts to parent
384 string tempFile = filename + toString(getpid()) + ".info.temp";
385 m->openOutputFile(tempFile, out);
387 //output totalBases and totalMatches
388 out << num << '\t' << totalBases << '\t' << totalMatches << endl << endl;
390 //output substitutionMatrix
391 for(int i = 0; i < substitutionMatrix.size(); i++) {
392 for (int j = 0; j < substitutionMatrix[i].size(); j++) {
393 out << substitutionMatrix[i][j] << '\t';
399 //output qScoreErrorMap
400 for (it = qScoreErrorMap.begin(); it != qScoreErrorMap.end(); it++) {
401 vector<int> thisScoreErrorMap = it->second;
402 out << it->first << '\t';
403 for (int i = 0; i < thisScoreErrorMap.size(); i++) {
404 out << thisScoreErrorMap[i] << '\t';
410 //output qualForwardMap
411 for(int i = 0; i < qualForwardMap.size(); i++) {
412 for (int j = 0; j < qualForwardMap[i].size(); j++) {
413 out << qualForwardMap[i][j] << '\t';
419 //output qualReverseMap
420 for(int i = 0; i < qualReverseMap.size(); i++) {
421 for (int j = 0; j < qualReverseMap[i].size(); j++) {
422 out << qualReverseMap[i][j] << '\t';
429 //output errorForward
430 for (it = errorForward.begin(); it != errorForward.end(); it++) {
431 vector<int> thisErrorForward = it->second;
432 out << it->first << '\t';
433 for (int i = 0; i < thisErrorForward.size(); i++) {
434 out << thisErrorForward[i] << '\t';
440 //output errorReverse
441 for (it = errorReverse.begin(); it != errorReverse.end(); it++) {
442 vector<int> thisErrorReverse = it->second;
443 out << it->first << '\t';
444 for (int i = 0; i < thisErrorReverse.size(); i++) {
445 out << thisErrorReverse[i] << '\t';
451 //output misMatchCounts
452 out << misMatchCounts.size() << endl;
453 for (int j = 0; j < misMatchCounts.size(); j++) {
454 out << misMatchCounts[j] << '\t';
459 //output megaAlignVector
460 for (int j = 0; j < megaAlignVector.size(); j++) {
461 out << megaAlignVector[j] << endl;
469 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
470 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
476 num = driver(filename, qFileName, rFileName, summaryFileName, errorOutputFileName, chimeraOutputFileName, lines[0], qLines[0], rLines[0]);
478 //force parent to wait until all the processes are done
479 for (int i=0;i<processIDS.size();i++) {
480 int temp = processIDS[i];
485 for(int i=0;i<processIDS.size();i++){
487 m->mothurOut("Appending files from process " + toString(processIDS[i])); m->mothurOutEndLine();
489 m->appendFiles((summaryFileName + toString(processIDS[i]) + ".temp"), summaryFileName);
490 m->mothurRemove((summaryFileName + toString(processIDS[i]) + ".temp"));
491 m->appendFiles((errorOutputFileName + toString(processIDS[i]) + ".temp"), errorOutputFileName);
492 m->mothurRemove((errorOutputFileName + toString(processIDS[i]) + ".temp"));
493 m->appendFiles((chimeraOutputFileName + toString(processIDS[i]) + ".temp"), chimeraOutputFileName);
494 m->mothurRemove((chimeraOutputFileName + toString(processIDS[i]) + ".temp"));
497 string tempFile = filename + toString(processIDS[i]) + ".info.temp";
498 m->openInputFile(tempFile, in);
500 //input totalBases and totalMatches
501 int tempBases, tempMatches, tempNumSeqs;
502 in >> tempNumSeqs >> tempBases >> tempMatches; m->gobble(in);
503 totalBases += tempBases; totalMatches += tempMatches; num += tempNumSeqs;
505 //input substitutionMatrix
507 for(int i = 0; i < substitutionMatrix.size(); i++) {
508 for (int j = 0; j < substitutionMatrix[i].size(); j++) {
509 in >> tempNum; substitutionMatrix[i][j] += tempNum;
515 //input qScoreErrorMap
517 for (int i = 0; i < qScoreErrorMap.size(); i++) {
519 vector<int> thisScoreErrorMap = qScoreErrorMap[first];
521 for (int i = 0; i < thisScoreErrorMap.size(); i++) {
522 in >> tempNum; thisScoreErrorMap[i] += tempNum;
524 qScoreErrorMap[first] = thisScoreErrorMap;
529 //input qualForwardMap
530 for(int i = 0; i < qualForwardMap.size(); i++) {
531 for (int j = 0; j < qualForwardMap[i].size(); j++) {
532 in >> tempNum; qualForwardMap[i][j] += tempNum;
538 //input qualReverseMap
539 for(int i = 0; i < qualReverseMap.size(); i++) {
540 for (int j = 0; j < qualReverseMap[i].size(); j++) {
541 in >> tempNum; qualReverseMap[i][j] += tempNum;
548 for (int i = 0; i < errorForward.size(); i++) {
550 vector<int> thisErrorForward = errorForward[first];
552 for (int i = 0; i < thisErrorForward.size(); i++) {
553 in >> tempNum; thisErrorForward[i] += tempNum;
555 errorForward[first] = thisErrorForward;
561 for (int i = 0; i < errorReverse.size(); i++) {
563 vector<int> thisErrorReverse = errorReverse[first];
565 for (int i = 0; i < thisErrorReverse.size(); i++) {
566 in >> tempNum; thisErrorReverse[i] += tempNum;
568 errorReverse[first] = thisErrorReverse;
573 //input misMatchCounts
575 in >> misMatchSize; m->gobble(in);
576 if (misMatchSize > misMatchCounts.size()) { misMatchCounts.resize(misMatchSize, 0); }
577 for (int j = 0; j < misMatchCounts.size(); j++) {
578 in >> tempNum; misMatchCounts[j] += tempNum;
582 //input megaAlignVector
584 for (int j = 0; j < megaAlignVector.size(); j++) {
585 thisLine = m->getline(in); m->gobble(in); megaAlignVector[j] += thisLine + '\n';
589 in.close(); m->mothurRemove(tempFile);
595 catch(exception& e) {
596 m->errorOut(e, "SeqErrorCommand", "createProcesses");
600 //**********************************************************************************************************************
601 int SeqErrorCommand::driver(string filename, string qFileName, string rFileName, string summaryFileName, string errorOutputFileName, string chimeraOutputFileName, linePair line, linePair qline, linePair rline) {
605 QualityScores quality;
607 misMatchCounts.resize(11, 0);
611 map<string, int>::iterator it;
612 qScoreErrorMap['m'].assign(41, 0);
613 qScoreErrorMap['s'].assign(41, 0);
614 qScoreErrorMap['i'].assign(41, 0);
615 qScoreErrorMap['a'].assign(41, 0);
617 errorForward['m'].assign(maxLength,0);
618 errorForward['s'].assign(maxLength,0);
619 errorForward['i'].assign(maxLength,0);
620 errorForward['d'].assign(maxLength,0);
621 errorForward['a'].assign(maxLength,0);
623 errorReverse['m'].assign(maxLength,0);
624 errorReverse['s'].assign(maxLength,0);
625 errorReverse['i'].assign(maxLength,0);
626 errorReverse['d'].assign(maxLength,0);
627 errorReverse['a'].assign(maxLength,0);
629 //open inputfiles and go to beginning place for this processor
631 m->openInputFile(filename, queryFile);
632 queryFile.seekg(line.start);
636 if(qFileName != "" && rFileName != ""){
637 m->openInputFile(qFileName, qualFile);
638 qualFile.seekg(qline.start);
641 if (rline.start == 0) { report = ReportFile(reportFile, rFileName); }
643 m->openInputFile(rFileName, reportFile);
644 reportFile.seekg(rline.start);
647 qualForwardMap.resize(maxLength);
648 qualReverseMap.resize(maxLength);
649 for(int i=0;i<maxLength;i++){
650 qualForwardMap[i].assign(41,0);
651 qualReverseMap[i].assign(41,0);
655 ofstream outChimeraReport;
656 m->openOutputFile(chimeraOutputFileName, outChimeraReport);
657 RefChimeraTest chimeraTest(referenceSeqs);
658 if (line.start == 0) { chimeraTest.printHeader(outChimeraReport); }
660 ofstream errorSummaryFile;
661 m->openOutputFile(summaryFileName, errorSummaryFile);
662 if (line.start == 0) { printErrorHeader(errorSummaryFile); }
664 ofstream errorSeqFile;
665 m->openOutputFile(errorOutputFileName, errorSeqFile);
667 megaAlignVector.resize(numRefs, "");
675 if (m->control_pressed) { queryFile.close(); if(qFileName != "" && rFileName != ""){ reportFile.close(); qualFile.close(); } outChimeraReport.close(); errorSummaryFile.close();errorSeqFile.close(); return 0; }
677 Sequence query(queryFile);
679 int numParentSeqs = chimeraTest.analyzeQuery(query.getName(), query.getAligned(), outChimeraReport);
680 int closestRefIndex = chimeraTest.getClosestRefIndex();
682 if(numParentSeqs > 1 && ignoreChimeras == 1) { ignoreSeq = 1; }
683 else { ignoreSeq = 0; }
685 Compare minCompare = getErrors(query, referenceSeqs[closestRefIndex]);
687 if(namesFileName != ""){
688 it = weights.find(query.getName());
689 minCompare.weight = it->second;
691 else{ minCompare.weight = 1; }
693 printErrorData(minCompare, numParentSeqs, errorSummaryFile, errorSeqFile);
697 for(int i=0;i<minCompare.sequence.length();i++){
698 char letter = minCompare.sequence[i];
701 errorForward[letter][i] += minCompare.weight;
702 errorReverse[letter][minCompare.total-i-1] += minCompare.weight;
707 if(qualFileName != "" && reportFileName != ""){
708 report = ReportFile(reportFile);
710 // int origLength = report.getQueryLength();
711 int startBase = report.getQueryStart();
712 int endBase = report.getQueryEnd();
714 quality = QualityScores(qualFile);
717 // cout << query.getName() << '\t';
719 quality.updateQScoreErrorMap(qScoreErrorMap, minCompare.sequence, startBase, endBase, minCompare.weight);
720 quality.updateForwardMap(qualForwardMap, startBase, endBase, minCompare.weight);
721 quality.updateReverseMap(qualReverseMap, startBase, endBase, minCompare.weight);
727 if(minCompare.errorRate < threshold && !ignoreSeq){
728 totalBases += (minCompare.total * minCompare.weight);
729 totalMatches += minCompare.matches * minCompare.weight;
730 if(minCompare.mismatches > maxMismatch){
731 maxMismatch = minCompare.mismatches;
732 misMatchCounts.resize(maxMismatch + 1, 0);
734 misMatchCounts[minCompare.mismatches] += minCompare.weight;
737 megaAlignVector[closestRefIndex] += query.getInlineSeq() + '\n';
742 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
743 unsigned long long pos = queryFile.tellg();
744 if ((pos == -1) || (pos >= line.end)) { break; }
746 if (queryFile.eof()) { break; }
749 if(index % 100 == 0){ m->mothurOut(toString(index) + '\n'); }
752 if(qFileName != "" && rFileName != ""){ reportFile.close(); qualFile.close(); }
753 errorSummaryFile.close();
754 errorSeqFile.close();
757 if(index % 100 != 0){ m->mothurOut(toString(index) + '\n'); }
761 catch(exception& e) {
762 m->errorOut(e, "SeqErrorCommand", "driver");
766 //***************************************************************************************************************
768 void SeqErrorCommand::getReferences(){
770 int numAmbigSeqs = 0;
773 int minEndPos = 100000;
775 if (referenceFileName == "saved") {
776 int start = time(NULL);
777 m->mothurOutEndLine(); m->mothurOut("Using sequences from " + rdb->getSavedReference() + " that are saved in memory."); m->mothurOutEndLine();
779 for (int i = 0; i < rdb->referenceSeqs.size(); i++) {
780 int numAmbigs = rdb->referenceSeqs[i].getAmbigBases();
781 if(numAmbigs > 0){ numAmbigSeqs++; }
783 // int startPos = rdb->referenceSeqs[i].getStartPos();
784 // if(startPos > maxStartPos) { maxStartPos = startPos; }
786 // int endPos = rdb->referenceSeqs[i].getEndPos();
787 // if(endPos < minEndPos) { minEndPos = endPos; }
789 referenceSeqs.push_back(rdb->referenceSeqs[i]);
791 referenceFileName = rdb->getSavedReference();
793 m->mothurOut("It took " + toString(time(NULL) - start) + " to load " + toString(rdb->referenceSeqs.size()) + " sequences.");m->mothurOutEndLine();
796 int start = time(NULL);
798 ifstream referenceFile;
799 m->openInputFile(referenceFileName, referenceFile);
801 while(referenceFile){
802 Sequence currentSeq(referenceFile);
803 int numAmbigs = currentSeq.getAmbigBases();
804 if(numAmbigs > 0){ numAmbigSeqs++; }
806 // int startPos = currentSeq.getStartPos();
807 // if(startPos > maxStartPos) { maxStartPos = startPos; }
809 // int endPos = currentSeq.getEndPos();
810 // if(endPos < minEndPos) { minEndPos = endPos; }
811 referenceSeqs.push_back(currentSeq);
813 if (rdb->save) { rdb->referenceSeqs.push_back(currentSeq); }
815 m->gobble(referenceFile);
817 referenceFile.close();
819 m->mothurOut("It took " + toString(time(NULL) - start) + " to read " + toString(referenceSeqs.size()) + " sequences.");m->mothurOutEndLine();
822 numRefs = referenceSeqs.size();
824 for(int i=0;i<numRefs;i++){
825 referenceSeqs[i].padToPos(maxStartPos);
826 referenceSeqs[i].padFromPos(minEndPos);
829 if(numAmbigSeqs != 0){
830 m->mothurOut("Warning: " + toString(numAmbigSeqs) + " reference sequences have ambiguous bases, these bases will be ignored\n");
834 catch(exception& e) {
835 m->errorOut(e, "SeqErrorCommand", "getReferences");
840 //***************************************************************************************************************
842 Compare SeqErrorCommand::getErrors(Sequence query, Sequence reference){
844 if(query.getAlignLength() != reference.getAlignLength()){
845 m->mothurOut("Warning: " + toString(query.getName()) + " and " + toString(reference.getName()) + " are different lengths\n");
847 int alignLength = query.getAlignLength();
849 string q = query.getAligned();
850 string r = reference.getAligned();
855 for(int i=0;i<alignLength;i++){
856 // cout << r[i] << '\t' << q[i] << '\t';
857 if(q[i] != '.' && r[i] != '.' && (q[i] != '-' || r[i] != '-')){ // no missing data and no double gaps
862 if(r[i] == 'A'){ errors.AA++; errors.matches++; errors.sequence += 'm'; }
863 if(r[i] == 'T'){ errors.AT++; errors.sequence += 's'; }
864 if(r[i] == 'G'){ errors.AG++; errors.sequence += 's'; }
865 if(r[i] == 'C'){ errors.AC++; errors.sequence += 's'; }
866 if(r[i] == '-'){ errors.Ai++; errors.sequence += 'i'; }
868 else if(q[i] == 'T'){
869 if(r[i] == 'A'){ errors.TA++; errors.sequence += 's'; }
870 if(r[i] == 'T'){ errors.TT++; errors.matches++; errors.sequence += 'm'; }
871 if(r[i] == 'G'){ errors.TG++; errors.sequence += 's'; }
872 if(r[i] == 'C'){ errors.TC++; errors.sequence += 's'; }
873 if(r[i] == '-'){ errors.Ti++; errors.sequence += 'i'; }
875 else if(q[i] == 'G'){
876 if(r[i] == 'A'){ errors.GA++; errors.sequence += 's'; }
877 if(r[i] == 'T'){ errors.GT++; errors.sequence += 's'; }
878 if(r[i] == 'G'){ errors.GG++; errors.matches++; errors.sequence += 'm'; }
879 if(r[i] == 'C'){ errors.GC++; errors.sequence += 's'; }
880 if(r[i] == '-'){ errors.Gi++; errors.sequence += 'i'; }
882 else if(q[i] == 'C'){
883 if(r[i] == 'A'){ errors.CA++; errors.sequence += 's'; }
884 if(r[i] == 'T'){ errors.CT++; errors.sequence += 's'; }
885 if(r[i] == 'G'){ errors.CG++; errors.sequence += 's'; }
886 if(r[i] == 'C'){ errors.CC++; errors.matches++; errors.sequence += 'm'; }
887 if(r[i] == '-'){ errors.Ci++; errors.sequence += 'i'; }
889 else if(q[i] == 'N'){
890 if(r[i] == 'A'){ errors.NA++; errors.sequence += 'a'; }
891 if(r[i] == 'T'){ errors.NT++; errors.sequence += 'a'; }
892 if(r[i] == 'G'){ errors.NG++; errors.sequence += 'a'; }
893 if(r[i] == 'C'){ errors.NC++; errors.sequence += 'a'; }
894 if(r[i] == '-'){ errors.Ni++; errors.sequence += 'a'; }
896 else if(q[i] == '-' && r[i] != '-'){
897 if(r[i] == 'A'){ errors.dA++; errors.sequence += 'd'; }
898 if(r[i] == 'T'){ errors.dT++; errors.sequence += 'd'; }
899 if(r[i] == 'G'){ errors.dG++; errors.sequence += 'd'; }
900 if(r[i] == 'C'){ errors.dC++; errors.sequence += 'd'; }
907 errors.sequence += 'd'; errors.total++;
910 errors.sequence += 'r';
915 else if(q[i] == '.' && r[i] != '.'){ // reference extends beyond query
916 if(started == 1){ break; }
918 else if(q[i] != '.' && r[i] == '.'){ // query extends beyond reference
919 if(started == 1){ break; }
921 else if(q[i] == '.' && r[i] == '.'){ // both are missing data
922 if(started == 1){ break; }
924 // cout << errors.sequence[errors.sequence.length()-1] << endl;
926 // cout << errors.sequence << endl;
927 errors.mismatches = errors.total-errors.matches;
928 errors.errorRate = (double)(errors.total-errors.matches) / (double)errors.total;
929 errors.queryName = query.getName();
930 errors.refName = reference.getName();
934 catch(exception& e) {
935 m->errorOut(e, "SeqErrorCommand", "getErrors");
940 //***************************************************************************************************************
942 map<string, int> SeqErrorCommand::getWeights(){
944 m->openInputFile(namesFileName, nameFile);
947 string redundantSeqs;
948 map<string, int> nameCountMap;
951 nameFile >> seqName >> redundantSeqs;
952 nameCountMap[seqName] = m->getNumNames(redundantSeqs);
959 //***************************************************************************************************************
961 void SeqErrorCommand::printErrorHeader(ofstream& errorSummaryFile){
963 errorSummaryFile << "query\treference\tweight\t";
964 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";
965 errorSummaryFile << "insertions\tdeletions\tsubstitutions\tambig\tmatches\tmismatches\ttotal\terror\tnumparents\n";
967 errorSummaryFile << setprecision(6);
968 errorSummaryFile.setf(ios::fixed);
970 catch(exception& e) {
971 m->errorOut(e, "SeqErrorCommand", "printErrorHeader");
976 //***************************************************************************************************************
978 void SeqErrorCommand::printErrorData(Compare error, int numParentSeqs, ofstream& errorSummaryFile, ofstream& errorSeqFile){
981 errorSummaryFile << error.queryName << '\t' << error.refName << '\t' << error.weight << '\t';
982 errorSummaryFile << error.AA << '\t' << error.AT << '\t' << error.AG << '\t' << error.AC << '\t';
983 errorSummaryFile << error.TA << '\t' << error.TT << '\t' << error.TG << '\t' << error.TC << '\t';
984 errorSummaryFile << error.GA << '\t' << error.GT << '\t' << error.GG << '\t' << error.GC << '\t';
985 errorSummaryFile << error.CA << '\t' << error.CT << '\t' << error.CG << '\t' << error.CC << '\t';
986 errorSummaryFile << error.NA << '\t' << error.NT << '\t' << error.NG << '\t' << error.NC << '\t';
987 errorSummaryFile << error.Ai << '\t' << error.Ti << '\t' << error.Gi << '\t' << error.Ci << '\t' << error.Ni << '\t';
988 errorSummaryFile << error.dA << '\t' << error.dT << '\t' << error.dG << '\t' << error.dC << '\t';
990 errorSummaryFile << error.Ai + error.Ti + error.Gi + error.Ci << '\t'; //insertions
991 errorSummaryFile << error.dA + error.dT + error.dG + error.dC << '\t'; //deletions
992 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
993 errorSummaryFile << error.NA + error.NT + error.NG + error.NC + error.Ni << '\t'; //ambiguities
994 errorSummaryFile << error.matches << '\t' << error.mismatches << '\t' << error.total << '\t' << error.errorRate << '\t' << numParentSeqs << endl;
996 errorSeqFile << '>' << error.queryName << "\tref:" << error.refName << '\n' << error.sequence << endl;
998 int a=0; int t=1; int g=2; int c=3;
1001 if(numParentSeqs == 1 || ignoreChimeras == 0){
1002 substitutionMatrix[a][a] += error.weight * error.AA;
1003 substitutionMatrix[a][t] += error.weight * error.TA;
1004 substitutionMatrix[a][g] += error.weight * error.GA;
1005 substitutionMatrix[a][c] += error.weight * error.CA;
1006 substitutionMatrix[a][gap] += error.weight * error.dA;
1007 substitutionMatrix[a][n] += error.weight * error.NA;
1009 substitutionMatrix[t][a] += error.weight * error.AT;
1010 substitutionMatrix[t][t] += error.weight * error.TT;
1011 substitutionMatrix[t][g] += error.weight * error.GT;
1012 substitutionMatrix[t][c] += error.weight * error.CT;
1013 substitutionMatrix[t][gap] += error.weight * error.dT;
1014 substitutionMatrix[t][n] += error.weight * error.NT;
1016 substitutionMatrix[g][a] += error.weight * error.AG;
1017 substitutionMatrix[g][t] += error.weight * error.TG;
1018 substitutionMatrix[g][g] += error.weight * error.GG;
1019 substitutionMatrix[g][c] += error.weight * error.CG;
1020 substitutionMatrix[g][gap] += error.weight * error.dG;
1021 substitutionMatrix[g][n] += error.weight * error.NG;
1023 substitutionMatrix[c][a] += error.weight * error.AC;
1024 substitutionMatrix[c][t] += error.weight * error.TC;
1025 substitutionMatrix[c][g] += error.weight * error.GC;
1026 substitutionMatrix[c][c] += error.weight * error.CC;
1027 substitutionMatrix[c][gap] += error.weight * error.dC;
1028 substitutionMatrix[c][n] += error.weight * error.NC;
1030 substitutionMatrix[gap][a] += error.weight * error.Ai;
1031 substitutionMatrix[gap][t] += error.weight * error.Ti;
1032 substitutionMatrix[gap][g] += error.weight * error.Gi;
1033 substitutionMatrix[gap][c] += error.weight * error.Ci;
1034 substitutionMatrix[gap][n] += error.weight * error.Ni;
1037 catch(exception& e) {
1038 m->errorOut(e, "SeqErrorCommand", "printErrorData");
1043 //***************************************************************************************************************
1045 void SeqErrorCommand::printSubMatrix(){
1047 string subMatrixFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.matrix";
1048 ofstream subMatrixFile;
1049 m->openOutputFile(subMatrixFileName, subMatrixFile);
1050 outputNames.push_back(subMatrixFileName); outputTypes["error.matrix"].push_back(subMatrixFileName);
1051 vector<string> bases(6);
1058 vector<int> refSums(5,1);
1060 for(int i=0;i<5;i++){
1061 subMatrixFile << "\tr" << bases[i];
1063 for(int j=0;j<6;j++){
1064 refSums[i] += substitutionMatrix[i][j];
1067 subMatrixFile << endl;
1069 for(int i=0;i<6;i++){
1070 subMatrixFile << 'q' << bases[i];
1071 for(int j=0;j<5;j++){
1072 subMatrixFile << '\t' << substitutionMatrix[j][i];
1074 subMatrixFile << endl;
1077 subMatrixFile << "total";
1078 for(int i=0;i<5;i++){
1079 subMatrixFile << '\t' << refSums[i];
1081 subMatrixFile << endl;
1082 subMatrixFile.close();
1084 catch(exception& e) {
1085 m->errorOut(e, "SeqErrorCommand", "printSubMatrix");
1089 //***************************************************************************************************************
1091 void SeqErrorCommand::printErrorFRFile(map<char, vector<int> > errorForward, map<char, vector<int> > errorReverse){
1093 string errorForwardFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.seq.forward";
1094 ofstream errorForwardFile;
1095 m->openOutputFile(errorForwardFileName, errorForwardFile);
1096 outputNames.push_back(errorForwardFileName); outputTypes["error.forward"].push_back(errorForwardFileName);
1098 errorForwardFile << "position\ttotalseqs\tmatch\tsubstitution\tinsertion\tdeletion\tambiguous" << endl;
1099 for(int i=0;i<maxLength;i++){
1100 float match = (float)errorForward['m'][i];
1101 float subst = (float)errorForward['s'][i];
1102 float insert = (float)errorForward['i'][i];
1103 float del = (float)errorForward['d'][i];
1104 float amb = (float)errorForward['a'][i];
1105 float total = match + subst + insert + del + amb;
1106 if(total == 0){ break; }
1107 errorForwardFile << i+1 << '\t' << total << '\t' << match/total << '\t' << subst/total << '\t' << insert/total << '\t' << del/total << '\t' << amb/total << endl;
1109 errorForwardFile.close();
1111 string errorReverseFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.seq.reverse";
1112 ofstream errorReverseFile;
1113 m->openOutputFile(errorReverseFileName, errorReverseFile);
1114 outputNames.push_back(errorReverseFileName); outputTypes["error.reverse"].push_back(errorReverseFileName);
1116 errorReverseFile << "position\ttotalseqs\tmatch\tsubstitution\tinsertion\tdeletion\tambiguous" << endl;
1117 for(int i=0;i<maxLength;i++){
1118 float match = (float)errorReverse['m'][i];
1119 float subst = (float)errorReverse['s'][i];
1120 float insert = (float)errorReverse['i'][i];
1121 float del = (float)errorReverse['d'][i];
1122 float amb = (float)errorReverse['a'][i];
1123 float total = match + subst + insert + del + amb;
1124 if(total == 0){ break; }
1125 errorReverseFile << i+1 << '\t' << total << '\t' << match/total << '\t' << subst/total << '\t' << insert/total << '\t' << del/total << '\t' << amb/total << endl;
1127 errorReverseFile.close();
1129 catch(exception& e) {
1130 m->errorOut(e, "SeqErrorCommand", "printErrorFRFile");
1135 //***************************************************************************************************************
1137 void SeqErrorCommand::printErrorQuality(map<char, vector<int> > qScoreErrorMap){
1140 string errorQualityFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.quality";
1141 ofstream errorQualityFile;
1142 m->openOutputFile(errorQualityFileName, errorQualityFile);
1143 outputNames.push_back(errorQualityFileName); outputTypes["error.quality"].push_back(errorQualityFileName);
1145 errorQualityFile << "qscore\tmatches\tsubstitutions\tinsertions\tambiguous" << endl;
1146 for(int i=0;i<41;i++){
1147 errorQualityFile << i << '\t' << qScoreErrorMap['m'][i] << '\t' << qScoreErrorMap['s'][i] << '\t' << qScoreErrorMap['i'][i] << '\t'<< qScoreErrorMap['a'][i] << endl;
1149 errorQualityFile.close();
1151 catch(exception& e) {
1152 m->errorOut(e, "SeqErrorCommand", "printErrorFRFile");
1158 //***************************************************************************************************************
1160 void SeqErrorCommand::printQualityFR(vector<vector<int> > qualForwardMap, vector<vector<int> > qualReverseMap){
1164 int numColumns = qualForwardMap[0].size();
1166 for(int i=0;i<qualForwardMap.size();i++){
1167 for(int j=0;j<numColumns;j++){
1168 if(qualForwardMap[i][j] != 0){
1169 if(numRows < i) { numRows = i+20; }
1174 string qualityForwardFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.qual.forward";
1175 ofstream qualityForwardFile;
1176 m->openOutputFile(qualityForwardFileName, qualityForwardFile);
1177 outputNames.push_back(qualityForwardFileName); outputTypes["error.qual.forward"].push_back(qualityForwardFileName);
1179 for(int i=0;i<numColumns;i++){ qualityForwardFile << '\t' << i; } qualityForwardFile << endl;
1181 for(int i=0;i<numRows;i++){
1182 qualityForwardFile << i+1;
1183 for(int j=0;j<numColumns;j++){
1184 qualityForwardFile << '\t' << qualForwardMap[i][j];
1187 qualityForwardFile << endl;
1189 qualityForwardFile.close();
1192 string qualityReverseFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.qual.reverse";
1193 ofstream qualityReverseFile;
1194 m->openOutputFile(qualityReverseFileName, qualityReverseFile);
1195 outputNames.push_back(qualityReverseFileName); outputTypes["error.qual.reverse"].push_back(qualityReverseFileName);
1197 for(int i=0;i<numColumns;i++){ qualityReverseFile << '\t' << i; } qualityReverseFile << endl;
1198 for(int i=0;i<numRows;i++){
1200 qualityReverseFile << i+1;
1201 for(int j=0;j<numColumns;j++){
1202 qualityReverseFile << '\t' << qualReverseMap[i][j];
1204 qualityReverseFile << endl;
1206 qualityReverseFile.close();
1208 catch(exception& e) {
1209 m->errorOut(e, "SeqErrorCommand", "printErrorFRFile");
1214 /**************************************************************************************************/
1216 int SeqErrorCommand::setLines(string filename, string qfilename, string rfilename, vector<unsigned long long>& fastaFilePos, vector<unsigned long long>& qfileFilePos, vector<unsigned long long>& rfileFilePos) {
1218 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
1219 //set file positions for fasta file
1220 fastaFilePos = m->divideFile(filename, processors);
1222 if (qfilename == "") { return processors; }
1224 //get name of first sequence in each chunk
1225 map<string, int> firstSeqNames;
1226 for (int i = 0; i < (fastaFilePos.size()-1); i++) {
1228 m->openInputFile(filename, in);
1229 in.seekg(fastaFilePos[i]);
1232 firstSeqNames[temp.getName()] = i;
1237 //make copy to use below
1238 map<string, int> firstSeqNamesReport = firstSeqNames;
1240 //seach for filePos of each first name in the qfile and save in qfileFilePos
1242 m->openInputFile(qfilename, inQual);
1245 while(!inQual.eof()){
1246 input = m->getline(inQual);
1248 if (input.length() != 0) {
1249 if(input[0] == '>'){ //this is a sequence name line
1250 istringstream nameStream(input);
1252 string sname = ""; nameStream >> sname;
1253 sname = sname.substr(1);
1255 map<string, int>::iterator it = firstSeqNames.find(sname);
1257 if(it != firstSeqNames.end()) { //this is the start of a new chunk
1258 unsigned long long pos = inQual.tellg();
1259 qfileFilePos.push_back(pos - input.length() - 1);
1260 firstSeqNames.erase(it);
1265 if (firstSeqNames.size() == 0) { break; }
1269 if (firstSeqNames.size() != 0) {
1270 for (map<string, int>::iterator it = firstSeqNames.begin(); it != firstSeqNames.end(); it++) {
1271 m->mothurOut(it->first + " is in your fasta file and not in your quality file, aborting."); m->mothurOutEndLine();
1273 m->control_pressed = true;
1277 //get last file position of qfile
1279 unsigned long long size;
1281 //get num bytes in file
1282 pFile = fopen (qfilename.c_str(),"rb");
1283 if (pFile==NULL) perror ("Error opening file");
1285 fseek (pFile, 0, SEEK_END);
1290 qfileFilePos.push_back(size);
1292 //seach for filePos of each first name in the rfile and save in rfileFilePos
1295 m->openInputFile(rfilename, inR);
1297 //read column headers
1298 for (int i = 0; i < 16; i++) {
1299 if (!inR.eof()) { inR >> junk; }
1305 if (m->control_pressed) { inR.close(); return processors; }
1307 input = m->getline(inR);
1309 if (input.length() != 0) {
1311 istringstream nameStream(input);
1312 string sname = ""; nameStream >> sname;
1314 map<string, int>::iterator it = firstSeqNamesReport.find(sname);
1316 if(it != firstSeqNamesReport.end()) { //this is the start of a new chunk
1317 unsigned long long pos = inR.tellg();
1318 rfileFilePos.push_back(pos - input.length() - 1);
1319 firstSeqNamesReport.erase(it);
1323 if (firstSeqNamesReport.size() == 0) { break; }
1328 if (firstSeqNamesReport.size() != 0) {
1329 for (map<string, int>::iterator it = firstSeqNamesReport.begin(); it != firstSeqNamesReport.end(); it++) {
1330 m->mothurOut(it->first + " is in your fasta file and not in your report file, aborting."); m->mothurOutEndLine();
1332 m->control_pressed = true;
1336 //get last file position of qfile
1338 unsigned long long sizeR;
1340 //get num bytes in file
1341 rFile = fopen (rfilename.c_str(),"rb");
1342 if (rFile==NULL) perror ("Error opening file");
1344 fseek (rFile, 0, SEEK_END);
1345 sizeR=ftell (rFile);
1349 rfileFilePos.push_back(sizeR);
1355 fastaFilePos.push_back(0); qfileFilePos.push_back(0);
1356 //get last file position of fastafile
1358 unsigned long long size;
1360 //get num bytes in file
1361 pFile = fopen (filename.c_str(),"rb");
1362 if (pFile==NULL) perror ("Error opening file");
1364 fseek (pFile, 0, SEEK_END);
1368 fastaFilePos.push_back(size);
1370 //get last file position of fastafile
1373 //get num bytes in file
1374 qFile = fopen (qfilename.c_str(),"rb");
1375 if (qFile==NULL) perror ("Error opening file");
1377 fseek (qFile, 0, SEEK_END);
1381 qfileFilePos.push_back(size);
1387 catch(exception& e) {
1388 m->errorOut(e, "SeqErrorCommand", "setLines");
1392 //***************************************************************************************************************