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