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