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