dnagent_domain/
translation.rs

1//! Translation of DNA to protein and open-reading-frame discovery.
2//!
3//! Genetic codes are the pinned NCBI tables in [`crate::genetic_codes`]. Translation is
4//! IUPAC-aware: a codon translates to an amino acid only when every expansion of its
5//! ambiguity codes agrees (otherwise `X`); a stop requires every expansion to be a stop;
6//! starts must be exact. Coordinates are zero-based reference positions on the forward
7//! strand. Multipart features are spliced in source order; reverse-strand features are
8//! the reverse complement of that splice (GenBank `complement(join(...))`). Nothing here
9//! infers genes, splicing, frames beyond the declared inputs, or experimental validity.
10
11use crate::genetic_codes::{GeneticCodeData, TABLES};
12use crate::{DnaSeq, Feature, Strand, Topology};
13use serde::Serialize;
14use thiserror::Error;
15
16pub use crate::genetic_codes::{GC_PRT_SHA256, GC_PRT_VERSION};
17
18/// Errors raised by translation requests.
19#[derive(Debug, Clone, PartialEq, Eq, Error)]
20pub enum TranslationError {
21    #[error("unknown NCBI genetic code {0}; valid tables are listed by `dnagent translate --help`")]
22    UnknownTable(u32),
23    #[error("feature {0:?} has unknown strand; translation needs a declared strand")]
24    UnknownStrand(String),
25    #[error("codon_start must be 1, 2 or 3, got {0:?}")]
26    InvalidCodonStart(String),
27    #[error("transl_table qualifier {0:?} is not an NCBI table number")]
28    InvalidTableQualifier(String),
29    #[error("frame offset must be 0, 1 or 2, got {0}")]
30    InvalidFrame(usize),
31    #[error("invalid range {start}..{end} for a {topology:?} molecule of length {length}")]
32    InvalidRange {
33        start: usize,
34        end: usize,
35        length: usize,
36        topology: Topology,
37    },
38    #[error("nothing to translate: fewer than three coding bases")]
39    TooShort,
40    #[error("minimum ORF length must be at least 1 codon")]
41    InvalidMinimum,
42}
43
44/// Which codons may open an ORF.
45#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize)]
46#[serde(rename_all = "snake_case")]
47pub enum StartPolicy {
48    /// Only ATG.
49    AtgOnly,
50    /// Every start codon of the selected table (e.g. GTG/TTG in table 11).
51    TableStarts,
52}
53
54/// Coding direction relative to the forward reference.
55#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize)]
56#[serde(rename_all = "snake_case")]
57pub enum CodingStrand {
58    Forward,
59    Reverse,
60}
61
62/// A pinned NCBI translation table.
63#[derive(Debug, Clone, Copy)]
64pub struct GeneticCode(&'static GeneticCodeData);
65
66fn base_index(base: u8) -> usize {
67    match base {
68        b'T' => 0,
69        b'C' => 1,
70        b'A' => 2,
71        b'G' => 3,
72        _ => unreachable!("expansions are ACGT"),
73    }
74}
75
76fn expansions(base: u8) -> &'static [u8] {
77    match base {
78        b'A' => b"A",
79        b'C' => b"C",
80        b'G' => b"G",
81        b'T' => b"T",
82        b'R' => b"AG",
83        b'Y' => b"CT",
84        b'S' => b"CG",
85        b'W' => b"AT",
86        b'K' => b"GT",
87        b'M' => b"AC",
88        b'B' => b"CGT",
89        b'D' => b"AGT",
90        b'H' => b"ACT",
91        b'V' => b"ACG",
92        b'N' => b"ACGT",
93        _ => unreachable!("DnaSeq validates IUPAC DNA"),
94    }
95}
96
97fn complement(base: u8) -> u8 {
98    match base {
99        b'A' => b'T',
100        b'T' => b'A',
101        b'C' => b'G',
102        b'G' => b'C',
103        b'R' => b'Y',
104        b'Y' => b'R',
105        b'S' => b'S',
106        b'W' => b'W',
107        b'K' => b'M',
108        b'M' => b'K',
109        b'B' => b'V',
110        b'V' => b'B',
111        b'D' => b'H',
112        b'H' => b'D',
113        b'N' => b'N',
114        _ => unreachable!("DnaSeq validates IUPAC DNA"),
115    }
116}
117
118impl GeneticCode {
119    /// Look up an NCBI table by number.
120    pub fn ncbi(id: u32) -> Result<Self, TranslationError> {
121        TABLES
122            .iter()
123            .find(|table| u32::from(table.id) == id)
124            .map(Self)
125            .ok_or(TranslationError::UnknownTable(id))
126    }
127
128    /// NCBI table 1.
129    #[must_use]
130    pub fn standard() -> Self {
131        Self(&TABLES[0])
132    }
133
134    #[must_use]
135    pub fn id(self) -> u8 {
136        self.0.id
137    }
138
139    #[must_use]
140    pub fn name(self) -> &'static str {
141        self.0.name
142    }
143
144    /// Every available table as `(id, name)`.
145    pub fn all() -> impl Iterator<Item = (u8, &'static str)> {
146        TABLES.iter().map(|table| (table.id, table.name))
147    }
148
149    fn each_expansion(codon: [u8; 3], mut visit: impl FnMut(usize)) {
150        for &a in expansions(codon[0]) {
151            for &b in expansions(codon[1]) {
152                for &c in expansions(codon[2]) {
153                    visit(16 * base_index(a) + 4 * base_index(b) + base_index(c));
154                }
155            }
156        }
157    }
158
159    /// Amino acid for an IUPAC codon: agreed by every expansion, else `X`. Stops are `*`.
160    #[must_use]
161    pub fn translate_codon(self, codon: [u8; 3]) -> char {
162        let mut agreed: Option<u8> = None;
163        let mut conflict = false;
164        Self::each_expansion(codon, |index| {
165            let amino = self.0.amino_acids[index];
166            match agreed {
167                None => agreed = Some(amino),
168                Some(previous) if previous != amino => conflict = true,
169                Some(_) => {}
170            }
171        });
172        if conflict {
173            'X'
174        } else {
175            char::from(agreed.expect("at least one expansion"))
176        }
177    }
178
179    /// True when every expansion is a stop codon.
180    #[must_use]
181    pub fn is_stop(self, codon: [u8; 3]) -> bool {
182        self.translate_codon(codon) == '*'
183    }
184
185    /// True for an exact (unambiguous) start codon under `policy`.
186    #[must_use]
187    pub fn is_start(self, codon: [u8; 3], policy: StartPolicy) -> bool {
188        if !codon.iter().all(|base| b"ACGT".contains(base)) {
189            return false;
190        }
191        match policy {
192            StartPolicy::AtgOnly => &codon == b"ATG",
193            StartPolicy::TableStarts => {
194                let index =
195                    16 * base_index(codon[0]) + 4 * base_index(codon[1]) + base_index(codon[2]);
196                self.0.starts[index] == b'M'
197            }
198        }
199    }
200}
201
202/// One translated codon with the reference positions of its three coding bases
203/// (in coding order, so they decrease on the reverse strand and may jump at joins).
204#[derive(Debug, Clone, PartialEq, Eq, Serialize)]
205pub struct Codon {
206    pub amino_acid: char,
207    /// Codon on the coding strand, 5′→3′.
208    pub codon: String,
209    pub positions: [usize; 3],
210}
211
212/// Translation of an explicit coding sequence.
213#[derive(Debug, Clone, PartialEq, Eq, Serialize)]
214pub struct Translation {
215    pub table: u8,
216    pub table_name: &'static str,
217    pub strand: CodingStrand,
218    /// 1-based offset into the coding sequence of the first translated base.
219    pub codon_start: u8,
220    /// One letter per codon; stops are `*` wherever they occur.
221    pub protein: String,
222    pub codons: Vec<Codon>,
223    /// Codon indices of stops before the final codon.
224    pub internal_stops: Vec<usize>,
225    pub terminal_stop: bool,
226    /// Coding bases left over after the last complete codon (0–2).
227    pub trailing_bases: usize,
228    /// Codons translated as `X` because their ambiguity codes disagree.
229    pub ambiguous_codons: usize,
230    pub starts_with_atg: bool,
231    /// A complete feature CDS (codon_start 1, table start codon, terminal stop) whose
232    /// first codon is an alternative start is translated with an initiator M
233    /// (NCBI/Biopython `cds=True` convention). Never applied to ranges.
234    pub initiator_as_methionine: bool,
235}
236
237type Coding = Vec<(u8, usize)>;
238
239fn translate_coding(
240    coding: &[(u8, usize)],
241    code: GeneticCode,
242    strand: CodingStrand,
243    codon_start: u8,
244) -> Result<Translation, TranslationError> {
245    let offset = usize::from(codon_start - 1);
246    let body = coding.get(offset..).unwrap_or_default();
247    if body.len() < 3 {
248        return Err(TranslationError::TooShort);
249    }
250    let codons: Vec<Codon> = body
251        .as_chunks::<3>()
252        .0
253        .iter()
254        .map(|chunk| {
255            let bases = [chunk[0].0, chunk[1].0, chunk[2].0];
256            Codon {
257                amino_acid: code.translate_codon(bases),
258                codon: String::from_utf8(bases.to_vec()).expect("ASCII bases"),
259                positions: [chunk[0].1, chunk[1].1, chunk[2].1],
260            }
261        })
262        .collect();
263    let protein: String = codons.iter().map(|codon| codon.amino_acid).collect();
264    let last = codons.len() - 1;
265    Ok(Translation {
266        table: code.id(),
267        table_name: code.name(),
268        strand,
269        codon_start,
270        internal_stops: codons
271            .iter()
272            .enumerate()
273            .filter(|(i, c)| *i < last && c.amino_acid == '*')
274            .map(|(i, _)| i)
275            .collect(),
276        terminal_stop: codons[last].amino_acid == '*',
277        trailing_bases: body.len() % 3,
278        ambiguous_codons: codons
279            .iter()
280            .filter(|codon| codon.amino_acid == 'X')
281            .count(),
282        starts_with_atg: codons[0].codon == "ATG",
283        initiator_as_methionine: false,
284        protein,
285        codons,
286    })
287}
288
289fn reverse_complement_coding(mut coding: Coding) -> Coding {
290    coding.reverse();
291    for (base, _) in &mut coding {
292        *base = complement(*base);
293    }
294    coding
295}
296
297/// Bases of `[start, start + length)` with reference positions, wrapping on circular molecules.
298fn forward_bases(sequence: &DnaSeq, start: usize, length: usize) -> Coding {
299    let bytes = sequence.as_str().as_bytes();
300    let total = bytes.len();
301    (0..length)
302        .map(|offset| {
303            let position = (start + offset) % total;
304            (bytes[position], position)
305        })
306        .collect()
307}
308
309/// Parsed `codon_start` (default 1) and `transl_table` (default: none) qualifiers.
310pub fn feature_translation_qualifiers(
311    feature: &Feature,
312) -> Result<(u8, Option<u32>), TranslationError> {
313    let mut codon_start = 1;
314    let mut table = None;
315    for qualifier in feature.qualifiers() {
316        let value = qualifier.value.as_deref().unwrap_or("").trim();
317        match qualifier.key.as_str() {
318            "codon_start" => {
319                codon_start = match value {
320                    "1" => 1,
321                    "2" => 2,
322                    "3" => 3,
323                    _ => return Err(TranslationError::InvalidCodonStart(value.to_owned())),
324                }
325            }
326            "transl_table" => {
327                table = Some(
328                    value
329                        .parse::<u32>()
330                        .map_err(|_| TranslationError::InvalidTableQualifier(value.to_owned()))?,
331                );
332            }
333            _ => {}
334        }
335    }
336    Ok((codon_start, table))
337}
338
339/// Translate a feature's spliced coding sequence. `table` overrides `transl_table`.
340pub fn translate_feature(
341    sequence: &DnaSeq,
342    feature: &Feature,
343    table: Option<u32>,
344) -> Result<Translation, TranslationError> {
345    let (codon_start, qualifier_table) = feature_translation_qualifiers(feature)?;
346    let code = GeneticCode::ncbi(table.or(qualifier_table).unwrap_or(1))?;
347    let mut coding: Coding = feature
348        .location()
349        .parts()
350        .iter()
351        .flat_map(|part| forward_bases(sequence, part.start().get(), part.length().get()))
352        .collect();
353    let strand = match feature.location().strand() {
354        Strand::Forward => CodingStrand::Forward,
355        Strand::Reverse => {
356            coding = reverse_complement_coding(coding);
357            CodingStrand::Reverse
358        }
359        Strand::Unknown => {
360            return Err(TranslationError::UnknownStrand(
361                feature.id().as_str().to_owned(),
362            ));
363        }
364    };
365    let mut translation = translate_coding(&coding, code, strand, codon_start)?;
366    let first = &translation.codons[0];
367    let bases: [u8; 3] = first.codon.as_bytes().try_into().expect("three bases");
368    if codon_start == 1
369        && translation.terminal_stop
370        && first.amino_acid != 'M'
371        && code.is_start(bases, StartPolicy::TableStarts)
372    {
373        translation.codons[0].amino_acid = 'M';
374        translation.protein.replace_range(0..1, "M");
375        translation.initiator_as_methionine = true;
376    }
377    Ok(translation)
378}
379
380/// Translate the forward-coordinate range `[start, end)` (on circular molecules `end <
381/// start` wraps through the origin) on either strand, skipping `frame` (0–2) bases
382/// from the 5′ end of the coding strand.
383pub fn translate_range(
384    sequence: &DnaSeq,
385    topology: Topology,
386    start: usize,
387    end: usize,
388    strand: CodingStrand,
389    frame: usize,
390    table: u32,
391) -> Result<Translation, TranslationError> {
392    let length = sequence.len();
393    let invalid = || TranslationError::InvalidRange {
394        start,
395        end,
396        length,
397        topology,
398    };
399    if frame > 2 {
400        return Err(TranslationError::InvalidFrame(frame));
401    }
402    if start >= length || end > length || end == start {
403        return Err(invalid());
404    }
405    let span = if end > start {
406        end - start
407    } else if topology == Topology::Circular {
408        length - start + end
409    } else {
410        return Err(invalid());
411    };
412    let code = GeneticCode::ncbi(table)?;
413    let mut coding = forward_bases(sequence, start, span);
414    if strand == CodingStrand::Reverse {
415        coding = reverse_complement_coding(coding);
416    }
417    translate_coding(
418        &coding,
419        code,
420        strand,
421        u8::try_from(frame + 1).expect("frame checked"),
422    )
423}
424
425/// Imported `translation` qualifier with SnapGene's per-part commas and whitespace removed.
426#[must_use]
427pub fn imported_translation(feature: &Feature) -> Option<String> {
428    feature
429        .qualifiers()
430        .iter()
431        .find(|qualifier| qualifier.key == "translation")
432        .and_then(|qualifier| qualifier.value.as_deref())
433        .map(|value| {
434            value
435                .chars()
436                .filter(|c| !c.is_whitespace() && *c != ',')
437                .collect()
438        })
439}
440
441/// Comparison between the computed protein and an imported translation.
442#[derive(Debug, Clone, PartialEq, Eq, Serialize)]
443pub struct ImportedComparison {
444    pub imported_length: usize,
445    pub matches: bool,
446    /// First differing residue index (0-based), if any.
447    pub first_difference: Option<usize>,
448}
449
450/// Compare with an imported translation, ignoring one terminal stop on either side.
451#[must_use]
452pub fn compare_imported(translation: &Translation, imported: &str) -> ImportedComparison {
453    let computed = translation
454        .protein
455        .strip_suffix('*')
456        .unwrap_or(&translation.protein);
457    let imported = imported.strip_suffix('*').unwrap_or(imported);
458    let first_difference = computed
459        .chars()
460        .zip(imported.chars())
461        .position(|(a, b)| a != b)
462        .or_else(|| (computed.len() != imported.len()).then(|| computed.len().min(imported.len())));
463    ImportedComparison {
464        imported_length: imported.chars().count(),
465        matches: first_difference.is_none(),
466        first_difference,
467    }
468}
469
470/// One complete open reading frame (start codon through stop codon).
471#[derive(Debug, Clone, PartialEq, Eq, Serialize)]
472pub struct Orf {
473    pub id: String,
474    pub strand: CodingStrand,
475    /// Lowest forward-reference position covered; with `length` it may wrap the origin.
476    pub start: usize,
477    /// Bases including the stop codon.
478    pub length: usize,
479    /// Offset of the 5′ codon base from the 5′ end of its strand, modulo 3
480    /// (`start % 3` forward), matching [`six_frames`] offsets.
481    pub frame: usize,
482    /// Amino acids, excluding the stop.
483    pub codons: usize,
484    pub protein: String,
485    pub start_codon: String,
486    pub stop_codon: String,
487    pub wraps_origin: bool,
488}
489
490/// Explicit ORF search parameters.
491#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize)]
492pub struct OrfOptions {
493    pub table: u8,
494    pub min_codons: usize,
495    pub starts: StartPolicy,
496}
497
498/// ORF search result with the effective parameters.
499#[derive(Debug, Clone, PartialEq, Eq, Serialize)]
500pub struct OrfScan {
501    pub options: OrfOptions,
502    pub table_name: &'static str,
503    pub topology: Topology,
504    pub orfs: Vec<Orf>,
505}
506
507/// Six-frame search for complete ORFs: the longest ORF per stop codon (first start after
508/// the previous in-frame stop), at least `min_codons` amino acids. On circular molecules
509/// ORFs may wrap the origin but never exceed one molecule length; frames without any
510/// stop are not reported.
511pub fn find_orfs(
512    sequence: &DnaSeq,
513    topology: Topology,
514    table: u32,
515    min_codons: usize,
516    starts: StartPolicy,
517) -> Result<OrfScan, TranslationError> {
518    if min_codons == 0 {
519        return Err(TranslationError::InvalidMinimum);
520    }
521    let code = GeneticCode::ncbi(table)?;
522    let length = sequence.len();
523    let circular = topology == Topology::Circular;
524    let forward: Coding = forward_bases(sequence, 0, length);
525    let reverse = reverse_complement_coding(forward.clone());
526    let mut found: Vec<Orf> = Vec::new();
527    for (strand, coding) in [
528        (CodingStrand::Forward, forward),
529        (CodingStrand::Reverse, reverse),
530    ] {
531        // Circular molecules are scanned over two turns so wrapping ORFs are complete.
532        let scan: Coding = if circular {
533            coding.iter().chain(coding.iter()).copied().collect()
534        } else {
535            coding
536        };
537        let mut best_per_stop: std::collections::BTreeMap<usize, Orf> =
538            std::collections::BTreeMap::new();
539        for frame in 0..3 {
540            let mut open: Option<usize> = None;
541            let mut index = frame;
542            while index + 3 <= scan.len() {
543                let codon = [scan[index].0, scan[index + 1].0, scan[index + 2].0];
544                if code.is_stop(codon) {
545                    if let Some(begin) = open.take() {
546                        let bases = index + 3 - begin;
547                        let amino_acids = bases / 3 - 1;
548                        if begin < length && bases <= length && amino_acids >= min_codons {
549                            let window = &scan[begin..index + 3];
550                            let protein: String = window
551                                .as_chunks::<3>()
552                                .0
553                                .iter()
554                                .take(amino_acids)
555                                .map(|c| code.translate_codon([c[0].0, c[1].0, c[2].0]))
556                                .collect();
557                            let positions: Vec<usize> = window.iter().map(|(_, p)| *p).collect();
558                            let start = match strand {
559                                CodingStrand::Forward => positions[0],
560                                CodingStrand::Reverse => *positions.last().expect("non-empty"),
561                            };
562                            let orf = Orf {
563                                id: String::new(),
564                                strand,
565                                start,
566                                length: bases,
567                                // Same offsets as `six_frames`: from the 5′ end of each strand.
568                                frame: match strand {
569                                    CodingStrand::Forward => start % 3,
570                                    CodingStrand::Reverse => {
571                                        (length - 1 - (start + bases - 1) % length) % 3
572                                    }
573                                },
574                                codons: amino_acids,
575                                protein,
576                                start_codon: window[..3]
577                                    .iter()
578                                    .map(|(b, _)| char::from(*b))
579                                    .collect(),
580                                stop_codon: window[bases - 3..]
581                                    .iter()
582                                    .map(|(b, _)| char::from(*b))
583                                    .collect(),
584                                wraps_origin: start + bases > length,
585                            };
586                            // Key by the stop's first reference base: one (longest) ORF per stop.
587                            let stop_key = window[bases - 3].1;
588                            let keep = best_per_stop
589                                .get(&stop_key)
590                                .is_none_or(|existing| existing.length < orf.length);
591                            if keep {
592                                best_per_stop.insert(stop_key, orf);
593                            }
594                        }
595                    }
596                } else if open.is_none() && code.is_start(codon, starts) {
597                    open = Some(index);
598                }
599                index += 3;
600            }
601        }
602        found.extend(best_per_stop.into_values());
603    }
604    found.sort_by_key(|orf| (orf.start, orf.strand == CodingStrand::Reverse, orf.length));
605    for (index, orf) in found.iter_mut().enumerate() {
606        orf.id = format!("orf-{:04}", index + 1);
607    }
608    Ok(OrfScan {
609        options: OrfOptions {
610            table: code.id(),
611            min_codons,
612            starts,
613        },
614        table_name: code.name(),
615        topology,
616        orfs: found,
617    })
618}
619
620/// One reading frame of the whole molecule for display: codon `i` covers forward
621/// positions `first + 3i .. first + 3i + 3` (forward) or `first - 3i - 2 ..= first - 3i`
622/// (reverse; `first` is the 5′ base on the reverse strand). Codons that would cross the
623/// origin are omitted.
624#[derive(Debug, Clone, PartialEq, Eq, Serialize)]
625pub struct FrameTranslation {
626    pub strand: CodingStrand,
627    pub offset: usize,
628    pub first: usize,
629    pub protein: String,
630}
631
632/// Six-frame translation for display, frames +1..+3 then -1..-3.
633pub fn six_frames(
634    sequence: &DnaSeq,
635    table: u32,
636) -> Result<Vec<FrameTranslation>, TranslationError> {
637    let code = GeneticCode::ncbi(table)?;
638    let bytes = sequence.as_str().as_bytes();
639    let length = bytes.len();
640    let mut frames = Vec::with_capacity(6);
641    for offset in 0..3 {
642        let protein = bytes
643            .get(offset..)
644            .unwrap_or_default()
645            .as_chunks::<3>()
646            .0
647            .iter()
648            .map(|c| code.translate_codon(*c))
649            .collect();
650        frames.push(FrameTranslation {
651            strand: CodingStrand::Forward,
652            offset,
653            first: offset,
654            protein,
655        });
656    }
657    for offset in 0..3 {
658        let top = length.saturating_sub(offset);
659        let protein = bytes[..top]
660            .rchunks_exact(3)
661            .map(|c| code.translate_codon([complement(c[2]), complement(c[1]), complement(c[0])]))
662            .collect();
663        frames.push(FrameTranslation {
664            strand: CodingStrand::Reverse,
665            offset,
666            first: top.saturating_sub(1),
667            protein,
668        });
669    }
670    Ok(frames)
671}
672
673#[cfg(test)]
674#[path = "translation_tests.rs"]
675mod tests;
dnagent_domain/translation.rs

dnagent_domain/
translation.rs
