rsomics-seq-grep 0.1.0

/// Case-preserving reverse complement over IUPAC nucleotide codes
/// (`seqkit`'s `DNAredundant` pairing table). Bytes outside the table —
/// protein letters, gaps, anything else — pass through unchanged, matching
/// `seqkit`'s behaviour of leaving unpaired letters as-is.
#[must_use]
pub fn revcomp(seq: &[u8]) -> Vec<u8> {
    seq.iter().rev().map(|&b| pair(b)).collect()
}

fn pair(b: u8) -> u8 {
    match b {
        b'A' => b'T',
        b'T' | b'U' => b'A',
        b'C' => b'G',
        b'G' => b'C',
        b'R' => b'Y',
        b'Y' => b'R',
        b'S' => b'S',
        b'W' => b'W',
        b'K' => b'M',
        b'M' => b'K',
        b'B' => b'V',
        b'V' => b'B',
        b'D' => b'H',
        b'H' => b'D',
        b'N' => b'N',
        b'a' => b't',
        b't' | b'u' => b'a',
        b'c' => b'g',
        b'g' => b'c',
        b'r' => b'y',
        b'y' => b'r',
        b's' => b's',
        b'w' => b'w',
        b'k' => b'm',
        b'm' => b'k',
        b'b' => b'v',
        b'v' => b'b',
        b'd' => b'h',
        b'h' => b'd',
        b'n' => b'n',
        other => other,
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn simple_dna() {
        assert_eq!(revcomp(b"ACGT"), b"ACGT");
        assert_eq!(revcomp(b"TTTTAAAACCCCGGGG"), b"CCCCGGGGTTTTAAAA");
    }

    #[test]
    fn case_preserved() {
        assert_eq!(revcomp(b"acgt"), b"acgt");
        assert_eq!(revcomp(b"AcGt"), b"aCgT");
    }

    #[test]
    fn ambiguity_codes() {
        assert_eq!(revcomp(b"RYSWKMBVDHN"), b"NDHBVKMWSRY");
    }
}