alg/sequence.hpp

   1 #ifndef _MUSSA_SEQUENCE_H_
   2 #define _MUSSA_SEQUENCE_H_
   3 //  This file is part of the Mussa source distribution.
   4 //  http://mussa.caltech.edu/
   5 //  Contact author: Tristan  De Buysscher, tristan@caltech.edu
   6
   7 // This program and all associated source code files are Copyright (C) 2005
   8 // the California Institute of Technology, Pasadena, CA, 91125 USA.  It is
   9 // under the GNU Public License; please see the included LICENSE.txt
  10 // file for more information, or contact Tristan directly.
  11
  12
  13 //                        ----------------------------------------
  14 //                           ---------- sequence.hh -----------
  15 //                        ----------------------------------------
  16
  17 #include <boost/filesystem/path.hpp>
  18 #include <boost/filesystem/fstream.hpp>
  19
  20 #include <boost/serialization/base_object.hpp>
  21 #include <boost/serialization/list.hpp>
  22 #include <boost/serialization/nvp.hpp>
  23 #include <boost/serialization/string.hpp>
  24 #include <boost/serialization/shared_ptr.hpp>
  25 #include <boost/serialization/utility.hpp>
  26 #include <boost/serialization/export.hpp>
  27
  28 #include <boost/shared_ptr.hpp>
  29
  30 #include <list>
  31 #include <string>
  32 #include <vector>
  33 #include <iostream>
  34
  35 // Sequence data class
  36
  37 //! Attach annotation information to a sequence track
  38 struct annot
  39 {
  40   annot();
  41   annot(int begin, int end, std::string type, std::string name);
  42   ~annot();
  43
  44   int begin;
  45   int end;
  46   std::string type;
  47   std::string name;
  48
  49   friend bool operator==(const annot& left, const annot& right);
  50 private:
  51   // boost::serialization support
  52   friend class boost::serialization::access;
  53   template<class Archive>
  54   void serialize(Archive& ar, const unsigned int /*version*/) {
  55     ar & BOOST_SERIALIZATION_NVP(begin);
  56     ar & BOOST_SERIALIZATION_NVP(end);
  57     ar & BOOST_SERIALIZATION_NVP(type);
  58     ar & BOOST_SERIALIZATION_NVP(name);
  59   }
  60 };
  61 BOOST_CLASS_EXPORT(annot);
  62
  63
  64 /* The way that motifs are found currently doesn't really
  65  * indicate that the match was a reverse compliment
  66  */
  67 struct motif : public annot
  68 {
  69   std::string sequence;
  70
  71   motif() :  annot(), sequence("") {};
  72   //! this constructor is for when we're adding motifs to our annotations
  73   motif(int begin, std::string motif);
  74   ~motif();
  75
  76   // boost::serialization support
  77 private:
  78   friend class boost::serialization::access;
  79   template<class Archive>
  80   void serialize(Archive& ar, const unsigned int /*version*/) {
  81     ar & BOOST_SERIALIZATION_BASE_OBJECT_NVP(annot);
  82     ar & BOOST_SERIALIZATION_NVP(sequence);
  83   }
  84 };
  85 BOOST_CLASS_EXPORT(motif);
  86
  87 //! sequence track for mussa.
  88 class Sequence
  89 {
  90 public:
  91   typedef std::string::difference_type difference_type;
  92   typedef std::string::iterator iterator;
  93   typedef std::string::const_iterator const_iterator;
  94   typedef std::string::reference reference;
  95   typedef std::string::const_reference const_reference;
  96   typedef std::string::size_type size_type;
  97   static const size_type npos = std::string::npos;
  98
  99   // some standard dna alphabets
 100   // Include nl (\012), and cr (\015) to make sequence parsing eol
 101   // convention independent.
 102
 103   static const std::string dna_alphabet;
 104   static const std::string rna_alphabet;
 105   //! this is the general iupac alphabet for nucleotides
 106   static const std::string nucleic_iupac_alphabet;
 107   //! the protein alphabet
 108   static const std::string protein_alphabet;
 109
 110   Sequence();
 111   ~Sequence();
 112   Sequence(const char* seq);
 113   Sequence(const std::string& seq);
 114   Sequence(const Sequence& seq);
 115   //! assignment to constant sequences
 116   Sequence &operator=(const Sequence&);
 117
 118   // dangerous since they create large copies
 119   //operator std::string();
 120   //operator std::string() const;
 121
 122   friend std::ostream& operator<<(std::ostream&, const Sequence&);
 123   friend bool operator<(const Sequence&, const Sequence&);
 124   friend bool operator==(const Sequence&, const Sequence&);
 125   const_reference operator[](size_type) const;
 126
 127   //! set sequence to a (sub)string containing nothing but AGCTN
 128   void set_filtered_sequence(const std::string& seq,
 129                              std::string::size_type start=0,
 130                              std::string::size_type count=0);
 131
 132   //! retrive element at specific position
 133   const_reference at(size_type n) const;
 134   //! clear the sequence and its annotations
 135   void clear();
 136   //! return c pointer to the sequence data
 137   const char *c_str() const;
 138   //! forward iterator
 139   const_iterator begin() const;
 140   //! last iterator
 141   const_iterator end() const;
 142   //! is our sequence empty?
 143   bool empty() const;
 144   //! how many base pairs are there in our sequence
 145   size_type size() const;
 146   //! alias for size (used by string)
 147   size_type length() const;
 148
 149   //! return a subsequence, copying over any appropriate annotation
 150   Sequence subseq(int start=0, int count = std::string::npos) const;
 151   //! return a reverse compliment
 152   std::string rev_comp() const;
 153
 154   //! set sequence (filtered)
 155   void set_sequence(const std::string &);
 156   //! get sequence
 157   std::string get_sequence() const;
 158   //! set species name
 159   void set_species(const std::string &);
 160   //! get species name
 161   std::string get_species() const;
 162   //! set the fasta header
 163   void set_fasta_header(std::string header);
 164   //! get the fasta header
 165   std::string get_fasta_header() const;
 166   //! get name (will return the first non-empty, of fasta_header, species)
 167   std::string get_name() const;
 168
 169   //! load sequence AGCT from fasta file
 170   //! \throw mussa_load_error
 171   //! \throw sequence_empty_error
 172   //! \throw sequence_empty_file_error
 173   void load_fasta(const boost::filesystem::path file_path, int seq_num=1,
 174                   int start_index=0, int end_index=0);
 175   //! load sequence from stream
 176   //! \throw mussa_load_error
 177   //! \throw sequence_empty_error
 178   //! \throw sequence_empty_file_error
 179   void load_fasta(std::iostream& file, int seq_num=1,
 180                   int start_index=0, int end_index=0);
 181   //! load sequence annotations
 182   //! \throws mussa_load_error
 183   void load_annot(const boost::filesystem::path file_path, int start_index, int end_index);
 184   //! load sequence annotations
 185   //! \throws mussa_load_error
 186   void load_annot(std::fstream& data_stream, int start_index, int end_index);
 187   bool parse_annot(std::string data, int start_index=0, int end_index=0);
 188   //! add an annotation to our list of annotations
 189   void add_annotation(const annot& a);
 190   const std::list<annot>& annotations() const;
 191   const std::list<motif>& motifs() const;
 192
 193   //! add a motif to our list of motifs
 194   //! \throws motif_normalize_error if there's something wrong with a_motif
 195   void add_motif(const Sequence& a_motif);
 196   //! clear our list of found motifs
 197   void clear_motifs();
 198   //! search a sequence for a_motif
 199   //! \throws motif_normalize_error if there's something wrong with a_motif
 200   std::vector<int> find_motif(const std::string& a_motif) const;
 201   //! search a sequence for a_motif
 202   //! \throws motif_normalize_error if there's something wrong with a_motif
 203   std::vector<int> find_motif(const Sequence& a_motif) const;
 204   //! convert IUPAC symbols to upperase
 205   //! \throws motif_normalize_error if there is an invalid symbol
 206   static std::string motif_normalize(const std::string& a_motif);
 207
 208   //! annotate the current sequence with other sequences
 209   void find_sequences(std::list<Sequence>::iterator start,
 210                       std::list<Sequence>::iterator end);
 211
 212   void save(boost::filesystem::fstream &save_file);
 213   void load_museq(boost::filesystem::path load_file_path, int seq_num);
 214
 215 private:
 216   boost::shared_ptr<const std::string> seq;
 217   std::string header;
 218   std::string species;
 219
 220   std::list<annot> annots;
 221   //! a seperate list for motifs since we're currently not saving them
 222   std::list<motif> motif_list;
 223
 224   void motif_scan(std::string a_motif, std::vector<int> * motif_match_starts) const;
 225   std::string rc_motif(std::string a_motif) const;
 226   //! look for a string sequence type and and it to an annotation list
 227   void add_string_annotation(std::string a_seq, std::string name);
 228
 229   // boost::serialization support
 230   friend class boost::serialization::access;
 231   template<class Archive>
 232   void serialize(Archive& ar, const unsigned int /*version*/) {
 233     ar & BOOST_SERIALIZATION_NVP(seq);
 234     ar & BOOST_SERIALIZATION_NVP(header);
 235     ar & BOOST_SERIALIZATION_NVP(species);
 236     ar & BOOST_SERIALIZATION_NVP(annots);
 237     ar & BOOST_SERIALIZATION_NVP(motif_list);
 238   }
 239 };
 240 //BOOST_CLASS_EXPORT(Sequence);
 241
 242 #endif