Introduction

bedpull extracts sequences from BAM, CRAM, or PAF alignment files using BED coordinates. It is designed for situations where standard coordinate-lifting tools (samtools faidx, bedtools getfasta, liftOver) give you the wrong answer because the region of interest contains structural variation relative to the reference.

The problem with reference-only extraction

BED coordinates are defined in reference space. A tool that naively slices a reference FASTA at those coordinates returns only the bases that exist in the reference. A tool that lifts BED intervals to query coordinates using only the alignment start/end positions makes the same mistake at the alignment level: it maps a reference interval to a query interval by applying a simple offset, which is correct for pure matches but wrong whenever the CIGAR string contains insertions.

An insertion in a CIGAR string means the read (or assembly) carries bases that have no reference coordinate at all. Those bases sit between two match blocks in read space, but they are invisible in reference space. Any tool that converts reference coordinates to read coordinates using only the alignment’s start position will produce a slice that skips over the inserted bases entirely.

bedpull avoids this by walking the CIGAR string one operation at a time. It tracks both the reference position and the read position simultaneously as it steps through each operation. When it reaches the start of a BED region it records the read position, continues walking (accumulating read positions through any insertions it encounters), and when it reaches the end of the BED region it records the read position again. The slice seq[read_start..read_end] includes every base (matched, mismatched, and inserted) that the alignment places between those two reference coordinates.

A concrete example: RFC1 VNTR

The RFC1 locus on chromosome 4 contains a variable-number tandem repeat (VNTR) that is relevant to spinocerebellar ataxia. A BED region spanning the annotated RFC1 repeat is 59 bp in reference coordinates (chr4:39318077-39318136 on hg38/hs1). A long-read assembly of a pathological haplotype aligns to this region with a 520 bp insertion inside those 59 reference bases; the VNTR expansion is carried entirely in the insertion block.

Without CIGAR-aware extraction you get 59 bp of flanking sequence and miss the expansion entirely. With bedpull, walking the CIGAR correctly produces a 579 bp extracted sequence (59 reference-matching bases plus 520 inserted bases). This is the actual unit of biological interest, and it cannot be obtained from any tool that works purely in reference coordinates.

When to use bedpull

Use bedpull when:

• Your target regions may contain insertions relative to the reference (VNTRs, STRs, mobile element insertions, structural variants).
• You are working with long-read sequencing data aligned to a reference and need the full read sequence over a locus, not just the reference-matching bases.
• You have an assembly-to-reference PAF alignment and need to extract the corresponding assembly sequence for a set of reference BED regions.

If your regions are simple SNP/indel loci with no large insertions, standard tools will give equivalent results. bedpull adds value specifically when the query sequence diverges structurally from the reference within the BED interval.

Installation

From crates.io

cargo install bedpull

Build from source

git clone https://github.com/Psy-Fer/bedpull
cd bedpull
cargo build --release
./target/release/bedpull --help

Static binary for HPC (musl)

Cluster environments often run older glibc versions that are incompatible with binaries built on a modern workstation. Building with the x86_64-unknown-linux-musl target produces a fully static binary with no shared library dependencies.

rustup target add x86_64-unknown-linux-musl
RUSTFLAGS='-C link-arg=-s' cargo build --release --target x86_64-unknown-linux-musl

The binary is at target/x86_64-unknown-linux-musl/release/bedpull. Copy it to your HPC login node or scratch space; it requires no runtime dependencies.

Requirements

• Rust 1.85 or later (edition 2024; install via rustup).
• A coordinate-sorted BAM for BAM mode, or a CRAM (optionally with --reference for reference-compressed CRAMs) for CRAM mode. If the .bai/.crai index is missing, bedpull builds it automatically the first time you run against that file, so samtools is not required.
• A FASTA for PAF mode query extraction (--query_ref) or CRAM reference decoding (--reference). If the .fai index is missing, bedpull builds it automatically too.
• A PAF file for PAF mode; the PAF must include the cg:Z: CIGAR tag (produced by minimap2 -c or --cs=long).

