source: trunk/yat/omic/DNA.cc @ 2342

Last change on this file since 2342 was 2342, checked in by Peter, 12 years ago

adding missing files. refs #377

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 3.3 KB
Line 
1// $Id: DNA.cc 2342 2010-10-17 03:23:48Z peter $
2
3/*
4  Copyright (C) 2010 Peter Johansson
5
6  This file is part of the yat library, http://dev.thep.lu.se/yat
7
8  The yat library is free software; you can redistribute it and/or
9  modify it under the terms of the GNU General Public License as
10  published by the Free Software Foundation; either version 3 of the
11  License, or (at your option) any later version.
12
13  The yat library 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 GNU
16  General Public License for more details.
17
18  You should have received a copy of the GNU General Public License
19  along with yat. If not, see <http://www.gnu.org/licenses/>.
20*/
21
22#include "DNA.h"
23
24#include <cassert>
25#include <ostream>
26#include <stdexcept>
27
28namespace theplu {
29namespace yat {
30namespace omic {
31
32  std::map<char, unsigned short> DNA::char2code_;
33  std::vector<char> DNA::code2char_;
34
35  DNA::DNA(void)
36    : code_(0) 
37  {
38    assert(code_<16);
39  }
40
41
42  DNA::DNA(char c)
43  {
44    init();
45    std::map<char, unsigned short>::const_iterator i = char2code_.find(c);
46    if (i==char2code_.end()) {
47      std::string msg("DNA('");
48      msg += c;
49      msg += "')";
50      throw std::invalid_argument(msg);
51    }
52    code_ = i->second;
53    assert(code_<16);
54  }
55
56
57  DNA::DNA(const DNA& other)
58    : code_(other.code_)
59  {
60    assert(other.code_<16);
61    assert(code_<16);
62  }
63
64
65  void DNA::init(void) const
66  {
67    if (!char2code_.empty())
68      return;
69
70    char2code_[' '] = 0;
71
72    char2code_['A'] = 1;
73    char2code_['C'] = 2;
74    char2code_['G'] = 4;
75    char2code_['T'] = 8;
76   
77    char2code_['M'] = 1+2;
78    char2code_['R'] = 1+4;
79    char2code_['W'] = 1+8;
80    char2code_['S'] = 2+4;
81    char2code_['Y'] = 2+8;
82    char2code_['K'] = 4+8;
83   
84    char2code_['B'] = 2+4+8;
85    char2code_['D'] = 1+4+8;
86    char2code_['H'] = 1+2+8;
87    char2code_['V'] = 1+2+4;
88   
89    char2code_['N'] = 1+2+4+8;
90
91    code2char_.resize(char2code_.size());
92    for (std::map<char, unsigned short>::const_iterator it=char2code_.begin();
93         it!=char2code_.end(); ++it) {
94      code2char_[it->second] = it->first;
95    }
96  }
97
98
99  DNA DNA::complement(void) const
100  {
101    DNA result;
102    assert(result.code_==0);
103    // FIXME, should probably hard-code this or some magic seq of binary ops
104    if ((*this & DNA('A')) != DNA(' '))
105      result |= DNA('T');
106    if ((*this & DNA('T')) != DNA(' '))
107      result |= DNA('A');
108    if ((*this & DNA('C')) != DNA(' '))
109      result |= DNA('G');
110    if ((*this & DNA('G')) != DNA(' '))
111      result |= DNA('C');
112    assert(result.code_<16);
113    return result;
114  }
115
116
117  std::string expand(const DNA& dna)
118  {
119    std::string result;
120    std::string str="ACGT";
121    for (size_t i=0; i<str.size(); ++i)
122      if ((dna & DNA(str[i])) != DNA(' '))
123        result.push_back(str[i]);
124    return result;
125  }
126
127
128  char DNA::get(void) const
129  {
130    assert(code_ < 16);
131    init();
132    assert(code_ < code2char_.size());
133    return code2char_[code_];
134  }
135   
136
137  bool operator==(const DNA& lhs, const DNA& rhs)
138  {
139    return lhs.code_ == rhs.code_;
140  }
141
142
143  std::ostream& operator<<(std::ostream& os, const DNA& dna)
144  {
145    os.put(dna.get());
146    return os;
147  }
148
149
150  DNA& DNA::operator&=(const DNA& other)
151  {
152    code_ &= other.code_;
153    assert(code_<16);
154    return *this;
155  }
156
157
158  DNA& DNA::operator|=(const DNA& other)
159  {
160    code_ |= other.code_;
161    assert(code_<16);
162    return *this;
163  }
164
165
166  DNA& DNA::operator^=(const DNA& other)
167  {
168    code_ ^= other.code_;
169    assert(code_<16);
170    return *this;
171  }
172}}}
Note: See TracBrowser for help on using the repository browser.