implement alphabet for sequence
[mussa.git] / alg / sequence.hpp
index 5010729243fc8b5caae76424419789db81c0b866..dc9d594a1f572753f73ced6a5de056db5b6c0209 100644 (file)
@@ -31,6 +31,8 @@
 
 #include <iostream>
 
+#include "alg/alphabet.hpp"
+
 // Sequence data class
 
 //! Attach annotation information to a sequence track
@@ -87,6 +89,11 @@ BOOST_CLASS_EXPORT(motif);
 //! functions need the serialization support to be in-class.
 class seq_string : public std::string
 {
+public:
+  typedef std::string::iterator iterator;
+  typedef std::string::reverse_iterator reverse_iterator;
+  typedef std::string::const_iterator const_iterator;
+  typedef std::string::const_reverse_iterator const_reverse_iterator;
 private:
   friend class boost::serialization::access;
   template<class Archive>
@@ -94,7 +101,7 @@ private:
     //ar & BOOST_SERIALIZATION_BASE_OBJECT_NVP(std::string);
     ar & boost::serialization::make_nvp("bases",
           boost::serialization::base_object<std::string>(*this)
-         );
+    );
   }
 };
 
@@ -102,43 +109,39 @@ private:
 class Sequence 
 {
 public:
+  typedef std::string::value_type value_type;
   typedef std::string::difference_type difference_type;
   typedef std::string::iterator iterator;
+  typedef std::string::reverse_iterator reverse_iterator;
   typedef std::string::const_iterator const_iterator;
+  typedef std::string::const_reverse_iterator const_reverse_iterator;
   typedef std::string::reference reference;
   typedef std::string::const_reference const_reference;
   typedef std::string::size_type size_type;
   static const size_type npos = std::string::npos;
   enum strand_type { UnknownStrand, PlusStrand, MinusStrand, BothStrand };
-
-  // some standard dna alphabets 
-  // Include nl (\012), and cr (\015) to make sequence parsing eol 
-  // convention independent.
-
-  static const std::string dna_alphabet;
-  static const std::string rna_alphabet;
-  //! this is the general iupac alphabet for nucleotides
-  static const std::string nucleic_iupac_alphabet;
-  //! the protein alphabet
-  static const std::string protein_alphabet;
-
-  Sequence();
-  ~Sequence();
-  Sequence(const char* seq);
-  Sequence(const std::string& seq);
+  enum alphabet_ref { reduced_dna_alphabet, reduced_rna_alphabet, reduced_nucleic_alphabet, 
+                      nucleic_alphabet, protein_alphabet };
+                      
+  Sequence(alphabet_ref a = reduced_nucleic_alphabet);
+  Sequence(const char* seq, alphabet_ref a = reduced_nucleic_alphabet);
+  Sequence(const std::string& seq, alphabet_ref a = reduced_nucleic_alphabet);
   Sequence(const Sequence& seq);
+  ~Sequence();
   //! assignment to constant sequences
   Sequence &operator=(const Sequence&);
 
   friend std::ostream& operator<<(std::ostream&, const Sequence&);
   friend bool operator<(const Sequence&, const Sequence&);
   friend bool operator==(const Sequence&, const Sequence&);
+  friend bool operator!=(const Sequence&, const Sequence&);
   const_reference operator[](size_type) const;
 
   //! set sequence to a (sub)string containing nothing but AGCTN
-  void set_filtered_sequence(const std::string& seq, 
-                            size_type start=0,
-                            size_type count=npos,
+  void set_filtered_sequence(const std::string& seq,
+                             alphabet_ref a, 
+                             size_type start=0,
+                             size_type count=npos,
                              strand_type strand=UnknownStrand);
 
   //! retrive element at specific position
@@ -164,11 +167,13 @@ public:
 
   //! return a subsequence, copying over any appropriate annotation
   Sequence subseq(int start=0, int count = std::string::npos);
+  //! reverse a character
+  std::string create_reverse_map() const;
   //! return a reverse compliment (this needs to be improved?)
-  std::string rev_comp() const;
+  Sequence rev_comp() const;
 
   //! set sequence (filtered)
-  void set_sequence(const std::string &);
+  void set_sequence(const std::string &, alphabet_ref);
   //! get sequence
   std::string get_sequence() const;
   //! set species name
@@ -181,19 +186,32 @@ public:
   std::string get_fasta_header() const;
   //! get name (will return the first non-empty, of fasta_header, species)
   std::string get_name() const;
+  //! return a reference to whichever alphabet we're currently representing
+  const Alphabet& get_alphabet() const; 
+  //! return a reference to whichever alphabet we're currently representing
+  const Alphabet& get_alphabet(alphabet_ref) const; 
   
+  //! load sequence from fasta file using the sequences current alphabet
+  void load_fasta(const boost::filesystem::path file_path, int seq_num=1,
+                  int start_index=0, int end_index=0);
   //! load sequence AGCT from fasta file
   //! \throw mussa_load_error
   //! \throw sequence_empty_error
   //! \throw sequence_empty_file_error
-  void load_fasta(const boost::filesystem::path file_path, int seq_num=1, 
-                 int start_index=0, int end_index=0);
+  void load_fasta(const boost::filesystem::path file_path, 
+                  alphabet_ref a, 
+                  int seq_num=1, 
+                             int start_index=0, int end_index=0);
+  void load_fasta(std::iostream& file, 
+                  int seq_num=1, int start_index=0, int end_index=0);
   //! load sequence from stream
   //! \throw mussa_load_error
   //! \throw sequence_empty_error
   //! \throw sequence_empty_file_error
-  void load_fasta(std::iostream& file, int seq_num=1, 
-                 int start_index=0, int end_index=0);
+  void load_fasta(std::iostream& file, 
+                  alphabet_ref a,
+                  int seq_num=1, 
+                             int start_index=0, int end_index=0);
   //! load sequence annotations
   //! \throws mussa_load_error 
   void load_annot(const boost::filesystem::path file_path, int start_index, int end_index);
@@ -207,20 +225,12 @@ public:
   const std::list<motif>& motifs() const;
 
   //! add a motif to our list of motifs
-  //! \throws motif_normalize_error if there's something wrong with a_motif
   void add_motif(const Sequence& a_motif);
   //! clear our list of found motifs
   void clear_motifs();
   //! search a sequence for a_motif
   //! \throws motif_normalize_error if there's something wrong with a_motif
-  std::vector<int> find_motif(const std::string& a_motif) const;
-  //! search a sequence for a_motif
-  //! \throws motif_normalize_error if there's something wrong with a_motif
-  std::vector<int> find_motif(const Sequence& a_motif) const;
-  //! convert IUPAC symbols to upperase
-  //! \throws motif_normalize_error if there is an invalid symbol
-  static std::string motif_normalize(const std::string& a_motif);
-  
+  std::vector<int> find_motif(const Sequence& a_motif) const;  
   //! annotate the current sequence with other sequences
   void find_sequences(std::list<Sequence>::iterator start, 
                      std::list<Sequence>::iterator end);
@@ -233,6 +243,8 @@ private:
   Sequence *parent;
   //! hold a shared pointer to our sequence string
   boost::shared_ptr<seq_string> seq;
+  //! which alphabet we're using
+  alphabet_ref alphabet;
   //! start offset into the sequence
   size_type seq_start;
   //! number of basepairs of the shared sequence we represent
@@ -249,7 +261,7 @@ private:
   //! a seperate list for motifs since we're currently not saving them
   std::list<motif> motif_list;
 
-  void motif_scan(std::string a_motif, std::vector<int> * motif_match_starts) const;
+  void motif_scan(const Sequence& a_motif, std::vector<int> * motif_match_starts) const;
   std::string rc_motif(std::string a_motif) const;
   //! look for a string sequence type and and it to an annotation list
   void add_string_annotation(std::string a_seq, std::string name);
@@ -260,6 +272,7 @@ private:
   void serialize(Archive& ar, const unsigned int /*version*/) {
     ar & BOOST_SERIALIZATION_NVP(parent);
     ar & BOOST_SERIALIZATION_NVP(seq);
+    ar & BOOST_SERIALIZATION_NVP(alphabet);
     ar & BOOST_SERIALIZATION_NVP(seq_start);
     ar & BOOST_SERIALIZATION_NVP(seq_count);
     ar & BOOST_SERIALIZATION_NVP(strand);