X-Git-Url: https://git.donarmstrong.com/?p=mothur.git;a=blobdiff_plain;f=seqerrorcommand.cpp;h=5879241372dd3b60b34683683a0c1edbcbd64d57;hp=8fd6368154773bcc38e2b950d7b170c88be2186e;hb=615301e57c25e241356a9c2380648d117709458d;hpb=49d2b7459c5027557564b21e9487dadafbbbdc96

diff --git a/seqerrorcommand.cpp b/seqerrorcommand.cpp
index 8fd6368..5879241 100644
--- a/seqerrorcommand.cpp
+++ b/seqerrorcommand.cpp
@@ -11,23 +11,27 @@
 #include "reportfile.h"
 #include "qualityscores.h"
 #include "refchimeratest.h"
+#include "myPerseus.h"
 #include "filterseqscommand.h"
 
 
 //**********************************************************************************************************************
+
 vector<string> SeqErrorCommand::setParameters(){	
 	try {
-		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 pprocessors("processors", "Number", "", "1", "", "", "",false,false); parameters.push_back(pprocessors);
-		CommandParameter psave("save", "Boolean", "", "F", "", "", "",false,false); parameters.push_back(psave);
-		CommandParameter pinputdir("inputdir", "String", "", "", "", "", "",false,false); parameters.push_back(pinputdir);
-		CommandParameter poutputdir("outputdir", "String", "", "", "", "", "",false,false); parameters.push_back(poutputdir);
+		CommandParameter pquery("fasta", "InputTypes", "", "", "none", "none", "none","errorType",false,true,true); parameters.push_back(pquery);
+		CommandParameter preference("reference", "InputTypes", "", "", "none", "none", "none","",false,true,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", "", "", "namecount", "none", "none","",false,false,true); parameters.push_back(pname);
+        CommandParameter pcount("count", "InputTypes", "", "", "namecount", "none", "none","",false,false,true); parameters.push_back(pcount);
+		CommandParameter pignorechimeras("ignorechimeras", "Boolean", "", "T", "", "", "","",false,false); parameters.push_back(pignorechimeras);
+		CommandParameter pthreshold("threshold", "Number", "", "1.0", "", "", "","",false,false); parameters.push_back(pthreshold);
+		CommandParameter paligned("aligned", "Boolean", "", "T", "", "", "","",false,false); parameters.push_back(paligned);
+		CommandParameter pprocessors("processors", "Number", "", "1", "", "", "","",false,false,true); parameters.push_back(pprocessors);
+		CommandParameter psave("save", "Boolean", "", "F", "", "", "","",false,false); parameters.push_back(psave);
+		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);		}
@@ -38,7 +42,9 @@ vector<string> SeqErrorCommand::setParameters(){
 		exit(1);
 	}
 }
+
 //**********************************************************************************************************************