BAM Mode

BAM mode queries an indexed BAM file for reads that overlap each BED region, walks each read’s CIGAR string to calculate exact read-coordinate boundaries for the region, and extracts the corresponding subsequence (and optionally the quality string).

Basic usage

bedpull \
    --bam reads.bam \
    --bed regions.bed \
    --output extracted.fasta

The BAM must be coordinate-sorted. If the companion .bam.bai index is missing, bedpull builds it automatically the first time you run against that file, so samtools is not required. The output is FASTA by default; use --fastq to include quality scores.

Omit --output (or pass -) to write to stdout:

bedpull --bam reads.bam --bed regions.bed | seqkit stats

FASTQ output

bedpull \
    --bam reads.bam \
    --bed regions.bed \
    --output extracted.fastq \
    --fastq

Filtering reads

Mapping quality

bedpull --bam reads.bam --bed regions.bed --output out.fasta \
    --min_mapq 20

Reads with MAPQ below 20 are skipped before any CIGAR processing. The default is 0 (all reads pass).

Secondary and supplementary alignments

By default only primary alignments are returned. Include secondary or supplementary alignments explicitly:

bedpull --bam reads.bam --bed regions.bed --output out.fasta \
    --include_secondary \
    --include_supplementary

Quality of the extracted region

Filter on the mean Phred score of the extracted slice. This reads the quality scores from the BAM record itself, so it applies regardless of --fastq; the effect is just most visible when you’re keeping the quality string in the output:

bedpull --bam reads.bam --bed regions.bed --output out.fastq \
    --fastq \
    --min_region_quality 20

The mean Phred score is computed in error-probability space (not a simple arithmetic mean of the Phred values), which is the statistically correct approach.

Partial overlaps

By default only reads whose alignment spans the entire BED region are returned (both the left and right boundaries are crossed). To include reads that overlap only part of the region:

bedpull --bam reads.bam --bed regions.bed --output out.fasta \
    --partial

Partial reads are clipped to whatever portion of the region they cover. Reads that start after region_start are extracted from read position 0; reads that end before region_end are extracted to the end of the read. When a read doesn’t fully cover the requested (region ± --flanks) window, its header gets a |missing_left=Nbp and/or |missing_right=Nbp suffix recording how many bases were missed on each side; see Output header format.

Minimum partial coverage

By default --partial accepts a read that overlaps the requested window by even a single base. To require a minimum fraction of coverage instead, similar to liftOver’s -minMatch, add --min_partial_coverage:

bedpull --bam reads.bam --bed regions.bed --output out.fasta \
    --partial \
    --min_partial_coverage 0.9

A read is kept only if the reference span it actually covers (accounting for --flanks, if used) is at least 90% of the requested window; anything less is dropped. --min_partial_coverage requires --partial and must be between 0.0 (the default, meaning no filter) and 1.0.

Flanks

To extend the extraction window beyond the BED coordinates, for example to capture an insertion that sits just outside the annotated boundary, use --flanks, --lflank, or --rflank:

# Symmetric 500 bp flank on both sides
bedpull --bam reads.bam --bed regions.bed --output out.fasta \
    --flanks 500

# Asymmetric: 1000 bp on the left, 200 bp on the right
bedpull --bam reads.bam --bed regions.bed --output out.fasta \
    --lflank 1000 \
    --rflank 200

Per-side values (--lflank, --rflank) take precedence over --flanks for their respective side. Flanks are applied in reference space before the CIGAR walk, so any insertions within the extended window are captured.

Haplotype splitting

If reads carry the HP aux tag (set by whatshap haplotag or similar), you can split the output by haplotype:

bedpull --bam reads.bam --bed regions.bed --output out.fasta \
    --hap_split

This creates three output files alongside the specified --output path:

Deduplication

