SHORE API
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
nucleotide.hpp
Go to the documentation of this file.
1 
2 /*
3  * Copyright 2008,2009,2010,2011,2012 Stephan Ossowski, Korbinian Schneeberger,
4  * Felix Ott, Joerg Hagmann, Alf Scotland, Sebastian Bender
5  *
6  * This file is part of SHORE.
7  *
8  * SHORE is free software: you can redistribute it and/or modify
9  * it under the terms of the GNU General Public License as published by
10  * the Free Software Foundation, either version 3 of the License, or
11  * (at your option) any later version.
12  *
13  * SHORE is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  * GNU General Public License for more details.
17  *
18  * You should have received a copy of the GNU General Public License
19  * along with SHORE. If not, see <http://www.gnu.org/licenses/>.
20  */
21 
25 
26 #ifndef SHORE_GENOME_NUCLEOTIDE_HPP__
27 #define SHORE_GENOME_NUCLEOTIDE_HPP__
28 
29 #include <cctype>
30 
31 #include <iostream>
32 #include <iterator>
33 #include <string>
34 
35 #include <boost/utility/enable_if.hpp>
36 #include <boost/type_traits.hpp>
37 
38 #include "shore/base/memops.hpp"
39 #include "shore/base/stringops.hpp"
40 
41 namespace shore {
42 
44 class nuc
45 {
46  public:
47 
49  // bitmask IUPAC abbrev.
50  enum base
51  { // GTCA
52  NONE=0, // 0000 -
53  ADENINE=1, // 0001 A
54  CYTOSINE=2, // 0010 C
55  AMINO=3, // 0011 M
56  THYMINE=4, // 0100 T
57  WEAK=5, // 0101 W
58  PYRIMIDINE=6, // 0110 Y
59  NOT_GUANINE=7, // 0111 H
60  GUANINE=8, // 1000 G
61  PURINE=9, // 1001 R
62  STRONG=10, // 1010 S
63  NOT_THYMINE=11, // 1011 V
64  KETO=12, // 1100 K
65  NOT_CYTOSINE=13,// 1101 D
66  NOT_ADENINE=14, // 1110 B
67  ANY=15 // 1111 N
68  };
69 
71  enum packed_base
72  {
73  PCK_ADENINE=0, // 0000
74  PCK_CYTOSINE=1, // 0001
75  PCK_GUANINE=2, // 0010
76  PCK_THYMINE=3, // 0011
77  PCK_ANY=4, // 0100
78  PCK_NONE=5, // 0101
79  PCK_PURINE=6, // 0110
80  PCK_PYRIMIDINE=7, // 0111
81  PCK_KETO=8, // 1000
82  PCK_AMINO=9, // 1001
83  PCK_WEAK=10, // 1010
84  PCK_STRONG=11, // 1011
85  PCK_NOT_THYMINE=12, // 1100
86  PCK_NOT_GUANINE=13, // 1101
87  PCK_NOT_CYTOSINE=14,// 1110
88  PCK_NOT_ADENINE=15 // 1111
89  };
90 
91  private:
92 
93  static const char*const m_abbreviations;
94  static const std::string m_basenameTable[16];
95 
96  static const packed_base PACK_LUT[16];
97  static const base UNPACK_LUT[16];
98 
100  static const int m_gctable[16];
101 
103  nuc();
104 
105  class symtab
106  {
107  public:
108 
109  base LUT[512];
110 
111  symtab()
112  {
113  for(int i=0;i<256;++i)
114  {
115  const char abbr=std::toupper(char(i));
116  switch(abbr)
117  {
118  case 'A':
119  LUT[i]=nuc::ADENINE; break;
120  case 'C':
121  LUT[i]=nuc::CYTOSINE; break;
122  case 'G':
123  LUT[i]=nuc::GUANINE; break;
124  case 'T':
125  case 'U':
126  LUT[i]=nuc::THYMINE; break;
127  case 'R':
128  LUT[i]=nuc::PURINE; break;
129  case 'Y':
130  LUT[i]=nuc::PYRIMIDINE; break;
131  case 'W':
132  LUT[i]=nuc::WEAK; break;
133  case 'S':
134  LUT[i]=nuc::STRONG; break;
135  case 'M':
136  LUT[i]=nuc::AMINO; break;
137  case 'K':
138  LUT[i]=nuc::KETO; break;
139  case 'H':
140  LUT[i]=nuc::NOT_GUANINE; break;
141  case 'B':
142  LUT[i]=nuc::NOT_ADENINE; break;
143  case 'V':
144  LUT[i]=nuc::NOT_THYMINE; break;
145  case 'D':
146  LUT[i]=nuc::NOT_CYTOSINE; break;
147  case 'N':
148  case '.':
149  LUT[i]=nuc::ANY; break;
150  case '-':
151  default:
152  LUT[i]=nuc::NONE;
153  }
154  }
155  }
156  };
157 
158  static const symtab m_symtab;
159 
160  public:
161 
163  static char enc(const base b)
164  {
165  return m_abbreviations[b];
166  }
167 
169  static std::string s_enc(const base b)
170  {
171  return m_basenameTable[b];
172  }
173 
175  static char rna_enc(const base b)
176  {
177  return (b==THYMINE)?'U':enc(b);
178  }
179 
181  static std::string rna_s_enc(const base b)
182  {
183  return (b==THYMINE)?"uracil":s_enc(b);
184  }
185 
187  static base comp(const base b)
188  {
189  // rotl 2 bits
190  const unsigned char low=b>>2;
191  const unsigned char high=(b<<2)&12;
192  return base(high|low);
193  }
194 
196  static base km_class(const base b)
197  {
198  return base((b&AMINO)?(b^AMINO):(b^KETO));
199  }
200 
202  static base ry_class(const base b)
203  {
204  return base((b&PURINE)?(b^PURINE):(b^PYRIMIDINE));
205  }
206 
210  static base tr(const base b,const base replacewhat,const base replacewith)
211  {
212  if((b&replacewhat)==0)
213  return b;
214  return static_cast<base>((b&(~replacewhat))|replacewith);
215  }
216 
221  static base dec(const unsigned char abbr)
222  {
223  return m_symtab.LUT[abbr];
224  }
225 
227  static base s_dec(const std::string& name)
228  {
229  return static_cast<base>(
230  shore::match_arg_idx<std::runtime_error>(name,
231  m_basenameTable,
232  shore::end(m_basenameTable)));
233  }
234 
236  static int amb(const base b)
237  {
238  return (b&1)+((b>>1)&1)+((b>>2)&1)+((b>>3)&1);
239  }
240 
247  static int gc(const base b)
248  {
249  return m_gctable[b];
250  }
251 
252  static packed_base pack(const base b)
253  {
254  return PACK_LUT[b];
255  }
256 
257  static base unpack(const packed_base b)
258  {
259  return UNPACK_LUT[b];
260  }
261 
263  template<typename Iter>
264  static bool is_nucleic_acid(Iter beg,Iter end,
265  typename boost::enable_if<boost::is_same<typename std::iterator_traits<Iter>::value_type,
266  char> >::type *dummy=0)
267  {
268  while(beg!=end)
269  if(nuc::amb(nuc::dec(*beg++))==0)
270  return false;
271  return true;
272  }
273 
276  template<typename Iter>
277  static bool is_plain_nucleic_acid(Iter beg,Iter end,
278  typename boost::enable_if<boost::is_same<typename std::iterator_traits<Iter>::value_type,
279  char> >::type *dummy=0)
280  {
281  while(beg!=end)
282  if(nuc::amb(nuc::dec(*beg++))!=1)
283  return false;
284  return true;
285  }
286 
288  template<typename Iter>
289  static bool is_rc_palindrome(Iter b,Iter e,
290  typename boost::enable_if<boost::is_same<typename std::iterator_traits<Iter>::value_type,
291  base> >::type *dummy=0)
292  {
293  if((e-b)&1)
294  return false; // must have an even number of bases
295  while(b!=e)
296  {
297  --e;
298  if((*b)!=comp(*e))
299  return false;
300  ++b;
301  }
302  return true;
303  }
304 
306  template<typename Iter>
307  static bool is_rc_palindrome(Iter b,Iter e,
308  typename boost::enable_if<boost::is_same<typename std::iterator_traits<Iter>::value_type,
309  packed_base> >::type *dummy=0)
310  {
311  if((e-b)&1)
312  return false; // must have an even number of bases
313  while(b!=e)
314  {
315  --e;
316  if(unpack(*b)!=comp(unpack(*e)))
317  return false;
318  ++b;
319  }
320  return true;
321  }
322 
324  template<typename Iter>
325  static bool is_rc_palindrome(Iter b,Iter e,
326  typename boost::enable_if<boost::is_same<typename std::iterator_traits<Iter>::value_type,
327  char> >::type *dummy=0)
328  {
329  if((e-b)&1)
330  return false; // must have an even number of bases
331  while(b!=e)
332  {
333  --e;
334  if(dec(*b)!=comp(dec(*e)))
335  return false;
336  ++b;
337  }
338  return true;
339  }
340 
342  static void rcomp(std::string& seq)
343  {
344  if(seq.empty())
345  return;
346 
347  base aux;
348  for(long i=0,j=seq.size()-1;i<=j;++i,--j)
349  {
350  aux=comp(dec(seq[i]));
351  seq[i]=enc(comp(dec(seq[j])));
352  seq[j]=enc(aux);
353  }
354  }
355 
356  template<typename Iter>
357  static void to_string(Iter beg,Iter end,std::string& res)
358  {
359  res.clear();
360  for(Iter i=beg;i!=end;++i)
361  {
362  const typename std::iterator_traits<Iter>::value_type v=*i;
363  res.push_back(enc(static_cast<base>(v)));
364  }
365  }
366 
367  template<typename Iter>
368  static std::string to_string(Iter beg,Iter end)
369  {
370  std::string res;
371  to_string(beg,end,res);
372  return res;
373  }
374 
375  template<typename Iter>
376  static void unpack_to_string(Iter beg,Iter end,std::string& res)
377  {
378  res.clear();
379  for(Iter i=beg;i!=end;++i)
380  {
381  const typename std::iterator_traits<Iter>::value_type v=*i;
382  res.push_back(enc(unpack(static_cast<packed_base>(v))));
383  }
384  }
385 
386  template<typename Iter>
387  static Iter unpack_to_string(
388  Iter beg,Iter end,std::string& res,size_t nmax)
389  {
390  res.clear();
391  Iter i;
392  for(i=beg;(i!=end)&&(nmax>0);++i)
393  {
394  const typename std::iterator_traits<Iter>::value_type v=*i;
395  res.push_back(enc(unpack(static_cast<packed_base>(v))));
396  --nmax;
397  }
398  return i;
399  }
400 
401  template<typename Iter>
402  static std::string unpack_to_string(Iter beg,Iter end)
403  {
404  std::string res;
405  unpack_to_string(beg,end,res);
406  return res;
407  }
408 };
409 
410 inline
411 std::ostream &operator<<(std::ostream &os,const nuc::base b)
412 {
413  os<<nuc::enc(b);
414  return os;
415 }
416 
417 } // namespace
418 
419 #endif // SHORE_GENOME_NUCLEOTIDE_HPP__
420