Annotate insertions at internal exon junctions as coding, not intergenic

A VCF insertion places its bases after the anchor, so an insertion anchored at
the last base of an internal coding exon lands at the exon-exon junction inside
the spliced CDS. The half-open exon interval used for overlap and per-segment
anchoring (< segment.end) excluded that anchor, so such insertions were dropped
from the transcript model and reported as intergenic.

Add Gene::insertion_at_internal_junction, which recognises an anchor on an
internal exon's last base (using transcript order: not the last segment on the
plus strand, not the first on the minus strand, so a CDS-terminal anchor stays
UTR). Use it in gene_overlaps_variant, variant_overlaps_coding_model, and the
per-segment insertion anchoring so the inserted bases are placed at the junction
of the reconstructed CDS and the protein effect is computed. The change is
additive: in-exon variants and every existing classification are unchanged.

Adds plus- and minus-strand regression tests (the minus-strand inserted bases
are reverse-complemented at the junction). 258 -> 260 lib tests; clippy, fmt and
the 36-scenario suite pass.

Author: Paururo

src/io/annotation.rs

@@ -7,7 +7,7 @@
 
 use super::vcf::VcfPosition;
 use crate::error::AppResult;
-use crate::variants::Gene;
+use crate::variants::{AlleleComponentKind, Gene};
 use std::collections::HashMap;
 
 mod cds_model;
