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