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
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.
13 // ----------------------------------------
14 // ---------- sequence.hh -----------
15 // ----------------------------------------
17 #include <boost/filesystem/path.hpp>
18 #include <boost/filesystem/fstream.hpp>
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>
30 #include <boost/shared_ptr.hpp>
31 #include <boost/enable_shared_from_this.hpp>
35 #include "alphabet.hpp"
37 #include "seq_span.hpp"
39 // Sequence data class
41 //! Attach annotation information to a sequence track
45 annot(int begin, int end, std::string type, std::string name);
53 friend bool operator==(const annot& left, const annot& right);
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);
65 BOOST_CLASS_EXPORT(annot);
68 /* The way that motifs are found currently doesn't really
69 * indicate that the match was a reverse compliment
71 struct motif : public annot
75 motif() : annot(), sequence("") {};
76 //! this constructor is for when we're adding motifs to our annotations
77 motif(int begin, std::string motif);
80 // boost::serialization support
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);
89 BOOST_CLASS_EXPORT(motif);
92 typedef boost::shared_ptr<Sequence> SequenceRef;
94 //! sequence track for mussa.
98 typedef SeqString::value_type value_type;
99 typedef SeqString::difference_type difference_type;
100 typedef SeqString::iterator iterator;
101 typedef SeqString::reverse_iterator reverse_iterator;
102 typedef SeqString::const_iterator const_iterator;
103 typedef SeqString::const_reverse_iterator const_reverse_iterator;
104 typedef SeqString::reference reference;
105 typedef SeqString::const_reference const_reference;
106 typedef SeqString::size_type size_type;
107 static const size_type npos = SeqString::npos;
109 typedef std::list<motif> MotifList;
110 typedef boost::shared_ptr<MotifList> MotifListRef;
112 Sequence(AlphabetRef a = reduced_dna_alphabet);
113 Sequence(const char* seq,
114 AlphabetRef a = reduced_dna_alphabet,
115 SeqSpan::strand_type strand = SeqSpan::PlusStrand);
116 Sequence(const std::string& seq,
117 AlphabetRef a = reduced_dna_alphabet,
118 SeqSpan::strand_type strand = SeqSpan::PlusStrand);
119 //! make a new sequence, with the same SeqSpan
120 Sequence(const Sequence& seq);
121 //! make a new sequence, with the same SeqSpan
122 Sequence(const Sequence *);
123 //! Make a new sequence using a copy of SeqSpan
124 Sequence(const SequenceRef);
125 Sequence(const SeqSpanRef&);
127 //! assignment to constant sequences
128 Sequence &operator=(const Sequence&);
130 friend std::ostream& operator<<(std::ostream&, const Sequence&);
131 friend bool operator<(const Sequence&, const Sequence&);
132 friend bool operator==(const Sequence&, const Sequence&);
133 friend bool operator!=(const Sequence&, const Sequence&);
134 const_reference operator[](size_type) const;
136 //! set sequence to a (sub)string containing nothing but AGCTN
137 void set_filtered_sequence(const std::string& seq,
138 AlphabetRef a=reduced_dna_alphabet,
140 size_type count=npos,
141 SeqSpan::strand_type strand=SeqSpan::PlusStrand);
143 //! retrive element at specific position
144 const_reference at(size_type i) const { return seq->at(i); }
145 //! clear the sequence and its annotations
147 //! return a non-null terminated c pointer to the sequence data
148 const char *data() const { return seq->data(); }
150 const_iterator begin() const { return seq->begin(); }
152 const_iterator end() const { return seq->end(); }
153 //! is our sequence empty?
154 bool empty() const { return (seq) ? seq->empty() : true ; }
156 size_type find_first_not_of(const std::string& q, size_type index=0) { return seq->find_first_not_of(q, index); }
157 //! how many base pairs are there in our sequence
158 size_type size() const { return (seq) ? seq->size() : 0; }
159 //! alias for size (used by string)
160 size_type length() const { return size(); }
162 const_reverse_iterator rbegin() const { return seq->rbegin(); }
163 //! reverse end iterator
164 const_reverse_iterator rend() const { return seq->rend(); }
165 //! is our sequence empty?
166 //! start position relative to "base" sequence
167 size_type start() const { return seq->parentStart(); }
168 //! one past the last position relative to "base" sequence
169 size_type stop() const { return seq->parentStop(); }
171 //! return a subsequence, copying over any appropriate annotation
172 Sequence subseq(size_type start=0,
173 size_type count = npos,
174 SeqSpan::strand_type strand = SeqSpan::SameStrand) const;
175 //! reverse a character
176 std::string create_reverse_map() const;
177 //! return a reverse compliment (this needs to be improved?)
178 Sequence rev_comp() const;
180 //! set sequence (filtered)
181 void set_sequence(const std::string &, AlphabetRef);
183 std::string get_sequence() const;
185 void set_species(const std::string &);
187 std::string get_species() const;
188 //! set the fasta header
189 void set_fasta_header(std::string header);
190 //! get the fasta header
191 std::string get_fasta_header() const;
192 //! get name (will return the first non-empty, of fasta_header, species)
193 std::string get_name() const;
194 //! return a reference to whichever alphabet we're currently representing
195 const Alphabet& get_alphabet() const;
197 //! load sequence from fasta file using the sequences current alphabet
198 void load_fasta(const boost::filesystem::path file_path, int seq_num=1,
199 int start_index=0, int end_index=0);
200 //! load sequence AGCT from fasta file
201 //! \throw mussa_load_error
202 //! \throw sequence_empty_error
203 //! \throw sequence_empty_file_error
204 void load_fasta(const boost::filesystem::path file_path,
207 int start_index=0, int end_index=0);
208 void load_fasta(std::istream& file,
209 int seq_num=1, int start_index=0, int end_index=0);
210 //! load sequence from stream
211 //! \throw mussa_load_error
212 //! \throw sequence_empty_error
213 //! \throw sequence_empty_file_error
214 void load_fasta(std::istream& file,
217 int start_index=0, int end_index=0);
218 //! load sequence annotations
219 //! \throws mussa_load_error
220 void load_annot(const boost::filesystem::path file_path, int start_index, int end_index);
221 //! load sequence annotations
222 //! \throws mussa_load_error
223 void load_annot(std::fstream& data_stream, int start_index, int end_index);
224 //! parse annotation file
225 /*! \throws annotation_load_error
227 void parse_annot(std::string data, int start_index=0, int end_index=0);
228 //! add an annotation to our list of annotations
229 void add_annotation(const annot& a);
230 const std::list<annot>& annotations() const;
231 const MotifList& motifs() const;
233 //! add a motif to our list of motifs
234 void add_motif(const Sequence& a_motif);
235 //! clear our list of found motifs
237 //! search a sequence for a_motif
238 //! \throws motif_normalize_error if there's something wrong with a_motif
239 std::vector<int> find_motif(const Sequence& a_motif) const;
240 //! annotate the current sequence with other sequences
241 void find_sequences(std::list<Sequence>::iterator start,
242 std::list<Sequence>::iterator end);
244 void save(boost::filesystem::fstream &save_file);
245 void load_museq(boost::filesystem::path load_file_path, int seq_num);
254 //! store our oldstyle annotations
255 std::list<annot> annots;
256 //! a seperate list for motifs since we're currently not saving them
257 MotifListRef motif_list;
259 //! copy over all our annotation children
260 void copy_children(Sequence &, size_type start, size_type count) const;
262 void motif_scan(const Sequence& a_motif, std::vector<int> * motif_match_starts) const;
263 std::string rc_motif(std::string a_motif) const;
264 //! look for a string sequence type and and it to an annotation list
265 void add_string_annotation(std::string a_seq, std::string name);
267 // boost::serialization support
268 friend class boost::serialization::access;
269 template<class Archive>
270 void serialize(Archive& ar, const unsigned int /*version*/) {
271 ar & BOOST_SERIALIZATION_NVP(seq);
272 ar & BOOST_SERIALIZATION_NVP(header);
273 ar & BOOST_SERIALIZATION_NVP(species);
274 ar & BOOST_SERIALIZATION_NVP(annots);
275 ar & BOOST_SERIALIZATION_NVP(motif_list);
278 BOOST_CLASS_EXPORT(Sequence);