A read whose alignment spans two BED regions will appear in the output once per region by default, which can be useful when you want region-specific context. To emit each read only once across all regions:

bedpull --bam reads.bam --bed regions.bed --output out.fasta --dedup

Output header format

Each FASTA/FASTQ record has a header of the form:

>read_name|chr:region_start-region_end|region_name

When --hap_split is active and the read carries an HP tag, a haplotype suffix is appended:

>read_name|chr:region_start-region_end|region_name|h1

region_name comes from the 4th column of the BED file, or defaults to chr:start-end if absent.

With --partial, a read that doesn’t fully span the requested window gets an extra suffix:

>read_name|chr:region_start-region_end|region_name|missing_left=12bp|missing_right=8bp

Diagnostics

By default bedpull only prints the mode banner and a completion line. Pass --debug to see verbose per-region diagnostics: the parsed CLI args, each BED record as it’s read, and a banner for every region as it’s processed.

Example: RFC1 VNTR from a long-read BAM

bedpull \
    --bam HG002.GRCh38.bam \
    --bed rfc1_vntr.bed \
    --output rfc1_reads.fasta \
    --min_mapq 20 \
    --flanks 200

Where rfc1_vntr.bed contains:

chr4	39318077	39318136	RFC1_VNTR

CRAM Mode

CRAM mode works identically to BAM mode but reads from an indexed CRAM file. All filtering flags, output formats, haplotype splitting, flanks, and deduplication work the same way. The only difference is the input format and the optional --reference flag required for reference-compressed CRAMs.

Requirements

• A CRAM file. If the <file>.cram.crai index is missing, bedpull builds it automatically the first time you run against that file, so samtools is not required.
• For reference-compressed CRAMs: the reference FASTA used during compression. If its .fai index is missing, it’s auto-built the same way.
• For CRAMs with embedded sequences (created with samtools view -C from a FASTA source with no external reference, or with embed_ref=2): no reference is needed.

Basic usage

bedpull \
    --cram reads.cram \
    --bed regions.bed \
    --output extracted.fasta

Reference-compressed CRAMs

If the CRAM was compressed against an external reference (the default for most pipelines), provide the same reference FASTA so bedpull can decode the sequences:

bedpull \
    --cram reads.cram \
    --reference GRCh38.fa \
    --bed regions.bed \
    --output extracted.fasta

If the reference FASTA’s .fai index is missing, bedpull builds it automatically before extraction.

FASTQ output

bedpull \
    --cram reads.cram \
    --reference GRCh38.fa \
    --bed regions.bed \
    --output extracted.fastq \
    --fastq

Filtering reads

CRAM mode supports the same read filters as BAM mode:

# Mapping quality filter
bedpull --cram reads.cram --bed regions.bed --output out.fasta \
    --min_mapq 20

# Include secondary and supplementary alignments
bedpull --cram reads.cram --bed regions.bed --output out.fasta \
    --include_secondary \
    --include_supplementary

# Mean quality of the extracted region (works regardless of --fastq)
bedpull --cram reads.cram --bed regions.bed --output out.fastq \
    --fastq \
    --min_region_quality 20

# Include partial overlaps
bedpull --cram reads.cram --bed regions.bed --output out.fasta \
    --partial

# Include partial overlaps, but require at least 90% coverage of the requested window
# (see BAM mode's "Minimum partial coverage" section for details)
bedpull --cram reads.cram --bed regions.bed --output out.fasta \
    --partial \
    --min_partial_coverage 0.9

Flanks

bedpull --cram reads.cram --bed regions.bed --output out.fasta \
    --flanks 500

Haplotype splitting

If reads carry the HP aux tag, split output by haplotype:

bedpull \
    --cram reads.cram \
    --bed regions.bed \
    --output out.fasta \
    --hap_split

This creates out.h0.fasta, out.h1.fasta, and out.h2.fasta. --hap_split requires a file --output (not stdout).

Deduplication

bedpull --cram reads.cram --bed regions.bed --output out.fasta --dedup