+
 string SeqErrorCommand::getHelpString(){	
 	try {
 		string helpString = "";
@@ -47,7 +53,8 @@ string SeqErrorCommand::getHelpString(){
 		helpString += "The reference parameter...\n";
 		helpString += "The qfile parameter...\n";
 		helpString += "The report parameter...\n";
-		helpString += "The name parameter...\n";
+		helpString += "The name parameter allows you to provide a name file associated with the fasta file.\n";
+        helpString += "The count parameter allows you to provide a count file associated with the fasta file.\n";
 		helpString += "The ignorechimeras parameter...\n";
 		helpString += "The threshold parameter...\n";
 		helpString += "The processors parameter...\n";
@@ -62,37 +69,36 @@ string SeqErrorCommand::getHelpString(){
 		exit(1);
 	}
 }
+
 //**********************************************************************************************************************
-string SeqErrorCommand::getOutputFileNameTag(string type, string inputName=""){	
-	try {
-        string outputFileName = "";
-		map<string, vector<string> >::iterator it;
+
+string SeqErrorCommand::getOutputPattern(string type) {
+    try {
+        string pattern = "";
         
-        //is this a type this command creates
-        it = outputTypes.find(type);
-        if (it == outputTypes.end()) {  m->mothurOut("[ERROR]: this command doesn't create a " + type + " output file.\n"); }
-        else {
-            if (type == "errorsummary")            {   outputFileName =  "error.summary";   }
-            else if (type == "errorseq")            {   outputFileName =  "error.seq";   }
-            else if (type == "errorquality")            {   outputFileName =  "error.quality";   }
-            else if (type == "errorqualforward")            {   outputFileName =  "error.qual.forward";   }
-            else if (type == "errorqualreverse")            {   outputFileName =  "error.qual.reverse";   }
-            else if (type == "errorforward")            {   outputFileName =  "error.seq.forward";   }
-            else if (type == "errorreverse")            {   outputFileName =  "error.seq.reverse";   }
-            else if (type == "errorcount")            {   outputFileName =  "error.count";   }
-            else if (type == "errormatrix")            {   outputFileName =  "error.matrix";   }
-            else if (type == "errorchimera")            {   outputFileName =  "error.chimera";   }
-            else if (type == "errorref-query")            {   outputFileName =  "error.ref-query";   }
-            else { m->mothurOut("[ERROR]: No definition for type " + type + " output file tag.\n"); m->control_pressed = true;  }
-        }
-        return outputFileName;
-	}
-	catch(exception& e) {
-		m->errorOut(e, "SeqErrorCommand", "getOutputFileNameTag");
-		exit(1);
-	}
+        if (type == "errorsummary")            {   pattern = "[filename],error.summary";   }
+        else if (type == "errorseq")            {   pattern = "[filename],error.seq";   }
+        else if (type == "errorquality")            {   pattern = "[filename],error.quality";   }
+        else if (type == "errorqualforward")            {   pattern = "[filename],error.qual.forward";   }
+        else if (type == "errorqualreverse")            {   pattern = "[filename],error.qual.reverse";   }
+        else if (type == "errorforward")            {   pattern = "[filename],error.seq.forward";   }
+        else if (type == "errorreverse")            {   pattern = "[filename],error.seq.reverse";   }
+        else if (type == "errorcount")            {   pattern = "[filename],error.count";   }
+        else if (type == "errormatrix")            {   pattern = "[filename],error.matrix";   }
+        else if (type == "errorchimera")            {   pattern = "[filename],error.chimera";   }
+        else if (type == "errorref-query")            {   pattern = "[filename],error.ref-query";   }
+        else { m->mothurOut("[ERROR]: No definition for type " + type + " output pattern.\n"); m->control_pressed = true;  }
+        
+        return pattern;
+    }
+    catch(exception& e) {
+        m->errorOut(e, "SeqErrorCommand", "getOutputPattern");
+        exit(1);
+    }
 }
+
 //**********************************************************************************************************************
+
 SeqErrorCommand::SeqErrorCommand(){	
 	try {
 		abort = true; calledHelp = true; 
@@ -115,6 +121,7 @@ SeqErrorCommand::SeqErrorCommand(){
 		exit(1);
 	}
 }
+
 //***************************************************************************************************************
 
 SeqErrorCommand::SeqErrorCommand(string option)  {
@@ -185,6 +192,14 @@ SeqErrorCommand::SeqErrorCommand(string option)  {
 					//if the user has not given a path then, add inputdir. else leave path alone.
 					if (path == "") {	parameters["name"] = inputDir + it->second;		}
 				}
+                
+                it = parameters.find("count");
+				//user has given a names 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["count"] = inputDir + it->second;		}
+				}
 
 				it = parameters.find("qfile");
 				//user has given a quality score file
@@ -222,6 +237,12 @@ SeqErrorCommand::SeqErrorCommand(string option)  {
 			if(namesFileName == "not found"){	namesFileName = "";	}
 			else if (namesFileName == "not open") { namesFileName = ""; abort = true; }	
 			else { m->setNameFile(namesFileName); }
+            
+            //check for optional parameters
+			countfile = validParameter.validFile(parameters, "count", true);
+			if(countfile == "not found"){	countfile = "";	}
+			else if (countfile == "not open") { countfile = ""; abort = true; }
+			else { m->setCountTableFile(countfile); }
 			
 			qualFileName = validParameter.validFile(parameters, "qfile", true);
 			if(qualFileName == "not found"){	qualFileName = "";	}
@@ -232,22 +253,19 @@ SeqErrorCommand::SeqErrorCommand(string option)  {
 			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.");
-				m->mothurOutEndLine();
-				abort = true; 
-			}
-			
 			outputDir = validParameter.validFile(parameters, "outputdir", false);
 			if (outputDir == "not found"){	
 				outputDir = "";	
 				outputDir += m->hasPath(queryFileName); //if user entered a file with a path then preserve it	
 			}
 			
+            if ((countfile != "") && (namesFileName != "")) { m->mothurOut("You must enter ONLY ONE of the following: count or name."); m->mothurOutEndLine(); abort = true; }
+            
 			//check for optional parameter and set defaults
 			// ...at some point should added some additional type checking...
 			temp = validParameter.validFile(parameters, "threshold", false);	if (temp == "not found") { temp = "1.00"; }
 			m->mothurConvert(temp, threshold);  
+            
 			
 			temp = validParameter.validFile(parameters, "save", false);			if (temp == "not found"){	temp = "f";				}
 			save = m->isTrue(temp); 
@@ -274,6 +292,9 @@ SeqErrorCommand::SeqErrorCommand(string option)  {
 			temp = validParameter.validFile(parameters, "ignorechimeras", false);	if (temp == "not found") { temp = "T"; }
 			ignoreChimeras = m->isTrue(temp);
 			
+            temp = validParameter.validFile(parameters, "aligned", false);			if (temp == "not found"){	temp = "t";				}
+			aligned = m->isTrue(temp); 
+
 			temp = validParameter.validFile(parameters, "processors", false);	if (temp == "not found"){	temp = m->getProcessors();	}
 			m->setProcessors(temp);
 			m->mothurConvert(temp, processors); 
@@ -285,6 +306,22 @@ SeqErrorCommand::SeqErrorCommand(string option)  {
 				vector<string> files; files.push_back(queryFileName); 
 				parser.getNameFile(files);
 			}
+
+            if(aligned == true){
+                if((reportFileName != "" && qualFileName == "") || (reportFileName == "" && qualFileName != "")){
+                    m->mothurOut("if you use either a qual file or a report file, you have to have both.");
+                    m->mothurOutEndLine();
+                    abort = true; 
+                }
+			}
+            else{
+                if(reportFileName != ""){
+                    m->mothurOut("we are ignoring the report file if your sequences are not aligned.  we will check that the sequences in your fasta and and qual file are the same length.");
+                    m->mothurOutEndLine();
+                }
+            }
+            
+
 		}
 	}
 	catch(exception& e) {
@@ -292,6 +329,7 @@ SeqErrorCommand::SeqErrorCommand(string option)  {
 		exit(1);
 	}
 }
+
 //***************************************************************************************************************
 
 int SeqErrorCommand::execute(){
@@ -299,23 +337,30 @@ int SeqErrorCommand::execute(){
 		if (abort == true) { if (calledHelp) { return 0; }  return 2;	}
 		
 		int start = time(NULL);
-		maxLength = 2000;
+		maxLength = 5000;
 		totalBases = 0;
 		totalMatches = 0;
 		
         string fileNameRoot = outputDir + m->getRootName(m->getSimpleName(queryFileName));
-		string errorSummaryFileName = fileNameRoot + getOutputFileNameTag("errorsummary");
+        map<string, string> variables; 
+		variables["[filename]"] = fileNameRoot;
+		string errorSummaryFileName = getOutputFileName("errorsummary",variables);
 		outputNames.push_back(errorSummaryFileName); outputTypes["errorsummary"].push_back(errorSummaryFileName);
 			
-		string errorSeqFileName = fileNameRoot + getOutputFileNameTag("errorseq");
+		string errorSeqFileName = getOutputFileName("errorseq",variables);
 		outputNames.push_back(errorSeqFileName); outputTypes["errorseq"].push_back(errorSeqFileName);
 		
-		string errorChimeraFileName = fileNameRoot + getOutputFileNameTag("errorchimera");
+		string errorChimeraFileName = getOutputFileName("errorchimera",variables);
 		outputNames.push_back(errorChimeraFileName); outputTypes["errorchimera"].push_back(errorChimeraFileName);
 		
 		getReferences();	//read in reference sequences - make sure there's no ambiguous bases
 
-		if(namesFileName != ""){	weights = getWeights();	}
+		if(namesFileName != "")     {	weights = getWeights();         }
+        else if (countfile != "")   {
+            CountTable ct;
+            ct.readTable(countfile, false);
+            weights = ct.getNameMap();
+        }
 		
 		vector<unsigned long long> fastaFilePos;
 		vector<unsigned long long> qFilePos;
@@ -328,22 +373,23 @@ int SeqErrorCommand::execute(){
 		for (int i = 0; i < (fastaFilePos.size()-1); i++) {
 			lines.push_back(linePair(fastaFilePos[i], fastaFilePos[(i+1)]));
 			if (qualFileName != "") {  qLines.push_back(linePair(qFilePos[i], qFilePos[(i+1)]));  }
-			if (reportFileName != "") {  rLines.push_back(linePair(reportFilePos[i], reportFilePos[(i+1)]));  }
+			if (reportFileName != "" && aligned == true) {  rLines.push_back(linePair(reportFilePos[i], reportFilePos[(i+1)]));  }
 		}	
 		if(qualFileName == "")	{	qLines = lines;	rLines = lines; } //fills with duds
-		
+        if(aligned == false){   rLines = lines; }
 		int numSeqs = 0;
 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
 		if(processors == 1){
 			numSeqs = driver(queryFileName, qualFileName, reportFileName, errorSummaryFileName, errorSeqFileName, errorChimeraFileName, lines[0], qLines[0], rLines[0]);
+            
 		}else{
 			numSeqs = createProcesses(queryFileName, qualFileName, reportFileName, errorSummaryFileName, errorSeqFileName, errorChimeraFileName);
 		}	
 #else
 		numSeqs = driver(queryFileName, qualFileName, reportFileName, errorSummaryFileName, errorSeqFileName, errorChimeraFileName, lines[0], qLines[0], rLines[0]);
 #endif
-		
-		if(qualFileName != "" && reportFileName != ""){		
+
+		if(qualFileName != ""){		
 			printErrorQuality(qScoreErrorMap);
 			printQualityFR(qualForwardMap, qualReverseMap);
 		}
@@ -352,7 +398,7 @@ int SeqErrorCommand::execute(){
 		
 		if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; }
 
-		string errorCountFileName = fileNameRoot + getOutputFileNameTag("errorcount");
+		string errorCountFileName = getOutputFileName("errorcount",variables);
 		ofstream errorCountFile;
 		m->openOutputFile(errorCountFileName, errorCountFile);
 		outputNames.push_back(errorCountFileName);  outputTypes["errorcount"].push_back(errorCountFileName);
@@ -365,23 +411,32 @@ int SeqErrorCommand::execute(){
 		}
 		errorCountFile.close();
 		
-		if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; }
+//		if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; }
 
 		printSubMatrix();
-				
-		string megAlignmentFileName = fileNameRoot + getOutputFileNameTag("errorref-query");
+        
+		string megAlignmentFileName = getOutputFileName("errorref-query",variables);
 		ofstream megAlignmentFile;
 		m->openOutputFile(megAlignmentFileName, megAlignmentFile);
 		outputNames.push_back(megAlignmentFileName);  outputTypes["errorref-query"].push_back(megAlignmentFileName);
-		
+        
 		for(int i=0;i<numRefs;i++){
 			megAlignmentFile << referenceSeqs[i].getInlineSeq() << endl;
 			megAlignmentFile << megaAlignVector[i] << endl;
 		}
+        megAlignmentFile.close();
+        
 		
 		m->mothurOut("It took " + toString(time(NULL) - start) + " secs to check " + toString(numSeqs) + " sequences.");
 		m->mothurOutEndLine();
 		
+        //set fasta file as new current fastafile
+		string current = "";
+		itTypes = outputTypes.find("errorseq");
+		if (itTypes != outputTypes.end()) {
+			if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setFastaFile(current); }
+		}
+        
 		m->mothurOutEndLine();
 		m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 		for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
@@ -394,7 +449,9 @@ int SeqErrorCommand::execute(){
 		exit(1);
 	}
 }
+
 //**********************************************************************************************************************
+
 int SeqErrorCommand::createProcesses(string filename, string qFileName, string rFileName, string summaryFileName, string errorOutputFileName, string chimeraOutputFileName) {	
 	try {
 		int process = 1;
@@ -609,7 +666,7 @@ int SeqErrorCommand::createProcesses(string filename, string qFileName, string r
 			int misMatchSize;
 			in >> misMatchSize; m->gobble(in);
 			if (misMatchSize > misMatchCounts.size()) {	misMatchCounts.resize(misMatchSize, 0);	}
-			for (int j = 0; j < misMatchCounts.size(); j++) {
+			for (int j = 0; j < misMatchSize; j++) {
 				in >> tempNum; misMatchCounts[j] += tempNum;
 			}
 			m->gobble(in);
@@ -632,7 +689,9 @@ int SeqErrorCommand::createProcesses(string filename, string qFileName, string r
 		exit(1);
 	}
 }
+
 //**********************************************************************************************************************
+
 int SeqErrorCommand::driver(string filename, string qFileName, string rFileName, string summaryFileName, string errorOutputFileName, string chimeraOutputFileName, linePair line, linePair qline, linePair rline) {	
 	
 	try {
@@ -644,10 +703,10 @@ int SeqErrorCommand::driver(string filename, string qFileName, string rFileName,
 		int numSeqs = 0;
 		
 		map<string, int>::iterator it;
-		qScoreErrorMap['m'].assign(41, 0);
-		qScoreErrorMap['s'].assign(41, 0);
-		qScoreErrorMap['i'].assign(41, 0);
-		qScoreErrorMap['a'].assign(41, 0);
+		qScoreErrorMap['m'].assign(101, 0);
+		qScoreErrorMap['s'].assign(101, 0);
+		qScoreErrorMap['i'].assign(101, 0);
+		qScoreErrorMap['a'].assign(101, 0);
 		
 		errorForward['m'].assign(maxLength,0);
 		errorForward['s'].assign(maxLength,0);
@@ -664,11 +723,12 @@ int SeqErrorCommand::driver(string filename, string qFileName, string rFileName,
 		//open inputfiles and go to beginning place for this processor
 		ifstream queryFile;
 		m->openInputFile(filename, queryFile);
+        
 		queryFile.seekg(line.start);
 		
 		ifstream reportFile;
 		ifstream qualFile;
-		if(qFileName != "" && rFileName != ""){
+		if((qFileName != "" && rFileName != "" && aligned)){
 			m->openInputFile(qFileName, qualFile);
 			qualFile.seekg(qline.start);  
 			
@@ -682,15 +742,32 @@ int SeqErrorCommand::driver(string filename, string qFileName, string rFileName,
 			qualForwardMap.resize(maxLength);
 			qualReverseMap.resize(maxLength);
 			for(int i=0;i<maxLength;i++){
-				qualForwardMap[i].assign(41,0);
-				qualReverseMap[i].assign(41,0);
+				qualForwardMap[i].assign(101,0);
+				qualReverseMap[i].assign(101,0);
 			}	
 		}
-		
+		else if(qFileName != "" && !aligned){
+
+            m->openInputFile(qFileName, qualFile);
+			qualFile.seekg(qline.start);  
+			
+			qualForwardMap.resize(maxLength);
+			qualReverseMap.resize(maxLength);
+			for(int i=0;i<maxLength;i++){
+				qualForwardMap[i].assign(101,0);
+				qualReverseMap[i].assign(101,0);
+			}	
+        }
+        
 		ofstream outChimeraReport;
 		m->openOutputFile(chimeraOutputFileName, outChimeraReport);
-		RefChimeraTest chimeraTest(referenceSeqs);
-		if (line.start == 0) { chimeraTest.printHeader(outChimeraReport); }
+		
+        
+        RefChimeraTest chimeraTest;
+        
+        chimeraTest = RefChimeraTest(referenceSeqs, aligned);
+        if (line.start == 0) { chimeraTest.printHeader(outChimeraReport); }        
+        
 		
 		ofstream errorSummaryFile;
 		m->openOutputFile(summaryFileName, errorSummaryFile);
@@ -699,48 +776,58 @@ int SeqErrorCommand::driver(string filename, string qFileName, string rFileName,
 		ofstream errorSeqFile;
 		m->openOutputFile(errorOutputFileName, errorSeqFile);
 		
-		megaAlignVector.resize(numRefs, "");
+		megaAlignVector.assign(numRefs, "");
 		
 		int index = 0;
 		bool ignoreSeq = 0;
 		
 		bool moreSeqs = 1;
 		while (moreSeqs) {
-			
-			if (m->control_pressed) { queryFile.close(); if(qFileName != "" && rFileName != ""){  reportFile.close(); qualFile.close(); } outChimeraReport.close(); errorSummaryFile.close();errorSeqFile.close(); return 0; }
-			
+						
 			Sequence query(queryFile);
-			
-			int numParentSeqs = chimeraTest.analyzeQuery(query.getName(), query.getAligned(), outChimeraReport);
-			int closestRefIndex = chimeraTest.getClosestRefIndex();
-			
+			Sequence reference;
+            int numParentSeqs = -1;
+            int closestRefIndex = -1;
+                        
+            string querySeq = query.getAligned();
+            if (!aligned) {  querySeq = query.getUnaligned();  }
+            
+            numParentSeqs = chimeraTest.analyzeQuery(query.getName(), querySeq, outChimeraReport);
+            
+            closestRefIndex = chimeraTest.getClosestRefIndex();
+            
+            reference = referenceSeqs[closestRefIndex];
+            
+            reference.setAligned(chimeraTest.getClosestRefAlignment());
+            query.setAligned(chimeraTest.getQueryAlignment());
+            
 			if(numParentSeqs > 1 && ignoreChimeras == 1)	{	ignoreSeq = 1;	}
-			else											{	ignoreSeq = 0;	}
-			
+            else											{	ignoreSeq = 0;	}
+
 			Compare minCompare;
-            getErrors(query, referenceSeqs[closestRefIndex], minCompare);
+            
+            getErrors(query, reference, minCompare);
 			
-			if(namesFileName != ""){
+			if((namesFileName != "") || (countfile != "")){
 				it = weights.find(query.getName());
 				minCompare.weight = it->second;
 			}
 			else{	minCompare.weight = 1;	}
 			
+            
 			printErrorData(minCompare, numParentSeqs, errorSummaryFile, errorSeqFile);
 			
 			if(!ignoreSeq){
-				
 				for(int i=0;i<minCompare.sequence.length();i++){
 					char letter = minCompare.sequence[i];
-					
 					if(letter != 'r'){
 						errorForward[letter][i] += minCompare.weight;
-						errorReverse[letter][minCompare.total-i-1] += minCompare.weight;				
+						errorReverse[letter][minCompare.total-i-1] += minCompare.weight;	
 					}
-				}
+				}                
 			}
 			
-			if(qualFileName != "" && reportFileName != ""){
+			if(aligned && qualFileName != "" && reportFileName != ""){
 				report = ReportFile(reportFile);
 				
 				//				int origLength = report.getQueryLength();
@@ -750,17 +837,29 @@ int SeqErrorCommand::driver(string filename, string qFileName, string rFileName,
 				quality = QualityScores(qualFile);
 				
 				if(!ignoreSeq){
-//					cout << query.getName() << '\t';
-					
 					quality.updateQScoreErrorMap(qScoreErrorMap, minCompare.sequence, startBase, endBase, minCompare.weight);
 					quality.updateForwardMap(qualForwardMap, startBase, endBase, minCompare.weight);
 					quality.updateReverseMap(qualReverseMap, startBase, endBase, minCompare.weight);
+				}
+			}
+            else if(aligned == false && qualFileName != ""){
 
-//					cout << endl;
+                quality = QualityScores(qualFile);
+                int qualityLength = quality.getLength();
+                
+                if(qualityLength != query.getNumBases()){   cout << "warning - quality and fasta sequence files do not match at " << query.getName() << '\t' << qualityLength <<'\t' << query.getNumBases() << endl;   }
+                
+                int startBase = 1;
+                int endBase = qualityLength;
+
+                if(!ignoreSeq){
+					quality.updateQScoreErrorMap(qScoreErrorMap, minCompare.sequence, startBase, endBase, minCompare.weight);
+					quality.updateForwardMap(qualForwardMap, startBase, endBase, minCompare.weight);
+					quality.updateReverseMap(qualReverseMap, startBase, endBase, minCompare.weight);
 				}
-			}			
-			
-			if(minCompare.errorRate < threshold && !ignoreSeq){
+            }
+            
+			if(minCompare.errorRate <= threshold && !ignoreSeq){                
 				totalBases += (minCompare.total * minCompare.weight);
 				totalMatches += minCompare.matches * minCompare.weight;
 				if(minCompare.mismatches > maxMismatch){
@@ -782,15 +881,18 @@ int SeqErrorCommand::driver(string filename, string qFileName, string rFileName,
 				if (queryFile.eof()) { break; }
 			#endif
 			
-			if(index % 100 == 0){	m->mothurOut(toString(index) + '\n');	}
+			if(index % 100 == 0){	m->mothurOutJustToScreen(toString(index)+"\n");	 }
 		}
 		queryFile.close();
-		if(qFileName != "" && rFileName != ""){  reportFile.close(); qualFile.close(); }
-		errorSummaryFile.close();	
+		outChimeraReport.close();
+        errorSummaryFile.close();	
 		errorSeqFile.close();
-		
+        
+        if(qFileName != "" && rFileName != "")      {   reportFile.close(); qualFile.close();   }
+		else if(qFileName != "" && aligned == false){   qualFile.close();                       }
+        
 		//report progress
-		if(index % 100 != 0){	m->mothurOut(toString(index) + '\n');	}
+		m->mothurOutJustToScreen(toString(index)+"\n");	
 		
 		return index;
 	}
@@ -799,6 +901,7 @@ int SeqErrorCommand::driver(string filename, string qFileName, string rFileName,
 		exit(1);
 	}
 }
+
 //***************************************************************************************************************
 
 void SeqErrorCommand::getReferences(){
@@ -821,12 +924,16 @@ void SeqErrorCommand::getReferences(){
 				//
 				//			int endPos = rdb->referenceSeqs[i].getEndPos();
 				//			if(endPos < minEndPos)		{	minEndPos = endPos;		}				
+				if (rdb->referenceSeqs[i].getNumBases() == 0) {
+                    m->mothurOut("[WARNING]: " + rdb->referenceSeqs[i].getName() + " is blank, ignoring.");m->mothurOutEndLine(); 
+                }else {
+                    referenceSeqs.push_back(rdb->referenceSeqs[i]);
+                }
 				
-				referenceSeqs.push_back(rdb->referenceSeqs[i]);
 			}
 			referenceFileName = rdb->getSavedReference();
 			
-			m->mothurOut("It took " + toString(time(NULL) - start) + " to load " + toString(rdb->referenceSeqs.size()) + " sequences.");m->mothurOutEndLine();  
+			m->mothurOut("It took " + toString(time(NULL) - start) + " to load " + toString(referenceSeqs.size()) + " sequences.");m->mothurOutEndLine();  
 		
 		}else {
 			int start = time(NULL);
@@ -844,9 +951,12 @@ void SeqErrorCommand::getReferences(){
 	//
 	//			int endPos = currentSeq.getEndPos();
 	//			if(endPos < minEndPos)		{	minEndPos = endPos;		}
-				referenceSeqs.push_back(currentSeq);
-				
-				if (rdb->save) { rdb->referenceSeqs.push_back(currentSeq); }
+                if (currentSeq.getNumBases() == 0) {
+                    m->mothurOut("[WARNING]: " + currentSeq.getName() + " is blank, ignoring.");m->mothurOutEndLine(); 
+                }else {
+                    referenceSeqs.push_back(currentSeq);
+                    if (rdb->save) { rdb->referenceSeqs.push_back(currentSeq); }
+                }
 					
 				m->gobble(referenceFile);
 			}
@@ -860,7 +970,7 @@ void SeqErrorCommand::getReferences(){
 		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");
@@ -887,9 +997,10 @@ int SeqErrorCommand::getErrors(Sequence query, Sequence reference, Compare& erro
 
 		int started = 0;
 		//Compare errors;
-
+        
+        errors.sequence = "";
 		for(int i=0;i<alignLength;i++){
-//			cout << r[i] << '\t' << q[i] << '\t';
+
 			if(q[i] != '.' && r[i] != '.' && (q[i] != '-' || r[i] != '-')){			//	no missing data and no double gaps
 				if(r[i] != 'N'){
 					started = 1;
@@ -946,8 +1057,7 @@ int SeqErrorCommand::getErrors(Sequence query, Sequence reference, Compare& erro
 						errors.sequence += 'r';
 					}
 				}
-			}
-				
+			}            
 			else if(q[i] == '.' && r[i] != '.'){		//	reference extends beyond query
 				if(started == 1){	break;	}
 			}
@@ -957,15 +1067,15 @@ int SeqErrorCommand::getErrors(Sequence query, Sequence reference, Compare& erro
 			else if(q[i] == '.' && r[i] == '.'){		//	both are missing data
 				if(started == 1){	break;	}			
 			}
-//			cout << errors.sequence[errors.sequence.length()-1] << endl;
 		}
-//		cout << errors.sequence << endl;
+
 		errors.mismatches = errors.total-errors.matches;
-		errors.errorRate = (double)(errors.total-errors.matches) / (double)errors.total;
+        if(errors.total != 0){  errors.errorRate = (double)(errors.total-errors.matches) / (double)errors.total;    }
+        else{   errors.errorRate = 0;   }
+
 		errors.queryName = query.getName();
 		errors.refName = reference.getName();
-		
-		//return errors;
+
         return 0;
 	}
 	catch(exception& e) {
@@ -989,10 +1099,12 @@ map<string, int> SeqErrorCommand::getWeights(){
 		nameCountMap[seqName] = m->getNumNames(redundantSeqs); 
 		m->gobble(nameFile);
 	}
+    
+    nameFile.close();
+    
 	return nameCountMap;
 }
 
-
 //***************************************************************************************************************
 
 void SeqErrorCommand::printErrorHeader(ofstream& errorSummaryFile){
@@ -1082,7 +1194,9 @@ void SeqErrorCommand::printErrorData(Compare error, int numParentSeqs, ofstream&
 void SeqErrorCommand::printSubMatrix(){
 	try {
         string fileNameRoot = outputDir + m->getRootName(m->getSimpleName(queryFileName));
-		string subMatrixFileName = fileNameRoot + getOutputFileNameTag("errormatrix");
+        map<string, string> variables; 
+		variables["[filename]"] = fileNameRoot;
+		string subMatrixFileName = getOutputFileName("errormatrix",variables);
 		ofstream subMatrixFile;
 		m->openOutputFile(subMatrixFileName, subMatrixFile);
 		outputNames.push_back(subMatrixFileName);  outputTypes["errormatrix"].push_back(subMatrixFileName);
@@ -1124,12 +1238,15 @@ void SeqErrorCommand::printSubMatrix(){
 		exit(1);
 	}
 }
+
 //***************************************************************************************************************
 
 void SeqErrorCommand::printErrorFRFile(map<char, vector<int> > errorForward, map<char, vector<int> > errorReverse){
 	try{
         string fileNameRoot = outputDir + m->getRootName(m->getSimpleName(queryFileName));
-		string errorForwardFileName = fileNameRoot + getOutputFileNameTag("errorforward");
+        map<string, string> variables; 
+		variables["[filename]"] = fileNameRoot;
+		string errorForwardFileName = getOutputFileName("errorforward",variables);
 		ofstream errorForwardFile;
 		m->openOutputFile(errorForwardFileName, errorForwardFile);
 		outputNames.push_back(errorForwardFileName);  outputTypes["errorforward"].push_back(errorForwardFileName);
@@ -1147,7 +1264,7 @@ void SeqErrorCommand::printErrorFRFile(map<char, vector<int> > errorForward, map
 		}
 		errorForwardFile.close();
 
-		string errorReverseFileName = fileNameRoot + getOutputFileNameTag("errorreverse");
+		string errorReverseFileName = getOutputFileName("errorreverse",variables);
 		ofstream errorReverseFile;
 		m->openOutputFile(errorReverseFileName, errorReverseFile);
 		outputNames.push_back(errorReverseFileName);  outputTypes["errorreverse"].push_back(errorReverseFileName);
@@ -1176,24 +1293,25 @@ void SeqErrorCommand::printErrorFRFile(map<char, vector<int> > errorForward, map
 void SeqErrorCommand::printErrorQuality(map<char, vector<int> > qScoreErrorMap){
 	try{
         string fileNameRoot = outputDir + m->getRootName(m->getSimpleName(queryFileName));
-		string errorQualityFileName = fileNameRoot + getOutputFileNameTag("errorquality");
+        map<string, string> variables; 
+		variables["[filename]"] = fileNameRoot;
+		string errorQualityFileName = getOutputFileName("errorquality",variables);
 		ofstream errorQualityFile;
 		m->openOutputFile(errorQualityFileName, errorQualityFile);
 		outputNames.push_back(errorQualityFileName);  outputTypes["errorquality"].push_back(errorQualityFileName);
 
 		errorQualityFile << "qscore\tmatches\tsubstitutions\tinsertions\tambiguous" << endl;
-		for(int i=0;i<41;i++){
+		for(int i=0;i<101;i++){
 			errorQualityFile << i << '\t' << qScoreErrorMap['m'][i] << '\t' << qScoreErrorMap['s'][i] << '\t' << qScoreErrorMap['i'][i] << '\t'<< qScoreErrorMap['a'][i] << endl;
 		}
 		errorQualityFile.close();
 	}
 	catch(exception& e) {
-		m->errorOut(e, "SeqErrorCommand", "printErrorFRFile");
+		m->errorOut(e, "SeqErrorCommand", "printErrorQuality");
 		exit(1);
 	}
 }
 
-
 //***************************************************************************************************************
 
 void SeqErrorCommand::printQualityFR(vector<vector<int> > qualForwardMap, vector<vector<int> > qualReverseMap){
@@ -1210,7 +1328,9 @@ void SeqErrorCommand::printQualityFR(vector<vector<int> > qualForwardMap, vector
 			}
 		}
         string fileNameRoot = outputDir + m->getRootName(m->getSimpleName(queryFileName));
-		string qualityForwardFileName = fileNameRoot + getOutputFileNameTag("errorqualforward");
+        map<string, string> variables; 
+		variables["[filename]"] = fileNameRoot;
+		string qualityForwardFileName = getOutputFileName("errorqualforward",variables);
 		ofstream qualityForwardFile;
 		m->openOutputFile(qualityForwardFileName, qualityForwardFile);
 		outputNames.push_back(qualityForwardFileName);  outputTypes["errorqualforward"].push_back(qualityForwardFileName);
@@ -1228,7 +1348,7 @@ void SeqErrorCommand::printQualityFR(vector<vector<int> > qualForwardMap, vector
 		qualityForwardFile.close();
 
 		
-		string qualityReverseFileName = fileNameRoot + getOutputFileNameTag("errorqualreverse");
+		string qualityReverseFileName = getOutputFileName("errorqualreverse",variables);
 		ofstream qualityReverseFile;
 		m->openOutputFile(qualityReverseFileName, qualityReverseFile);
 		outputNames.push_back(qualityReverseFileName);  outputTypes["errorqualreverse"].push_back(qualityReverseFileName);
@@ -1245,11 +1365,12 @@ void SeqErrorCommand::printQualityFR(vector<vector<int> > qualForwardMap, vector
 		qualityReverseFile.close();
 	}
 	catch(exception& e) {
-		m->errorOut(e, "SeqErrorCommand", "printErrorFRFile");
+		m->errorOut(e, "SeqErrorCommand", "printQualityFR");
 		exit(1);
 	}
 	
 }
+
 /**************************************************************************************************/
 
 int SeqErrorCommand::setLines(string filename, string qfilename, string rfilename, vector<unsigned long long>& fastaFilePos, vector<unsigned long long>& qfileFilePos, vector<unsigned long long>& rfileFilePos) {
@@ -1290,6 +1411,8 @@ int SeqErrorCommand::setLines(string filename, string qfilename, string rfilenam
 					
 					string sname = "";  nameStream >> sname;
 					sname = sname.substr(1);
+                    
+                    m->checkName(sname);
 					
 					map<string, int>::iterator it = firstSeqNames.find(sname);
 					
@@ -1328,65 +1451,67 @@ int SeqErrorCommand::setLines(string filename, string qfilename, string rfilenam
 		
 		qfileFilePos.push_back(size);
 		
-		//seach for filePos of each first name in the rfile and save in rfileFilePos
-		string junk;
-		ifstream inR;
-		m->openInputFile(rfilename, inR);
-		
-		//read column headers
-		for (int i = 0; i < 16; i++) {  
-			if (!inR.eof())	{	inR >> junk;	}
-			else			{	break;			}
-		}
-		
-		while(!inR.eof()){
-			
-			if (m->control_pressed) { inR.close();  return processors; }
-			
-			input = m->getline(inR);	
-			
-			if (input.length() != 0) {
-				
-				istringstream nameStream(input);
-				string sname = "";  nameStream >> sname;
-				
-				map<string, int>::iterator it = firstSeqNamesReport.find(sname);
-			
-				if(it != firstSeqNamesReport.end()) { //this is the start of a new chunk
-					unsigned long long pos = inR.tellg(); 
-					rfileFilePos.push_back(pos - input.length() - 1);	
-					firstSeqNamesReport.erase(it);
-				}
-			}
-			
-			if (firstSeqNamesReport.size() == 0) { break; }
-			m->gobble(inR);
-		}
-		inR.close();
-		
-		if (firstSeqNamesReport.size() != 0) { 
-			for (map<string, int>::iterator it = firstSeqNamesReport.begin(); it != firstSeqNamesReport.end(); it++) {
-				m->mothurOut(it->first + " is in your fasta file and not in your report file, aborting."); m->mothurOutEndLine();
-			}
-			m->control_pressed = true;
-			return processors;
-		}
-		
-		//get last file position of qfile
-		FILE * rFile;
-		unsigned long long sizeR;
-		
-		//get num bytes in file
-		rFile = fopen (rfilename.c_str(),"rb");
-		if (rFile==NULL) perror ("Error opening file");
-		else{
-			fseek (rFile, 0, SEEK_END);
-			sizeR=ftell (rFile);
-			fclose (rFile);
+        if(aligned){
+            //seach for filePos of each first name in the rfile and save in rfileFilePos
+            string junk;
+            ifstream inR;
+
+            m->openInputFile(rfilename, inR);
+
+            //read column headers
+            for (int i = 0; i < 16; i++) {  
+                if (!inR.eof())	{	inR >> junk;	}
+                else			{	break;			}
+            }
+            
+            while(!inR.eof()){
+                
+                input = m->getline(inR);	
+                
+                if (input.length() != 0) {
+                    
+                    istringstream nameStream(input);
+                    string sname = "";  nameStream >> sname;
+                    
+                    m->checkName(sname);
+                    
+                    map<string, int>::iterator it = firstSeqNamesReport.find(sname);
+                
+                    if(it != firstSeqNamesReport.end()) { //this is the start of a new chunk
+                        unsigned long long pos = inR.tellg(); 
+                        rfileFilePos.push_back(pos - input.length() - 1);	
+                        firstSeqNamesReport.erase(it);
+                    }
+                }
+                
+                if (firstSeqNamesReport.size() == 0) { break; }
+                m->gobble(inR);
+            }
+            inR.close();
+            
+            if (firstSeqNamesReport.size() != 0) { 
+                for (map<string, int>::iterator it = firstSeqNamesReport.begin(); it != firstSeqNamesReport.end(); it++) {
+                    m->mothurOut(it->first + " is in your fasta file and not in your report file, aborting."); m->mothurOutEndLine();
+                }
+                m->control_pressed = true;
+                return processors;
+            }
+            
+            //get last file position of qfile
+            FILE * rFile;
+            unsigned long long sizeR;
+            
+            //get num bytes in file
+            rFile = fopen (rfilename.c_str(),"rb");
+            if (rFile==NULL) perror ("Error opening file");
+            else{
+                fseek (rFile, 0, SEEK_END);
+                sizeR=ftell (rFile);
+                fclose (rFile);
+            }
+            
+            rfileFilePos.push_back(sizeR);
 		}
-		
-		rfileFilePos.push_back(sizeR);
-		
 		return processors;
 		
 #else
@@ -1428,4 +1553,5 @@ int SeqErrorCommand::setLines(string filename, string qfilename, string rfilenam
 		exit(1);
 	}
 }
+
 //***************************************************************************************************************