Skip to content

Commit ba9a03f

Browse files
committed
Fix iVar deletion convention: POS is the anchor before the gap
iVar reports a deletion at the reference base immediately BEFORE the deleted sequence: POS is that anchor, REF is the base at POS, and ALT='-<bases>' lists the bases deleted starting at POS+1 (andersen-lab/ivar issue #86 and the iVar manual; this is also what nf-core/viralrecon's ivar_variants_to_vcf.py assumes). normalize_ivar_allele instead treated POS as the first deleted base: it read the deleted bases from reference[POS..POS+len] and anchored at POS-1. On real iVar output this rejected every deletion as a FASTA mismatch (or, in a homopolymer, silently shifted it one base left). Read the deleted bases from reference[POS+1..] and anchor at POS, producing the VCF allele (POS, REF=anchor+deleted, ALT=anchor) - symmetric with the insertion branch, which already anchors at POS. Updates the unit test and the iVar deletion scenario to real iVar inputs; the emitted DEL:31:G frameshift row is unchanged, only the input notation (POS 30 anchor instead of 31). 252 lib + 37 integration + 9 doc tests; clippy -D warnings, fmt, 36 scenarios pass.
1 parent 455baff commit ba9a03f

2 files changed

Lines changed: 27 additions & 26 deletions

File tree

src/io/ivar.rs

Lines changed: 22 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -185,34 +185,30 @@ fn normalize_ivar_allele(
185185
}
186186
validate_vcf_allele(deleted_raw, record_idx, chrom, pos, "ALT")?;
187187
let deleted_len = deleted_raw.chars().count();
188-
let deleted_from_ref = reference_slice(references, chrom, pos, deleted_len, record_idx)?;
189-
190-
// The deletion is anchored using reference coordinates, so the bases
191-
// iVar reports as deleted must match the reference at POS..POS+len. A
192-
// mismatch means the FASTA does not correspond to the iVar run (or the
188+
// iVar reports a deletion at the anchor base immediately BEFORE the gap:
189+
// POS is that anchor, REF is the base at POS, and ALT='-<bases>' lists the
190+
// bases deleted starting at POS+1 (see andersen-lab/ivar issue #86 and the
191+
// iVar manual). Convert to a VCF-anchored allele: REF = anchor base plus
192+
// the deleted bases at POS+1.., ALT = the anchor base alone.
193+
let deleted_from_ref =
194+
reference_slice(references, chrom, pos + 1, deleted_len, record_idx)?;
195+
196+
// The bases iVar reports as deleted must match the reference at POS+1.. .
197+
// A mismatch means the FASTA does not correspond to the iVar run (or the
193198
// coordinates are off) and would otherwise silently mis-place the
194199
// deletion; fail fast instead.
195200
if !deleted_raw.eq_ignore_ascii_case(deleted_from_ref) {
196201
return Err(format!(
197202
"iVar record {record_idx}: deletion at {chrom}:{pos} lists deleted bases \
198-
'{deleted_raw}' but the reference reads '{deleted_from_ref}' there. \
199-
Check that the FASTA matches the iVar run."
203+
'{deleted_raw}' but the reference reads '{deleted_from_ref}' at {chrom}:{}. \
204+
Check that the FASTA matches the iVar run.",
205+
pos + 1
200206
)
201207
.into());
202208
}
203209

204-
if pos > 1 {
205-
let anchor = reference_base(references, chrom, pos - 1, record_idx)?;
206-
return Ok(Some((
207-
pos - 1,
208-
format!("{anchor}{deleted_from_ref}"),
209-
anchor,
210-
)));
211-
}
212-
213-
let after_pos = pos + deleted_len;
214-
let anchor = reference_base(references, chrom, after_pos, record_idx)?;
215-
return Ok(Some((pos, format!("{deleted_from_ref}{anchor}"), anchor)));
210+
let anchor = reference_base(references, chrom, pos, record_idx)?;
211+
return Ok(Some((pos, format!("{anchor}{deleted_from_ref}"), anchor)));
216212
}
217213

218214
Ok(Some((pos, ref_allele.to_string(), alt_allele.to_string())))
@@ -401,17 +397,20 @@ chr1\t5\tA\tA\t5\t5\t0.5\t10\tTRUE\n",
401397

402398
#[test]
403399
fn test_load_ivar_tsv_converts_deletion_to_vcf_allele() {
400+
// Real iVar anchors a deletion at the base BEFORE the gap: POS=3 is the
401+
// anchor (G), and ALT='-TA' lists the bases deleted at positions 4-5 of
402+
// reference "ACGTACGT". This converts to the VCF allele (3, GTA, G).
404403
let path = write_temp_ivar(
405404
"REGION\tPOS\tREF\tALT\tREF_DP\tALT_DP\tALT_FREQ\tTOTAL_DP\tPASS\n\
406-
chr1\t3\tG\t-GT\t5\t5\t0.5\t10\tTRUE\n",
405+
chr1\t3\tG\t-TA\t5\t5\t0.5\t10\tTRUE\n",
407406
);
408407
let references = ReferenceMap::from([("chr1".to_string(), "ACGTACGT".to_string())]);
409408
let parsed = load_ivar_tsv(&path, &references).unwrap();
410409
let chr1 = parsed.get("chr1").unwrap();
411410
assert_eq!(chr1.len(), 1);
412-
assert_eq!(chr1[0].position, 2);
413-
assert_eq!(chr1[0].ref_allele, "CGT");
414-
assert_eq!(chr1[0].alt_allele, "C");
411+
assert_eq!(chr1[0].position, 3);
412+
assert_eq!(chr1[0].ref_allele, "GTA");
413+
assert_eq!(chr1[0].alt_allele, "G");
415414
let _ = std::fs::remove_file(path);
416415
}
417416
}

tests/scenarios/scenarios.py

Lines changed: 5 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1223,10 +1223,12 @@
12231223

12241224

12251225
# 29. iVar TSV input - deletion (notacion -SEQ)
1226-
# iVar: POS=31, REF=G, ALT=-G -> borrar pos 31 (G). Anchor en pos 30.
1226+
# iVar ancla la delecion en la base ANTERIOR al hueco: POS=30 (ancla, T),
1227+
# ALT=-G borra la G en pos 31 (primera base del codon 11). Convierte a la
1228+
# variante VCF (30, TG, T) -> DEL:31:G, frameshift desde el codon 11.
12271229
scenario_ivar_deletion = Scenario(
12281230
name="29_ivar_tsv_deletion",
1229-
description="iVar TSV: delecion frameshift -G en pos 31 (notacion -SEQ)",
1231+
description="iVar TSV: delecion frameshift -G (ancla POS=30, borra pos 31; notacion -SEQ)",
12301232
variants=[],
12311233
reads=[
12321234
ReadGroup(
@@ -1238,7 +1240,7 @@
12381240
),
12391241
],
12401242
ivar_records=[
1241-
IvarRecord(pos=31, ref="G", alt="-G", total_dp=20, alt_dp=20, alt_freq=1.0),
1243+
IvarRecord(pos=30, ref="T", alt="-G", total_dp=20, alt_dp=20, alt_freq=1.0),
12421244
],
12431245
expected=[
12441246
ExpectedRow(

0 commit comments

Comments
 (0)