Bio Alignment Indexing

Version Compatibility

Reference examples tested with: pysam 0.22+, samtools 1.19+

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

Which Index for Which Genome

For polyploid plants and salamander-scale genomes, increase the bin shift:

bashCopy code
# 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

bashCopy code
samtools index input.bam
# Creates input.bam.bai

Create CSI Index

bashCopy code
samtools index -c input.bam
# Creates input.bam.csi

Specify Output Name

bashCopy code
samtools index input.bam output.bai

Multi-threaded Indexing

bashCopy code
samtools index -@ 4 input.bam

Index CRAM

bashCopy code
samtools index input.cram
# Creates input.cram.crai

Index Requirements

Indexing requires coordinate-sorted files:

bashCopy code
# 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

bashCopy code
# Requires index file present
samtools view input.bam chr1:1000000-2000000

Multiple Regions

bashCopy code
samtools view input.bam chr1:1000-2000 chr2:3000-4000

Regions from BED File

bashCopy code
samtools view -L regions.bed input.bam

pysam Python Alternative

Create Index

pythonCopy code
import pysam

pysam.index('input.bam')
# Creates input.bam.bai

Create CSI Index

pythonCopy code
# 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

pythonCopy code
with pysam.AlignmentFile('input.bam', 'rb') as bam:
    # fetch() requires index
    for read in bam.fetch('chr1', 1000000, 2000000):
        print(read.query_name)

Check if Indexed

pythonCopy code
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')

Fetch Multiple Regions

pythonCopy code
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')

Count Reads in Region

pythonCopy code
with pysam.AlignmentFile('input.bam', 'rb') as bam:
    count = bam.count('chr1', 1000000, 2000000)
    print(f'Reads in region: {count}')

Get Reads Covering Position

pythonCopy code
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:

Copy code
input.bam.bai   # Standard location
input.bai       # Alternative location

For CRAM:

Copy code
input.cram.crai

idxstats - Index Statistics

Get Per-Chromosome Counts

bashCopy code
samtools idxstats input.bam

Output format:

Copy code
chr1    248956422    5000000    0
chr2    242193529    4500000    0
*       0            0          10000

Columns: reference name, length, mapped reads, unmapped reads.

What "mapped" Actually Counts (Caveat)

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:

bashCopy code
samtools view -c -F 2304 input.bam chr1   # primary only

Cross-check unmapped consistency (a senior sanity check):

bashCopy code
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)

Sum Total Mapped Reads

bashCopy code
samtools idxstats input.bam | awk '{sum += $3} END {print sum}'

pysam idxstats

pythonCopy code
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:

bashCopy code
samtools faidx reference.fa
# Creates reference.fa.fai

# Fetch region from indexed FASTA
samtools faidx reference.fa chr1:1000-2000

pysam FastaFile

pythonCopy code
with pysam.FastaFile('reference.fa') as ref:
    seq = ref.fetch('chr1', 1000, 2000)
    print(seq)

Quick Reference

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:

bashCopy code
if [ input.bam -nt input.bam.bai ]; then
    echo "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:

bashCopy code
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

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()