BioUnfold.in
All protocols

Extract reads from a genomic region in a BAM

Pull every read overlapping a genomic locus out of an indexed BAM into a smaller region BAM with samtools view, then count reads in the interval.

SSSudipta SardarUpdated July 21, 2026

Prerequisites

  • samtools installed (v1.x, htslib)
  • A coordinate-sorted, indexed BAM (e.g. sorted.bam)
  • Basic familiarity with genomic coordinates (contig:start-end)

1Index the BAM

Region queries need a coordinate index (.bai) sitting next to the BAM so samtools can seek straight to the locus. This writes sorted.bam.bai and prints nothing on success; if it errors, the BAM is almost certainly not coordinate-sorted.

bash
samtools index sorted.bam

Expected output

$ samtools index sorted.bam
$

2Extract the region to a new BAM

The region string chr1:1000000-2000000 selects every read overlapping that interval, and -b keeps the output in BAM (drop it for SAM text). Output is redirected to region.bam, so a clean run prints nothing.

bash
samtools view -b sorted.bam chr1:1000000-2000000 > region.bam

Expected output

$ samtools view -b sorted.bam chr1:1000000-2000000 > region.bam
$ ls -lh region.bam
-rw-r--r--  1 you  staff   3.1M Jul 21 14:32 region.bam

3Count reads in the region (optional check)

To sanity-check the interval before extracting — or when you only need the number rather than the reads themselves — swap -b for -c and samtools prints just the count of overlapping reads.

bash
samtools view -c sorted.bam chr1:1000000-2000000

Expected output

$ samtools view -c sorted.bam chr1:1000000-2000000
42317

Troubleshooting

[main_samview] region "chr1:1000000-2000000" specifies an unknown reference name. Continue anyway.

The contig name in your region string must match the BAM header exactly — many references use "1" rather than "chr1" (or vice versa). List the exact names with `samtools view -H sorted.bam | grep '@SQ'` and copy one into your query.

[E::idx_find_and_load] Could not retrieve index file for 'sorted.bam' — fail to open index.

The BAM has no .bai index, or the file isn't coordinate-sorted (samtools index refuses to index an unsorted BAM and exits with an error rather than producing a usable .bai). Coordinate-sort into a NEW file — never overwrite the input in place — with `samtools sort sorted.bam -o sorted.coord.bam`, then `samtools index sorted.coord.bam` and run your region query against sorted.coord.bam.

References