Output header format

Identical to BAM mode:

>read_name|chr:region_start-region_end|region_name

When --hap_split is active and the read carries an HP tag:

>read_name|chr:region_start-region_end|region_name|h1

With --partial, a read that doesn’t fully span the requested window gets a |missing_left=Nbp/|missing_right=Nbp suffix; see the BAM mode page for details.

Diagnostics

By default bedpull only prints the mode banner and a completion line. Pass --debug to see verbose per-region diagnostics: the parsed CLI args, each BED record as it’s read, and a banner for every region as it’s processed.

stdout

Omit --output (or pass -) to write to stdout:

bedpull --cram reads.cram --reference ref.fa --bed regions.bed | seqkit stats

Converting a BAM to an embedded-reference CRAM

If you do not have the original reference available at runtime, create a self-contained CRAM with embedded sequences:

samtools view -C --no-PG reads.bam -o reads.cram

The samtools index step is optional; bedpull builds the .crai itself on first use if it’s missing.

The resulting CRAM stores the decoded sequences internally; no --reference flag is needed.

PAF Mode

PAF mode extracts sequences from a query FASTA using coordinates derived from a PAF (Pairwise mApping Format) alignment file. The typical use case is assembly-to-reference mapping: you have aligned an assembly to a reference with minimap2, you have a set of reference BED regions, and you want to pull out the corresponding assembly sequence, including any insertions the assembly carries within those regions.

Requirements

• A PAF file with the cg:Z: CIGAR tag. Produce this with minimap2 -c or minimap2 --cs=long.
• The query FASTA (the assembly). If its .fai index is missing, bedpull builds it automatically the first time you run against that file, so samtools is not required.

--partial (see BAM mode) is BAM/CRAM-only and errors if combined with --paf; PAF mode always extracts whatever portion of the region an overlapping alignment covers (see Alignment coverage warnings below).

Basic usage

bedpull \
    --paf assembly_to_ref.paf \
    --query_ref assembly.fasta \
    --bed regions.bed \
    --output extracted.fasta

Output is always FASTA in PAF mode (quality scores are not available from a FASTA). Omit --output (or pass -) to write to stdout.

The PAF index

Scanning a multi-gigabyte PAF file for every BED region would be slow. bedpull builds a byte-offset index the first time you run on a PAF file and saves it alongside the PAF as <paf_file>.idx. On subsequent runs the index is loaded instead of rebuilt.

Index format

The index is a plain TSV with four columns:

<chrom>  <byte_offset>  <target_start>  <target_end>

Each row records where one PAF record lives on disk (byte_offset is the byte position of the start of that line) plus its target chromosome coordinates so that overlap queries can be answered without seeking to the record.

Index lifecycle

If you modify the PAF file, delete the .idx file to force a rebuild.

Disabling the index

bedpull \
    --paf assembly_to_ref.paf \
    --query_ref assembly.fasta \
    --bed regions.bed \
    --output extracted.fasta \
    --use_paf_index false

Flanks

As with BAM mode, --flanks, --lflank, and --rflank expand the extraction window in reference coordinates before the CIGAR walk:

bedpull \
    --paf assembly_to_ref.paf \
    --query_ref assembly.fasta \
    --bed regions.bed \
    --output extracted.fasta \
    --flanks 500

The effective window is clamped to the actual extent of each overlapping alignment, so a flank that extends beyond the alignment boundary does not cause an error; it is silently truncated.

BED liftover output

PAF mode can simultaneously write a BED6 file containing the lifted-over query coordinates for every extracted region. This is useful for the two-step liftover workflow: extract assembly coordinates from the PAF, then use those coordinates to query reads aligned to the assembly.

bedpull \
    --paf assembly_to_ref.paf \
    --query_ref assembly.fasta \
    --bed ref_regions.bed \
    --output extracted.fasta \
    --bed_out assembly_regions.bed

The output BED6 has the form:

query_name  query_start  query_end  bed_name  0  strand

