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   //! find first
 160   size_type find_first_not_of(const std::string&, size_type index=0);
 161   //! how many base pairs are there in our sequence
 162   size_type size() const;
 163   //! alias for size (used by string)
 164   size_type length() const;
 165   //! reverse iterator
 166   const_reverse_iterator rbegin() const;
 167   //! reverse end iterator
 168   const_reverse_iterator rend() const;
 169   //! is our sequence empty?
 170   //! start position relative to "base" sequence
 171   size_type start() const;
 172   //! one past the last position relative to "base" sequence
 173   size_type stop() const;
 174
 175   //! return a subsequence, copying over any appropriate annotation
 176   Sequence subseq(int start=0, int count = std::string::npos);
 177   //! reverse a character
 178   std::string create_reverse_map() const;
 179   //! return a reverse compliment (this needs to be improved?)
 180   Sequence rev_comp() const;
 181
 182   //! set sequence (filtered)
 183   void set_sequence(const std::string &, alphabet_ref);
 184   //! get sequence
 185   std::string get_sequence() const;
 186   //! set species name
 187   void set_species(const std::string &);
 188   //! get species name
 189   std::string get_species() const;
 190   //! set the fasta header
 191   void set_fasta_header(std::string header);
 192   //! get the fasta header
 193   std::string get_fasta_header() const;
 194   //! get name (will return the first non-empty, of fasta_header, species)
 195   std::string get_name() const;
 196   //! return a reference to whichever alphabet we're currently representing
 197   const Alphabet& get_alphabet() const;
 198   //! return a reference to whichever alphabet we're currently representing
 199   const Alphabet& get_alphabet(alphabet_ref) const;
 200
 201   //! load sequence from fasta file using the sequences current alphabet
 202   void load_fasta(const boost::filesystem::path file_path, int seq_num=1,
 203                   int start_index=0, int end_index=0);
 204   //! load sequence AGCT from fasta file
 205   //! \throw mussa_load_error
 206   //! \throw sequence_empty_error
 207   //! \throw sequence_empty_file_error
 208   void load_fasta(const boost::filesystem::path file_path,
 209                   alphabet_ref a,
 210                   int seq_num=1,
 211                               int start_index=0, int end_index=0);
 212   void load_fasta(std::istream& file,
 213                   int seq_num=1, int start_index=0, int end_index=0);
 214   //! load sequence from stream
 215   //! \throw mussa_load_error
 216   //! \throw sequence_empty_error
 217   //! \throw sequence_empty_file_error
 218   void load_fasta(std::istream& file,
 219                   alphabet_ref a,
 220                   int seq_num=1,
 221                               int start_index=0, int end_index=0);
 222   //! load sequence annotations
 223   //! \throws mussa_load_error
 224   void load_annot(const boost::filesystem::path file_path, int start_index, int end_index);
 225   //! load sequence annotations
 226   //! \throws mussa_load_error
 227   void load_annot(std::fstream& data_stream, int start_index, int end_index);
 228   //! parse annotation file
 229   /*! \throws annotation_load_error
 230    */
 231   void parse_annot(std::string data, int start_index=0, int end_index=0);
 232   //! add an annotation to our list of annotations
 233   void add_annotation(const annot& a);
 234   const std::list<annot>& annotations() const;
 235   const std::list<motif>& motifs() const;
 236
 237   //! add a motif to our list of motifs
 238   void add_motif(const Sequence& a_motif);
 239   //! clear our list of found motifs
 240   void clear_motifs();
 241   //! search a sequence for a_motif
 242   //! \throws motif_normalize_error if there's something wrong with a_motif
 243   std::vector<int> find_motif(const Sequence& a_motif) const;
 244   //! annotate the current sequence with other sequences
 245   void find_sequences(std::list<Sequence>::iterator start,
 246                       std::list<Sequence>::iterator end);
 247
 248   void save(boost::filesystem::fstream &save_file);
 249   void load_museq(boost::filesystem::path load_file_path, int seq_num);
 250
 251 private:
 252   //! parent sequence
 253   Sequence *parent;
 254   //! hold a shared pointer to our sequence string
 255   boost::shared_ptr<seq_string> seq;
 256   //! which alphabet we're using
 257   alphabet_ref alphabet;
 258   //! start offset into the sequence
 259   size_type seq_start;
 260   //! number of basepairs of the shared sequence we represent
 261   size_type seq_count;
 262   //! strand orientation
 263   strand_type strand;
 264   //! fasta header
 265   std::string header;
 266   //! species name
 267   std::string species;
 268
 269   //! store our oldstyle annotations
 270   std::list<annot> annots;
 271   //! a seperate list for motifs since we're currently not saving them
 272   std::list<motif> motif_list;
 273
 274   void motif_scan(const Sequence& a_motif, std::vector<int> * motif_match_starts) const;
 275   std::string rc_motif(std::string a_motif) const;
 276   //! look for a string sequence type and and it to an annotation list
 277   void add_string_annotation(std::string a_seq, std::string name);
 278
 279   // boost::serialization support
 280   friend class boost::serialization::access;
 281   template<class Archive>
 282   void serialize(Archive& ar, const unsigned int /*version*/) {
 283     ar & BOOST_SERIALIZATION_NVP(parent);
 284     ar & BOOST_SERIALIZATION_NVP(seq);
 285     ar & BOOST_SERIALIZATION_NVP(alphabet);
 286     ar & BOOST_SERIALIZATION_NVP(seq_start);
 287     ar & BOOST_SERIALIZATION_NVP(seq_count);
 288     ar & BOOST_SERIALIZATION_NVP(strand);
 289     ar & BOOST_SERIALIZATION_NVP(header);
 290     ar & BOOST_SERIALIZATION_NVP(species);
 291     ar & BOOST_SERIALIZATION_NVP(annots);
 292     ar & BOOST_SERIALIZATION_NVP(motif_list);
 293   }
 294 };
 295 BOOST_CLASS_EXPORT(Sequence);
 296
 297 #endif