fix: VNTR coverage calculation bug and penalty factor recalibration (#85, #86) - #87
Conversation
…ions (#85) Without -a, samtools depth only reports positions with >=1 read, inflating the mean coverage. This aligns with the sister function calculate_vntr_coverage() which already uses -a. Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
After downsampling, re-measure VNTR/target coverage and log the result. Emit a warning if the actual coverage deviates more than 20% from target. Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
…86) Empirically re-derived from 1,043 CerKiD Berlin Twist v2 exomes. Median implied penalty from observed VNTR/flanking ratio 3.66 vs simulated base ratio 9.32. Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
|
Codecov Report✅ All modified and coverable lines are covered by tests. 📢 Thoughts on this report? Let us know! |
- Track post-validation depth file in intermediate_files for cleanup - Add test for post-downsampling validation in non_vntr mode - Fix doc: "Coverage ratio" → "VNTR:Flanking ratio" in expected behavior table Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
There was a problem hiding this comment.
Pull request overview
This PR addresses inflated VNTR coverage measurements during simulation by ensuring zero-depth bases are included in mean coverage calculations, adds a post-downsampling coverage validation step, and recalibrates the VNTR capture efficiency penalty factor based on a larger empirical cohort.
Changes:
- Fix mean coverage calculation to include zero-coverage positions via
samtools depth -a. - Add post-downsampling coverage re-measurement + deviation logging/warning.
- Update default VNTR penalty factor from 0.375 → 0.39 and refresh tests/docs/config defaults accordingly.
Reviewed changes
Copilot reviewed 9 out of 9 changed files in this pull request and generated 3 comments.
Show a summary per file
| File | Description |
|---|---|
muc_one_up/read_simulator/utils/samtools.py |
Updates calculate_mean_coverage() to call samtools depth -a. |
muc_one_up/read_simulator/stages/alignment.py |
Adds post-downsampling validation coverage measurement and warning on large deviation. |
muc_one_up/read_simulator/vntr_efficiency.py |
Updates default penalty factor and inline documentation. |
muc_one_up/config.py |
Updates config schema default penalty factor. |
config.json |
Updates shipped/default config penalty factor value. |
docs/guides/vntr-capture-efficiency.md |
Rewrites/recalibrates documentation for penalty factor and expected behavior. |
tests/test_samtools.py |
Adds regression test asserting -a is used and zeros contribute to the mean. |
tests/read_simulator/test_alignment_stage.py |
Updates mocks/tests for the additional post-downsampling coverage measurement. |
tests/test_vntr_efficiency.py |
Updates expected default penalty factor and related assertions. |
Comments suppressed due to low confidence (1)
docs/guides/vntr-capture-efficiency.md:126
- The log example still says
Coverage ratio: 3.30 (expected ~2.7-3.5)which appears to be the older expected range from the 0.375 calibration. With the updated penalty factor (0.39) and earlier “Expected (penalty 0.39) ~3.6” statement, this expected range should be updated to avoid conflicting guidance.
10. Applying VNTR capture efficiency bias
Downsampling VNTR reads by 0.39...
VNTR efficiency bias applied successfully
VNTR coverage: 177.3x
Non-VNTR coverage: 53.7x
Coverage ratio: 3.30 (expected ~2.7-3.5) ✓
</details>
---
💡 <a href="/berntpopp/MucOneUp/new/main?filename=.github/instructions/*.instructions.md" class="Link--inTextBlock" target="_blank" rel="noopener noreferrer">Add Copilot custom instructions</a> for smarter, more guided reviews. <a href="https://docs.github.com/en/copilot/customizing-copilot/adding-repository-custom-instructions-for-github-copilot" class="Link--inTextBlock" target="_blank" rel="noopener noreferrer">Learn how to get started</a>.
|
|
||
| try: | ||
| result = run_command( | ||
| ["samtools", "depth", "-b", str(region_bed), str(bam_file)], | ||
| ["samtools", "depth", "-a", "-b", str(region_bed), str(bam_file)], | ||
| capture=True, | ||
| ) | ||
|
|
There was a problem hiding this comment.
calculate_mean_coverage() now uses samtools depth -a and captures stdout into memory. Combined with later materializing all depth values into a list, this can become very memory/time intensive for large BEDs (e.g., capture target BEDs), since -a emits one line per base. Consider avoiding capture=True + list materialization (stream and keep running totals), or avoid -a by dividing summed depth by total BED length to account for zero-depth bases.
| ### Statistical Results | ||
|
|
||
| | Metric | Value | 95% CI | | ||
| |--------|-------|--------| | ||
| | **Penalty factor** | **0.375** | [0.357, 0.395] | | ||
| | **Statistical significance** | p < 10⁻³⁰ | Highly significant | | ||
| | **Effect size** | ~3.0x coverage reduction | Strong effect | | ||
| | **Penalty factor** | **0.39** | [0.20, 1.28] | | ||
| | **Cohort size** | 1,043 samples | All Twist Bioscience v2 | |
There was a problem hiding this comment.
The reported 95% CI for the penalty factor includes values > 1.0 ([0.20, 1.28]), but the model/config schema constrain penalty_factor to ≤ 1.0. Either cap/transform the CI to the valid parameter range or clarify that this CI is for an unbounded implied penalty (and describe how it maps to the bounded parameter actually used).
| | Metric | No Bias | With Bias (0.39) | | ||
| |--------|---------|-------------------| | ||
| | VNTR coverage | 200x | 75x | | ||
| | VNTR coverage | 200x | 78x | | ||
| | Flanking coverage | 200x | 200x | | ||
| | Coverage ratio | 1.0 | ~3.0 | | ||
| | VNTR:Flanking ratio | 1.0 | ~2.6 | | ||
| | Reads retained | 100% | ~52% | |
There was a problem hiding this comment.
In the “Expected Behavior” table, the listed “VNTR:Flanking ratio” (~2.6) is inconsistent with the preceding rows: if VNTR coverage drops from 200x to 78x while flanking remains 200x, the ratio would be ~0.39, not >1. Please reconcile the example numbers and ensure the ratio definition matches the earlier “VNTR/Flanking ratio” sections (where ratios are >1).
- Stream depth output with running totals instead of materializing list (avoids memory overhead for large BED regions with -a flag) - Clarify 95% CI: use 5th-95th percentile [0.20, 0.72] for implied penalty; add note explaining outlier values >1.0 - Fix expected behavior table: show realistic VNTR enrichment (~9.3x over flanking) instead of misleading 1:1 ratio Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
Summary
-aflag bug incalculate_mean_coverage()— was excluding zero-coverage positions from samtools depth, inflating reported VNTR coverage stats (VNTR coverage calculation excludes uncovered bases, causing incorrect downsampling #85)Empirical Data
Test plan
Closes #85, closes #86
🤖 Generated with Claude Code