--bed_out can be combined with --output independently. Pass - to write the BED to stdout (but not at the same time as --output -).

Haplotype splitting

If PAF records carry the hp:i: optional tag (set by assemblers or phasing tools that emit phased contigs), you can split the output by haplotype:

bedpull \
    --paf assembly_to_ref.paf \
    --query_ref assembly.fasta \
    --bed regions.bed \
    --output out.fasta \
    --hap_split

This creates three output files:

Deduplication

A query contig that aligns to two BED regions will appear in the output once per region by default. To emit each contig only once across all regions:

bedpull \
    --paf assembly_to_ref.paf \
    --query_ref assembly.fasta \
    --bed regions.bed \
    --output out.fasta \
    --dedup

Output header format

Each FASTA record produced in PAF mode has a header of the form:

>contig_name|chr:ref_start-ref_end|bed_name|contig_name:query_start-query_end|strand

When the alignment carries an hp:i: tag, a haplotype suffix is appended:

>contig_name|chr:ref_start-ref_end|bed_name|contig_name:query_start-query_end|strand|h1

The query coordinates (query_start, query_end) are the 0-based half-open positions within the query contig that were extracted. These are the same values written to --bed_out.

Alignment coverage warnings

If an overlapping alignment does not fully span the requested BED region, bedpull emits a warning but still extracts whatever portion of the region the alignment covers:

Warning: Alignment starts after region start, may be incomplete
Warning: Alignment ends before region end, may be incomplete

These warnings are always printed. A per-alignment Query coords: contig:start-end (strand X) line is also available but only shown with --debug (see below), since it’s noisy for PAF files with many overlapping alignments.

Cross-record stitching

A single PAF record only ever covers a contiguous block of one alignment. A large structural variant, most often a big novel insertion with no homologous target sequence, frequently causes the aligner to emit two or more separate, chained records instead of one record with a large indel operation. Without help, a BED region spanning that boundary is only ever partially extracted from whichever single record happens to overlap it (or produces two separate partial hits, one per record).

--stitch_records looks for a chain of records that share a query contig and strand, are contiguous in target space (within --max_stitch_gap bp), and together span the requested region, then extracts one sequence across the whole chain in a single slice instead of walking each record’s CIGAR separately. Any “gap” in query space between chain members (the inserted sequence itself, which has no aligned counterpart in either record) is spliced in directly from the raw query FASTA, since that gap is the structural variant the aligner split around.

bedpull --paf assembly.paf --query_ref assembly.fa --bed regions.bed \
    --stitch_records --max_stitch_gap 10000 --output output.fasta

--max_stitch_gap (default 10000) bounds only the target-side gap between consecutive chain members. The reconstructed query-side splice (the insertion itself) can be much larger and isn’t separately bounded. Off by default; when a single record already fully covers a region, stitching is skipped even if enabled, since there’s nothing to bridge.

Diagnostics

By default bedpull only prints index-build/load status, the mode banner, and a completion line. Pass --debug to see verbose per-region diagnostics: the parsed CLI args, each BED record as it’s read, a banner for every region as it’s processed, the overlapping-alignment count per region, and per-alignment query coordinates. With --stitch_records, a successful stitch also prints Stitched N chained record(s): contig:start-end (strand X).

Example: extracting RFC1 from a paternal assembly

minimap2 -a -c chm13v2.0.fa HG002_paternal.fasta > hg002pat_to_chm13.paf

bedpull \
    --paf hg002pat_to_chm13.paf \
    --query_ref HG002_paternal.fasta \
    --bed rfc1_vntr.bed \
    --output rfc1_paternal.fasta \
    --bed_out rfc1_paternal_coords.bed

The first run builds hg002pat_to_chm13.paf.idx. The second run loads it in milliseconds. rfc1_paternal.fasta contains the paternal assembly sequence over the RFC1 region, including inserted bases invisible to reference-only tools. rfc1_paternal_coords.bed contains the corresponding query coordinates, ready to use as input to a subsequent BAM-mode run against reads aligned to the paternal assembly.

