argenus/
snp.rs

1//! SNP Verification Module
2//!
3//! Verifies point mutations in SNP-based ARG genes after targeted assembly.
4//! Parses SNP information from gene names and checks if the mutation exists
5//! in the assembled contig sequence.
6//!
7//! # SNP Gene Name Format
8//! SNP-based ARG names follow the pattern: `gene_WtPosVar`
9//! - `gene`: Base gene name (e.g., "rpoB", "gyrA")
10//! - `Wt`: Wild-type amino acid (single letter)
11//! - `Pos`: Position in the protein sequence
12//! - `Var`: Variant (mutant) amino acid
13//!
14//! Examples: `rpoB_V146F`, `gyrA_S83L`, `parC_S80I`
15
16use rustc_hash::FxHashMap;
17use std::sync::LazyLock;
18
19// ============================================================================
20// SNP Status
21// ============================================================================
22
23/// Result of SNP verification.
24#[derive(Debug, Clone, PartialEq)]
25pub enum SnpStatus {
26    /// Not a SNP-based gene (acquired gene)
27    NotApplicable,
28    /// SNP confirmed - mutant allele present
29    Confirmed,
30    /// Wild-type allele found - no resistance mutation
31    WildType,
32    /// Different amino acid found at position
33    NovelVariant(char),
34    /// SNP position not covered by contig
35    NotCovered,
36    /// Could not verify (alignment issue, etc.)
37    Unverified(String),
38}
39
40impl std::fmt::Display for SnpStatus {
41    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
42        match self {
43            SnpStatus::NotApplicable => write!(f, "Acquired"),
44            SnpStatus::Confirmed => write!(f, "Confirmed"),
45            SnpStatus::WildType => write!(f, "WildType"),
46            SnpStatus::NovelVariant(aa) => write!(f, "Novel({})", aa),
47            SnpStatus::NotCovered => write!(f, "NotCovered"),
48            SnpStatus::Unverified(reason) => write!(f, "Unverified({})", reason),
49        }
50    }
51}
52
53// ============================================================================
54// SNP Info
55// ============================================================================
56
57/// Parsed SNP information from gene name.
58#[derive(Debug, Clone)]
59pub struct SnpInfo {
60    /// Base gene name (e.g., "rpoB")
61    #[allow(dead_code)]
62    pub gene: String,
63    /// Wild-type amino acid
64    pub wildtype: char,
65    /// Position in protein sequence (1-based)
66    pub position: usize,
67    /// Mutant amino acid
68    pub mutant: char,
69}
70
71impl SnpInfo {
72    /// Calculates the nucleotide position range for this SNP (0-based).
73    /// Returns (start, end) for the codon.
74    pub fn nucleotide_range(&self) -> (usize, usize) {
75        let start = (self.position - 1) * 3;
76        let end = start + 3;
77        (start, end)
78    }
79}
80
81// ============================================================================
82// Codon Table
83// ============================================================================
84
85/// Standard genetic code codon table.
86static CODON_TABLE: LazyLock<FxHashMap<&'static str, char>> = LazyLock::new(|| {
87    let mut table = FxHashMap::default();
88    // Phenylalanine (F)
89    table.insert("TTT", 'F'); table.insert("TTC", 'F');
90    // Leucine (L)
91    table.insert("TTA", 'L'); table.insert("TTG", 'L');
92    table.insert("CTT", 'L'); table.insert("CTC", 'L');
93    table.insert("CTA", 'L'); table.insert("CTG", 'L');
94    // Isoleucine (I)
95    table.insert("ATT", 'I'); table.insert("ATC", 'I'); table.insert("ATA", 'I');
96    // Methionine (M) - Start
97    table.insert("ATG", 'M');
98    // Valine (V)
99    table.insert("GTT", 'V'); table.insert("GTC", 'V');
100    table.insert("GTA", 'V'); table.insert("GTG", 'V');
101    // Serine (S)
102    table.insert("TCT", 'S'); table.insert("TCC", 'S');
103    table.insert("TCA", 'S'); table.insert("TCG", 'S');
104    table.insert("AGT", 'S'); table.insert("AGC", 'S');
105    // Proline (P)
106    table.insert("CCT", 'P'); table.insert("CCC", 'P');
107    table.insert("CCA", 'P'); table.insert("CCG", 'P');
108    // Threonine (T)
109    table.insert("ACT", 'T'); table.insert("ACC", 'T');
110    table.insert("ACA", 'T'); table.insert("ACG", 'T');
111    // Alanine (A)
112    table.insert("GCT", 'A'); table.insert("GCC", 'A');
113    table.insert("GCA", 'A'); table.insert("GCG", 'A');
114    // Tyrosine (Y)
115    table.insert("TAT", 'Y'); table.insert("TAC", 'Y');
116    // Stop codons (*)
117    table.insert("TAA", '*'); table.insert("TAG", '*'); table.insert("TGA", '*');
118    // Histidine (H)
119    table.insert("CAT", 'H'); table.insert("CAC", 'H');
120    // Glutamine (Q)
121    table.insert("CAA", 'Q'); table.insert("CAG", 'Q');
122    // Asparagine (N)
123    table.insert("AAT", 'N'); table.insert("AAC", 'N');
124    // Lysine (K)
125    table.insert("AAA", 'K'); table.insert("AAG", 'K');
126    // Aspartic acid (D)
127    table.insert("GAT", 'D'); table.insert("GAC", 'D');
128    // Glutamic acid (E)
129    table.insert("GAA", 'E'); table.insert("GAG", 'E');
130    // Cysteine (C)
131    table.insert("TGT", 'C'); table.insert("TGC", 'C');
132    // Tryptophan (W)
133    table.insert("TGG", 'W');
134    // Arginine (R)
135    table.insert("CGT", 'R'); table.insert("CGC", 'R');
136    table.insert("CGA", 'R'); table.insert("CGG", 'R');
137    table.insert("AGA", 'R'); table.insert("AGG", 'R');
138    // Glycine (G)
139    table.insert("GGT", 'G'); table.insert("GGC", 'G');
140    table.insert("GGA", 'G'); table.insert("GGG", 'G');
141    table
142});
143
144/// Translates a codon (3 nucleotides) to amino acid.
145pub fn translate_codon(codon: &str) -> Option<char> {
146    let upper = codon.to_uppercase();
147    CODON_TABLE.get(upper.as_str()).copied()
148}
149
150// ============================================================================
151// SNP Parsing
152// ============================================================================
153
154/// Checks if a gene name represents a SNP-based ARG.
155///
156/// SNP genes have patterns like: gene_X###Y where X and Y are amino acids
157/// and ### is a number.
158#[allow(dead_code)]
159pub fn is_snp_gene(gene_name: &str) -> bool {
160    parse_snp_info(gene_name).is_some()
161}
162
163/// Parses SNP information from a gene name.
164///
165/// # Examples
166/// ```
167/// use argenus::snp::parse_snp_info;
168///
169/// let info = parse_snp_info("rpoB_V146F").unwrap();
170/// assert_eq!(info.gene, "rpoB");
171/// assert_eq!(info.wildtype, 'V');
172/// assert_eq!(info.position, 146);
173/// assert_eq!(info.mutant, 'F');
174/// ```
175pub fn parse_snp_info(gene_name: &str) -> Option<SnpInfo> {
176    // Handle full AMR database format: gene_X###Y|CLASS|SUBCLASS|...
177    let name = gene_name.split('|').next().unwrap_or(gene_name);
178
179    // Find the last underscore that precedes a SNP pattern
180    let parts: Vec<&str> = name.rsplitn(2, '_').collect();
181    if parts.len() != 2 {
182        return None;
183    }
184
185    let snp_part = parts[0];  // e.g., "V146F"
186    let gene = parts[1];      // e.g., "rpoB"
187
188    // SNP pattern: single letter + digits + single letter
189    // Minimum length: 3 (e.g., "A1B")
190    if snp_part.len() < 3 {
191        return None;
192    }
193
194    let chars: Vec<char> = snp_part.chars().collect();
195
196    // First char must be amino acid letter
197    let wildtype = chars[0].to_ascii_uppercase();
198    if !is_amino_acid(wildtype) {
199        return None;
200    }
201
202    // Last char must be amino acid letter
203    let mutant = chars[chars.len() - 1].to_ascii_uppercase();
204    if !is_amino_acid(mutant) {
205        return None;
206    }
207
208    // Middle part must be digits (position)
209    let pos_str: String = chars[1..chars.len()-1].iter().collect();
210    let position: usize = pos_str.parse().ok()?;
211
212    if position == 0 {
213        return None;
214    }
215
216    Some(SnpInfo {
217        gene: gene.to_string(),
218        wildtype,
219        position,
220        mutant,
221    })
222}
223
224/// Checks if a character is a valid amino acid single-letter code.
225fn is_amino_acid(c: char) -> bool {
226    matches!(c, 'A' | 'C' | 'D' | 'E' | 'F' | 'G' | 'H' | 'I' | 'K' | 'L' |
227                'M' | 'N' | 'P' | 'Q' | 'R' | 'S' | 'T' | 'V' | 'W' | 'Y')
228}
229
230// ============================================================================
231// SNP Verification
232// ============================================================================
233
234/// Verifies SNP presence in a contig aligned to an ARG reference.
235///
236/// # Arguments
237/// * `contig_seq` - The assembled contig sequence
238/// * `gene_name` - Full ARG gene name (e.g., "rpoB_V146F|RIFAMYCIN|...")
239/// * `ref_start` - Reference (ARG) start position in alignment (0-based)
240/// * `ref_end` - Reference (ARG) end position in alignment
241/// * `contig_start` - Contig start position in alignment (0-based)
242/// * `contig_end` - Contig end position in alignment
243/// * `strand` - Alignment strand ('+' or '-')
244///
245/// # Returns
246/// SNP verification status
247pub fn verify_snp(
248    contig_seq: &str,
249    gene_name: &str,
250    ref_start: usize,
251    ref_end: usize,
252    contig_start: usize,
253    contig_end: usize,
254    strand: char,
255) -> SnpStatus {
256    // Parse SNP info from gene name
257    let snp_info = match parse_snp_info(gene_name) {
258        Some(info) => info,
259        None => return SnpStatus::NotApplicable,
260    };
261
262    // Calculate nucleotide position of SNP in reference
263    let (snp_nt_start, snp_nt_end) = snp_info.nucleotide_range();
264
265    // Check if SNP position is covered by alignment
266    if snp_nt_start < ref_start || snp_nt_end > ref_end {
267        return SnpStatus::NotCovered;
268    }
269
270    // Calculate corresponding position in contig
271    let offset_in_ref = snp_nt_start - ref_start;
272    let contig_snp_start = if strand == '+' {
273        contig_start + offset_in_ref
274    } else {
275        // For reverse strand, positions are mirrored
276        contig_end - offset_in_ref - 3
277    };
278    let contig_snp_end = contig_snp_start + 3;
279
280    // Check bounds
281    if contig_snp_end > contig_seq.len() {
282        return SnpStatus::NotCovered;
283    }
284
285    // Extract codon from contig
286    let codon_seq = &contig_seq[contig_snp_start..contig_snp_end];
287
288    // Handle reverse strand - need reverse complement
289    let codon = if strand == '-' {
290        reverse_complement(codon_seq)
291    } else {
292        codon_seq.to_uppercase()
293    };
294
295    // Translate codon
296    let amino_acid = match translate_codon(&codon) {
297        Some(aa) => aa,
298        None => return SnpStatus::Unverified(format!("invalid_codon:{}", codon)),
299    };
300
301    // Compare with expected
302    if amino_acid == snp_info.mutant {
303        SnpStatus::Confirmed
304    } else if amino_acid == snp_info.wildtype {
305        SnpStatus::WildType
306    } else {
307        SnpStatus::NovelVariant(amino_acid)
308    }
309}
310
311/// Computes the reverse complement of a DNA sequence.
312fn reverse_complement(seq: &str) -> String {
313    seq.chars()
314        .rev()
315        .map(|c| match c.to_ascii_uppercase() {
316            'A' => 'T',
317            'T' => 'A',
318            'G' => 'C',
319            'C' => 'G',
320            _ => 'N',
321        })
322        .collect()
323}
324
325// ============================================================================
326// Tests
327// ============================================================================
328
329#[cfg(test)]
330mod tests {
331    use super::*;
332
333    #[test]
334    fn test_parse_snp_info() {
335        // Basic SNP pattern
336        let info = parse_snp_info("rpoB_V146F").unwrap();
337        assert_eq!(info.gene, "rpoB");
338        assert_eq!(info.wildtype, 'V');
339        assert_eq!(info.position, 146);
340        assert_eq!(info.mutant, 'F');
341
342        // With full AMR format
343        let info = parse_snp_info("gyrA_S83L|QUINOLONE|QUINOLONE|J01M").unwrap();
344        assert_eq!(info.gene, "gyrA");
345        assert_eq!(info.wildtype, 'S');
346        assert_eq!(info.position, 83);
347        assert_eq!(info.mutant, 'L');
348
349        // Single digit position
350        let info = parse_snp_info("test_A1G").unwrap();
351        assert_eq!(info.position, 1);
352
353        // Non-SNP gene
354        assert!(parse_snp_info("blaTEM-1").is_none());
355        assert!(parse_snp_info("tet(M)").is_none());
356    }
357
358    #[test]
359    fn test_is_snp_gene() {
360        assert!(is_snp_gene("rpoB_V146F"));
361        assert!(is_snp_gene("gyrA_S83L|QUINOLONE"));
362        assert!(!is_snp_gene("blaTEM-1"));
363        assert!(!is_snp_gene("tet(M)"));
364    }
365
366    #[test]
367    fn test_translate_codon() {
368        assert_eq!(translate_codon("ATG"), Some('M'));
369        assert_eq!(translate_codon("TTT"), Some('F'));
370        assert_eq!(translate_codon("TTC"), Some('F'));
371        assert_eq!(translate_codon("GTT"), Some('V'));
372        assert_eq!(translate_codon("TAA"), Some('*'));
373        assert_eq!(translate_codon("XXX"), None);
374    }
375
376    #[test]
377    fn test_nucleotide_range() {
378        let info = SnpInfo {
379            gene: "test".to_string(),
380            wildtype: 'V',
381            position: 146,
382            mutant: 'F',
383        };
384        // Position 146 (1-based) = nucleotides 435-437 (0-based)
385        assert_eq!(info.nucleotide_range(), (435, 438));
386
387        let info2 = SnpInfo {
388            gene: "test".to_string(),
389            wildtype: 'A',
390            position: 1,
391            mutant: 'B',
392        };
393        assert_eq!(info2.nucleotide_range(), (0, 3));
394    }
395
396    #[test]
397    fn test_verify_snp() {
398        // Create a mock contig with known codon at position
399        // Position 2 (aa) = nucleotides 3-5 (0-based)
400        // TTT = F (Phe), GTT = V (Val)
401
402        // Test confirmed (mutant found)
403        let contig = "ATGTTTGGG"; // Met-Phe-Gly
404        let status = verify_snp(
405            contig,
406            "test_V2F",  // V at position 2 mutated to F
407            0, 9,        // ref covers full sequence
408            0, 9,        // contig covers full sequence
409            '+',
410        );
411        assert_eq!(status, SnpStatus::Confirmed);
412
413        // Test wildtype found
414        let contig = "ATGGTTGGG"; // Met-Val-Gly (wildtype V at position 2)
415        let status = verify_snp(
416            contig,
417            "test_V2F",
418            0, 9,
419            0, 9,
420            '+',
421        );
422        assert_eq!(status, SnpStatus::WildType);
423
424        // Test novel variant
425        let contig = "ATGCTTGGG"; // Met-Leu-Gly (neither V nor F)
426        let status = verify_snp(
427            contig,
428            "test_V2F",
429            0, 9,
430            0, 9,
431            '+',
432        );
433        assert_eq!(status, SnpStatus::NovelVariant('L'));
434    }
435
436    #[test]
437    fn test_snp_status_display() {
438        assert_eq!(format!("{}", SnpStatus::Confirmed), "Confirmed");
439        assert_eq!(format!("{}", SnpStatus::WildType), "WildType");
440        assert_eq!(format!("{}", SnpStatus::NovelVariant('L')), "Novel(L)");
441        assert_eq!(format!("{}", SnpStatus::NotApplicable), "Acquired");
442    }
443}
argenus/snp.rs

argenus/
snp.rs
