/samtools-bam-processing

name: "samtools-bam-processing"
description: "CLI toolkit for SAM/BAM/CRAM: sort, index, convert, filter, QC alignments. Core commands: view, sort, index, flagstat, stats, depth, markdup, merge. Required between alignment and variant/peak calling. Use pysam for Python-native BAM access; deeptools for normalized coverage tracks."
license: "MIT"

samtools — SAM/BAM/CRAM Alignment Toolkit

Overview

samtools is the standard command-line toolkit for processing sequence alignment files in SAM, BAM, and CRAM formats. It handles the complete alignment file lifecycle: format conversion, coordinate sorting, index creation, quality control statistics, read filtering, duplicate marking, and multi-file merging. samtools is a near-universal component of NGS pipelines between alignment (STAR, BWA) and downstream analysis (variant calling, peak calling, coverage).

When to Use

• Sorting BAM files by coordinate after alignment (required before indexing)
• Indexing sorted BAM files for random access and region queries
• Converting between SAM, BAM, and CRAM formats to save storage
• Generating alignment QC metrics: mapping rates, insert sizes, per-chromosome stats
• Filtering reads by mapping quality, FLAG bits, or genomic regions
• Marking or removing PCR duplicates before variant calling
• Merging multiple BAM files from different lanes or samples
• Calculating per-base depth or coverage breadth for target regions
• Use `pysam` instead for Python-native BAM manipulation in custom scripts
• Use `deeptools bamCoverage` instead when you need normalized bigWig coverage tracks
• Use `mosdepth` instead for whole-genome per-base depth (faster, parallelized)

Prerequisites

• **Installation**: samtools 1.17+ recommended
• **Input requirements**: SAM/BAM/CRAM files; CRAM requires FASTA reference
• **Companion tools**: `samtools faidx` for FASTA indexing; `samtools sort` before `samtools index`

> **Check before installing**: The tool may already be available in the current environment (e.g., inside a `pixi` / `conda` env). Run `command -v samtools` first and skip the install commands below if it returns a path. When running inside a pixi project, invoke the tool via `pixi run samtools` rather than bare `samtools`.

# Bioconda (recommended)
conda install -c bioconda samtools

# Homebrew (macOS)
brew install samtools

# Verify
samtools --version | head -1

Pre-flight Interview

Settle these with the user before writing any analysis code.

decisions:
  - id: D1
    param: mappingQualityFloor
    kind: required
    source: user
    ask: "Below what mapping confidence should reads be dropped? Everything counted downstream inherits this cut."
    default: "0 - keep every alignment, including ambiguously placed ones"

  - id: D2
    param: flagFilter
    kind: required
    source: user
    ask: "Which read classes should be excluded - unmapped, secondary, supplementary, failed-QC, duplicates?"
    default: "keep all classes"

  - id: D3
    param: duplicateHandling
    kind: required
    source: user
    ask: "Should duplicates be marked and left in place, or physically removed from the file?"
    default: "marked, not removed - downstream tools can then choose"

  - id: D4
    param: sortOrder
    kind: derived
    source: upstream
    ask: "Does the next step need coordinate order, or name order for mate pairing?"
    default: "coordinate"

  - id: D5
    param: outputFormat
    kind: optional
    source: user
    ask: "Write BAM, or CRAM against a reference to save space?"
    default: "BAM"

  - id: D6
    param: opticalDuplicateDistance
    kind: optional_conditional
    source: data
    ask: "Should optical duplicates be distinguished from PCR duplicates, using this platform's pixel distance?"
    default: "not distinguished"

  - id: D7
    param: threadsAndMemory
    kind: never_ask
    source: data
    reason: "Compression threads and per-thread sort memory affect runtime and peak RAM, not the records"
    default: "available cores, 768M per thread"

D1 and D2 look like plumbing and are not. A MAPQ floor silently removes multi-mapping regions - paralogs, recent duplications, repeat-adjacent genes - from everything computed afterwards, and a run that keeps secondary alignments counts the same fragment several times. Neither is visible in the output file.

Quick Start

# Typical post-alignment workflow: sort → index → QC
samtools sort -@ 8 -o sorted.bam input.bam
samtools index sorted.bam
samtools flagstat sorted.bam

Core API

Module 1: BAM/SAM I/O and Format Conversion

Convert between SAM/BAM/CRAM formats and extract subsets.

# SAM → BAM (saves ~75% disk space)
samtools view -b -h input.sam -o output.bam

# BAM → CRAM (saves additional 40-50%)
samtools view -C -T reference.fa input.bam -o output.cram

# Filter: mapping quality ≥20, exclude unmapped (-F 4)
samtools view -q 20 -F 4 input.bam -o filtered.bam

# Extract specific region (requires index)
samtools view -h sorted.bam "chr1:1000000-2000000" -o region.bam

# Count reads matching filter
samtools view -c -F 4 input.bam
# Output: 45231923 (number of mapped reads)

# Extract reads as FASTQ (for realignment or de novo assembly)
samtools fastq -@ 4 -1 R1.fastq.gz -2 R2.fastq.gz -0 unpaired.fastq.gz input.bam

# Extract reads as FASTA
samtools fasta input.bam > reads.fasta

# Filter by read group
samtools view -r SAMPLE_001 multi_rg.bam -o sample001.bam

Module 2: Sorting and Indexing

Organize BAM files for efficient random access.

# Sort by coordinate (required before indexing)
samtools sort -@ 8 -m 2G input.bam -o sorted.bam

# Sort by read name (required for fixmate/markdup)
samtools sort -n -@ 8 input.bam -o namesorted.bam

# Index sorted BAM (creates sorted.bam.bai)
samtools index sorted.bam

# For chromosomes > 512 Mbp: use CSI index instead
samtools index -c sorted.bam

# Group reads by name (fast, for fixmate — no full sort needed)
samtools col