align.check command added

diff --git a/commandfactory.cpp b/commandfactory.cpp
index fce55bbe9e856634a25b8141a3715965e88251b0..dd5b36317deae9007ccf8951f8133f573e75e316 100644 (file)
--- 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 48ea1549aa78537966745dbe11fbc9d1afbf8f51..7156082f8b1b43228a62c6cf93abec9f2d900109 100644 (file)
--- 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 d94cc95582f8d1482233a3cac0aee4b418dcc925..a82cf4e5edd89a5536cfbecafb5def8f6a7e3ed4 100644 (file)
--- a/secondarystructurecommand.h
+++ b/secondarystructurecommand.h
@@ -44,6 +44,7 @@ class AlignCheckCommand : public Command {
                int seqLength;
                
                void readMap();
+               statData getStats(string sequence);
 };
 
 /**************************************************************************************************/