v0.3.0 update to fix read cuts boundary bug and add soft-clip fields

.gitignore

@@ -30,6 +30,7 @@ TODO.md
 *.bam
 /docs/book/
 *.paf.idx
+bedpull_wishlist.md
 
 # Benchmark harness: local config (absolute paths, machine-specific) and
 # generated output/cache are not checked in — only the scripts/templates are.

CHANGELOG.md

@@ -7,6 +7,46 @@ This project uses [Semantic Versioning](https://semver.org/).
 
 ---
 
+## [0.3.0] — 2026-07-21
+
+Coordinate-math correctness pass, driven by a differential test against
+[bladerunner](https://github.com/Psy-Fer/bladerunner) and aimed at letting bladerunner link
+`bedpull` as a library.
+
+### Fixed
+
+- **Boundary-coincidence bug in `get_read_cuts`** — a read whose alignment began (or ended)
+  exactly at a region boundary was mis-sliced or silently dropped. The `start > 0` sentinel
+  conflated "start not found yet" with "start found at read offset 0". Replaced with an explicit
+  `found_start` flag plus an op-entry guard, and the op-level `== ref_end` shortcut no longer
+  skips a mid-op `ref_start` crossing. This affected BAM/CRAM extraction (the PAF path already
+  routed around it via `read_pos_at_ref`).
+
+- **Spurious `missing_left=1bp` on fully-spanning BAM/CRAM reads** — `get_bam_reads` /
+  `get_cram_reads` returned `ref_start`/`ref_end` in the 1-based CIGAR-walk frame while the
+  header suffix compared them against 0-based BED coordinates, so every spanning read was
+  mislabelled as missing 1bp on the left. The returned coordinates are now normalised to
+  0-based (a pre-existing bug, also present in v0.2.0; extraction/sequences were never affected).
+
+### Changed
+
+- **`get_read_cuts` signature** is now
+  `get_read_cuts(cigar_ops, align_start, align_end, region_start, region_end)`. It takes the
+  alignment end and clamps its fire-on boundaries to the alignment span internally
+  (`ref_start = max(region_start, align_start)`, `ref_end = min(region_end, align_end)`), so a
+  partially-overlapping read yields a real slice instead of a `read_end == 0` sentinel.
+  `region_start`/`region_end` are the desired (flank-expanded) window; the caller no longer
+  pre-clamps. **`align_end` is exclusive (one past the last reference base)** — callers using
+  noodles' inclusive `alignment_end()` must add 1.
+- **`ReadCuts` gains `softclip_lead_start` and `softclip_trail_end`** (read offsets of the
+  leading/trailing soft-clip runs; equal to `read_start`/`read_end` when no extension is
+  available). Layout is now field-identical to bladerunner's `ReadCuts`.
+- **noodles bumped `0.110` → `0.111`** to match bladerunner across the crate boundary.
+
+### Added
+
+- Regression battery for the boundary-coincidence cases and the soft-clip extension fields.
+
 ## [0.2.0] — 2026-07-16
 
 ### Added

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "bedpull"
-version = "0.2.0"
+version = "0.3.0"
 edition = "2024"
 authors = ["James M. Ferguson <j.ferguson@garvan.org.au>"]
 repository = "https://github.com/Psy-Fer/bedpull"
@@ -22,7 +22,7 @@ tempfile = "3"
 
 [dependencies]
 clap = { version = "4", features = ["derive", "cargo", "wrap_help"] }
-noodles = { version = "0.110.0", features = ["core", "bgzf", "bam", "sam", "bed", "fasta", "fastq", "cram"] }
+noodles = { version = "0.111.0", features = ["core", "bgzf", "bam", "sam", "bed", "fasta", "fastq", "cram"] }
 itertools = "0.11"
 anyhow = "1.0"
 csv = "1.3"

docs/src/concepts.md

@@ -30,33 +30,39 @@ To include them, you must walk the CIGAR and count read positions and reference
 
 This is exactly what `get_read_cuts` in `src/utils.rs` does.
 
-## The mixed 0-based / 1-based coordinate convention
-
-`get_read_cuts` takes three coordinate arguments with different bases:
-
-- `align_start`: **1-based**
-- `region_start`: **0-based** (from BED)
-- `region_end`: **0-based** (from BED)
-
-This is intentional and non-obvious. Here is why.
-
-### BAM coordinates are 1-based in noodles
-
-The noodles library represents BAM alignment positions as `Position`, which is a 1-based type. When bedpull reads `record.alignment_start()` and converts it with `usize::from(...)`, the result is a 1-based integer. This value is passed directly to `get_read_cuts` as `align_start`.
-
-### BED coordinates are 0-based
-
-The BED format is 0-based, half-open: a BED record with `start=100, end=200` refers to reference bases 100–199 (0-based). bedpull reads BED coordinates directly from the file and passes them to `get_read_cuts` as `region_start` and `region_end`.
-
-### The offset cancels correctly
-
-Inside `get_read_cuts`, the walk begins with `ref_pos = align_start` (1-based). The trigger condition is `ref_pos == region_start` or `ref_pos == region_end` (0-based). These conditions fire when the 1-based running position equals the 0-based BED coordinate. Because BED `start` is 0-based and the alignment start is 1-based, the trigger fires exactly one position earlier in 0-based space than it would if both were 1-based, which is exactly the behaviour needed: BED `start=100` means "the first base of the region is at 0-based position 100", and the 1-based ref_pos of `100` is that same base.
-
-In short: the mixed convention is the natural consequence of BED being 0-based and noodles BAM positions being 1-based, and the arithmetic works out correctly because the two offsets cancel. **Do not normalise `align_start` to 0-based before calling `get_read_cuts`.**
+## Coordinate frames
+
+`get_read_cuts` takes `align_start`, `align_end`, `region_start`, and `region_end`. It is
+**frame-agnostic**: all four must be in the *same* coordinate frame. The walk is relative —
+`ref_pos` is initialised to `align_start` and advances one step per reference-consuming base — so
+only the *difference* between the boundaries affects the returned read offsets. The returned
+`ref_start`/`ref_end` come back in whatever frame you passed in. `align_end` is **exclusive**:
+one position past the last reference base.
+
+### How BAM/CRAM mode keeps the frame consistent
+
+noodles represents BAM/CRAM alignment positions as 1-based `Position` values. `get_bam_reads` /
+`get_cram_reads` pass `alignment_start()` as `align_start` and `alignment_end() + 1` as the
+exclusive `align_end`. The region bounds come from the noodles `Region` built by `read_bed`, which
+stores each BED coordinate shifted by `+1` — a workaround for `Position` being `NonZeroUsize`
+(which cannot represent a BED `start` of `0`). That same `+1` shift puts the region into the same
+1-based frame as the alignment positions, so the walk is internally consistent. Before returning,
+both functions normalise `ref_start`/`ref_end` back to 0-based, so callers — and the
+`|missing_left`/`|missing_right` header suffix — see plain BED coordinates. (Skipping that
+normalisation is what caused the historical `missing_left=1bp`-on-every-read bug.)
+
+Do not put `align_start` and the region bounds in different frames — mixing them shifts every
+boundary by one.
 
 ### PAF coordinates
 
-PAF `target_start` is 0-based. However, bedpull passes it to `get_read_cuts` as `align_start` (the 1-based slot). In practice this means the PAF-mode walk starts at one position "before" the true alignment start relative to what the BAM walk does. The PAF integration (`get_paf_reads` in `src/reads.rs`) was developed and validated against the RFC1 test case, and the extracted lengths are correct (`579 bp = 59 ref bp + 520 inserted bp`). The test in `src/main.rs::tests::get_read_cuts_rfc1_captures_insertion` pins this behaviour at the CIGAR-cut level. Any change to the coordinate handling must keep that test passing.
+PAF-mode extraction does **not** go through `get_read_cuts`. Because a single BED window's two
+edges can fall on two *different* chained PAF records, `get_paf_reads` resolves each edge
+independently with `read_pos_at_ref`, a single-boundary reference→read mapper. PAF `target_start`
+is 0-based and is used consistently with 0-based region bounds throughout that path. It is
+validated against the RFC1 test case (`579 bp = 59 ref bp + 520 inserted bp`):
+`src/main.rs::tests::get_read_cuts_rfc1_captures_insertion` pins the CIGAR-cut math and the
+`get_paf_reads` tests in `src/reads.rs` pin the end-to-end extraction.
 
 Note that the CIGAR-cut math and the final extracted FASTA sequence are two separate steps:
 `get_paf_reads` converts the 0-based half-open cut coordinates into query-contig coordinates and
@@ -68,23 +74,32 @@ silently shifts the extracted sequence by one base without changing the *reporte
 
 Soft-clipped bases (`S`) are present in the read sequence but are not aligned to the reference. Like insertions, they advance the read position without advancing the reference position. `get_read_cuts` handles them the same way as `I` operations, so soft-clipped bases at the start of a read shift the read-position cursor before the first aligned base.
 
+`ReadCuts` also reports `softclip_lead_start` and `softclip_trail_end` — the read offsets of any leading/trailing soft-clip runs. A caller can use these to optionally extend an extraction into clipped "expansion" sequence that sits just outside the aligned window (bladerunner uses this for flanking-expansion detection). When there is no such extension they equal `read_start`/`read_end` respectively.
+
 ## Partial overlaps
 
-When a read's alignment starts after `region_start`, the `region_start` trigger in `get_read_cuts` never fires (because `ref_pos` is initialised to `align_start`, which is already past `region_start`). If `region_end` is subsequently reached, the position is stored in `read_start` (because `start == 0` at that point, indicating the start sentinel has not been set). `read_end` remains `0`.
+`get_read_cuts` clamps its fire-on boundaries to the alignment's own span:
+`ref_start = max(region_start, align_start)` and `ref_end = min(region_end, align_end)`. A read
+that only partially overlaps the window therefore still produces a real slice — `ref_start` /
+`ref_end` report the covered sub-span — instead of a sentinel. `read_end == 0` now means exactly
+one thing: the alignment never reached the window at all.
+
+The start boundary is tracked with an explicit `found_start` flag (not a `start > 0` sentinel) and
+is recorded at op entry when `ref_pos == ref_start`. This is what makes `read_start == 0` a valid
+result: a read whose alignment begins exactly at the window start correctly yields `read_start = 0`
+instead of being dropped (the boundary-coincidence bug fixed in 0.3.0).
 
-`resolve_cuts` in `src/reads.rs` (a private helper shared by `get_bam_reads` and
-`get_cram_reads`) detects this sentinel pattern when `BamConfig::partial` is `true` and remaps
-the fields into the actual `(read_start, read_end, ref_start, ref_end)` covered by the returned
-slice:
+`resolve_cuts` in `src/reads.rs` (a private helper shared by `get_bam_reads` and `get_cram_reads`)
+then applies the spanning policy:
 
-- If `align_start > region_start` (left-partial or contained): the read is extracted from position `0` to `read_cuts.read_start` (where the `region_end` was reached).
-- If `read_cuts.read_end == 0` and partial mode is on (right-partial): the read is extracted from `read_cuts.read_start` to the end of the read.
-- If `read_cuts.read_end == 0` and partial mode is off: the read is skipped.
+- `read_end == 0` (no overlap): the read is skipped.
+- Non-partial mode and the alignment does not fully cover the region (i.e.
+  `align_start <= region_start && align_end >= region_end` is false): the read is skipped.
+- Otherwise: the clamped `(read_start, read_end, ref_start, ref_end)` is returned unchanged.
 
-The corrected `ref_start`/`ref_end` this produces are what drive the `missing_left=Nbp` /
-`missing_right=Nbp` header annotations for `--partial` reads that don't fully cover the
-requested window.
+The resulting `ref_start`/`ref_end` (normalised to 0-based) drive the `missing_left=Nbp` /
+`missing_right=Nbp` header annotations for `--partial` reads that don't fully cover the requested
+window.
 
-The unit tests for `get_read_cuts`'s underlying cut semantics are in `src/utils.rs` under the
-`// --- partial overlap ---` section; the tests for `resolve_cuts`'s field-remapping on top of
-that are in `src/reads.rs` (`resolve_cuts_*`).
+The unit tests for `get_read_cuts`'s cut semantics are in `src/utils.rs`; the tests for
+`resolve_cuts`'s spanning policy are in `src/reads.rs` (`resolve_cuts_*`).

docs/src/library.md

@@ -6,7 +6,7 @@ bedpull is published as a Rust library crate as well as a binary. You can use it
 
 ```toml
 [dependencies]
-bedpull = "0.2"
+bedpull = "0.3"
 ```
 
 ## Import paths
@@ -88,13 +88,22 @@ logic with data from another source.
 use anyhow::Result;
 use bedpull::{ToCigarOps, get_read_cuts};
 
-fn show_cuts(cigar: &str, align_start: usize, region_start: usize, region_end: usize) -> Result<()> {
+fn show_cuts(
+    cigar: &str,
+    align_start: usize,
+    align_end: usize,
+    region_start: usize,
+    region_end: usize,
+) -> Result<()> {
     // Parse a raw CIGAR string from a PAF cg:Z: tag or any str source.
     let ops = cigar.to_cigar_ops()?;
 
-    // align_start is 1-based (as from a BAM record or PAF target_start field).
-    // region_start / region_end are 0-based (as from a BED file).
-    let cuts = get_read_cuts(&ops, align_start, region_start, region_end);
+    // align_start / align_end and region_start / region_end must all be in the same
+    // coordinate frame (the walk is relative). align_end is exclusive — one past the
+    // last reference base. region_start / region_end are the desired window (expand
+    // them yourself if you want flanks); get_read_cuts clamps its fire-on boundaries
+    // to the alignment span internally.
+    let cuts = get_read_cuts(&ops, align_start, align_end, region_start, region_end);
 
     println!(
         "read slice: [{}..{}]  ({} bases)",
@@ -112,12 +121,12 @@ fn show_cuts(cigar: &str, align_start: usize, region_start: usize, region_end: u
 
 fn main() -> Result<()> {
     // A read aligned at ref position 1, with a 5-base insertion inside the region.
-    // CIGAR: 3M5I4M
-    // Region: ref bases 3 to 7 (0-based BED coordinates)
+    // CIGAR: 3M5I4M consumes 7 reference bases, so align_end = 1 + 7 = 8.
+    // Region: ref bases 3 to 7.
     //
     // Expected: the insertion is captured, so the slice is 9 bases
     // (1 match + 5 inserted + 3 match), not 4 reference bases.
-    show_cuts("3M5I4M", 1, 3, 7)?;
+    show_cuts("3M5I4M", 1, 8, 3, 7)?;
     Ok(())
 }
 ```
@@ -138,7 +147,7 @@ cost in PAF-mode extraction once you're processing more than a handful of alignm
 
 ```rust
 use anyhow::Result;
-use bedpull::{PafIndex, get_paf_reads, write_fasta_record};
+use bedpull::{PafIndex, StitchConfig, get_paf_reads, write_fasta_record};
 use noodles::fasta;
 use std::fs::File;
 use std::io::{self, BufReader};
@@ -177,6 +186,7 @@ fn extract_from_paf(
         region_end,
         0, // lflank
         0, // rflank
+        StitchConfig::default(), // no cross-record stitching
         false, // debug
     )?;
 
@@ -197,12 +207,17 @@ and call `get_paf_reads` once per region, exactly what `bedpull`'s own CLI does
 
 ## Notes on the coordinate convention
 
-`get_read_cuts` uses a **mixed coordinate system** that matches the conventions of its two callers:
-
-- `align_start` is **1-based** (noodles `Position` from BAM; PAF `target_start` is 0-based but is used with the CIGAR walk in a way that matches 1-based BAM behaviour; see the [Concepts](concepts.md) chapter for details).
-- `region_start` and `region_end` are **0-based** (directly from BED).
-
-Do not normalise both to the same base before calling `get_read_cuts`. The test suite in `src/utils.rs` documents the exact expected behaviour for all edge cases.
+`get_read_cuts` is **frame-agnostic**: `align_start`, `align_end`, `region_start`, and
+`region_end` must all be in the *same* coordinate frame, because the walk is relative
+(`ref_pos` starts at `align_start` and advances). Only the *difference* between the boundaries
+determines the returned read offsets; the returned `ref_start`/`ref_end` inherit whichever frame
+you passed in. `align_end` is **exclusive** (one past the last reference base) — with noodles'
+1-based inclusive `alignment_end()`, add 1.
+
+In BAM/CRAM mode, `get_bam_reads`/`get_cram_reads` supply all four in noodles' 1-based frame and
+then normalise the returned `ref_start`/`ref_end` back to 0-based (BED) before handing them to you.
+See the [Concepts](concepts.md) chapter for the full derivation. The test suite in `src/utils.rs`
+documents the exact expected behaviour for all edge cases.
 
 `extract_from_fasta_coords`/`extract_from_fasta_coords_reader` take `start`/`end` as **0-based
 half-open**, the same convention as everywhere else in the crate, and convert internally to

src/main.rs

@@ -854,10 +854,33 @@ mod tests {
 
     #[test]
     fn get_read_cuts_rfc1_captures_insertion() {
+        use noodles::sam::alignment::record::cigar::op::Kind;
         let r = read_paf_record_at_offset(PAF_PATH, 0).unwrap();
         let cigar_str = r.cigar.as_deref().expect("no CIGAR");
         let ops = cigar_str.to_cigar_ops().expect("valid CIGAR from PAF file");
-        let cuts = get_read_cuts(&ops, RFC1_TARGET_START, RFC1_REGION_START, RFC1_REGION_END);
+        // align_end = target_start + reference bases consumed (M/=/X/D/N).
+        let ref_len: usize = ops
+            .iter()
+            .filter(|op| {
+                matches!(
+                    op.kind,
+                    Kind::Match
+                        | Kind::SequenceMatch
+                        | Kind::SequenceMismatch
+                        | Kind::Deletion
+                        | Kind::Skip
+                )
+            })
+            .map(|op| op.len)
+            .sum();
+        let align_end = RFC1_TARGET_START + ref_len;
+        let cuts = get_read_cuts(
+            &ops,
+            RFC1_TARGET_START,
+            align_end,
+            RFC1_REGION_START,
+            RFC1_REGION_END,
+        );
         let extracted_len = cuts.read_end - cuts.read_start;
         assert_eq!(
             extracted_len, RFC1_EXPECTED_BP,