X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=seqerrorcommand.cpp;h=392f5877d9099308593d0cba26a8ff7baddd1d3b;hb=55ec7cde88d5512e177fe9488d5ee13793853bad;hp=00f78626b6270455f5a8f35464923fe72e9ff527;hpb=c69e2e9749626cfbf1d6cb0125ae94f869e00b18;p=mothur.git

diff --git a/seqerrorcommand.cpp b/seqerrorcommand.cpp
index 00f7862..392f587 100644
--- a/seqerrorcommand.cpp
+++ b/seqerrorcommand.cpp
@@ -13,14 +13,39 @@
 #include "refchimeratest.h"
 
 //**********************************************************************************************************************
-vector<string> SeqErrorCommand::getValidParameters(){	
+vector<string> SeqErrorCommand::setParameters(){	
 	try {
-		string Array[] =  {"query", "reference", "name", "qfile", "report", "threshold", "inputdir", "outputdir"};
-		vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
+		CommandParameter pquery("fasta", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(pquery);
+		CommandParameter preference("reference", "InputTypes", "", "", "none", "none", "none",false,true); parameters.push_back(preference);
+		CommandParameter pqfile("qfile", "InputTypes", "", "", "none", "none", "QualReport",false,false); parameters.push_back(pqfile);
+		CommandParameter preport("report", "InputTypes", "", "", "none", "none", "QualReport",false,false); parameters.push_back(preport);
+		CommandParameter pname("name", "InputTypes", "", "", "none", "none", "none",false,false); parameters.push_back(pname);
+		CommandParameter pignorechimeras("ignorechimeras", "Boolean", "", "T", "", "", "",false,false); parameters.push_back(pignorechimeras);
+		CommandParameter pthreshold("threshold", "Number", "", "1.0", "", "", "",false,false); parameters.push_back(pthreshold);
+		CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
+		CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
+		
+		vector<string> myArray;
+		for (int i = 0; i < parameters.size(); i++) {	myArray.push_back(parameters[i].name);		}
 		return myArray;
 	}
 	catch(exception& e) {
-		m->errorOut(e, "SeqErrorCommand", "getValidParameters");
+		m->errorOut(e, "SeqErrorCommand", "setParameters");
+		exit(1);
+	}
+}
+//**********************************************************************************************************************
+string SeqErrorCommand::getHelpString(){	
+	try {
+		string helpString = "";
+		helpString += "The seq.error command reads a query alignment file and a reference alignment file and creates .....\n";
+		helpString += "Example seq.error(...).\n";
+		helpString += "Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFasta).\n";
+		helpString += "For more details please check out the wiki http://www.mothur.org/wiki/seq.error .\n";
+		return helpString;
+	}
+	catch(exception& e) {
+		m->errorOut(e, "SeqErrorCommand", "getHelpString");
 		exit(1);
 	}
 }
@@ -44,29 +69,6 @@ SeqErrorCommand::SeqErrorCommand(){
 		exit(1);
 	}
 }
-//**********************************************************************************************************************
-vector<string> SeqErrorCommand::getRequiredParameters(){	
-	try {
-		string Array[] =  {"query","reference"};
-		vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
-		return myArray;
-	}
-	catch(exception& e) {
-		m->errorOut(e, "SeqErrorCommand", "getRequiredParameters");
-		exit(1);
-	}
-}
-//**********************************************************************************************************************
-vector<string> SeqErrorCommand::getRequiredFiles(){	
-	try {
-		vector<string> myArray;
-		return myArray;
-	}
-	catch(exception& e) {
-		m->errorOut(e, "SeqErrorCommand", "getRequiredFiles");
-		exit(1);
-	}
-}
 //***************************************************************************************************************
 
 SeqErrorCommand::SeqErrorCommand(string option)  {
@@ -79,11 +81,7 @@ SeqErrorCommand::SeqErrorCommand(string option)  {
 		
 		else {
 			string temp;
-			
-			//valid paramters for this command
-			string AlignArray[] =  {"query", "reference", "name", "qfile", "report", "threshold", "inputdir", "ignorechimeras", "outputdir"};
-			
-			vector<string> myArray (AlignArray, AlignArray+(sizeof(AlignArray)/sizeof(string)));
+			vector<string> myArray = setParameters();
 			
 			OptionParser parser(option);
 			map<string,string> parameters = parser.getParameters();
@@ -114,12 +112,12 @@ SeqErrorCommand::SeqErrorCommand(string option)  {
 			if (inputDir == "not found"){	inputDir = "";		}
 			else {
 				string path;
-				it = parameters.find("query");
+				it = parameters.find("fasta");
 				//user has given a template file
 				if(it != parameters.end()){ 
 					path = m->hasPath(it->second);
 					//if the user has not given a path then, add inputdir. else leave path alone.
-					if (path == "") {	parameters["query"] = inputDir + it->second;		}
+					if (path == "") {	parameters["fasta"] = inputDir + it->second;		}
 				}
 				
 				it = parameters.find("reference");
@@ -156,8 +154,12 @@ SeqErrorCommand::SeqErrorCommand(string option)  {
 				
 			}
 			//check for required parameters
-			queryFileName = validParameter.validFile(parameters, "query", true);
-			if (queryFileName == "not found") { m->mothurOut("query is a required parameter for the seq.error command."); m->mothurOutEndLine(); abort = true; }
+			queryFileName = validParameter.validFile(parameters, "fasta", true);
+			if (queryFileName == "not found") { 
+				queryFileName = m->getFastaFile(); 
+				if (queryFileName != "") { m->mothurOut("Using " + queryFileName + " as input file for the fasta parameter."); m->mothurOutEndLine(); }
+				else { 	m->mothurOut("You have no current fasta file and the fasta parameter is required."); m->mothurOutEndLine(); abort = true; }
+			}
 			else if (queryFileName == "not open") { abort = true; }	
 			
 			referenceFileName = validParameter.validFile(parameters, "reference", true);
@@ -168,12 +170,15 @@ SeqErrorCommand::SeqErrorCommand(string option)  {
 			//check for optional parameters
 			namesFileName = validParameter.validFile(parameters, "name", true);
 			if(namesFileName == "not found"){	namesFileName = "";	}
+			else if (namesFileName == "not open") { namesFileName = ""; abort = true; }	
 			
 			qualFileName = validParameter.validFile(parameters, "qfile", true);
 			if(qualFileName == "not found"){	qualFileName = "";	}
+			else if (qualFileName == "not open") { qualFileName = ""; abort = true; }	
 
 			reportFileName = validParameter.validFile(parameters, "report", true);
 			if(reportFileName == "not found"){	reportFileName = "";	}
+			else if (reportFileName == "not open") { reportFileName = ""; abort = true; }	
 			
 			if((reportFileName != "" && qualFileName == "") || (reportFileName == "" && qualFileName != "")){
 				m->mothurOut("if you use either a qual file or a report file, you have to have both.");
@@ -196,7 +201,7 @@ SeqErrorCommand::SeqErrorCommand(string option)  {
 			convert(temp, ignoreChimeras);  
 
 			substitutionMatrix.resize(6);
-			for(int i=0;i<6;i++){	substitutionMatrix[i].assign(6,0);	}
+			for(int i=0;i<6;i++){	substitutionMatrix[i].resize(6,0);	}
 		}
 	}
 	catch(exception& e) {
@@ -204,34 +209,14 @@ SeqErrorCommand::SeqErrorCommand(string option)  {
 		exit(1);
 	}
 }
-
-//**********************************************************************************************************************
-
-void SeqErrorCommand::help(){
-	try {
-		m->mothurOut("The seq.error command reads a query alignment file and a reference alignment file and creates .....\n");
-		m->mothurOut("Example seq.error(...).\n");
-		m->mothurOut("Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFasta).\n");
-		m->mothurOut("For more details please check out the wiki http://www.mothur.org/wiki/seq.error .\n\n");
-	}
-	catch(exception& e) {
-		m->errorOut(e, "SeqErrorCommand", "help");
-		exit(1);
-	}
-}
-
-//***************************************************************************************************************
-
-SeqErrorCommand::~SeqErrorCommand(){
-
-}
-
 //***************************************************************************************************************
 
 int SeqErrorCommand::execute(){
 	try{
 		if (abort == true) { if (calledHelp) { return 0; }  return 2;	}
 
+		maxLength = 2000;
+		
 		string errorSummaryFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.summary";
 		m->openOutputFile(errorSummaryFileName, errorSummaryFile);
 		outputNames.push_back(errorSummaryFileName); outputTypes["error.summary"].push_back(errorSummaryFileName);
@@ -261,11 +246,11 @@ int SeqErrorCommand::execute(){
 			m->openInputFile(qualFileName, qualFile);
 			report = ReportFile(reportFile, reportFileName);
 			
-			qualForwardMap.resize(1000);
-			qualReverseMap.resize(1000);
-			for(int i=0;i<1000;i++){
-				qualForwardMap[i].assign(100,0);
-				qualReverseMap[i].assign(100,0);
+			qualForwardMap.resize(maxLength);
+			qualReverseMap.resize(maxLength);
+			for(int i=0;i<maxLength;i++){
+				qualForwardMap[i].assign(41,0);
+				qualReverseMap[i].assign(41,0);
 			}				
 		}
 		
@@ -284,66 +269,69 @@ int SeqErrorCommand::execute(){
 		qScoreErrorMap['a'].assign(41, 0);
 		
 		map<char, vector<int> > errorForward;
-		errorForward['m'].assign(1000,0);
-		errorForward['s'].assign(1000,0);
-		errorForward['i'].assign(1000,0);
-		errorForward['d'].assign(1000,0);
-		errorForward['a'].assign(1000,0);
+		errorForward['m'].assign(maxLength,0);
+		errorForward['s'].assign(maxLength,0);
+		errorForward['i'].assign(maxLength,0);
+		errorForward['d'].assign(maxLength,0);
+		errorForward['a'].assign(maxLength,0);
 		
 		map<char, vector<int> > errorReverse;
-		errorReverse['m'].assign(1000,0);
-		errorReverse['s'].assign(1000,0);
-		errorReverse['i'].assign(1000,0);
-		errorReverse['d'].assign(1000,0);
-		errorReverse['a'].assign(1000,0);	
+		errorReverse['m'].assign(maxLength,0);
+		errorReverse['s'].assign(maxLength,0);
+		errorReverse['i'].assign(maxLength,0);
+		errorReverse['d'].assign(maxLength,0);
+		errorReverse['a'].assign(maxLength,0);	
 		
 
 		string errorChimeraFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.chimera";
 		RefChimeraTest chimeraTest(referenceSeqs, errorChimeraFileName);
 		outputNames.push_back(errorChimeraFileName); outputTypes["error.chimera"].push_back(errorChimeraFileName);
 		
+		vector<string> megaAlignVector(numRefs, "");
+
 		int index = 0;
 		bool ignoreSeq = 0;
 		
 		while(queryFile){
-			
+
 			if (m->control_pressed) { errorSummaryFile.close();	errorSeqFile.close(); for (int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); } return 0; }
 		
 			Sequence query(queryFile);
-						
+			
 			int numParentSeqs = chimeraTest.analyzeQuery(query.getName(), query.getAligned());
 			int closestRefIndex = chimeraTest.getClosestRefIndex();
 
 			if(numParentSeqs > 1 && ignoreChimeras == 1)	{	ignoreSeq = 1;	}
 			else											{	ignoreSeq = 0;	}
 
-
 			Compare minCompare = getErrors(query, referenceSeqs[closestRefIndex]);
 			
 			if(namesFileName != ""){
 				it = weights.find(query.getName());
 				minCompare.weight = it->second;
 			}
-			else	{	minCompare.weight = 1;	}
+			else{	minCompare.weight = 1;	}
 
 			printErrorData(minCompare, numParentSeqs);
-
+		
 			if(!ignoreSeq){
-				for(int i=0;i<minCompare.total;i++){
+
+				for(int i=0;i<minCompare.sequence.length();i++){
 					char letter = minCompare.sequence[i];
+
 					errorForward[letter][i] += minCompare.weight;
 					errorReverse[letter][minCompare.total-i-1] += minCompare.weight;				
 				}
 			}
-			
+
 			if(qualFileName != "" && reportFileName != ""){
 				report = ReportFile(reportFile);
 				
-				int origLength = report.getQueryLength();
+//				int origLength = report.getQueryLength();
 				int startBase = report.getQueryStart();
 				int endBase = report.getQueryEnd();
 
-				quality = QualityScores(qualFile, origLength);
+				quality = QualityScores(qualFile);
 
 				if(!ignoreSeq){
 					quality.updateQScoreErrorMap(qScoreErrorMap, minCompare.sequence, startBase, endBase, minCompare.weight);
@@ -351,7 +339,7 @@ int SeqErrorCommand::execute(){
 					quality.updateReverseMap(qualReverseMap, startBase, endBase, minCompare.weight);
 				}
 			}			
-			
+
 			if(minCompare.errorRate < threshold && !ignoreSeq){
 				totalBases += (minCompare.total * minCompare.weight);
 				totalMatches += minCompare.matches * minCompare.weight;
@@ -361,10 +349,13 @@ int SeqErrorCommand::execute(){
 				}				
 				misMatchCounts[minCompare.mismatches] += minCompare.weight;
 				numSeqs++;
+				
+				megaAlignVector[closestRefIndex] += query.getInlineSeq() + '\n';
 			}
-			
+
 			index++;
-			if(index % 1000 == 0){	cout << index << endl;	}
+			
+			if(index % 1000 == 0){	m->mothurOut(toString(index) + '\n');	}
 		}
 		queryFile.close();
 		errorSummaryFile.close();	
@@ -383,7 +374,7 @@ int SeqErrorCommand::execute(){
 		ofstream errorCountFile;
 		m->openOutputFile(errorCountFileName, errorCountFile);
 		outputNames.push_back(errorCountFileName);  outputTypes["error.count"].push_back(errorCountFileName);
-		m->mothurOut("Overall error rate:\t" + toString((double)(totalBases - totalMatches) / (double)totalBases) + "\n\n");
+		m->mothurOut("Overall error rate:\t" + toString((double)(totalBases - totalMatches) / (double)totalBases) + "\n");
 		m->mothurOut("Errors\tSequences\n");
 		errorCountFile << "Errors\tSequences\n";		
 		for(int i=0;i<misMatchCounts.size();i++){
@@ -396,6 +387,17 @@ int SeqErrorCommand::execute(){
 
 		printSubMatrix();
 				
+		string megAlignmentFileName = queryFileName.substr(0,queryFileName.find_last_of('.')) + ".error.ref-query";
+		ofstream megAlignmentFile;
+		m->openOutputFile(megAlignmentFileName, megAlignmentFile);
+		outputNames.push_back(megAlignmentFileName);  outputTypes["error.ref-query"].push_back(megAlignmentFileName);
+		
+		for(int i=0;i<numRefs;i++){
+			megAlignmentFile << referenceSeqs[i].getInlineSeq() << endl;
+			megAlignmentFile << megaAlignVector[i] << endl;
+		}
+		
+		
 		m->mothurOutEndLine();
 		m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 		for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
@@ -417,20 +419,38 @@ void SeqErrorCommand::getReferences(){
 		ifstream referenceFile;
 		m->openInputFile(referenceFileName, referenceFile);
 		
+		int numAmbigSeqs = 0;
+		
+		int maxStartPos = 0;
+		int minEndPos = 100000;
+		
 		while(referenceFile){
 			Sequence currentSeq(referenceFile);
 			int numAmbigs = currentSeq.getAmbigBases();
+			if(numAmbigs > 0){	numAmbigSeqs++;	}
 			
-			if(numAmbigs != 0){
-				m->mothurOut("Warning: " + toString(currentSeq.getName()) + " has " + toString(numAmbigs) + " ambiguous bases, these bases will be removed\n");
-				currentSeq.removeAmbigBases();
-			}
+//			int startPos = currentSeq.getStartPos();
+//			if(startPos > maxStartPos)	{	maxStartPos = startPos;	}
+//
+//			int endPos = currentSeq.getEndPos();
+//			if(endPos < minEndPos)		{	minEndPos = endPos;		}
 			referenceSeqs.push_back(currentSeq);
+				
 			m->gobble(referenceFile);
 		}
+		referenceFile.close();
 		numRefs = referenceSeqs.size();
+
+		
+		for(int i=0;i<numRefs;i++){
+			referenceSeqs[i].padToPos(maxStartPos);
+			referenceSeqs[i].padFromPos(minEndPos);
+		}
+		
+		if(numAmbigSeqs != 0){
+			m->mothurOut("Warning: " + toString(numAmbigSeqs) + " reference sequences have ambiguous bases, these bases will be ignored\n");
+		}		
 		
-		referenceFile.close();
 	}
 	catch(exception& e) {
 		m->errorOut(e, "SeqErrorCommand", "getReferences");
@@ -454,7 +474,7 @@ Compare SeqErrorCommand::getErrors(Sequence query, Sequence reference){
 		Compare errors;
 
 		for(int i=0;i<alignLength;i++){
-			if(q[i] != '.' && r[i] != '.' && (q[i] != '-' || r[i] != '-')){			//	no missing data and no double gaps
+			if(r[i] != 'N' && q[i] != '.' && r[i] != '.' && (q[i] != '-' || r[i] != '-')){			//	no missing data and no double gaps
 				started = 1;
 				
 				if(q[i] == 'A'){
@@ -505,7 +525,6 @@ Compare SeqErrorCommand::getErrors(Sequence query, Sequence reference){
 				if(started == 1){	break;	}
 			}
 			else if(q[i] != '.' && r[i] == '.'){		//	query extends beyond reference
-				m->mothurOut("Warning: " + toString(query.getName()) + " extend beyond " + toString(reference.getName()) + ".  Ignoring the extra bases in the query\n");
 				if(started == 1){	break;	}
 			}
 			else if(q[i] == '.' && r[i] == '.'){		//	both are missing data
@@ -566,13 +585,14 @@ void SeqErrorCommand::printErrorHeader(){
 
 void SeqErrorCommand::printErrorData(Compare error, int numParentSeqs){
 	try {
+
 		errorSummaryFile << error.queryName << '\t' << error.refName << '\t' << error.weight << '\t';
 		errorSummaryFile << error.AA << '\t' << error.AT << '\t' << error.AG << '\t' << error.AC << '\t';
 		errorSummaryFile << error.TA << '\t' << error.TT << '\t' << error.TG << '\t' << error.TC << '\t';
 		errorSummaryFile << error.GA << '\t' << error.GT << '\t' << error.GG << '\t' << error.GC << '\t';
 		errorSummaryFile << error.CA << '\t' << error.CT << '\t' << error.CG << '\t' << error.CC << '\t';
 		errorSummaryFile << error.NA << '\t' << error.NT << '\t' << error.NG << '\t' << error.NC << '\t';
-		errorSummaryFile << error.Ai << '\t' << error.Ti << '\t' << error.Gi << '\t' << error.Ci << '\t' << error.Ni << '\t' ;
+		errorSummaryFile << error.Ai << '\t' << error.Ti << '\t' << error.Gi << '\t' << error.Ci << '\t' << error.Ni << '\t';
 		errorSummaryFile << error.dA << '\t' << error.dT << '\t' << error.dG << '\t' << error.dC << '\t';
 		
 		errorSummaryFile << error.Ai + error.Ti + error.Gi + error.Ci << '\t';			//insertions
@@ -580,12 +600,12 @@ void SeqErrorCommand::printErrorData(Compare error, int numParentSeqs){
 		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
 		errorSummaryFile << error.NA + error.NT + error.NG + error.NC + error.Ni << '\t';	//ambiguities
 		errorSummaryFile << error.matches << '\t' << error.mismatches << '\t' << error.total << '\t' << error.errorRate << '\t' << numParentSeqs << endl;
-		
+
 		errorSeqFile << '>' << error.queryName << "\tref:" << error.refName << '\n' << error.sequence << endl;
 		
-		
 		int a=0;		int t=1;		int g=2;		int c=3;
 		int gap=4;		int n=5;
+
 		if(numParentSeqs == 1 || ignoreChimeras == 0){
 			substitutionMatrix[a][a] += error.weight * error.AA;
 			substitutionMatrix[a][t] += error.weight * error.TA;
@@ -593,7 +613,7 @@ void SeqErrorCommand::printErrorData(Compare error, int numParentSeqs){
 			substitutionMatrix[a][c] += error.weight * error.CA;
 			substitutionMatrix[a][gap] += error.weight * error.dA;
 			substitutionMatrix[a][n] += error.weight * error.NA;
-
+			
 			substitutionMatrix[t][a] += error.weight * error.AT;
 			substitutionMatrix[t][t] += error.weight * error.TT;
 			substitutionMatrix[t][g] += error.weight * error.GT;
@@ -684,7 +704,7 @@ void SeqErrorCommand::printErrorFRFile(map<char, vector<int> > errorForward, map
 		outputNames.push_back(errorForwardFileName);  outputTypes["error.forward"].push_back(errorForwardFileName);
 
 		errorForwardFile << "position\ttotalseqs\tmatch\tsubstitution\tinsertion\tdeletion\tambiguous" << endl;
-		for(int i=0;i<1000;i++){
+		for(int i=0;i<maxLength;i++){
 			float match = (float)errorForward['m'][i];
 			float subst = (float)errorForward['s'][i];
 			float insert = (float)errorForward['i'][i];
@@ -702,7 +722,7 @@ void SeqErrorCommand::printErrorFRFile(map<char, vector<int> > errorForward, map
 		outputNames.push_back(errorReverseFileName);  outputTypes["error.reverse"].push_back(errorReverseFileName);
 
 		errorReverseFile << "position\ttotalseqs\tmatch\tsubstitution\tinsertion\tdeletion\tambiguous" << endl;
-		for(int i=0;i<1000;i++){
+		for(int i=0;i<maxLength;i++){
 			float match = (float)errorReverse['m'][i];
 			float subst = (float)errorReverse['s'][i];
 			float insert = (float)errorReverse['i'][i];
@@ -746,17 +766,15 @@ void SeqErrorCommand::printErrorQuality(map<char, vector<int> > qScoreErrorMap){
 //***************************************************************************************************************
 
 void SeqErrorCommand::printQualityFR(vector<vector<int> > qualForwardMap, vector<vector<int> > qualReverseMap){
-	try{
-
 
-		int lastRow = 0;
-		int lastColumn = 0;
+	try{
+		int numRows = 0;
+		int numColumns = qualForwardMap[0].size();
 
 		for(int i=0;i<qualForwardMap.size();i++){
-			for(int j=0;j<qualForwardMap[i].size();j++){
+			for(int j=0;j<numColumns;j++){
 				if(qualForwardMap[i][j] != 0){
-					if(lastRow < i)		{	lastRow = i+2;		}
-					if(lastColumn < j)	{	lastColumn = j+2;	}
+					if(numRows < i)		{	numRows = i+20;		}
 				}
 			}
 		}
@@ -766,11 +784,11 @@ void SeqErrorCommand::printQualityFR(vector<vector<int> > qualForwardMap, vector
 		m->openOutputFile(qualityForwardFileName, qualityForwardFile);
 		outputNames.push_back(qualityForwardFileName);  outputTypes["error.qual.forward"].push_back(qualityForwardFileName);
 
-		for(int i=0;i<lastColumn;i++){	qualityForwardFile << '\t' << i;	}	qualityForwardFile << endl;
+		for(int i=0;i<numColumns;i++){	qualityForwardFile << '\t' << i;	}	qualityForwardFile << endl;
 
-		for(int i=0;i<lastRow;i++){
+		for(int i=0;i<numRows;i++){
 			qualityForwardFile << i+1;
-			for(int j=0;j<lastColumn;j++){
+			for(int j=0;j<numColumns;j++){
 				qualityForwardFile << '\t' << qualForwardMap[i][j];
 			}
 
@@ -784,11 +802,11 @@ void SeqErrorCommand::printQualityFR(vector<vector<int> > qualForwardMap, vector
 		m->openOutputFile(qualityReverseFileName, qualityReverseFile);
 		outputNames.push_back(qualityReverseFileName);  outputTypes["error.qual.reverse"].push_back(qualityReverseFileName);
 		
-		for(int i=0;i<lastColumn;i++){	qualityReverseFile << '\t' << i;	}	qualityReverseFile << endl;
-		for(int i=0;i<lastRow;i++){
+		for(int i=0;i<numColumns;i++){	qualityReverseFile << '\t' << i;	}	qualityReverseFile << endl;
+		for(int i=0;i<numRows;i++){
 			
 			qualityReverseFile << i+1;
-			for(int j=0;j<lastColumn;j++){
+			for(int j=0;j<numColumns;j++){
 				qualityReverseFile << '\t' << qualReverseMap[i][j];
 			}
 			qualityReverseFile << endl;
@@ -799,7 +817,7 @@ void SeqErrorCommand::printQualityFR(vector<vector<int> > qualForwardMap, vector
 		m->errorOut(e, "SeqErrorCommand", "printErrorFRFile");
 		exit(1);
 	}
+	
 }
 
-
 //***************************************************************************************************************