Create and use BAI/CSI indices for BAM/CRAM files using samtools and pysam. Use when enabling random access to alignment files or fetching specific genomic regions.
Reference examples tested with: pysam 0.22+, samtools 1.19+
Before using code patterns, verify installed versions match. If versions differ:
pip show <package> then help(module.function) to check signatures<tool> --version then <tool> --help to confirm flagsIf code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.
Create indices for random access to alignment files using samtools and pysam.
"Index a BAM file" → Create a .bai/.csi index enabling random access to genomic regions.
samtools index file.bampysam.index('file.bam')| Index |
|---|
| Extension |
|---|
| Use Case |
|---|
| BAI | .bai | Standard BAM index, chromosomes < 512 Mbp |
| CSI | .csi | Large chromosomes, custom bin sizes |
| CRAI | .crai | CRAM index |
samtools index input.bam
# Creates input.bam.bai
samtools index -c input.bam
# Creates input.bam.csi
samtools index input.bam output.bai
samtools index -@ 4 input.bam
samtools index input.cram
# Creates input.cram.crai
Indexing requires coordinate-sorted files:
# Check sort order
samtools view -H input.bam | grep "^@HD"
# Should show SO:coordinate
# Sort if needed, then index
samtools sort -o sorted.bam input.bam
samtools index sorted.bam
Goal: Extract reads overlapping specific genomic coordinates from an indexed BAM.
Approach: With the index present, samtools view or pysam.fetch() can jump directly to the relevant file offset instead of scanning the entire file.
# Requires index file present
samtools view input.bam chr1:1000000-2000000
samtools view input.bam chr1:1000-2000 chr2:3000-4000
samtools view -L regions.bed input.bam
import pysam
pysam.index('input.bam')
# Creates input.bam.bai
pysam.index('input.bam', 'input.bam.csi', csi=True)
with pysam.AlignmentFile('input.bam', 'rb') as bam:
# fetch() requires index
for read in bam.fetch('chr1', 1000000, 2000000):
print(read.query_name)
import pysam
from pathlib import Path
def is_indexed(bam_path):
bam_path = Path(bam_path)
return (bam_path.with_suffix('.bam.bai').exists() or
Path(str(bam_path) + '.bai').exists() or
bam_path.with_suffix('.bam.csi').exists())
if not is_indexed('input.bam'):
pysam.index('input.bam')
regions = [('chr1', 1000, 2000), ('chr1', 5000, 6000), ('chr2', 1000, 2000)]
with pysam.AlignmentFile('input.bam', 'rb') as bam:
for chrom, start, end in regions:
count = sum(1 for _ in bam.fetch(chrom, start, end))
print(f'{chrom}:{start}-{end}: {count} reads')
with pysam.AlignmentFile('input.bam', 'rb') as bam:
count = bam.count('chr1', 1000000, 2000000)
print(f'Reads in region: {count}')
with pysam.AlignmentFile('input.bam', 'rb') as bam:
for read in bam.fetch('chr1', 1000000, 1000001):
if read.reference_start <= 1000000 < read.reference_end:
print(f'{read.query_name} covers position 1000000')
samtools looks for indices in two locations:
input.bam.bai # Standard location
input.bai # Alternative location
For CRAM:
input.cram.crai
samtools idxstats input.bam
Output format:
chr1 248956422 5000000 0
chr2 242193529 4500000 0
* 0 0 10000
Columns: reference name, length, mapped reads, unmapped reads
samtools idxstats input.bam | awk '{sum += $3} END {print sum}'
with pysam.AlignmentFile('input.bam', 'rb') as bam:
for stat in bam.get_index_statistics():
print(f'{stat.contig}: {stat.mapped} mapped, {stat.unmapped} unmapped')
Related but different - index reference FASTA for random access:
samtools faidx reference.fa
# Creates reference.fa.fai
# Fetch region from indexed FASTA
samtools faidx reference.fa chr1:1000-2000
with pysam.FastaFile('reference.fa') as ref:
seq = ref.fetch('chr1', 1000, 2000)
print(seq)
| Task | samtools | pysam |
|---|---|---|
| Create BAI | samtools index file.bam | pysam.index('file.bam') |
| Create CSI | samtools index -c file.bam | pysam.index('file.bam', csi=True) |
| Fetch region | samtools view file.bam chr1:1-1000 | bam.fetch('chr1', 0, 1000) |
| Count in region | samtools view -c file.bam chr1:1-1000 | bam.count('chr1', 0, 1000) |
| Index stats | samtools idxstats file.bam | bam.get_index_statistics() |
| Index FASTA | samtools faidx ref.fa | Automatic with FastaFile |
| Error | Cause | Solution |
|---|---|---|
random alignment retrieval only works for indexed BAM | Missing index | Run samtools index file.bam |
file is not sorted | Unsorted BAM | Sort first with samtools sort |
chromosome not found | Wrong chromosome name | Check names with samtools view -H |