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.
Install with Codex or Claude Copy this prompt, paste it into Codex, Claude, or another assistant, and let it review the skill page and install it for you.
A direct command skips the review prompt. Inspect the source before running it.
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.
Before using code patterns, verify installed versions match. If versions differ:
Python: pip show <package> then help(module.function) to check signatures
CLI: <tool> --version then <tool> --help to confirm flags
If code throws ImportError, AttributeError, or TypeError, introspect the installed
package and adapt the example to match the actual API rather than retrying.
Alignment Indexing
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.
CLI: samtools index file.bam
Python: pysam.index('file.bam')
Index Types
Index
Extension
Max contig
Bin shift
When required
BAI
.bai / .bam.bai
2^29 bp ≈ 537 Mbp
fixed (16 kb)
Default for human, mouse, fly, fish
CSI
.csi / .bam.csi
2^(min_shift + depth*3)
configurable via -m
Required for any contig >537 Mbp
CRAI
.crai / .cram.crai
chunk-based
n/a
CRAM only
TBI
.tbi
2^29-1
fixed
tabix VCF/BED -- same limit as BAI
Which Index for Which Genome
Genome
Largest contig
Index
GRCh38 / GRCh37 (human)
248 Mbp
BAI
GRCm39 (mouse)
195 Mbp
BAI
GRCz11 (zebrafish), TAIR10 (Arabidopsis)
78 Mbp / 30 Mbp
BAI
Wheat IWGSC (Triticum aestivum)
~830 Mbp (chr3B)
CSI
Pine, fir, axolotl, sugar pine
multi-Gbp
CSI with larger -m
Long-read assembly with very large contigs
varies
check cut -f2 ref.fa.fai | sort -nr | head -1
For polyploid plants and salamander-scale genomes, increase the bin shift:
# Default CSI matches BAI bin layout: 2^(14 + 5*3) = 2^29 bp ≈ 537 Mbp per contig
samtools index -c file.bam
# Larger min_shift for contigs >537 Mbp (wheat, axolotl, sugar pine)
samtools index -c -m 18 file.bam # 2^(18+15) = 2^33 = ~8.5 Gbp per contig
Index file precedence: htslib (and therefore the samtools CLI) tries .csi before .bai on auto-load, so when both exist the .csi is used -- a stale .bai left over after re-indexing to CSI is generally ignored. (Note: htsjdk/Java prefers .bai, the opposite order.) Removing the obsolete .bai still avoids confusion.
samtools index
Create BAI Index
samtools index input.bam
# Creates input.bam.bai
Create CSI Index
samtools index -c input.bam
# Creates input.bam.csi
Specify Output Name
samtools index input.bam output.bai
Multi-threaded Indexing
samtools index -@ 4 input.bam
Index CRAM
samtools index input.cram
# Creates input.cram.crai
Index Requirements
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
Using Indices for Region Access
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.
samtools view with Region
# Requires index file present
samtools view input.bam chr1:1000000-2000000
# pysam.index passes through to samtools index; pass the -c flag for CSI.
pysam.index('-c', 'input.bam')
# Produces input.bam.csi.
Fetch with Index
with pysam.AlignmentFile('input.bam', 'rb') as bam:
# fetch() requires indexfor read in bam.fetch('chr1', 1000000, 2000000):
print(read.query_name)
Check if Indexed
import pysam
from pathlib import Path
defis_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())
ifnot is_indexed('input.bam'):
pysam.index('input.bam')
Fetch Multiple Regions
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(1for _ in bam.fetch(chrom, start, end))
print(f'{chrom}:{start}-{end}: {count} reads')
Count Reads in Region
with pysam.AlignmentFile('input.bam', 'rb') as bam:
count = bam.count('chr1', 1000000, 2000000)
print(f'Reads in region: {count}')
Get Reads Covering Position
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')
Index File Locations
samtools looks for indices in two locations:
input.bam.bai # Standard location
input.bai # Alternative location
The mapped column counts every alignment record with that RNAME, including secondary AND supplementary. For long-read minimap2 output, where a single read can produce many supplementary chimeric alignments, idxstats overcounts input reads -- typically 1.5-3x.
For unique read counts, use primary-only:
samtools view -c -F 2304 input.bam chr1 # primary only
Cross-check unmapped consistency (a senior sanity check):
samtools idxstats file.bam | awk '{sum+=$4} END {print sum}'# idxstats unmapped (sum across all rows; PE orphans get a contig RNAME)
samtools view -c -f 4 -F 2304 file.bam # primary unmapped (should match)
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')
FASTA Index (faidx)
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
pysam FastaFile
with pysam.FastaFile('reference.fa') as ref:
seq = ref.fetch('chr1', 1000, 2000)
print(seq)
Quick Reference
Task
samtools
pysam
Create BAI
samtools index file.bam
pysam.index('file.bam')
Create CSI
samtools index -c file.bam
pysam.index('-c', 'file.bam')
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
Index Staleness
If the BAM was modified after indexing, the index points to wrong file offsets and region queries return wrong (or zero) reads. Quick check:
if [ input.bam -nt input.bam.bai ]; thenecho"Index older than BAM; re-indexing"
samtools index input.bam
fi
Contig-Naming Sanity Check
A leading cause of "my variant calling produced empty VCFs" tickets: querying chrM against a BAM that uses MT (or chr1 vs 1). Always inspect contig conventions before region queries:
samtools view -H input.bam | grep '^@SQ' | head -3
# Compare with reference dict:
samtools dict ref.fa | head -3
UCSC convention uses chr1/chrM; Ensembl/NCBI uses 1/MT. The two are not interchangeable; tools fail with "contig not found" or silently return zero reads.
Common Errors
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
Region query returns zero reads on a known-populated locus
Stale BAI / chr vs no-chr mismatch
Re-index; verify naming convention
BAI silently truncates reads on contigs >537 Mbp
Plant / amphibian / amplified genome
Use CSI: samtools index -c file.bam
Related Skills
sam-bam-basics - View and convert alignment files
alignment-sorting - Sort BAM files (required before indexing)
alignment-filtering - Filter by regions using index
bam-statistics - Use idxstats for quick counts
sequence-io/read-sequences - Index FASTA with SeqIO.index_db()
