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/export.hpp>
  22 #include <boost/serialization/list.hpp>
  23 #include <boost/serialization/nvp.hpp>
  24 #include <boost/serialization/string.hpp>
  25 #include <boost/serialization/shared_ptr.hpp>
  26 #include <boost/serialization/utility.hpp>
  27 #include <boost/serialization/version.hpp>
  28 #include <boost/serialization/vector.hpp>
  29
  30 #include <boost/shared_ptr.hpp>
  31
  32 #include <iostream>
  33
  34 #include "alg/alphabet.hpp"
  35
  36 // Sequence data class
  37
  38 //! Attach annotation information to a sequence track
  39 struct annot
  40 {
  41   annot();
  42   annot(int begin, int end, std::string type, std::string name);
  43   ~annot();
  44
  45   int begin;
  46   int end;
  47   std::string type;
  48   std::string name;
  49
  50   friend bool operator==(const annot& left, const annot& right);
  51 private:
  52   // boost::serialization support
  53   friend class boost::serialization::access;
  54   template<class Archive>
  55   void serialize(Archive& ar, const unsigned int /*version*/) {
  56     ar & BOOST_SERIALIZATION_NVP(begin);
  57     ar & BOOST_SERIALIZATION_NVP(end);
  58     ar & BOOST_SERIALIZATION_NVP(type);
  59     ar & BOOST_SERIALIZATION_NVP(name);
  60   }
  61 };
  62 BOOST_CLASS_EXPORT(annot);
  63
  64
  65 /* The way that motifs are found currently doesn't really
  66  * indicate that the match was a reverse compliment
  67  */
  68 struct motif : public annot
  69 {
  70   std::string sequence;
  71
  72   motif() :  annot(), sequence("") {};
  73   //! this constructor is for when we're adding motifs to our annotations
  74   motif(int begin, std::string motif);
  75   ~motif();
  76
  77   // boost::serialization support
  78 private:
  79   friend class boost::serialization::access;
  80   template<class Archive>
  81   void serialize(Archive& ar, const unsigned int /*version*/) {
  82     ar & BOOST_SERIALIZATION_BASE_OBJECT_NVP(annot);
  83     ar & BOOST_SERIALIZATION_NVP(sequence);
  84   }
  85 };
  86 BOOST_CLASS_EXPORT(motif);
  87
  88 //! the only purpose of this class is that the shared_ptr template
  89 //! functions need the serialization support to be in-class.
  90 class seq_string : public std::string
  91 {
  92 public:
  93   typedef std::string::iterator iterator;
  94   typedef std::string::reverse_iterator reverse_iterator;
  95   typedef std::string::const_iterator const_iterator;
  96   typedef std::string::const_reverse_iterator const_reverse_iterator;
  97 private:
  98   friend class boost::serialization::access;
  99   template<class Archive>
 100   void serialize(Archive& ar, const unsigned int /*version*/) {
 101     //ar & BOOST_SERIALIZATION_BASE_OBJECT_NVP(std::string);
 102     ar & boost::serialization::make_nvp("bases",
 103            boost::serialization::base_object<std::string>(*this)
 104     );
 105   }
 106 };
 107
 108 //! sequence track for mussa.
 109 class Sequence
 110 {
 111 public:
 112   typedef std::string::value_type value_type;
 113   typedef std::string::difference_type difference_type;
 114   typedef std::string::iterator iterator;
 115   typedef std::string::reverse_iterator reverse_iterator;
 116   typedef std::string::const_iterator const_iterator;
 117   typedef std::string::const_reverse_iterator const_reverse_iterator;
 118   typedef std::string::reference reference;
 119   typedef std::string::const_reference const_reference;
 120   typedef std::string::size_type size_type;
 121   static const size_type npos = std::string::npos;
 122   enum strand_type { UnknownStrand, PlusStrand, MinusStrand, BothStrand };
 123   enum alphabet_ref { reduced_dna_alphabet, reduced_rna_alphabet, reduced_nucleic_alphabet,
 124                       nucleic_alphabet, protein_alphabet };
 125
 126   Sequence(alphabet_ref a = reduced_nucleic_alphabet);
 127   Sequence(const char* seq, alphabet_ref a = reduced_nucleic_alphabet);
 128   Sequence(const std::string& seq, alphabet_ref a = reduced_nucleic_alphabet);
 129   Sequence(const Sequence& seq);
 130   ~Sequence();
 131   //! assignment to constant sequences
 132   Sequence &operator=(const Sequence&);
 133
 134   friend std::ostream& operator<<(std::ostream&, const Sequence&);
 135   friend bool operator<(const Sequence&, const Sequence&);
 136   friend bool operator==(const Sequence&, const Sequence&);
 137   friend bool operator!=(const Sequence&, const Sequence&);
 138   const_reference operator[](size_type) const;
 139
 140   //! set sequence to a (sub)string containing nothing but AGCTN
 141   void set_filtered_sequence(const std::string& seq,
 142                              alphabet_ref a,
 143                              size_type start=0,
 144                              size_type count=npos,
 145                              strand_type strand=UnknownStrand);
 146
 147   //! retrive element at specific position
 148   const_reference at(size_type n) const;
 149   //! clear the sequence and its annotations
 150   void clear();
 151   //! return c pointer to the sequence data
 152   const char *c_str() const;
 153   //! forward iterator
 154   const_iterator begin() const;
 155   //! last iterator
 156   const_iterator end() const;
 157   //! is our sequence empty?
 158   bool empty() const;
 159   //! how many base pairs are there in our sequence
 160   size_type size() const;
 161   //! alias for size (used by string)
 162   size_type length() const;
 163   //! start position relative to "base" sequence
 164   size_type start() const;
 165   //! one past the last position relative to "base" sequence
 166   size_type stop() const;
 167
 168   //! return a subsequence, copying over any appropriate annotation
 169   Sequence subseq(int start=0, int count = std::string::npos);
 170   //! reverse a character
 171   std::string create_reverse_map() const;
 172   //! return a reverse compliment (this needs to be improved?)
 173   Sequence rev_comp() const;
 174
 175   //! set sequence (filtered)
 176   void set_sequence(const std::string &, alphabet_ref);
 177   //! get sequence
 178   std::string get_sequence() const;
 179   //! set species name
 180   void set_species(const std::string &);
 181   //! get species name
 182   std::string get_species() const;
 183   //! set the fasta header
 184   void set_fasta_header(std::string header);
 185   //! get the fasta header
 186   std::string get_fasta_header() const;
 187   //! get name (will return the first non-empty, of fasta_header, species)
 188   std::string get_name() const;
 189   //! return a reference to whichever alphabet we're currently representing
 190   const Alphabet& get_alphabet() const;
 191   //! return a reference to whichever alphabet we're currently representing
 192   const Alphabet& get_alphabet(alphabet_ref) const;
 193
 194   //! load sequence from fasta file using the sequences current alphabet
 195   void load_fasta(const boost::filesystem::path file_path, int seq_num=1,
 196                   int start_index=0, int end_index=0);
 197   //! load sequence AGCT from fasta file
 198   //! \throw mussa_load_error
 199   //! \throw sequence_empty_error
 200   //! \throw sequence_empty_file_error
 201   void load_fasta(const boost::filesystem::path file_path,
 202                   alphabet_ref a,
 203                   int seq_num=1,
 204                               int start_index=0, int end_index=0);
 205   void load_fasta(std::iostream& file,
 206                   int seq_num=1, int start_index=0, int end_index=0);
 207   //! load sequence from stream
 208   //! \throw mussa_load_error
 209   //! \throw sequence_empty_error
 210   //! \throw sequence_empty_file_error
 211   void load_fasta(std::iostream& file,
 212                   alphabet_ref a,
 213                   int seq_num=1,
 214                               int start_index=0, int end_index=0);
 215   //! load sequence annotations
 216   //! \throws mussa_load_error
 217   void load_annot(const boost::filesystem::path file_path, int start_index, int end_index);
 218   //! load sequence annotations
 219   //! \throws mussa_load_error
 220   void load_annot(std::fstream& data_stream, int start_index, int end_index);
 221   //! parse annotation file
 222   /*! \throws annotation_load_error
 223    */
 224   void parse_annot(std::string data, int start_index=0, int end_index=0);
 225   //! add an annotation to our list of annotations
 226   void add_annotation(const annot& a);
 227   const std::list<annot>& annotations() const;
 228   const std::list<motif>& motifs() const;
 229
 230   //! add a motif to our list of motifs
 231   void add_motif(const Sequence& a_motif);
 232   //! clear our list of found motifs
 233   void clear_motifs();
 234   //! search a sequence for a_motif
 235   //! \throws motif_normalize_error if there's something wrong with a_motif
 236   std::vector<int> find_motif(const Sequence& a_motif) const;
 237   //! annotate the current sequence with other sequences
 238   void find_sequences(std::list<Sequence>::iterator start,
 239                       std::list<Sequence>::iterator end);
 240
 241   void save(boost::filesystem::fstream &save_file);
 242   void load_museq(boost::filesystem::path load_file_path, int seq_num);
 243
 244 private:
 245   //! parent sequence
 246   Sequence *parent;
 247   //! hold a shared pointer to our sequence string
 248   boost::shared_ptr<seq_string> seq;
 249   //! which alphabet we're using
 250   alphabet_ref alphabet;
 251   //! start offset into the sequence
 252   size_type seq_start;
 253   //! number of basepairs of the shared sequence we represent
 254   size_type seq_count;
 255   //! strand orientation
 256   strand_type strand;
 257   //! fasta header
 258   std::string header;
 259   //! species name
 260   std::string species;
 261
 262   //! store our oldstyle annotations
 263   std::list<annot> annots;
 264   //! a seperate list for motifs since we're currently not saving them
 265   std::list<motif> motif_list;
 266
 267   void motif_scan(const Sequence& a_motif, std::vector<int> * motif_match_starts) const;
 268   std::string rc_motif(std::string a_motif) const;
 269   //! look for a string sequence type and and it to an annotation list
 270   void add_string_annotation(std::string a_seq, std::string name);
 271
 272   // boost::serialization support
 273   friend class boost::serialization::access;
 274   template<class Archive>
 275   void serialize(Archive& ar, const unsigned int /*version*/) {
 276     ar & BOOST_SERIALIZATION_NVP(parent);
 277     ar & BOOST_SERIALIZATION_NVP(seq);
 278     ar & BOOST_SERIALIZATION_NVP(alphabet);
 279     ar & BOOST_SERIALIZATION_NVP(seq_start);
 280     ar & BOOST_SERIALIZATION_NVP(seq_count);
 281     ar & BOOST_SERIALIZATION_NVP(strand);
 282     ar & BOOST_SERIALIZATION_NVP(header);
 283     ar & BOOST_SERIALIZATION_NVP(species);
 284     ar & BOOST_SERIALIZATION_NVP(annots);
 285     ar & BOOST_SERIALIZATION_NVP(motif_list);
 286   }
 287 };
 288 BOOST_CLASS_EXPORT(Sequence);
 289
 290 #endif