1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
#[cfg(test)]
mod tests {
use super::hamming::*;
use super::*;
#[test]
fn test_hamming_str_dist() {
assert_eq!(hamming_str("ACTG", "ACTT"), 1);
assert_eq!(hamming_str("ACTG", "ACTTT"), 1);
}
#[test]
fn test_base_bits() {
let alpha = BaseBits::new(b"ACTG").unwrap();
let beta = BaseBits::new(b"ACTT").unwrap();
assert_eq!(hamming_dist(&alpha, &beta), 1);
}
#[test]
fn test_cases_bb_hamming() {
assert_eq!(hamming_dist(&BaseBits::new(b"ACTN").unwrap(), &BaseBits::new(b"ACTG").unwrap()), 1);
assert_eq!(hamming_dist(&BaseBits::new(b"ACT*").unwrap(), &BaseBits::new(b"ACTG").unwrap()), 0);
assert_eq!(hamming_dist(&BaseBits::new(b"ACT9").unwrap(), &BaseBits::new(b"ACTG").unwrap()), 1);
assert_eq!(hamming_dist(&BaseBits::new(b"ACTG").unwrap(), &BaseBits::new(b"ACTG").unwrap()), 0);
assert_eq!(hamming_dist(&BaseBits::new(b"GATACA").unwrap(), &BaseBits::new(b"GATACT").unwrap()), 1);
assert_eq!(hamming_dist(&BaseBits::new(b"CATACAGATACTTCCATAGCT").unwrap(),
&BaseBits::new(b"GATACAGATACTTCCATAGCA").unwrap()), 2);
}
#[test]
fn test_bb_to_string() {
let alpha = BaseBits::new(b"GCTAN").unwrap();
let beta = BaseBits::new(b"ACTG*").unwrap();
assert_eq!(alpha.to_string(), "GCTAN".to_string());
assert_eq!(beta.to_string(), "ACTG*".to_string());
let long = BaseBits::new(b"GATACAGATACAACNATAGCA").unwrap();
assert_eq!(long.to_string(), "GATACAGATACAACNATAGCA".to_string());
}
#[test]
fn test_encoding() {
let bb = BaseBits::new(b"ACTG").unwrap();
assert_eq!(bb.code, 0b000110101011);
}
}
use std::fmt;
use std::str;
use std::u64;
pub const ENCODING_DIST: u32 = 2;
pub const ENCODING_LENGTH: u32 = 3;
pub const CONTAINER_WIDTH: u32 = 64;
pub const MAX_BASES: usize = (CONTAINER_WIDTH / ENCODING_LENGTH) as usize;
pub const UNDETERMINED: u64 = 0b100;
pub const ANY: u64 = 0b111;
pub const MAX_VAL: u64 = u64::MAX;
struct Bases;
impl Bases {
const A: u64 = 0b000;
const C: u64 = 0b110;
const T: u64 = 0b101;
const G: u64 = 0b011;
const N: u64 = UNDETERMINED;
const STAR: u64 = ANY;
}
#[derive(Hash, PartialEq, Eq, Debug, Copy, Clone)]
pub struct BaseBits {
pub code: u64,
len: usize,
}
impl BaseBits {
pub fn new(seq: &[u8])-> Result<BaseBits, &'static str> {
let mut code: u64 = 0;
let len = seq.len();
if len > MAX_BASES {
return Err("Length of string to encode exceeds MAX_BASES");
}
for c in seq.iter() {
code = (code << ENCODING_LENGTH) | match c {
b'A' => Bases::A,
b'C' => Bases::C,
b'T' => Bases::T,
b'G' => Bases::G,
b'N' => Bases::N,
b'*' => Bases::STAR,
_ => Bases::N,
}
}
Ok(BaseBits{code, len})
}
pub fn decode(&self) -> Vec<u8> {
let mut s = Vec::new();
let mut code = self.code;
for _ in 0..self.len {
let base = extract_bits(code, ENCODING_LENGTH);
code = code >> ENCODING_LENGTH;
s.push(match base {
Bases::A => b'A',
Bases::C => b'C',
Bases::T => b'T',
Bases::G => b'G',
Bases::N => b'N',
Bases::STAR => b'*',
_ => b'N'
});
}
s.into_iter().rev().collect()
}
}
impl fmt::Display for BaseBits {
fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
write!(f, "{}", str::from_utf8(&self.decode()).unwrap())
}
}
#[inline]
pub fn hamming_dist(alpha: &BaseBits, beta: &BaseBits) -> u32 {
(alpha.code ^ beta.code).count_ones() / 2
}
#[inline]
fn extract_bits(n: u64, k: u32) -> u64 {
!(!0u64 << k) & n
}
pub mod hamming {
#[inline]
pub fn hamming_code(alpha: u64, beta: u64) -> u32 {
(alpha ^ beta).count_ones() / 2
}
pub fn hamming_str(alpha: &str, beta: &str) -> u32 {
let mut dist = 0;
for (a, b) in alpha.chars().zip(beta.chars()) {
if a != b {
dist += 1
}
}
dist
}
}