fix: two-step VNTR downsampling to match real exome coverage profiles

[berntpopp]

Summary

Fixes three bugs in downsample_mode: "vntr" that made it produce completely unrealistic coverage profiles:

1. ^region syntax broken -- downsample_bam() used samtools view ^chr1:... to extract non-VNTR reads, which silently returns zero reads. All flanking reads were dropped. Replaced with two-pass BED-based approach (-L/-U for extraction, -L/-s for downsampling).

2. Single-step VNTR-only downsampling -- Downsampled only the VNTR region, leaving flanking at 700x+ while reducing VNTR to 150x. Produced inverted VNTR:flanking ratio (0.5x vs real 1.4x). Replaced with two-step:

   • Step 1: downsample entire BAM so non-VNTR BED = target coverage
   • Step 2: downsample VNTR region to match empirical ratio (default 1.4, configurable via vntr_to_flanking_ratio)

3. Same seed for both steps -- samtools -s uses hash-based subsampling, so reads surviving step 1's lower fraction always survive step 2 at the same seed. Step 2 now uses seed+1.

Validation (1,043 real exomes vs 50 simulations)

Metric	Real (median)	Simulated (median)	Sim/Real	Verdict
VNTR mean	169x	204x	1.21	GOOD
VNTR median	106x	77x	0.73	GOOD
Non-VNTR BED	107x	150x	1.40	GOOD
VNTR:BED ratio	1.6x	1.4x	0.87	GOOD
Total reads	29,784	8,827	0.30	POOR*
VNTR % uncov	9.9%	24.4%	2.48	POOR*

*POOR metrics are pre-existing read generation characteristics, not downsampling issues.

Full validation report: .planning/reports/2026-04-06-vntr-downsampling-validation.md

Test plan

• 1,423 tests pass (89% coverage)
• ruff clean
• mypy clean
• 50 simulations (VNTR lengths 30-90) validated against 1,043 real exomes
• CI passes across Python 3.10-3.12

Closes #89

🤖 Generated with Claude Code

[codecov-commenter]

⚠️ Please install the to ensure uploads and comments are reliably processed by Codecov.

Codecov Report

❌ Patch coverage is 92.30769% with 4 lines in your changes missing coverage. Please review.

Files with missing lines	Patch %	Lines
muc_one_up/read_simulator/stages/alignment.py	92.00%	2 Missing and 2 partials ⚠️

📢 Thoughts on this report? Let us know!

[Copilot]

Pull request overview

This PR reworks Illumina downsampling behavior for downsample_mode: "vntr" to better match empirical exome VNTR vs flanking coverage patterns by using a two-step approach based on a non-VNTR target BED and a configurable VNTR:flanking ratio.

Changes:

• Implement two-step VNTR downsampling: (1) downsample entire BAM to hit target coverage over a non-VNTR BED, then (2) downsample VNTR region to a target VNTR:flanking ratio using a different seed.
• Update alignment-stage unit tests to reflect the new call sequence/requirements for VNTR downsampling.
• Add a validation report documenting real-vs-simulated coverage metrics.

Reviewed changes

Copilot reviewed 3 out of 3 changed files in this pull request and generated 4 comments.

File	Description
muc_one_up/read_simulator/stages/alignment.py	Implements the two-step VNTR downsampling flow, adds sample_target_bed requirement, and logs post-downsampling VNTR:flanking ratio.
tests/read_simulator/test_alignment_stage.py	Updates tests to match the new two-step VNTR logic and adjusted validation paths.
.planning/reports/2026-04-06-vntr-downsampling-validation.md	Adds validation results comparing simulated coverage distributions to a real exome cohort.

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

[Copilot]

downsample_bam() is called with a samtools-style region string (chr:start-end), but the current implementation in wrappers/samtools_coverage.py converts that string into a BED for samtools view -L. BED coordinates are 0-based, end-exclusive, so writing start/end directly will shift/truncate the region by 1bp. This can cause boundary reads to be incorrectly treated as VNTR vs non-VNTR during step 2. Consider fixing the region→BED conversion in downsample_bam() (or passing a correctly-formatted BED) so the extracted region exactly matches the VNTR coordinates used for coverage measurement (samtools depth -r).

[Copilot]

vntr_to_flanking_ratio is pulled from config and used to compute vntr_target without validation. If it is missing/typed as a string, <=0, or otherwise invalid, step 2 can behave unexpectedly (e.g., downsample VNTR to zero when ratio=0) or raise at runtime. Consider coercing to float and validating it is >0 (and possibly setting an upper bound) with a clear ConfigurationError if invalid.

Suggested change

	vntr_to_flanking_ratio: float = rs_config.get("vntr_to_flanking_ratio", 1.4)
	raw_vntr_to_flanking_ratio = rs_config.get("vntr_to_flanking_ratio", 1.4)
	try:
	vntr_to_flanking_ratio = float(raw_vntr_to_flanking_ratio)
	except (TypeError, ValueError) as exc:
	raise ConfigurationError(
	"Invalid read_simulator config: 'vntr_to_flanking_ratio' "
	f"must be a positive float, got {raw_vntr_to_flanking_ratio!r}"
	) from exc
	if (
	vntr_to_flanking_ratio <= 0
	or vntr_to_flanking_ratio != vntr_to_flanking_ratio
	or vntr_to_flanking_ratio in (float("inf"), float("-inf"))
	):
	raise ConfigurationError(
	"Invalid read_simulator config: 'vntr_to_flanking_ratio' "
	f"must be a positive finite float, got {raw_vntr_to_flanking_ratio!r}"
	)

[Copilot]

VNTR mode now requires sample_target_bed, but there is no unit test covering the new ConfigurationError path when downsample_mode: "vntr" is selected and sample_target_bed is missing. Adding a test for this would prevent regressions in config validation.

[Copilot]

This report includes absolute local filesystem paths (e.g., /home/..., /tmp/...). If this repository/branch is shared beyond the author’s machine, these paths can leak personal usernames or internal directory structure and reduce reproducibility. Consider replacing them with generalized placeholders or relative paths (and, if needed, describing the layout in words).

Suggested change

	- Real exome BAMs: `/home/bernt-popp/development/vntyper-analyses/results/screening/Bernt/data/cerkid-exome-chr1/`
	- Simulated BAMs: `/tmp/sim_batch/len_{30..90}/sim_{len}.001.simulated_reads_downsampled.bam`
	- Real exome BAMs: local CerKiD exome chr1 BAM directory (machine-specific path omitted)
	- Simulated BAMs: local simulation batch output pattern `len_{30..90}/sim_{len}.001.simulated_reads_downsampled.bam` (base directory omitted)