From c60a92a7d52608e351dca01909fe1d4c0bb06f46 Mon Sep 17 00:00:00 2001 From: Daniel Sprockett Date: Wed, 10 Jun 2026 23:55:27 -0400 Subject: [PATCH] Fix broken bwa hybrid alignment path; align hybrid docs with implementation MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The bwa hybrid path was broken: _merge_long_reads_bwa ran minimap2 through the _run helper, which captures and discards stdout, so the long-read SAM was never written and the subsequent open(long.sam) raised FileNotFoundError. Hybrid short+long samples under a bwa aligner therefore always failed. - Route the long-read alignment through _align_minimap2_preset, which redirects minimap2's SAM to long.sam, then append its records (header skipped) to the short SAM — producing one merged BAM per sample as intended. Dropped the unused aligner_bin/samtools_bin params; minimap2 binary is now a parameter. - Resolve the hybrid coverage design question: keep the merged-BAM total-depth approach (one BAM per sample, CoverM once; relative cross-sample signal preserved because every sample is processed identically). Reject the previously-documented read-count-weighted-mean merge, which needs separate per-technology CoverM passes for no benefit to relative coverage. - CLAUDE.md corrected: merging is at the alignment stage, not "in coverage.py weighted by read count"; the rejected alternative and the one real reason to revisit it (per-technology identity filtering) are noted. Tests: test_align_hybrid.py covers the merge (long records appended, header not duplicated) and is a regression guard for the discarded-stdout bug. 136 pass. --- CLAUDE.md | 4 +-- src/refrover/align.py | 25 ++++++++-------- tests/test_align_hybrid.py | 58 ++++++++++++++++++++++++++++++++++++++ 3 files changed, 73 insertions(+), 14 deletions(-) create mode 100644 tests/test_align_hybrid.py diff --git a/CLAUDE.md b/CLAUDE.md index 5cfc22e..325c44f 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -137,7 +137,7 @@ Rules: - `assembly` is required for every row - At least one of `r1` or `long_reads` required - `r2` requires `r1`; absent `r2` = single-end short reads -- Presence of both `r1` and `long_reads` = hybrid sample (aligns both, merges depth) +- Presence of both `r1` and `long_reads` = hybrid sample (aligns both technologies into one BAM; total depth — see Key design decisions) ### Prototype assignment output (written by `refrover select`) @@ -276,7 +276,7 @@ pytest tests/ # unit tests (all should pass without external tools **Adaptive k (not yet implemented)**: The current selectors all require an explicit `--k`. The intended design: for each query, compute the marginal gain in predicted differential signal as k increases, and stop at the elbow. This requires a signal prediction model — deferred to after the fixed-k selectors are benchmarked. -**Coverage merging for hybrid samples**: For samples with both short and long reads, align separately, then take the mean depth per contig weighted by read count. This is handled in `coverage.py`. +**Hybrid samples (short + long reads)**: Merging happens at the **alignment** stage, not in `coverage.py`. Short reads are aligned with bwa-mem2/bwa (or minimap2 `sr`) and long reads with minimap2 `map-ont`; both alignment sets are written into a single sorted BAM per sample (see `align._merge_long_reads_bwa` and the hybrid branch of `align._align_minimap2`). CoverM then runs once per sample, so a contig's reported depth is the total depth across both technologies. This keeps one BAM per sample and a clean sample→BAM mapping, and because every sample is processed identically the relative cross-sample signal that differential binning depends on is preserved. We deliberately do **not** compute per-technology depths and combine them by a read-count-weighted mean: that needs separate CoverM passes and a sample→(short_bam, long_bam) mapping for no benefit to *relative* coverage. The one real advantage of separate passes — applying technology-appropriate read-identity filters (short reads map at higher identity than ONT) — is noted as future work. **Idempotent stages**: Each stage writes its output to a predictable path under `--outdir`. Re-running a stage skips existing outputs (file-existence check). The `--force` flag disables this. diff --git a/src/refrover/align.py b/src/refrover/align.py index 71c2ecb..d8ab246 100644 --- a/src/refrover/align.py +++ b/src/refrover/align.py @@ -161,9 +161,7 @@ def _align_bwa( raise RuntimeError(f"bwa mem failed:\n{result.stderr[:500]}") if pd.notna(long_reads): - _merge_long_reads_bwa( - sample_row, ref_fa, sam_path, work_dir, aligner_bin, samtools_bin, threads - ) + _merge_long_reads_bwa(sample_row, ref_fa, sam_path, work_dir, threads) elif pd.notna(long_reads): # Long-reads only: fall through to minimap2 @@ -180,21 +178,24 @@ def _merge_long_reads_bwa( ref_fa: Path, short_sam: Path, work_dir: Path, - aligner_bin: str, - samtools_bin: str, threads: int, + minimap2_bin: str = "minimap2", ) -> None: """ - For hybrid samples: align long reads with minimap2, merge with short-read SAM. - Modifies short_sam in place (appends long-read alignments). + For hybrid samples under a bwa aligner: align the long reads with minimap2 + (bwa cannot map long reads) and append them to the short-read SAM, producing + one merged alignment per sample. Modifies short_sam in place. """ long_sam = work_dir / "long.sam" - _run( - ["minimap2", "-ax", "map-ont", "-t", str(threads), - str(ref_fa), str(sample_row["long_reads"])], - "minimap2 long reads", + # _align_minimap2_preset redirects minimap2's SAM to long_sam; _run does not + # (it captures stdout), which is why the long alignments must go through the + # preset helper here. + _align_minimap2_preset( + minimap2_bin, ref_fa, [str(sample_row["long_reads"])], long_sam, + preset="map-ont", threads=threads, ) - # Append long-read alignments (skip header lines that start with @) + # Append long-read alignments, skipping the long SAM's header (the short + # SAM already has a header for the same reference). with open(short_sam, "a") as fout, open(long_sam) as fin: for line in fin: if not line.startswith("@"): diff --git a/tests/test_align_hybrid.py b/tests/test_align_hybrid.py new file mode 100644 index 0000000..965a0d9 --- /dev/null +++ b/tests/test_align_hybrid.py @@ -0,0 +1,58 @@ +""" +Tests for hybrid (short + long read) alignment merging. + +Hybrid samples align short reads with bwa/minimap2 and long reads with minimap2, +then merge both into one SAM per sample. The long-read alignment must actually be +written to disk and appended (a prior bug discarded it via the capturing _run +helper, leaving the merge to fail on a missing file). +""" + +from pathlib import Path + +import pandas as pd + +import refrover.align as align_mod +from refrover.align import _merge_long_reads_bwa + + +def test_merge_long_reads_appends_long_alignments(tmp_path, monkeypatch): + short_sam = tmp_path / "aln.sam" + short_sam.write_text( + "@HD\tVN:1.6\n" + "@SQ\tSN:c1\tLN:100\n" + "shortread1\t0\tc1\t1\t60\t4M\t*\t0\t0\tACGT\tIIII\n" + ) + + def fake_preset(minimap2_bin, ref_fa, reads, out_sam, preset, threads): + # Stand in for minimap2: write a long-read SAM with its own header. + assert preset == "map-ont" + Path(out_sam).write_text( + "@HD\tVN:1.6\n" + "@SQ\tSN:c1\tLN:100\n" + "longread1\t0\tc1\t1\t60\t100M\t*\t0\t0\t*\t*\n" + ) + + monkeypatch.setattr(align_mod, "_align_minimap2_preset", fake_preset) + + row = pd.Series({"sample_id": "S1", "long_reads": str(tmp_path / "S1_ont.fq")}) + _merge_long_reads_bwa(row, tmp_path / "ref.fa", short_sam, tmp_path, threads=4) + + content = short_sam.read_text() + assert "shortread1" in content # original short alignment kept + assert "longread1" in content # long alignment appended (the fix) + assert content.count("@HD") == 1 # long header not duplicated into body + + +def test_merge_long_reads_does_not_raise_on_valid_inputs(tmp_path, monkeypatch): + """Regression: the old _run-based path left long.sam unwritten and raised.""" + short_sam = tmp_path / "aln.sam" + short_sam.write_text("@HD\tVN:1.6\nr1\t0\tc1\t1\t60\t4M\t*\t0\t0\tACGT\tIIII\n") + + monkeypatch.setattr( + align_mod, "_align_minimap2_preset", + lambda *a, **k: Path(a[3]).write_text("@HD\tVN:1.6\nlr\t0\tc1\t1\t60\t9M\t*\t0\t0\t*\t*\n"), + ) + + row = pd.Series({"sample_id": "S1", "long_reads": str(tmp_path / "x.fq")}) + _merge_long_reads_bwa(row, tmp_path / "ref.fa", short_sam, tmp_path, threads=1) + assert "lr" in short_sam.read_text()