nanomonsv norm

src/annotation/colorscb.rs

@@ -9,9 +9,13 @@ use crate::{
     annotation::{Annotation, Annotations},
     io::{readers::TabixReader, vcf::fetch_vcf_with_reader},
     positions::GenomeRange,
-    variant::vcf_variant::{Info, Infos, VcfVariant},
+    variant::vcf_variant::{AlterationCategory, Info, Infos, VcfVariant},
 };
 
+const SV_MATCH_WINDOW: u32 = 1_000;
+const SV_MIN_RECIPROCAL_OVERLAP: f32 = 0.85;
+const SV_MIN_SIZE_RATIO: f32 = 0.7;
+
 #[derive(Debug, Clone, Serialize, Deserialize, PartialEq, Encode, Decode)]
 pub struct ColorsDbEntry {
     pub ac: u32,
@@ -104,11 +108,7 @@ pub fn annotate_colorsdb(
                     }
                 }
 
-                let region = GenomeRange::new(
-                    &variant.position.contig(),
-                    variant.position.position,
-                    variant.position.position + 1,
-                );
+                let region = colorsdb_query_region(variant);
 
                 let colorsdb_records = match fetch_vcf_with_reader(reader, &region) {
                     Ok(v) => v,
@@ -123,7 +123,7 @@ pub fn annotate_colorsdb(
                 };
 
                 for dv in &colorsdb_records {
-                    if dv != variant {
+                    if !colorsdb_variant_matches(variant, dv) {
                         continue;
                     }
 
@@ -155,3 +155,409 @@ pub fn annotate_colorsdb(
 
     Ok(())
 }
+
+fn colorsdb_query_region(variant: &VcfVariant) -> GenomeRange {
+    let contig = variant.position.contig();
+    if variant.has_svtype() {
+        GenomeRange::new(
+            &contig,
+            variant.position.position.saturating_sub(SV_MATCH_WINDOW),
+            variant
+                .position
+                .position
+                .saturating_add(SV_MATCH_WINDOW)
+                .saturating_add(1),
+        )
+    } else {
+        GenomeRange::new(
+            &contig,
+            variant.position.position,
+            variant.position.position + 1,
+        )
+    }
+}
+
+fn colorsdb_variant_matches(query: &VcfVariant, colorsdb: &VcfVariant) -> bool {
+    if colorsdb == query {
+        return true;
+    }
+
+    if query.position.contig != colorsdb.position.contig || !query.has_svtype() {
+        return false;
+    }
+
+    let Some(query_svtype) = query.svtype() else {
+        return false;
+    };
+    if colorsdb.svtype() != Some(query_svtype) {
+        return false;
+    }
+
+    match query.alteration_category() {
+        AlterationCategory::DEL | AlterationCategory::DUP | AlterationCategory::INV => {
+            interval_variant_matches(query, colorsdb)
+        }
+        AlterationCategory::INS => insertion_variant_matches(query, colorsdb),
+        _ => false,
+    }
+}
+
+fn interval_variant_matches(query: &VcfVariant, colorsdb: &VcfVariant) -> bool {
+    let Some((query_start, query_end)) = variant_interval(query) else {
+        return false;
+    };
+    let Some((colorsdb_start, colorsdb_end)) = variant_interval(colorsdb) else {
+        return false;
+    };
+
+    if query_start.abs_diff(colorsdb_start) > SV_MATCH_WINDOW {
+        return false;
+    }
+
+    let query_len = query_end.saturating_sub(query_start).max(1);
+    let colorsdb_len = colorsdb_end.saturating_sub(colorsdb_start).max(1);
+    if size_ratio(query_len, colorsdb_len) < SV_MIN_SIZE_RATIO {
+        return false;
+    }
+
+    let overlap_start = query_start.max(colorsdb_start);
+    let overlap_end = query_end.min(colorsdb_end);
+    if overlap_end <= overlap_start {
+        return false;
+    }
+
+    let overlap = overlap_end - overlap_start;
+    overlap as f32 / query_len as f32 >= SV_MIN_RECIPROCAL_OVERLAP
+        && overlap as f32 / colorsdb_len as f32 >= SV_MIN_RECIPROCAL_OVERLAP
+}
+
+fn insertion_variant_matches(query: &VcfVariant, colorsdb: &VcfVariant) -> bool {
+    if query
+        .position
+        .position
+        .abs_diff(colorsdb.position.position)
+        > SV_MATCH_WINDOW
+    {
+        return false;
+    }
+
+    match (variant_len(query), variant_len(colorsdb)) {
+        (Some(query_len), Some(colorsdb_len)) => {
+            size_ratio(query_len.max(1), colorsdb_len.max(1)) >= SV_MIN_SIZE_RATIO
+        }
+        _ => true,
+    }
+}
+
+fn variant_interval(variant: &VcfVariant) -> Option<(u32, u32)> {
+    if let Some(deletion) = variant.deletion_desc() {
+        return Some((deletion.start, deletion.end));
+    }
+
+    let start = variant.position.position + 1;
+    variant_len(variant).map(|len| (start, start.saturating_add(len.max(1))))
+}
+
+fn variant_len(variant: &VcfVariant) -> Option<u32> {
+    if let Some(len) = variant.deletion_len() {
+        return Some(len);
+    }
+
+    variant
+        .infos
+        .0
+        .iter()
+        .find_map(|info| match info {
+            Info::SVLEN(len) => Some(len.unsigned_abs()),
+            Info::END(end) => Some(end.saturating_sub(variant.position.position + 1)),
+            Info::SVINSLEN(len) => Some(*len),
+            Info::INSLEN(len) => Some(len.unsigned_abs()),
+            _ => None,
+        })
+        .or_else(|| variant.inserted_seq().map(|seq| seq.len() as u32))
+}
+
+fn size_ratio(a: u32, b: u32) -> f32 {
+    a.min(b) as f32 / a.max(b) as f32
+}
+
+#[cfg(test)]
+mod tests {
+    use std::env;
+
+    use crate::{
+        annotation::{Annotation, Annotations},
+        config::Config,
+        io::{readers::TabixReader, vcf::fetch_vcf_with_reader},
+        positions::GenomeRange,
+        variant::{variant_collection::Variants, vcf_variant::VcfVariant},
+    };
+
+    use super::{annotate_colorsdb, colorsdb_variant_matches};
+
+    #[test]
+    fn matches_symbolic_deletion_to_sequence_resolved_colorsdb_deletion() -> anyhow::Result<()> {
+        let query: VcfVariant =
+            "chr1\t16762328\t.\tT\t<DEL>\t.\tPASS\tSVTYPE=DEL;END=16762585".parse()?;
+        let colorsdb: VcfVariant =
+            "chr1\t16762328\t.\tTTTTTTTTTTCTTTTTTCTTTTTTTGAGACGGAGTCTTGCTCTGTTGCTTAGGCTGGATTGCAGTGGCGCAATCTCGGCTCACGCCCGGCTTCCCAGCAACTTTGGGAGGCTGAGGCGGGCGGATCACGAGGTCAGGAGATCGAGACCATCCTGGCTAACACGGTGAAACCCCGTCTCTACTAAAAATACAAAAATTAGCCAGGCGTGGTGGCAGGTACCTGTAGTCCCAGCTACTCGGGAGGCTGAGGCAGGAGAATCACTTGAACCCGGGAGGCGGAGGTTGCAGTGAGCCAAGATCACGCCACTGCACTCCAGCCTGGGCGACAGAGTGAGACTCCATCTCAAAAAAAA\tT\t.\t.\tSVTYPE=DEL;SVLEN=-257;AC=1;AN=2;NS=1;AC_Het=1;AC_Hom=0;AC_Hemi=0;AF=0.5;HWE=1;ExcHet=1;nhomalt=0"
+                .parse()?;
+
+        assert!(colorsdb_variant_matches(&query, &colorsdb));
+        Ok(())
+    }
+
+    #[test]
+    fn matches_nearby_sequence_resolved_insertions_by_size() -> anyhow::Result<()> {
+        let query: VcfVariant =
+            "chr1\t2168842\t.\tG\tAACAGGACCCTGCAGCCCCGGGTGAGGATCATCAGACAGCCTGGA\t.\tPASS\tSVTYPE=INS;SVLEN=46"
+                .parse()?;
+        let colorsdb: VcfVariant =
+            "chr1\t2168842\t.\tG\tGAACAGGACCCTGCAGCCCCGGGTGAGGATCAGACAGCCTGG\t.\t.\tSVTYPE=INS;SVLEN=45;AC=1;AN=2;NS=1;AC_Het=1;AC_Hom=0;AC_Hemi=0;AF=0.5;HWE=1;ExcHet=1;nhomalt=0"
+                .parse()?;
+
+        assert!(colorsdb_variant_matches(&query, &colorsdb));
+        Ok(())
+    }
+
+    fn count_colorsdb_annotations(annotations: &Annotations) -> usize {
+        annotations
+            .store
+            .iter()
+            .map(|entry| {
+                entry
+                    .iter()
+                    .filter(|annotation| matches!(annotation, Annotation::CoLoRSdb(_)))
+                    .count()
+            })
+            .sum()
+    }
+
+    fn seed_existing_annotations(variants: &Variants) -> Annotations {
+        let annotations = Annotations::default();
+
+        for variant in &variants.data {
+            for vcf_variant in &variant.vcf_variants {
+                annotations.insert_update(vcf_variant.hash, &variant.annotations);
+            }
+        }
+
+        annotations
+    }
+
+    fn population_annotation_counts(annotations: &Annotations) -> (usize, usize) {
+        let mut gnomad = 0;
+        let mut dbsnp = 0;
+
+        for entry in annotations.store.iter() {
+            if entry.iter().any(|a| matches!(a, Annotation::GnomAD(_))) {
+                gnomad += 1;
+            }
+            if entry.iter().any(|a| matches!(a, Annotation::DbSnp(_, _))) {
+                dbsnp += 1;
+            }
+        }
+
+        (gnomad, dbsnp)
+    }
+
+    #[derive(Clone)]
+    struct AnnotatedVcfVariant {
+        variant: VcfVariant,
+        callers: Vec<Annotation>,
+    }
+
+    fn sv_vcf_variants(variants: &Variants) -> Vec<AnnotatedVcfVariant> {
+        variants
+            .data
+            .iter()
+            .flat_map(|variant| {
+                let callers = variant
+                    .annotations
+                    .iter()
+                    .filter(|annotation| matches!(annotation, Annotation::Callers(_, _)))
+                    .cloned()
+                    .collect::<Vec<_>>();
+
+                variant
+                    .vcf_variants
+                    .iter()
+                    .filter(|vcf_variant| vcf_variant.has_svtype())
+                    .cloned()
+                    .map(move |vcf_variant| AnnotatedVcfVariant {
+                        variant: vcf_variant,
+                        callers: callers.clone(),
+                    })
+            })
+            .collect()
+    }
+
+    fn short_allele(value: impl ToString) -> String {
+        let value = value.to_string();
+        const MAX_LEN: usize = 80;
+        if value.len() <= MAX_LEN {
+            value
+        } else {
+            format!("{}...({} bp)", &value[..MAX_LEN], value.len())
+        }
+    }
+
+    fn variant_summary(variant: &VcfVariant) -> String {
+        format!(
+            "{}:{} {}>{} {:?}",
+            variant.position.contig(),
+            variant.position.position + 1,
+            short_allele(&variant.reference),
+            short_allele(&variant.alternative),
+            variant.svtype()
+        )
+    }
+
+    fn caller_summary(callers: &[Annotation]) -> String {
+        if callers.is_empty() {
+            "unknown".to_string()
+        } else {
+            callers
+                .iter()
+                .map(ToString::to_string)
+                .collect::<Vec<_>>()
+                .join(", ")
+        }
+    }
+
+    fn print_colorsdb_fetch_diagnostics(
+        variants: &[AnnotatedVcfVariant],
+        colorsdb_path: &str,
+    ) -> anyhow::Result<()> {
+        let mut reader = TabixReader::open(colorsdb_path)?;
+        let mut exact_positions_with_records = 0;
+        let mut exact_records = 0;
+        let mut exact_matches = 0;
+        let mut same_svtype_at_exact_pos = 0;
+        let mut mismatch_examples = Vec::new();
+
+        for annotated_variant in variants {
+            let variant = &annotated_variant.variant;
+            let region = GenomeRange::new(
+                &variant.position.contig(),
+                variant.position.position,
+                variant.position.position + 1,
+            );
+            let records = fetch_vcf_with_reader(&mut reader, &region)?;
+            if records.is_empty() {
+                continue;
+            }
+
+            exact_positions_with_records += 1;
+            exact_records += records.len();
+
+            for colorsdb_variant in records {
+                if colorsdb_variant == *variant {
+                    exact_matches += 1;
+                } else if colorsdb_variant.svtype() == variant.svtype() {
+                    same_svtype_at_exact_pos += 1;
+                    if mismatch_examples.len() < 5 {
+                        mismatch_examples.push(format!(
+                            "caller={} | query={} | colorsdb={}",
+                            caller_summary(&annotated_variant.callers),
+                            variant_summary(variant),
+                            variant_summary(&colorsdb_variant)
+                        ));
+                    }
+                }
+            }
+        }
+
+        let mut reader = TabixReader::open(colorsdb_path)?;
+        let mut window_positions_with_records = 0;
+        let mut window_records = 0;
+        for annotated_variant in variants {
+            let variant = &annotated_variant.variant;
+            let region = GenomeRange::new(
+                &variant.position.contig(),
+                variant.position.position.saturating_sub(1_000),
+                variant.position.position.saturating_add(1_001),
+            );
+            let records = fetch_vcf_with_reader(&mut reader, &region)?;
+            if !records.is_empty() {
+                window_positions_with_records += 1;
+                window_records += records.len();
+            }
+        }
+
+        println!(
+            "CoLoRSdb exact-position fetch: {exact_positions_with_records}/{} query SV positions returned records; {exact_records} records fetched; {same_svtype_at_exact_pos} same-SVTYPE allele mismatches; {exact_matches} exact REF/ALT matches",
+            variants.len()
+        );
+        println!(
+            "CoLoRSdb +/-1kb fetch: {window_positions_with_records}/{} query SV positions returned records; {window_records} records fetched",
+            variants.len()
+        );
+
+        for example in mismatch_examples {
+            println!("CoLoRSdb exact-position mismatch example: {example}");
+        }
+
+        Ok(())
+    }
+
+    #[test]
+    fn annotates_sv_from_pandora_file_with_colorsdb() -> anyhow::Result<()> {
+        // let pandora_path = env::var("PANDORA_COLORSDB_TEST_PATH").map_err(|_| {
+        //     anyhow::anyhow!(
+        //         "set PANDORA_COLORSDB_TEST_PATH to a .pandora, .bit, or .json.gz variants file"
+        //     )
+        // })?;
+
+        let pandora_path = "/mnt/beegfs02/scratch/t_steimle/data/wgs/DUMCO/diag/DUMCO.pandora";
+
+        let colorsdb_path = "/home/t_steimle/ref/hs1/CoLoRSdb.CHM13.v1.2.0.pbsv.jasmine.vcf.gz";
+        // let colorsdb_small_path = "/home/t_steimle/ref/hs1/CoLoRSdb.CHM13.v1.0.0.deepvariant.glnexus.vcf.gz";
+        //
+        // let colorsdb_path = env::var("PANDORA_COLORSDB_TEST_COLORSDB")
+        //     .unwrap_or_else(|_| Config::default().colors_db);
+
+        let variants = Variants::load_from_path(&pandora_path)?;
+        let sv_variants = sv_vcf_variants(&variants);
+
+        println!(
+            "Loaded {} merged variants from {pandora_path}; {} underlying SV VCF records",
+            variants.len(),
+            sv_variants.len()
+        );
+
+        assert!(
+            !sv_variants.is_empty(),
+            "test file contains no VCF records with INFO/SVTYPE"
+        );
+
+        print_colorsdb_fetch_diagnostics(&sv_variants, &colorsdb_path)?;
+
+        let raw_sv_variants = sv_variants
+            .iter()
+            .map(|variant| variant.variant.clone())
+            .collect::<Vec<_>>();
+
+        let empty_annotations = Annotations::default();
+        annotate_colorsdb(&raw_sv_variants, &empty_annotations, &colorsdb_path)?;
+        let n_empty_store = count_colorsdb_annotations(&empty_annotations);
+        println!("CoLoRSdb annotations with empty annotation store: {n_empty_store}");
+
+        let seeded_annotations = seed_existing_annotations(&variants);
+        let (gnomad, dbsnp) = population_annotation_counts(&seeded_annotations);
+        println!(
+            "Seeded annotation store contains {gnomad} GnomAD and {dbsnp} dbSNP annotated VCF records"
+        );
+
+        annotate_colorsdb(&raw_sv_variants, &seeded_annotations, &colorsdb_path)?;
+        let n_seeded_store = count_colorsdb_annotations(&seeded_annotations);
+        println!("CoLoRSdb annotations with Pandora annotations seeded: {n_seeded_store}");
+
+        assert!(
+            n_empty_store > 0 || n_seeded_store > 0,
+            "expected at least one SV from {pandora_path} to receive a CoLoRSdb annotation"
+        );
+
+        Ok(())
+    }
+}

src/callers/nanomonsv.rs

@@ -89,14 +89,17 @@ use crate::{
     },
     config::Config,
     helpers::{is_file_older, remove_dir_if_exists},
-    io::vcf::read_vcf,
+    io::{
+        fasta::{open_indexed_fasta, reference_alternative_from_svlen},
+        vcf::read_vcf,
+    },
     locker::SampleLock,
     pipes::{Initialize, InitializeSolo, ShouldRun, Version},
     run,
     runners::Run,
     variant::{
         variant_collection::VariantCollection,
-        vcf_variant::{Label, Variants},
+        vcf_variant::{Label, SVType, Variants, VcfVariant},
     },
 };
 
@@ -407,6 +410,36 @@ impl Version for NanomonSV {
     }
 }
 
+fn nanomonsv_var_norm(
+    mut v: VcfVariant,
+    fa_reader: &mut noodles_fasta::io::IndexedReader<noodles_fasta::io::BufReader<std::fs::File>>,
+) -> VcfVariant {
+    match v.svtype() {
+        Some(SVType::DEL) => {
+            if let Some(svlen) = v.deletion_len() {
+                if let Ok(reference) =
+                    reference_alternative_from_svlen(fa_reader, &v.position, svlen as i32)
+                {
+                    v.alternative = v.reference.clone();
+                    v.reference = reference;
+                }
+            }
+        }
+
+        Some(SVType::INS) => {
+            if let Some(ins_seq) = v.inserted_seq() {
+                if let Ok(alternative) = format!("{}{}", v.reference, ins_seq).parse() {
+                    v.alternative = alternative;
+                }
+            }
+        }
+
+        _ => {}
+    }
+
+    v
+}
+
 impl Variants for NanomonSV {
     /// Loads and annotates the variants from the NanomonSV PASS VCF.
     ///
@@ -420,8 +453,17 @@ impl Variants for NanomonSV {
 
         info!("Loading variants from {caller}: {vcf_passed}");
 
-        let variants = read_vcf(&vcf_passed)
-            .map_err(|e| anyhow::anyhow!("Failed to read NanomonSV VCF {}.\n{e}", vcf_passed))?;
+        // let variants = read_vcf(&vcf_passed)
+        //     .map_err(|e| anyhow::anyhow!("Failed to read NanomonSV VCF {}.\n{e}", vcf_passed))?;
+
+        let reference = self.config.reference.to_string();
+        let path = PathBuf::from(reference);
+        let mut fa_reader = open_indexed_fasta(&path)?;
+        let variants: Vec<VcfVariant> = read_vcf(&vcf_passed)
+            .map_err(|e| anyhow::anyhow!("Failed to read NanomonSV VCF {vcf_passed}.\n{e}"))?
+            .into_iter()
+            .map(|v| nanomonsv_var_norm(v, &mut fa_reader))
+            .collect();
 
         variants.par_iter().for_each(|v| {
             annotations.insert_update(v.hash(), &add);
@@ -630,7 +672,16 @@ impl Variants for NanomonSVSolo {
             "Loading variant from NanomonSVSolo {} with annotations: {:?}",
             self.id, add
         );
-        let variants = read_vcf(&self.vcf_passed)?;
+        // let variants = read_vcf(&self.vcf_passed)?;
+
+        let reference = self.config.reference.to_string();
+        let path = PathBuf::from(reference);
+        let mut fa_reader = open_indexed_fasta(&path)?;
+        let variants: Vec<VcfVariant> = read_vcf(&self.vcf_passed)
+            .map_err(|e| anyhow::anyhow!("Failed to read NanomonSV VCF {}.\n{e}", self.vcf_passed))?
+            .into_iter()
+            .map(|v| nanomonsv_var_norm(v, &mut fa_reader))
+            .collect();
 
         variants.par_iter().for_each(|v| {
             // let var_cat = v.alteration_category();
@@ -1113,7 +1164,7 @@ pub fn patch_simple_repeat_filter_header<P: AsRef<Path>>(
 #[cfg(test)]
 mod tests {
     use super::*;
-    use crate::helpers::test_init;
+    use crate::{helpers::test_init, variant::vcf_variant::SVType};
 
     #[test]
     fn nanomonsv_version() -> anyhow::Result<()> {
@@ -1132,4 +1183,22 @@ mod tests {
         caller.run()?;
         Ok(())
     }
+
+    #[test]
+    fn nanomonsv_load() -> anyhow::Result<()> {
+        test_init();
+        let config = Config::default();
+
+        let caller = NanomonSV::initialize("DUMCO", &config)?;
+        let annotations = Annotations::default();
+        let v = caller.variants(&annotations)?;
+
+        v.variants.iter().for_each(|vv| {
+            if let Some(SVType::INS) = vv.svtype() {
+                println!("{vv:?}")
+            }
+        });
+
+        Ok(())
+    }
 }

src/callers/savana.rs

@@ -413,7 +413,7 @@ impl Variants for Savana {
             .map(|mut v| {
                 if let Some(crate::variant::vcf_variant::SVType::INS) = v.svtype() {
                     if let Some(sequence) = get_longest_sequence(&insertion_store, &v.id) {
-                        if let Ok(parsed) = sequence.parse() {
+                        if let Ok(parsed) = format!("{}{}", v.reference, sequence).parse() {
                             v.alternative = parsed;
                             v.infos.0.push(Info::INSLEN(sequence.len() as i32));
                             v.infos.0.push(Info::INSSEQ(sequence));

src/io/fasta.rs

@@ -15,6 +15,8 @@ use std::{
 use anyhow::Context;
 use noodles_fasta::io::{BufReader as NoodlesBufReader, IndexedReader};
 
+use crate::{positions::GenomePosition, variant::vcf_variant::ReferenceAlternative};
+
 /// Open an indexed FASTA file (requires a `.fai` index alongside the path).
 ///
 /// # Errors
@@ -93,6 +95,106 @@ pub fn sequence_range(
     Ok(String::from_utf8(record.sequence().as_ref().to_vec())?.to_uppercase())
 }
 
+/// Fetch a `ReferenceAlternative` allele from FASTA using an SV anchor position and `SVLEN`.
+///
+/// `position` is the 0-based VCF anchor position used throughout the codebase.
+///
+/// - For negative `SVLEN` values (`DEL`), the alternative allele is the anchor
+///   base at `position`.
+/// - For positive `SVLEN` values (`INS`/`DUP`-like sequence-resolved alleles),
+///   this returns the anchored reference span `[position, position + SVLEN]`.
+///
+/// This mirrors VCF's anchored allele representation, where indels include the
+/// base immediately before the changed sequence.
+pub fn reference_alternative_from_svlen(
+    fasta_reader: &mut IndexedReader<noodles_fasta::io::BufReader<File>>,
+    position: &GenomePosition,
+    svlen: i32,
+) -> anyhow::Result<ReferenceAlternative> {
+    anyhow::ensure!(svlen != 0, "SVLEN must be non-zero");
+
+    let start0 = usize::try_from(position.position).context("GenomePosition overflows usize")?;
+    let len = usize::try_from(svlen.unsigned_abs()).context("SVLEN overflows usize")?;
+    let end0_inclusive = if svlen < 0 {
+        start0
+    } else {
+        start0
+            .checked_add(len)
+            .context("SVLEN FASTA range end coordinate overflow")?
+    };
+
+    sequence_range(fasta_reader, &position.contig(), start0, end0_inclusive)?
+    .parse()
+    .context("failed to parse FASTA sequence as ReferenceAlternative")
+}
+
+/// Fetch the VCF reference allele for an insertion from FASTA.
+///
+/// Insertions are represented with an anchor base in VCF, so the reference
+/// allele is just the base at `position`. `svlen` is validated for consistency
+/// but cannot identify the inserted sequence by itself.
+pub fn insertion_reference_from_svlen(
+    fasta_reader: &mut IndexedReader<noodles_fasta::io::BufReader<File>>,
+    position: &GenomePosition,
+    svlen: i32,
+) -> anyhow::Result<ReferenceAlternative> {
+    anyhow::ensure!(svlen > 0, "insertion SVLEN must be positive, got {svlen}");
+
+    let pos0 = usize::try_from(position.position).context("GenomePosition overflows usize")?;
+    sequence_range(fasta_reader, &position.contig(), pos0, pos0)?
+        .parse()
+        .context("failed to parse FASTA anchor base as ReferenceAlternative")
+}
+
+/// Build an insertion alternative allele from FASTA anchor base plus inserted bases.
+///
+/// This is the correct helper for ordinary insertions where the inserted
+/// sequence comes from the caller (`SVINSSEQ`, sequence-resolved ALT, etc.).
+pub fn insertion_alternative_from_sequence(
+    fasta_reader: &mut IndexedReader<noodles_fasta::io::BufReader<File>>,
+    position: &GenomePosition,
+    inserted_sequence: &str,
+) -> anyhow::Result<ReferenceAlternative> {
+    anyhow::ensure!(
+        !inserted_sequence.is_empty(),
+        "inserted sequence must not be empty"
+    );
+
+    let pos0 = usize::try_from(position.position).context("GenomePosition overflows usize")?;
+    let anchor = sequence_range(fasta_reader, &position.contig(), pos0, pos0)?;
+    format!("{anchor}{}", inserted_sequence.to_uppercase())
+        .parse()
+        .context("failed to parse insertion alternative as ReferenceAlternative")
+}
+
+/// Build a reference-derived insertion alternative from FASTA using `SVLEN`.
+///
+/// This is only valid for tandem-duplication-like insertions where the inserted
+/// bases are copied from the reference immediately after the anchor. For novel
+/// inserted sequence, use [`insertion_alternative_from_sequence`] instead.
+pub fn reference_derived_insertion_alternative_from_svlen(
+    fasta_reader: &mut IndexedReader<noodles_fasta::io::BufReader<File>>,
+    position: &GenomePosition,
+    svlen: i32,
+) -> anyhow::Result<ReferenceAlternative> {
+    anyhow::ensure!(svlen > 0, "insertion SVLEN must be positive, got {svlen}");
+
+    let start0 = usize::try_from(position.position).context("GenomePosition overflows usize")?;
+    let len = usize::try_from(svlen).context("SVLEN overflows usize")?;
+    let end0_inclusive = start0
+        .checked_add(len)
+        .context("SVLEN FASTA range end coordinate overflow")?;
+
+    sequence_range(
+        fasta_reader,
+        &position.contig(),
+        start0,
+        end0_inclusive,
+    )?
+    .parse()
+    .context("failed to parse FASTA sequence as ReferenceAlternative")
+}
+
 /// A single-contig FASTA file produced by [`split_fasta`].
 pub struct ContigFasta {
     pub name: String,