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() {
        // Test N encoding
        assert_eq!(hamming_dist(&BaseBits::new(b"ACTN").unwrap(), &BaseBits::new(b"ACTG").unwrap()), 1);
        // Test * encoding
        assert_eq!(hamming_dist(&BaseBits::new(b"ACT*").unwrap(), &BaseBits::new(b"ACTG").unwrap()), 0);
        // Test that unkown chars treated like Ns
        assert_eq!(hamming_dist(&BaseBits::new(b"ACT9").unwrap(), &BaseBits::new(b"ACTG").unwrap()), 1);
        // Test regular equality
        assert_eq!(hamming_dist(&BaseBits::new(b"ACTG").unwrap(), &BaseBits::new(b"ACTG").unwrap()), 0);

        // Test Other string
        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);
        // This one should fail since it overflows and wraps around. needs thought
        //assert_eq!(hamming_dist(&BaseBits::new("TATACAGATACTTCCATAGCATC"),
                                //&BaseBits::new("GATACAGATACAACNATAGCATT")), 4);
    }

    #[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);
    }
}

/// Encode a DNA string of up to 21 bases as a u64 for fast hamming distance calculations.
/// TODO: Add a bump to use u128 or maybe bigint if 21 chars is not enough. 
/// TODO: Add equalities and hash function stuff so this type can be used in data structures
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(Copy, Clone)]
#[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
}

/// Extract 'k' bits from the end of a u64 integer
#[inline]
fn extract_bits(n: u64, k: u32) -> u64 {
    !(!0u64 << k) & n
}


// Hamming distance functions that don't depend on BaseBits types
pub mod hamming {
    #[inline]
    pub fn hamming_code(alpha: u64, beta: u64) -> u32 {
        //let x = alpha ^ beta;
        //x.count_ones() / 2
        (alpha ^ beta).count_ones() / 2
    }

    pub fn hamming_str(alpha: &str, beta: &str) -> u32 {
        // skip length check. will default to up to length of alpha
        let mut dist = 0;
        for (a, b) in alpha.chars().zip(beta.chars()) {
            if a != b {
                dist += 1
            }
        }
        dist
    }
}


// dist (pop_count in the paper)
// https://github.com/Daniel-Liu-c0deb0t/UMICollapse/blob/master/src/umicollapse/util/Read.java