CLI Reference

Synopsis

bedpull [OPTIONS] --bed <BED>

At least one input mode must be specified: --bam, --cram, or --paf + --query_ref. --bam and --cram are mutually exclusive. --paf cannot be combined with --bam or --cram.

Input flags

Output flags

Filtering flags

Extraction flags

Flank precedence

When both --flanks and a per-side flag are set, the per-side value takes precedence:

effective_lflank = lflank  (if lflank > 0)  else  flanks
effective_rflank = rflank  (if rflank > 0)  else  flanks

Output header formats

BAM and CRAM

>read_name|chr:region_start-region_end|region_name

With --hap_split and a non-zero HP tag:

>read_name|chr:region_start-region_end|region_name|h1

With --partial and a read that doesn’t fully span the requested (region ± flank) window:

>read_name|chr:region_start-region_end|region_name|missing_left=12bp|missing_right=8bp

PAF

>contig_name|chr:ref_start-ref_end|bed_name|contig_name:query_start-query_end|strand

With --hap_split and a non-zero hp:i: tag:

>contig_name|chr:ref_start-ref_end|bed_name|contig_name:query_start-query_end|strand|h1

region_name / bed_name comes from the 4th column of the BED file, or chr:start-end if the column is absent.

Unmapped regions

Pass --unmapped <file> to get a report of every input BED region that produced zero final output rows, across all three modes. The format mirrors liftOver’s own unmapped file: a #reason comment line followed by the input BED record (chrom, start, end, name):

#no overlapping reads found
chr4	1000	2000	NOMATCH
#44 candidate read(s) found but all were filtered out (--min_mapq/--include_secondary/--include_supplementary/--partial/--min_partial_coverage/--min_region_quality)
chr4	39318077	39318136	RFC1
#40 matching read(s) found but all were already emitted for another region (--dedup)
chr4	39318077	39318136	RFC1_b
#region skipped (chromosome name contains '#')
chr4	1000	2000	COMMENTED

Reasons distinguish four cases:

• No overlapping reads/alignments at all: nothing in the BAM/CRAM/PAF overlapped the region.
• Candidates existed but all were filtered out: one or more of --min_mapq, --include_secondary/--include_supplementary, --partial/--min_partial_coverage, or --min_region_quality excluded every candidate. The exact count of raw candidates is included, but not which specific filter caused each exclusion.
• All matches already emitted for another region: only possible with --dedup; the reads or alignments existed and passed every filter, but were all skipped because they’d already been written out for an earlier region.
• Region explicitly skipped: the region’s chromosome name contains #.

Without --unmapped, bedpull still prints an equivalent one-line notice to stderr for each empty region, but doesn’t distinguish these cases or write a structured file.

Exit codes

Diagnostics

bedpull writes progress and diagnostic messages to stderr so that stdout (or the output file) contains only sequence data.

