Compute and output Grantham distance, SO terms/impact and MNV shift

Wire the biologist annotations into the indel engine's SNP/MNV path:
compute the Grantham distance and MNV-vs-SNV consequence shift in
gene_path::process_codon, map every variant to a Sequence Ontology term
and impact (so_consequence, covering the engine's StartLost and
FrameshiftDownstreamOfStop change types), and emit them as new TSV
columns (SO Term, Impact, Grantham, MNV Consequence Shift) and VCF INFO
fields (SO, IMPACT, GD, MNVSHIFT) with headers. Output matches the
standalone 1.1.5 feature branch on the example dataset.

Author: Paururo

src/output/common.rs

@@ -34,9 +34,45 @@ fn is_reserved_info_key(key: &str) -> bool {
             | "ERR"
             | "EDP"
             | "EFREQ"
+            | "SO"
+            | "IMPACT"
+            | "GD"
+            | "MNVSHIFT"
     )
 }
 
+/// Map a variant to a Sequence Ontology consequence term and its impact level
+/// (`HIGH` / `MODERATE` / `LOW` / `MODIFIER`), following SnpEff/VEP conventions.
+pub(crate) fn so_consequence(variant: &VariantInfo) -> (&'static str, &'static str) {
+    use crate::variants::ChangeType;
+    if variant.gene == "intergenic" {
+        return ("intergenic_variant", "MODIFIER");
+    }
+    match variant.change_type {
+        ChangeType::Synonymous => ("synonymous_variant", "LOW"),
+        ChangeType::NonSynonymous => ("missense_variant", "MODERATE"),
+        ChangeType::StartLost => ("start_lost", "HIGH"),
+        ChangeType::StopGained | ChangeType::FrameshiftStopGained => ("stop_gained", "HIGH"),
+        ChangeType::StopLost | ChangeType::FrameshiftStopLost => ("stop_lost", "HIGH"),
+        ChangeType::FrameshiftIndel
+        | ChangeType::FrameshiftSynonymous
+        | ChangeType::FrameshiftNonSynonymous
+        | ChangeType::FrameshiftUnknown => ("frameshift_variant", "HIGH"),
+        // A frameshift after the stop codon does not alter the protein.
+        ChangeType::FrameshiftDownstreamOfStop => ("coding_sequence_variant", "MODIFIER"),
+        ChangeType::InFrameIndel => {
+            let ref_len = variant.ref_bases.first().map_or(0, |b| b.len());
+            let alt_len = variant.base_changes.first().map_or(0, |b| b.len());
+            if alt_len > ref_len {
+                ("inframe_insertion", "MODERATE")
+            } else {
+                ("inframe_deletion", "MODERATE")
+            }
+        }
+        ChangeType::IndelOverlap | ChangeType::Unknown => ("coding_sequence_variant", "MODIFIER"),
+    }
+}
+
 #[derive(Default)]
 struct InfoBuilder {
     fields: Vec<(String, String)>,
@@ -192,6 +228,15 @@ pub(crate) fn build_info_string(
     }
     builder.push_text("CT", &variant.change_type.to_string());
     builder.push_text("TYPE", variant_type);
+    let (so_term, impact) = so_consequence(variant);
+    builder.push_text("SO", so_term);
+    builder.push_text("IMPACT", impact);
+    if let Some(gd) = variant.annotations.grantham {
+        builder.push("GD", gd);
+    }
+    if variant.annotations.consequence_shift != crate::variants::ConsequenceShift::NotApplicable {
+        builder.push_text("MNVSHIFT", variant.annotations.consequence_shift.as_str());
+    }
 
     if let Some((sr, srf, srr)) = snp_metrics {
         builder.push("SR", sr);
@@ -534,6 +579,22 @@ pub(crate) fn write_info_header(
         writer,
         "##INFO=<ID=CT,Number=1,Type=String,Description=\"Change type\">"
     )?;
+    writeln!(
+        writer,
+        "##INFO=<ID=SO,Number=1,Type=String,Description=\"Sequence Ontology consequence term\">"
+    )?;
+    writeln!(
+        writer,
+        "##INFO=<ID=IMPACT,Number=1,Type=String,Description=\"Predicted impact (HIGH, MODERATE, LOW, MODIFIER)\">"
+    )?;
+    writeln!(
+        writer,
+        "##INFO=<ID=GD,Number=1,Type=Integer,Description=\"Grantham distance of the (combined) missense change\">"
+    )?;
+    writeln!(
+        writer,
+        "##INFO=<ID=MNVSHIFT,Number=1,Type=String,Description=\"Combined MNV consequence vs. individual SNVs (MNV-gained / MNV-masked / Concordant)\">"
+    )?;
     writeln!(
         writer,
         "##INFO=<ID=EC,Number=1,Type=String,Description=\"Canonical allele event class derived from REF/ALT\">"

src/output/tsv.rs

@@ -15,6 +15,22 @@ fn is_intergenic(variant: &VariantInfo) -> bool {
     variant.gene == "intergenic"
 }
 
+/// Trailing annotation cells appended to every TSV row: SO consequence term,
+/// impact, Grantham distance (with category), and MNV-vs-SNV consequence shift.
+fn annotation_cells(variant: &VariantInfo) -> [String; 4] {
+    let grantham = match variant.annotations.grantham {
+        Some(d) => format!("{d} ({})", crate::utils::grantham_category(d)),
+        None => "-".to_string(),
+    };
+    let (so_term, impact) = super::common::so_consequence(variant);
+    [
+        so_term.to_string(),
+        impact.to_string(),
+        grantham,
+        variant.annotations.consequence_shift.to_string(),
+    ]
+}
+
 /// Render a "Local …" column. Falls back to the protein-wide column when the
 /// local vector is empty (e.g. records produced before this field existed and
 /// then deserialized via `#[serde(default)]`).
@@ -431,6 +447,10 @@ impl TsvWriter {
                 "Event Reverse Reads",
                 "Event Depth",
                 "Event Frequency",
+                "SO Term",
+                "Impact",
+                "Grantham",
+                "MNV Consequence Shift",
             ]
         } else {
             vec![
@@ -450,6 +470,10 @@ impl TsvWriter {
                 "MNV Codon",
                 "Event Class",
                 "Event Components",
+                "SO Term",
+                "Impact",
+                "Grantham",
+                "MNV Consequence Shift",
             ]
         };
         writer.write_record(&header)?;
@@ -474,10 +498,12 @@ impl TsvWriter {
                 if !passes_filters(variant, self.filters)? {
                     continue;
                 }
-                let row = build_tsv_row_with_reads(variant)?;
+                let mut row = build_tsv_row_with_reads(variant)?;
+                row.extend(annotation_cells(variant));
                 self.writer.write_record(&row)?;
             } else {
-                let row = build_tsv_row_without_reads(variant)?;
+                let mut row = build_tsv_row_without_reads(variant)?;
+                row.extend(annotation_cells(variant));
                 self.writer.write_record(&row)?;
             }
         }

src/variants/codon/gene_path.rs

@@ -1,7 +1,46 @@
 //! Genomic-coordinate codon construction, processing, and codon-bounds math.
 
-use crate::utils::{determine_change_type, iupac_aa, reverse_complement};
-use crate::variants::{ChangeType, CodonInfo, Gene, Snp, Strand, VariantInfo, VariantType};
+use crate::utils::{determine_change_type, grantham_distance, iupac_aa, reverse_complement};
+
+/// Coarse severity of an amino-acid change: stop (3) > missense (2) >
+/// synonymous (1) > ambiguous (0).
+fn aa_severity(orig: char, mutated: char) -> u8 {
+    if orig == 'X' || mutated == 'X' {
+        0
+    } else if orig == '*' || mutated == '*' {
+        3
+    } else if orig == mutated {
+        1
+    } else {
+        2
+    }
+}
+
+/// Compare the combined MNV consequence against its individual SNVs.
+fn compute_consequence_shift(
+    orig: char,
+    combined_mut: char,
+    single_muts: &[char],
+) -> ConsequenceShift {
+    if single_muts.len() <= 1 || orig == 'X' {
+        return ConsequenceShift::NotApplicable;
+    }
+    let combined = aa_severity(orig, combined_mut);
+    let max_individual = single_muts
+        .iter()
+        .map(|&m| aa_severity(orig, m))
+        .max()
+        .unwrap_or(0);
+    match combined.cmp(&max_individual) {
+        std::cmp::Ordering::Greater => ConsequenceShift::Gained,
+        std::cmp::Ordering::Less => ConsequenceShift::Masked,
+        std::cmp::Ordering::Equal => ConsequenceShift::Concordant,
+    }
+}
+use crate::variants::{
+    ChangeType, CodonInfo, ConsequenceShift, Gene, Snp, Strand, VariantAnnotations, VariantInfo,
+    VariantType,
+};
 
 use super::grouping::{collect_all_f64, collect_all_usize, merge_original_info};
 
@@ -123,7 +162,9 @@ pub fn process_codon(
     let combined_aa_local = iupac_aa(&combined_change_local);
     let change_type = ChangeType::from_label(&determine_change_type(&combined_change));
 
-    let snp_changes: Vec<String> = codon_info
+    // Single-residue result of each SNV on its own (reused for the per-SNP
+    // columns and the MNV-vs-SNV consequence comparison).
+    let single_aas: Vec<char> = codon_info
         .codon_list
         .iter()
         .map(|snp| {
@@ -132,23 +173,32 @@ pub fn process_codon(
                 Strand::Minus => reverse_complement(&single_codon),
                 Strand::Plus => single_codon,
             };
-            let single_aa = genetic_code.translate_seq(single.as_bytes());
-            iupac_aa(&format!("{orig_aa}{aa_pos}{single_aa}"))
+            genetic_code
+                .translate_seq(single.as_bytes())
+                .chars()
+                .next()
+                .unwrap_or('X')
         })
         .collect();
-    let snp_changes_local: Vec<String> = codon_info
-        .codon_list
+    let snp_changes: Vec<String> = single_aas
         .iter()
-        .map(|snp| {
-            let single_codon = construct_codon(&codon_info, &[snp]);
-            let single = match strand {
-                Strand::Minus => reverse_complement(&single_codon),
-                Strand::Plus => single_codon,
-            };
-            let single_aa = genetic_code.translate_seq(single.as_bytes());
-            iupac_aa(&format!("{orig_aa}{local_aa_pos}{single_aa}"))
-        })
+        .map(|aa| iupac_aa(&format!("{orig_aa}{aa_pos}{aa}")))
         .collect();
+    let snp_changes_local: Vec<String> = single_aas
+        .iter()
+        .map(|aa| iupac_aa(&format!("{orig_aa}{local_aa_pos}{aa}")))
+        .collect();
+
+    let orig_c = orig_aa.chars().next().unwrap_or('X');
+    let mut_c = mut_aa.chars().next().unwrap_or('X');
+    let annotations = VariantAnnotations {
+        grantham: if change_type == ChangeType::NonSynonymous {
+            grantham_distance(orig_c, mut_c)
+        } else {
+            None
+        },
+        consequence_shift: compute_consequence_shift(orig_c, mut_c, &single_aas),
+    };
 
     VariantInfo {
         chrom: chrom.to_string(),
@@ -202,7 +252,7 @@ pub fn process_codon(
             .iter()
             .map(|s| format!("SNV:{}:{}>{}", s.position, s.ref_base, s.base))
             .collect(),
-        annotations: crate::variants::VariantAnnotations::default(),
+        annotations,
     }
 }