norm: realign large biallelic indels in linear time

[davmlaw]

Human PR summary

Faster left-align algorithm for <DEL> and <DUP>. 100mb example goes from 65 mins to 1.5 seconds

This is a C implementation of existing code - which was recently used to perform HGVS normalization when I discovered the faster algorithm, I had the agent investigate what other tools did, run benchmarks and evaluate whether bcftools would gain from it as well. I had the agent make a few minor edits based on my feedback, it looks good and tests pass but I haven't written C for ~20 years

This has been previously reported in #2500 - 23 kb DRAGEN <DUP> in #2500, which led to --no-realign=NUM

AI summary

Motivation

Since #1919 and #2145, bcftools norm -f left-aligns symbolic <DEL> and <DUP> alleles by expanding them to sequence, which made realignment reachable for CNV-scale events. Left-aligning in realign_left() pads every allele with aln_win (100) reference bases and trims again until the shift stops, and each padding step memmoves the whole long allele, so a tandem duplication of length L costs about L/100 copies of L bytes. A 10 Mb <DUP> takes 40 s and a 100 Mb one 65 minutes (measured below). That is what a user ran into with a 23 kb DRAGEN <DUP> in #2500, which led to --no-realign=NUM, and it is the realignment behind the SV reports in #2330. #2601/#2602 extend <DUP> realignment further, so more large events will take this path.

Change

Shifting a biallelic insertion or deletion one base to the left is possible iff the reference base to the left equals the last base of the inserted or deleted sequence X, and the step rotates X right by one base. So the shift can be counted with a single cyclic comparison and X rotated once, instead of rediscovering it 100 bases at a time. realign_left() now hands the biallelic case to a new realign_left_indel() the first time an allele is trimmed away (which is the point at which the record is known to be a plain insertion or deletion). Multiallelic records keep using the generic loop. Right-alignment only appends to the alleles and already runs in linear time, so it is unchanged.

The function comment carries a worked example. The same approach is used by roll_left() in the biocommons bioutils package (bioutils.normalize, 2019) and is the basis of the ambiguity-region expansion in NCBI's SPDI notation (Holmes et al., Bioinformatics 2020, doi:10.1093/bioinformatics/btz856); the pad-and-trim form is the one described for vt normalize (Tan, Abecasis & Kang, Bioinformatics 2015, doi:10.1093/bioinformatics/btv112).

The output is identical to the generic loop, including the case of the bases when the reference is soft-masked: bases rotated in from the reference keep the reference case, the rest of X stays uppercase.

Benchmarks

Single-record VCFs on a synthetic reference with planted tandem duplications, bcftools norm -f ref.fa one.vcf, best of 3 with a warm cache, Linux x86-64. "floor" is the same command with --no-realign, i.e. parsing and REF checking only.

event	size	before	after	floor
literal dup, left	1 Mb	0.26 s	0.016 s	0.008 s
literal dup, left	10 Mb	40.0 s	0.16 s	0.05 s
literal del, left	10 Mb	40.2 s	0.17 s	0.04 s
<DEL>, left	10 Mb	40.1 s	0.087 s	0.001 s
<DUP>, left	10 Mb	40.1 s	0.091 s	0.001 s
<DUP>, left	100 Mb	3931 s	1.56 s
literal dup, right	10 Mb	0.083 s	0.086 s	0.047 s
200k indels of 1-30 bp, left		0.234 s	0.218 s
200k indels of 1-30 bp, right		0.356 s	0.367 s

Peak RSS for the 100 Mb <DUP> was 200 MB before and 300 MB after, both dominated by the allele strings.

Checks

• test/test.pl -f test_vcf_norm passes (170 tests, one new).

• Output compared byte for byte with the previous code on ~270k random records: insertions, deletions, insertions copying nearby sequence, delins, SNVs, multiallelic records, symbolic <DEL>/<DUP>, positions at both ends of the contig, soft-masked references with IUPAC codes and N, left- and right-aligned. All identical.

• Output also compared byte for byte with the previous code on real VCFs against the Ensembl GRCh38 toplevel reference, all identical:

  input	records	realigned to a new POS
  ClinVar GRCh38 (clinvar.vcf.gz), left-aligned	4,438,773	0
  ClinVar GRCh38, --right-align	4,438,773	223,946
  GIAB HG001 GRCh38 v4.2.1 benchmark calls, left-aligned	3,893,341	39
  gnomAD v4.0 SV, chr21+chr22 <DEL>/<DUP> with |SVLEN| ≤ 500 kb, left-aligned (-c s, as gnomAD writes REF=N)	51,754	13,395

  ClinVar and the GIAB calls are already left-aligned, so the random records and the gnomAD SVs are what exercise the new path.

• An AddressSanitizer/LeakSanitizer build runs the norm test suite, 68k of the random records and a 1 Mb tandem duplication with no reports.

• New test norm.realign-indel covers a rotated shift that picks up lowercase bases, the first-base anchor rule, a shift stopped by an N, symbolic alleles and a multiallelic record.

Assisted-by: Claude:claude-fable-5-1