From cc1838f288a7712d44929fca5175ce89f3a58f12 Mon Sep 17 00:00:00 2001
From: westcott <westcott>
Date: Mon, 21 Sep 2009 13:19:02 +0000
Subject: [PATCH] align.check command added

---
 commandfactory.cpp            |   3 +
 secondarystructurecommand.cpp | 104 +++++++++++++++++++++++++++++++---
 secondarystructurecommand.h   |   1 +
 3 files changed, 99 insertions(+), 9 deletions(-)

diff --git a/commandfactory.cpp b/commandfactory.cpp
index fce55bb..dd5b363 100644
--- a/commandfactory.cpp
+++ b/commandfactory.cpp
@@ -54,6 +54,7 @@
 #include "getseqscommand.h"
 #include "removeseqscommand.h"
 #include "systemcommand.h"
+#include "secondarystructurecommand.h"
 
 /***********************************************************/
 
@@ -106,6 +107,7 @@ CommandFactory::CommandFactory(){
 	commands["get.seqs"]			= "get.seqs";
 	commands["remove.seqs"]			= "get.seqs";
 	commands["system"]				= "system";
+	commands["align.check"]			= "align.check";
 	commands["quit"]				= "quit"; 
 
 }
@@ -168,6 +170,7 @@ Command* CommandFactory::getCommand(string commandName, string optionString){
 		else if(commandName == "remove.seqs")			{	command = new RemoveSeqsCommand(optionString);			}
 		else if(commandName == "merge.files")			{	command = new MergeFileCommand(optionString);			}
 		else if(commandName == "system")				{	command = new SystemCommand(optionString);				}
+		else if(commandName == "align.check")			{	command = new AlignCheckCommand(optionString);			}
 		else											{	command = new NoCommand(optionString);					}
 
 		return command;
diff --git a/secondarystructurecommand.cpp b/secondarystructurecommand.cpp
index 48ea154..7156082 100644
--- a/secondarystructurecommand.cpp
+++ b/secondarystructurecommand.cpp
@@ -55,12 +55,12 @@ AlignCheckCommand::AlignCheckCommand(string option){
 
 void AlignCheckCommand::help(){
 	try {
-		//mothurOut("The remove.seqs command reads an .accnos file and one of the following file types: fasta, name, group or alignreport file.\n");
-		//mothurOut("It outputs a file containing the sequences NOT in the .accnos file.\n");
-		//mothurOut("The remove.seqs command parameters are accnos, fasta, name, group and alignreport.  You must provide accnos and one of the other parameters.\n");
-		//mothurOut("The remove.seqs command should be in the following format: remove.seqs(accnos=yourAccnos, fasta=yourFasta).\n");
-		//mothurOut("Example remove.seqs(accnos=amazon.accnos, fasta=amazon.fasta).\n");
-		//mothurOut("Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFasta).\n\n");
+		mothurOut("The align.check command reads a fasta file and map file.\n");
+		mothurOut("It outputs a file containing the secondary structure matches in the .align.check file.\n");
+		mothurOut("The align.check command parameters are fasta and map, both are required.\n");
+		mothurOut("The align.check command should be in the following format: align.check(fasta=yourFasta, map=yourMap).\n");
+		mothurOut("Example align.check(map=silva.ss.map, fasta=amazon.fasta).\n");
+		mothurOut("Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFasta).\n\n");
 	}
 	catch(exception& e) {
 		errorOut(e, "AlignCheckCommand", "help");
@@ -78,7 +78,28 @@ int AlignCheckCommand::execute(){
 		//get secondary structure info.
 		readMap();
 		
-	
+		ifstream in;
+		openInputFile(fastafile, in);
+		
+		ofstream out;
+		string outfile = getRootName(fastafile) + "align.check";
+		openOutputFile(outfile, out);
+		
+		out << "name" << '\t' << "pound" << '\t' << "dash" << '\t' << "plus" << '\t' << "equal" << '\t';
+		out << "loop" << '\t' << "tilde" << '\t' << "total" << endl;
+
+		
+		while(!in.eof()){
+			
+			Sequence seq(in);
+			statData data = getStats(seq.getAligned());
+			
+			out << seq.getName() << '\t' << data.pound << '\t' << data.dash << '\t' << data.plus << '\t' << data.equal << '\t';
+			out << data.loop << '\t' << data.tilde << '\t' << data.total << endl;
+		}
+
+		in.close();
+		out.close();
 		
 		return 0;		
 	}
@@ -104,7 +125,7 @@ void AlignCheckCommand::readMap(){
 			gobble(in);
 		}
 		in.close();
-		
+
 		seqLength = structMap.size();
 		
 		
@@ -120,9 +141,74 @@ void AlignCheckCommand::readMap(){
 		
 	}
 	catch(exception& e) {
-		errorOut(e, "AlignCheckCommand", "readFasta");
+		errorOut(e, "AlignCheckCommand", "readMap");
 		exit(1);
 	}
 }
+/**************************************************************************************************/
+
+statData AlignCheckCommand::getStats(string sequence){
+	try {
+	
+		statData data;
+		sequence = "*" + sequence; // need to pad the sequence so we can index it by 1
+		
+		int seqLength = sequence.length();
+		for(int i=1;i<seqLength;i++){
+			if(structMap[i] != 0){
+				if(sequence[i] == 'A'){
+					if(sequence[structMap[i]] == 'T')		{	data.tilde++;	}
+					else if(sequence[structMap[i]] == 'A')	{	data.pound++;	}
+					else if(sequence[structMap[i]] == 'G')	{	data.equal++;	}
+					else if(sequence[structMap[i]] == 'C')	{	data.pound++;	}
+					else if(sequence[structMap[i]] == '-')	{	data.pound++;	}
+					data.total++;
+				}
+				else if(sequence[i] == 'T'){
+					if(sequence[structMap[i]] == 'T')		{	data.plus++;	}
+					else if(sequence[structMap[i]] == 'A')	{	data.tilde++;	}
+					else if(sequence[structMap[i]] == 'G')	{	data.dash++;	}
+					else if(sequence[structMap[i]] == 'C')	{	data.pound++;	}
+					else if(sequence[structMap[i]] == '-')	{	data.pound++;	}
+					data.total++;
+				}
+				else if(sequence[i] == 'G'){
+					if(sequence[structMap[i]] == 'T')		{	data.dash++;	}
+					else if(sequence[structMap[i]] == 'A')	{	data.equal++;	}
+					else if(sequence[structMap[i]] == 'G')	{	data.pound++;	}
+					else if(sequence[structMap[i]] == 'C')	{	data.tilde++;	}
+					else if(sequence[structMap[i]] == '-')	{	data.pound++;	}
+					data.total++;
+				}
+				else if(sequence[i] == 'C'){
+					if(sequence[structMap[i]] == 'T')		{	data.pound++;	}
+					else if(sequence[structMap[i]] == 'A')	{	data.pound++;	}
+					else if(sequence[structMap[i]] == 'G')	{	data.tilde++;	}
+					else if(sequence[structMap[i]] == 'C')	{	data.pound++;	}
+					else if(sequence[structMap[i]] == '-')	{	data.pound++;	}
+					data.total++;
+				}
+				else if(sequence[i] == '-'){
+					if(sequence[structMap[i]] == 'T')		{	data.pound++;	data.total++;	}
+					else if(sequence[structMap[i]] == 'A')	{	data.pound++;	data.total++;	}
+					else if(sequence[structMap[i]] == 'G')	{	data.pound++;	data.total++;	}
+					else if(sequence[structMap[i]] == 'C')	{	data.pound++;	data.total++;	}
+					else if(sequence[structMap[i]] == '-')	{		/*donothing*/				}
+				}			
+			}
+			else if(isalnum(sequence[i])){
+				data.loop++;
+				data.total++;
+			}
+		}
+		return data;
+		
+	}
+	catch(exception& e) {
+		errorOut(e, "AlignCheckCommand", "getStats");
+		exit(1);
+	}
+}
+
 
 //**********************************************************************************************************************
diff --git a/secondarystructurecommand.h b/secondarystructurecommand.h
index d94cc95..a82cf4e 100644
--- a/secondarystructurecommand.h
+++ b/secondarystructurecommand.h
@@ -44,6 +44,7 @@ class AlignCheckCommand : public Command {
 		int seqLength;
 		
 		void readMap();
+		statData getStats(string sequence);
 };
 
 /**************************************************************************************************/
-- 
2.39.2