Give every variant its H37Rv coordinate and numbering, on any reference

[Paururo]

What changed and why

On a cohort mapped against references other than H37Rv, the H37Rv coordinate and numbering the pipeline gives each variant were wrong or missing. On the 185-sample cohort mapped against two assemblies (E1ASM0035, L7; E1ASM0057, A4), the report showed a wrong H37Rv coordinate for 15,653 of E1ASM0057's 16,296 variant sites (97%). Five faults:

• Several references. The variant positions of every reference were lifted once, from the first reference's sequence, and the map was applied by position alone. LIFT_VARIANTS now runs once per reference, each from its own sequence and for the positions on its own contigs, and the report looks each variant up by contig and position.
• Several contigs. Only the first contig of a FASTA was read, so on a draft assembly a variant of any other contig was lifted with the first contig's sequence. Every contig is now read and chained on its own, in whichever orientation it lies in H37Rv, and the map names both contigs and the strand. The blind-spot mask reaches every contig too.
• The H37Rv numbering (--annotate_canonical, on in the TB profile). Each VCF was annotated in place against the H37Rv snpEff database: another reference's position was read as an H37Rv one, and since the database calls the chromosome Chromosome, nothing was annotated. The unlifted position was then shown as the H37Rv coordinate, over the lift. bin/lift_vcf.py now moves each SNP to its H37Rv position first (H37Rv's base as REF, the allele complemented where the reference lies reversed, the mapping coordinate kept in INFO/OPOS), under the database's chromosome name (--canonical_chrom).
• Resistance tags in the dynamics panel. A variant's HGVS change (p.Ile66Met) was compared with the catalogue's short form (I66M), so no trajectory was ever tagged with its catalogue entry. Both are now read in one form, in H37Rv's gene and numbering first, and a stop or frameshift matches the gene's LoF entry.
• Positions the lift could not interpolate (an indel between two anchors, or no shared unique k-mer for longer than ~2k) were dropped: 4 to 5% of the variant sites. With --align-gaps the stretch between the two anchors is aligned instead; a position inside a stretch H37Rv lacks is reported as such and shown as not in H37Rv.

Checked against minimap2's alignment of each reference to H37Rv:

	Before	Now
Positions of 27 resistance genes	99.95 to 100% same coordinate as minimap2	100%
The report's variant sites, E1ASM0057	0.3% (lifted from the other reference)	99.8%
The report's variant sites, E1ASM0035	96.4%	99.8%
20,000 random positions, either reference	98%	99.8%
E1ASM0035 cut into 10 contigs, 6 reversed	49% wrong	99.8%

None of the placements that differ fits H37Rv worse than minimap2's coordinate: they are ties in identical repeat units, bases at the junction of a structural difference, or places where minimap2 aligned another copy. About 0.2% are left without a coordinate, in PE_PGRS-type repeats where a position could sit on either side of an indel.

Outputs and defaults

• LIFT_VARIANTS runs once per reference instead of once per run: about 25 s and 2.5 GB for a 4.4 Mb reference, on 1 CPU (was 4).
• --annotate_canonical now runs the liftover as well, since it annotates each SNP at its lifted position.
• *.canonical.ann.vcf.gz holds SNPs at their H37Rv position (CHROM Chromosome), each with its mapping coordinate in INFO/OPOS. A SNP whose allele is H37Rv's own base is left out, since in H37Rv numbering nothing changed.
• New parameter --canonical_chrom (default Chromosome, the name the bundled H37Rv snpEff database uses).
• pathotypr_liftover.py lift --global-chain writes a contig-aware map (src_contig src_pos tgt_contig tgt_pos strand, . as the target of a position H37Rv lacks), and has new options --align-gaps, --align-min-identity and --max-align-gap. --contig and --source-contig default to the FASTA's own names. A single-contig reference gets the same map as before.
• The report's SNP matrix and dynamics carry the H37Rv gene beside the H37Rv amino-acid change.

How it was verified

tests/run_tests.sh unit   # 2,176 passed, 39 skipped
tests/run_tests.sh js     # 116 passed
tests/run_tests.sh lint   # clean
mkdocs build --strict     # passes

The coordinates were checked against minimap2 on the two real references as in the table, and rpoB S450L, katG S315T and atpE I66M written in A4 coordinates come out of lift_vcf.py at 761155, 2155168 and 1461242 with H37Rv's base.

Not run here: Nextflow, so the new wiring is left to CI's stub run, and snpEff, so the annotation of the lifted records against the H37Rv database needs a run on the cluster.

Checklist

• tests/run_tests.sh passes locally (unit, js and lint; the pipeline leg needs Nextflow, which is not installed here)
• ruff check . is clean, and I did not restyle code the change does not touch
• New or changed behaviour has a test, and tests pin what the code does today
• A new process in modules/ has a stub: block (not applicable: no new process)
• tests/data/ was not edited by hand (fixtures come from make_test_data.py)
• User-visible changes have a CHANGELOG.md entry and, where relevant, a docs update