Samtools Bam Processing

# Bioconda (recommended)
conda install -c bioconda samtools

# Homebrew (macOS)
brew install samtools

# Verify
samtools --version | head -1

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"

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

# 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

# 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 collate -o collated.bam input.bam

# Quick summary: total, mapped, paired, properly paired
samtools flagstat sorted.bam
# Example output:
# 50000000 + 0 in total (QC-passed reads + QC-failed reads)
# 48523111 + 0 mapped (97.05% : N/A)
# 50000000 + 0 paired in sequencing
# 48490234 + 0 properly paired (96.98% : N/A)

# Per-chromosome mapped/unmapped read counts
samtools idxstats sorted.bam
# chr1  248956422  12345678  0
# chr2  242193529  11234567  0

# Comprehensive stats (insert sizes, GC content, base quality)
samtools stats -r reference.fa sorted.bam > full_stats.txt
grep "^SN" full_stats.txt | cut -f2,3  # Summary Numbers only

# Coverage report (min/max/mean per region/chromosome)
samtools coverage sorted.bam

# Per-base read depth for specific regions
samtools depth -b target_regions.bed sorted.bam > depth.txt
# Output: chr  pos  depth (e.g., chr1  1000  45)

# Statistics split by read group
samtools stats -S RG sorted.bam > per_rg_stats.txt

# FLAG reference — common masks:
# 1    = paired         4  = unmapped
# 2    = proper pair    8  = mate unmapped
# 16   = reverse strand 64 = R1 (first in pair)
# 128  = R2             256= secondary alignment
# 1024 = PCR duplicate  2048= supplementary

# Extract properly paired, mapped reads (FLAG 2 set, 4 unset)
samtools view -f 2 -F 4 sorted.bam -o proper_pairs.bam

# Extract R1 reads only
samtools view -f 64 sorted.bam -o R1.bam

# Remove secondary and supplementary alignments
samtools view -F 2304 sorted.bam -o primary.bam

# Extract reads from BED file regions
samtools view -L regions.bed -b sorted.bam -o regions.bam

# Full duplicate marking workflow (collate → fixmate → sort → markdup)
samtools collate -@ 8 -o collated.bam input.bam
samtools fixmate -m -@ 8 collated.bam fixmated.bam
samtools sort -@ 8 -o sorted.bam fixmated.bam
samtools markdup -@ 8 sorted.bam marked.bam
samtools index marked.bam

# Check duplication rate
samtools flagstat marked.bam | grep "duplicates"
# Output: 2345678 + 0 duplicates (4.83%)

# NovaSeq optical duplicate detection (2500 pixel distance)
samtools markdup -d 2500 sorted.bam marked_novaseq.bam

# Remove duplicates instead of marking
samtools markdup -r sorted.bam deduped.bam

# Get duplication stats without writing output
samtools markdup -s sorted.bam /dev/null

# Merge multiple BAM files (all must be sorted)
samtools merge -@ 8 merged.bam lane1.bam lane2.bam lane3.bam

# Merge files listed in a text file (one per line)
samtools merge -b bam_list.txt -@ 8 merged.bam

# Merge with read group tags from filenames
samtools merge -r merged.bam sample1.bam sample2.bam

# Extract specific chromosome region from merged output
samtools view -h merged.bam chr1 -b -o chr1.bam

#!/bin/bash
SAMPLE="sample_001"
REF="reference.fa"
THREADS=8

# 1. Sort and index (aligner often outputs unsorted SAM/BAM)
samtools sort -@ $THREADS -o ${SAMPLE}.sorted.bam ${SAMPLE}.bam
samtools index ${SAMPLE}.sorted.bam

# 2. QC metrics
samtools flagstat ${SAMPLE}.sorted.bam > ${SAMPLE}.flagstat.txt
samtools stats -r $REF ${SAMPLE}.sorted.bam > ${SAMPLE}.stats.txt
samtools coverage ${SAMPLE}.sorted.bam > ${SAMPLE}.coverage.txt

# 3. Per-chromosome stats
samtools idxstats ${SAMPLE}.sorted.bam > ${SAMPLE}.idxstats.txt

echo "QC complete: $(grep 'mapped (' ${SAMPLE}.flagstat.txt | head -1)"

#!/bin/bash
INPUT="aligned.bam"
FINAL="deduped.bam"
THREADS=8

# Collate → fixmate → sort → markdup
samtools collate -@ $THREADS -o collated.bam $INPUT
samtools fixmate -m -@ $THREADS collated.bam fixmated.bam
samtools sort -@ $THREADS -o sorted.bam fixmated.bam
samtools markdup -@ $THREADS -s sorted.bam $FINAL

# Clean up intermediates
rm collated.bam fixmated.bam sorted.bam

# Index and verify
samtools index $FINAL
samtools flagstat $FINAL | grep "duplic"
# Expected: 3-15% duplicates (WGS); 10-30% for amplicon