Fix broken bwa hybrid alignment; resolve hybrid coverage design

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.
 

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("@"):

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()