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 then help(module.function) to check signatures
• CLI: --version then --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-1 = ~536 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 >536 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 avg | 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 ≈ 512 Mbp per contig
samtools index -c file.bam

# Larger min_shift for contigs >512 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 trap: when both .bai and .csi exist, samtools uses .bai. After re-indexing to CSI for a long contig, delete the old .bai or operations fail confusingly.

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

Multiple Regions

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

Regions from BED File

samtools view -L regions.bed input.bam

pysam Python Alternative

Create Index

import pysam

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

Create CSI Index

# 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 index
    for read in bam.fetch('chr1', 1000000, 2000000):
        print(read.query_name)

Check if Indexed

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

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

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

## Source & license

This open-source skill is cataloged on AgentStack and links to its original source — we do not rehost the code.

- **Author:** [GPTomics](https://github.com/GPTomics)
- **Source:** [GPTomics/bioSkills](https://github.com/GPTomics/bioSkills)
- **License:** MIT

Install and usage instructions live in the source repository linked above.