By default this is minimal: a mode banner, index-build/load status (PAF mode), and a completion line, plus warnings that indicate a real data issue (an unexpected HP tag value, a PAF alignment that doesn’t fully cover a region, a region skipped because it contains #, or a region with no overlapping reads/alignments). These warnings are always shown, regardless of --debug.

Pass --debug to additionally see: the fully parsed Opts struct, each BED record as it’s parsed, a banner for every region as it’s processed, and (PAF mode) the overlapping-alignment count per region plus per-alignment query coordinates.

Input validation

Before extraction starts, bedpull validates the input files and argument combinations and fails fast with an actionable message rather than partway through processing. This includes: file existence, output/--bed_out directory existence and writability, mutually exclusive modes (--bam+--cram, alignment modes+--paf), --paf/--query_ref must be used together, --fastq/--min_region_quality/--min_mapq/--include_secondary/--include_supplementary/--partial requiring --bam/--cram, --min_partial_coverage requiring --partial and being within [0.0, 1.0], and --hap_split requiring a file --output (not stdout).

Library Usage

bedpull is published as a Rust library crate as well as a binary. You can use it programmatically to integrate CIGAR-aware extraction into your own pipelines.

Adding to Cargo.toml

[dependencies]
bedpull = "0.3"

Import paths

All public items are re-exported from the crate root:

#![allow(unused)]
fn main() {
use bedpull::{
    // CIGAR types
    CigarOp, CigarOps, ToCigarOps,
    // PAF index and record access
    PafIndex, PafIndexEntry, PafRecord, read_paf_record_at_offset, read_paf_record_from_reader,
    // BAM/CRAM/PAF extraction
    BamConfig, BamRead, PafRead, StitchConfig, get_bam_reads, get_cram_reads, get_paf_reads,
    // Core utilities
    ReadCuts, get_read_cuts, calculate_qscore, revcomp,
    read_bed, write_fasta_record, write_fastq_record,
    extract_from_fasta_coords, extract_from_fasta_coords_reader,
};
}

read_paf_record_at_offset and extract_from_fasta_coords each open their file fresh on every call, which is fine for a one-off lookup. read_paf_record_from_reader and extract_from_fasta_coords_reader take an already-open reader instead, so you can open the PAF/FASTA once and reuse it across many records; see Example 3.

Example 1: BAM extraction with BamConfig

#![allow(unused)]
fn main() {
use anyhow::Result;
use bedpull::{BamConfig, get_bam_reads};
use noodles::bam;
use noodles::core::Region;

fn extract_region(bam_path: &str, region_str: &str) -> Result<()> {
    let config = BamConfig {
        min_mapq: 20,
        include_secondary: false,
        include_supplementary: false,
        partial: false,
        min_region_quality: 0.0,
    };

    let mut reader = bam::io::indexed_reader::Builder::default()
        .build_from_path(bam_path)?;
    let header = reader.read_header()?;

    let region: Region = region_str.parse()?;
    let query = reader.query(&header, &region)?;

    // No flanks; extract exactly the BED region.
    let reads = get_bam_reads(&config, query, &region, 0, 0)?;

    for (name, seq, _qual, ref_start, ref_end, hap) in reads {
        let seq_str = String::from_utf8(seq)?;
        println!(
            ">{} ref={}:{}-{} hp={}",
            name, region_str, ref_start, ref_end, hap
        );
        println!("{}", seq_str);
    }

    Ok(())
}
}

get_bam_reads returns a Vec<BamRead>, where each element is a tuple (name, sequence, quality_string, ref_start, ref_end, haplotype). The sequence field contains only the bases that fall within the requested region after the CIGAR walk: inserted bases are included, deleted bases are absent.

Example 2: CIGAR coordinate math with get_read_cuts

This example shows how to use get_read_cuts directly to calculate the read slice for a known CIGAR string. This is useful if you are integrating bedpull’s coordinate logic with data from another source.

use anyhow::Result;
use bedpull::{ToCigarOps, get_read_cuts};

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 / 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)",
        cuts.read_start,
        cuts.read_end,
        cuts.read_end - cuts.read_start,
    );
    println!(
        "reference covered: [{}-{}]",
        cuts.ref_start, cuts.ref_end,
    );

    Ok(())
}

fn main() -> Result<()> {
    // A read aligned at ref position 1, with a 5-base insertion inside the region.
    // 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, 8, 3, 7)?;
    Ok(())
}

Output:

read slice: [2..11]  (9 bases)
reference covered: [3-7]

Example 3: PAF-based extraction

get_paf_reads is the high-level entry point: it drives the CIGAR walk, the FASTA extraction, and the minus-strand reverse-complement for you. It takes an already-open PAF reader and an already-open indexed FASTA reader so you can call it once per region (or even across many regions) without reopening either file each time. Reopening per record is the single biggest cost in PAF-mode extraction once you’re processing more than a handful of alignments.