@@ -52,9 +52,20 @@ pub fn gene_overlaps_variant(gene: &Gene, variant: &VcfPosition) -> bool {
     if gene.cds_segments.is_empty() {
         return variant.overlaps_interval(gene.start, gene.end);
     }
-    gene.cds_segments
+    if gene
+        .cds_segments
         .iter()
         .any(|segment| variant.overlaps_interval(segment.start, segment.end))
+    {
+        return true;
+    }
+    // An insertion anchored at an internal exon-exon junction is coding even
+    // though overlaps_interval (half-open) excludes the exon's last base; keep it
+    // out of the intergenic set so it is annotated on the spliced CDS instead.
+    variant.event().components.iter().any(|component| {
+        matches!(component.kind, AlleleComponentKind::Insertion)
+            && gene.insertion_at_internal_junction(component.position)
+    })
 }
 
 fn gene_has_variant(gene: &Gene, snp_list: &[VcfPosition]) -> bool {

src/variants/codon/allele_apply.rs

@@ -57,7 +57,7 @@ pub(super) fn apply_allele_to_transcript(
                     }
                 }
                 AlleleComponentKind::Insertion => {
-                    if insertion_anchor_in_segment(segment, component.position) {
+                    if insertion_anchor_in_segment(gene, segment, component.position) {
                         match insertions.get(&component.position) {
                             Some(existing)
                                 if !existing.eq_ignore_ascii_case(&component.alt_allele) =>

src/variants/codon/tests.rs

@@ -476,6 +476,94 @@ fn test_minus_strand_legacy_deletion_reports_consistent_codons() {
     assert_eq!(variants[0].mnv_codon.as_deref(), Some("GGG"));
 }
 
+#[test]
+fn test_plus_strand_insertion_at_internal_exon_junction_is_coding() {
+    // Regression: an insertion anchored at the last base of an internal exon
+    // lands at the exon-exon junction inside the spliced CDS, so it must be
+    // annotated as a coding indel, not dropped (and later reported intergenic).
+    // exon1 = [1,3] "ATG", exon2 = [10,15] "AAATTT"; CDS = "ATGAAATTT"
+    // (Met-Lys-Phe). Inserting GGG after pos 3 -> spliced CDS "ATGGGGAAATTT".
+    let gene = Gene {
+        name: "tx_cds".to_string(),
+        start: 1,
+        end: 15,
+        strand: Strand::Plus,
+        phase: 0,
+        protein_offset: 0,
+        transcript_id: Some("tx1".to_string()),
+        cds_segments: vec![
+            CdsSegment { start: 1, end: 3 },
+            CdsSegment { start: 10, end: 15 },
+        ],
+    };
+    let reference = crate::io::Reference {
+        sequence: "ATGCCCCCCAAATTT",
+    };
+    let variants = crate::variants::get_mnv_variants_for_gene(
+        &gene,
+        &[crate::io::VcfPosition {
+            position: 3,
+            ref_allele: "G".to_string(),
+            alt_allele: "GGGG".to_string(),
+            original_dp: None,
+            original_freq: None,
+            original_info: None,
+        }],
+        &reference,
+        "chr1",
+        crate::genetic_code::GeneticCode::default(),
+    );
+    assert_eq!(variants.len(), 1);
+    assert_ne!(variants[0].gene, "intergenic");
+    assert_eq!(variants[0].variant_type, VariantType::Indel);
+    assert_eq!(variants[0].change_type, ChangeType::InFrameIndel);
+    assert_eq!(variants[0].event_components, vec!["INS:3:+GGG".to_string()]);
+}
+
+#[test]
+fn test_minus_strand_insertion_at_internal_exon_junction_is_coding() {
+    // Same junction case on the minus strand. cds_segments are in transcript
+    // order (descending genomic): exon0 = [10,15], exon1 = [1,3]. CDS =
+    // RC("TTTCAT") + RC("AAA") = "ATGAAATTT". An insertion of GGG after genomic
+    // pos 3 (exon1's genomic end = transcript-5' boundary, an internal junction)
+    // becomes RC("GGG")="CCC" at the exon0|exon1 junction -> "ATGAAACCCTTT".
+    let gene = Gene {
+        name: "tx_cds".to_string(),
+        start: 1,
+        end: 15,
+        strand: Strand::Minus,
+        phase: 0,
+        protein_offset: 0,
+        transcript_id: Some("tx1".to_string()),
+        cds_segments: vec![
+            CdsSegment { start: 10, end: 15 },
+            CdsSegment { start: 1, end: 3 },
+        ],
+    };
+    let reference = crate::io::Reference {
+        sequence: "AAACCCCCCTTTCAT",
+    };
+    let variants = crate::variants::get_mnv_variants_for_gene(
+        &gene,
+        &[crate::io::VcfPosition {
+            position: 3,
+            ref_allele: "A".to_string(),
+            alt_allele: "AGGG".to_string(),
+            original_dp: None,
+            original_freq: None,
+            original_info: None,
+        }],
+        &reference,
+        "chr1",
+        crate::genetic_code::GeneticCode::default(),
+    );
+    assert_eq!(variants.len(), 1);
+    assert_ne!(variants[0].gene, "intergenic");
+    assert_eq!(variants[0].variant_type, VariantType::Indel);
+    assert_eq!(variants[0].change_type, ChangeType::InFrameIndel);
+    assert_eq!(variants[0].event_components, vec!["INS:3:+GGG".to_string()]);
+}
+
 // H5: an in-frame insertion that introduces a stop codon must be classified
 // as Stop gained rather than the generic In-frame Indel. Inserting TAA after
 // the start codon of ATGAAATTT yields ATG-TAA-AAA-TTT (Met-*-Lys-Phe).

src/variants/codon/transcript_model.rs

@@ -56,26 +56,41 @@ pub(super) fn transcript_offset_for_position(gene: &Gene, position: usize) -> Op
     None
 }
 
-/// Whether an insertion anchored at `position` falls strictly inside this CDS
-/// segment. The upper bound is exclusive (`< segment.end`) on purpose: a VCF
-/// insertion places the inserted bases *after* the anchor, so an anchor at the
-/// segment's last base sits on the exon|intron boundary, where the inserted
-/// sequence is genomically ambiguous between "end of this exon" and "start of
-/// the intron" (spliced out). Such boundary insertions are intentionally left
-/// out of the spliced-CDS model rather than assigned a possibly-wrong codon;
-/// they surface as intergenic instead. Anchors inside an exon are unaffected.
-pub(super) fn insertion_anchor_in_segment(segment: &CdsSegment, position: usize) -> bool {
-    position >= segment.start && position < segment.end
+/// Whether an insertion anchored at `position` belongs to this CDS `segment`
+/// when reconstructing the spliced transcript. True when the anchor is strictly
+/// inside the exon, or when it sits exactly on the exon's last base and that base
+/// is an internal exon-exon junction: a VCF insertion places its bases *after*
+/// the anchor, so an anchor at an internal exon's last base lands the inserted
+/// sequence at the junction with the next exon, inside the spliced CDS. An anchor
+/// at the CDS terminus is UTR and is excluded (see
+/// [`Gene::insertion_at_internal_junction`]).
+pub(super) fn insertion_anchor_in_segment(
+    gene: &Gene,
+    segment: &CdsSegment,
+    position: usize,
+) -> bool {
+    (position >= segment.start && position < segment.end)
+        || (position == segment.end && gene.insertion_at_internal_junction(position))
 }
 
 pub(super) fn variant_overlaps_coding_model(gene: &Gene, variant: &VcfPosition) -> bool {
     if !has_transcript_cds_model(gene) {
         return variant.overlaps_interval(gene.start, gene.end);
     }
 
-    gene.cds_segments
+    if gene
+        .cds_segments
         .iter()
         .any(|segment| variant.overlaps_interval(segment.start, segment.end))
+    {
+        return true;
+    }
+    // An insertion anchored at an internal exon-exon junction is coding even
+    // though overlaps_interval (half-open) excludes the exon's last base.
+    variant.event().components.iter().any(|component| {
+        matches!(component.kind, AlleleComponentKind::Insertion)
+            && gene.insertion_at_internal_junction(component.position)
+    })
 }
 
 pub(super) fn variant_touched_transcript_offsets(gene: &Gene, variant: &VcfPosition) -> Vec<usize> {
@@ -93,7 +108,7 @@ pub(super) fn variant_touched_transcript_offsets(gene: &Gene, variant: &VcfPosit
             }
             AlleleComponentKind::Insertion => {
                 for segment in &gene.cds_segments {
-                    if insertion_anchor_in_segment(segment, component.position) {
+                    if insertion_anchor_in_segment(gene, segment, component.position) {
                         if let Some(offset) =
                             transcript_offset_for_position(gene, component.position)
                         {

src/variants/types.rs

@@ -195,6 +195,28 @@ pub struct Gene {
     pub cds_segments: Vec<CdsSegment>,
 }
 
+impl Gene {
+    /// Whether an insertion anchored at genomic `position` lands at an internal
+    /// exon-exon junction of the spliced CDS — i.e. immediately after the last
+    /// base of an internal coding exon. The inserted bases then fall between two
+    /// exons in the spliced transcript, so the insertion is coding, even though
+    /// the exon's terminal base is excluded from the half-open exon interval used
+    /// elsewhere. An anchor at the CDS terminus (after the last exon on the plus
+    /// strand, before the first on the minus strand) is UTR, not an internal
+    /// junction. `cds_segments` is in transcript order, so the terminal exon is
+    /// the last entry on the plus strand and the first on the minus strand.
+    pub(crate) fn insertion_at_internal_junction(&self, position: usize) -> bool {
+        let segment_count = self.cds_segments.len();
+        self.cds_segments.iter().enumerate().any(|(idx, segment)| {
+            position == segment.end
+                && match self.strand {
+                    Strand::Plus => idx + 1 < segment_count,
+                    Strand::Minus => idx > 0,
+                }
+        })
+    }
+}
+
 #[derive(Debug, Serialize, Deserialize, Clone)]
 pub struct Snp {
     pub index: usize,