remove commented out string conversion
[mussa.git] / 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
113   // some standard dna alphabets 
114   // Include nl (\012), and cr (\015) to make sequence parsing eol 
115   // convention independent.
116
117   static const std::string dna_alphabet;
118   static const std::string rna_alphabet;
119   //! this is the general iupac alphabet for nucleotides
120   static const std::string nucleic_iupac_alphabet;
121   //! the protein alphabet
122   static const std::string protein_alphabet;
123
124   Sequence();
125   ~Sequence();
126   Sequence(const char* seq);
127   Sequence(const std::string& seq);
128   Sequence(const Sequence& seq);
129   //! assignment to constant sequences
130   Sequence &operator=(const Sequence&);
131
132   friend std::ostream& operator<<(std::ostream&, const Sequence&);
133   friend bool operator<(const Sequence&, const Sequence&);
134   friend bool operator==(const Sequence&, const Sequence&);
135   const_reference operator[](size_type) const;
136
137   //! set sequence to a (sub)string containing nothing but AGCTN
138   void set_filtered_sequence(const std::string& seq, 
139                              std::string::size_type start=0, 
140                              std::string::size_type count=0);
141
142   //! retrive element at specific position
143   const_reference at(size_type n) const;
144   //! clear the sequence and its annotations
145   void clear();
146   //! return c pointer to the sequence data
147   const char *c_str() const;
148   //! forward iterator
149   const_iterator begin() const;
150   //! last iterator
151   const_iterator end() const;
152   //! is our sequence empty?
153   bool empty() const;
154   //! how many base pairs are there in our sequence
155   size_type size() const;
156   //! alias for size (used by string)
157   size_type length() const;
158
159   //! return a subsequence, copying over any appropriate annotation
160   Sequence subseq(int start=0, int count = std::string::npos) const;
161   //! return a reverse compliment
162   std::string rev_comp() const;
163
164   //! set sequence (filtered)
165   void set_sequence(const std::string &);
166   //! get sequence
167   std::string get_sequence() const;
168   //! set species name
169   void set_species(const std::string &);
170   //! get species name
171   std::string get_species() const;
172   //! set the fasta header
173   void set_fasta_header(std::string header);
174   //! get the fasta header
175   std::string get_fasta_header() const;
176   //! get name (will return the first non-empty, of fasta_header, species)
177   std::string get_name() const;
178   
179   //! load sequence AGCT from fasta file
180   //! \throw mussa_load_error
181   //! \throw sequence_empty_error
182   //! \throw sequence_empty_file_error
183   void load_fasta(const boost::filesystem::path file_path, int seq_num=1, 
184                   int start_index=0, int end_index=0);
185   //! load sequence from stream
186   //! \throw mussa_load_error
187   //! \throw sequence_empty_error
188   //! \throw sequence_empty_file_error
189   void load_fasta(std::iostream& file, int seq_num=1, 
190                   int start_index=0, int end_index=0);
191   //! load sequence annotations
192   //! \throws mussa_load_error 
193   void load_annot(const boost::filesystem::path file_path, int start_index, int end_index);
194   //! load sequence annotations
195   //! \throws mussa_load_error 
196   void load_annot(std::fstream& data_stream, int start_index, int end_index);
197   bool parse_annot(std::string data, int start_index=0, int end_index=0);
198   //! add an annotation to our list of annotations
199   void add_annotation(const annot& a);
200   const std::list<annot>& annotations() const;
201   const std::list<motif>& motifs() const;
202
203   //! add a motif to our list of motifs
204   //! \throws motif_normalize_error if there's something wrong with a_motif
205   void add_motif(const Sequence& a_motif);
206   //! clear our list of found motifs
207   void clear_motifs();
208   //! search a sequence for a_motif
209   //! \throws motif_normalize_error if there's something wrong with a_motif
210   std::vector<int> find_motif(const std::string& a_motif) const;
211   //! search a sequence for a_motif
212   //! \throws motif_normalize_error if there's something wrong with a_motif
213   std::vector<int> find_motif(const Sequence& a_motif) const;
214   //! convert IUPAC symbols to upperase
215   //! \throws motif_normalize_error if there is an invalid symbol
216   static std::string motif_normalize(const std::string& a_motif);
217   
218   //! annotate the current sequence with other sequences
219   void find_sequences(std::list<Sequence>::iterator start, 
220                       std::list<Sequence>::iterator end);
221   
222   void save(boost::filesystem::fstream &save_file);
223   void load_museq(boost::filesystem::path load_file_path, int seq_num); 
224   
225 private:
226   boost::shared_ptr<seq_string> seq;
227   std::string header;
228   //! species name
229   std::string species;
230
231   //! store our oldstyle annotations
232   std::list<annot> annots;
233   //! a seperate list for motifs since we're currently not saving them
234   std::list<motif> motif_list;
235
236   void motif_scan(std::string a_motif, std::vector<int> * motif_match_starts) const;
237   std::string rc_motif(std::string a_motif) const;
238   //! look for a string sequence type and and it to an annotation list
239   void add_string_annotation(std::string a_seq, std::string name);
240   
241   // boost::serialization support
242   friend class boost::serialization::access;
243   template<class Archive>
244   void serialize(Archive& ar, const unsigned int /*version*/) {
245     ar & BOOST_SERIALIZATION_NVP(seq);
246     ar & BOOST_SERIALIZATION_NVP(header);
247     ar & BOOST_SERIALIZATION_NVP(species);
248     ar & BOOST_SERIALIZATION_NVP(annots);
249     ar & BOOST_SERIALIZATION_NVP(motif_list);
250   }
251 };
252 BOOST_CLASS_EXPORT(Sequence);
253
254 #endif