#![allow(unused)]
fn main() {
use anyhow::Result;
use bedpull::{PafIndex, StitchConfig, get_paf_reads, write_fasta_record};
use noodles::fasta;
use std::fs::File;
use std::io::{self, BufReader};

fn extract_from_paf(
    paf_path: &str,
    query_fasta: &str,
    chrom: &str,
    region_start: usize,
    region_end: usize,
) -> Result<()> {
    // Build or load the byte-offset index.
    let index_path = format!("{}.idx", paf_path);
    let index = if std::path::Path::new(&index_path).exists() {
        PafIndex::load(&index_path)?
    } else {
        let idx = PafIndex::build(paf_path)?;
        idx.save(&index_path)?;
        idx
    };

    // Open both files once; reuse them across every region you process.
    let mut paf_reader = BufReader::new(File::open(paf_path)?);
    let mut fasta_reader = fasta::io::indexed_reader::Builder::default()
        .build_from_path(query_fasta)?;

    let entries = index.query(chrom, region_start, region_end);
    let stdout = io::stdout();
    let mut out = stdout.lock();

    let reads = get_paf_reads(
        &mut paf_reader,
        &mut fasta_reader,
        &entries,
        region_start,
        region_end,
        0, // lflank
        0, // rflank
        StitchConfig::default(), // no cross-record stitching
        false, // debug
    )?;

    for (sequence, query_name, query_start, query_end, strand, _hap) in reads {
        let header = format!(
            "{}|{}:{}-{}|{}:{}-{}|{}",
            query_name, chrom, region_start, region_end, query_name, query_start, query_end, strand
        );
        write_fasta_record(&mut out, &header, &sequence)?;
    }

    Ok(())
}
}

To process many regions, open paf_reader and fasta_reader once outside a loop over regions and call get_paf_reads once per region, exactly what bedpull’s own CLI does internally.

Notes on the coordinate convention

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 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 noodles’ 1-based-inclusive region syntax. Pass read_cuts/get_paf_reads output straight through without adjusting it yourself.

Concepts

This chapter explains the technical foundations that underpin bedpull’s extraction logic. It is aimed at readers who want to understand why the tool behaves the way it does, or who are integrating the library into their own code.

CIGAR strings

A CIGAR string encodes how an aligned sequence maps to a reference. It is a series of (length, operation) pairs. The operations relevant to coordinate extraction are:

The key asymmetry: insertions advance the read position without advancing the reference position. Inserted bases exist in the read (or assembly) but have no corresponding reference coordinate. This is why reference-coordinate slicing misses them.

Why insertions have no reference coordinate

Reference coordinates number the bases of the reference sequence. An insertion is a sequence of bases present in the read that does not correspond to any reference base. In the CIGAR string, an I operation of length n means: “the next n bases in the read have no reference counterpart.”

If you convert a BED interval [start, end) to read coordinates by taking align_start + (start - align_start) and align_start + (end - align_start), which is what a simple offset-based tool does, you are only counting reference positions. You skip over inserted bases because they are not counted in the reference arithmetic.

To include them, you must walk the CIGAR and count read positions and reference positions simultaneously. When you encounter an I operation, you increment only the read position counter. The read position counter therefore “runs ahead” of the reference position counter, and the difference between them at any point is the total number of inserted bases seen so far. When you eventually reach region_end in reference coordinates, the read position cursor is already past all insertions that occurred between region_start and region_end, and they are included in the slice.

This is exactly what get_read_cuts in src/utils.rs does.

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-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 extracts via extract_from_fasta_coords_reader, which itself converts them into noodles’ 1-based-inclusive region syntax internally. Both conversions matter: a bug in either one silently shifts the extracted sequence by one base without changing the reported length.

Soft clips

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

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) then applies the spanning policy:

• 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 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 cut semantics are in src/utils.rs; the tests for resolve_cuts’s spanning policy are in src/reads.rs (resolve_cuts_*).