samtools-bam-processing

samtools — SAM/BAM/CRAM 比对工具包

概述

samtools 是用于处理 SAM、BAM 和 CRAM 格式序列比对文件的标准命令行工具包。它涵盖比对文件的完整生命周期：格式转换、坐标排序、创建索引、质量控制统计、读段筛选、重复标记和多文件合并。samtools 几乎是 NGS 工作流中的通用组件，连接比对（STAR、BWA）和下游分析（变异检测、峰值检测、覆盖度）。

适用场景

• 比对后按坐标对 BAM 文件排序（创建索引前必须进行）
• 为已排序的 BAM 文件创建索引，以便随机访问和区域查询
• 在 SAM、BAM 和 CRAM 格式之间转换，以节省存储空间
• 生成比对 QC 指标：比对率、插入片段长度、各染色体统计信息
• 按比对质量、FLAG 位或基因组区域筛选读段
• 在变异检测前标记或移除 PCR 重复
• 合并来自不同测序通道或样本的多个 BAM 文件
• 计算目标区域的逐碱基深度或覆盖广度
• 在自定义脚本中进行原生 Python BAM 操作时，改用 pysam
• 需要标准化 bigWig 覆盖度轨迹时，改用 deeptools bamCoverage
• 需要全基因组逐碱基深度（速度更快、支持并行）时，改用 mosdepth

先决条件

• 安装：建议使用 samtools 1.17+
• 输入要求：SAM/BAM/CRAM 文件；CRAM 需要 FASTA 参考序列
• 配套工具：使用 samtools faidx 为 FASTA 创建索引；使用 samtools index 前先运行 samtools sort

安装前请检查：当前环境中可能已提供此工具（例如在 pixi / conda 环境中）。先运行 command -v samtools，如果返回路径，则跳过下面的安装命令。在 pixi 项目中运行时，请通过 pixi run samtools 调用工具，而不要直接使用 samtools。

# Bioconda (recommended)
conda install -c bioconda samtools

# Homebrew (macOS)
brew install samtools

# Verify
samtools --version | head -1

快速开始

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

核心 API

模块 1：BAM/SAM 输入输出与格式转换

在 SAM/BAM/CRAM 格式之间转换并提取子集。

# 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

模块 2：排序和创建索引

整理 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

模块 3：质量控制和统计

生成比对 QC 指标和覆盖度报告。

# 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

模块 4：读段筛选和 FLAG 操作

使用 SAM FLAG 位筛选特定子集的读段。

# 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

模块 5：重复处理

在变异检测前标记或移除 PCR 重复。

# 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

模块 6：多文件操作和区域分析

合并 BAM 文件并进行区域级分析。

# 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

核心概念

SAM FLAG 位

FLAG 将读段属性编码为位值之和。常见筛选模式：

常见筛选条件	-f（要求）	-F（排除）	选取内容
已比对读段	—	4	所有已比对读段
正确配对	2	—	正确配对且两端均已比对
唯一主比对	—	2308	无次要比对/补充比对/重复
仅 R1	64	—	每对中的第一条读段
未比对	4	—	比对失败

CRAM 与 BAM 与 SAM

格式	大小	速度	要求
SAM	约为 BAM 的 10 倍	I/O 较慢	无
BAM	1 倍	快	随机访问需要 .bai 索引
CRAM	约为 BAM 的 0.6 倍	稍慢	参考 FASTA + 索引

长期存储使用 CRAM；正在分析时使用 BAM。

常见工作流

工作流 1：比对后 QC 和准备

目标：将比对器输出转换为可用于分析的 BAM，并生成 QC 指标。

#!/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)"

工作流 2：完整的重复标记流程

目标：为 GATK 或其他要求输入文件去重的变异检测器准备 BAM。

#!/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

关键参数

参数	命令	默认值	范围/选项	作用
-@	Most	0	1–N cores	额外的压缩/I/O 线程
-m	sort	768M	例如：2G、4G	排序时每个线程使用的内存
-q	view	0	0–60	最低比对质量过滤值
-f	view	0	FLAG bits	包含设置了所有这些位的 reads
-F	view	0	FLAG bits	排除设置了任意这些位的 reads
-b	view	—	flag	输出 BAM 格式
-C	view	—	flag	输出 CRAM（需要 -T）
-T	view	—	FASTA path	CRAM 输出的参考序列
-d	markdup	0	0–2500	光学重复像素距离
-r	markdup	—	flag	移除重复项（而非仅标记）
-n	sort	—	flag	按 read 名称而非位置排序
-c	index	—	flag	创建 CSI 索引（chr > 512 Mb 时需要）

最佳实践

1. 始终先排序再建立索引：samtools index 要求输入按坐标排序。尝试为未排序的 BAM 建立索引会失败或产生错误结果。

2. 所有生产运行都使用 -@：大多数 samtools 命令受 I/O 限制。添加 -@ 8 可在压缩/解压时获得接近线性的加速，且开销极小。

3. 任何分析前先运行 flagstat：samtools flagstat 只需几秒即可运行，并能在下游步骤耗费时间前发现比对失败（比对率低、双端测序比例异常）。

4. 使用 collate → fixmate → sort → markdup 流程：在未运行 fixmate 的情况下，直接对按坐标排序的 BAM 运行 samtools markdup 会导致重复检测结果不正确。fixmate -m 添加的 mate 信息至关重要。

5. 归档时优先使用 CRAM：与 BAM 相比，CRAM 可减少 40-50% 的存储空间，且不会造成损失。请始终将参考 FASTA 与 CRAM 文件一起保存。

6. 针对性分析使用 -L bed_file：将 samtools view 限定在 BED 定义的目标区域（WES 捕获区域、扩增子），可大幅减少下游步骤的 I/O。

常用命令示例

示例：对多个样本批量运行 Flagstat

# Process all BAM files in directory
for bam in *.sorted.bam; do
    echo "=== $bam ==="
    samtools flagstat $bam | grep -E "mapped|properly paired|duplicates"
done

示例：提取未比对 reads 用于从头组装

# Pull both unmapped reads (useful for pathogen detection)
samtools view -f 4 -b input.bam -o unmapped.bam
samtools fastq -@ 4 -1 unmapped_R1.fastq -2 unmapped_R2.fastq unmapped.bam
echo "Unmapped pairs ready for de novo assembly"

示例：将 BAM 下采样至目标覆盖度

# Estimate current depth, then subsample to ~30×
TOTAL=$(samtools flagstat input.bam | grep "mapped (" | head -1 | awk '{print $1}')
GENOME_SIZE=3100000000  # hg38
READ_LEN=150
CURRENT_COV=$(echo "scale=1; $TOTAL * $READ_LEN / $GENOME_SIZE" | bc)
TARGET_FRAC=$(echo "scale=3; 30 / $CURRENT_COV" | bc)
echo "Current: ${CURRENT_COV}×; subsample fraction: $TARGET_FRAC"
samtools view -b -s $TARGET_FRAC input.bam -o downsampled.bam
samtools index downsampled.bam

故障排除

问题	原因	解决方案
[bam_index_build2] fail to index	BAM 未按坐标排序	先排序：samtools sort -o sorted.bam input.bam
BAI index too large for chromosome	染色体 > 512 Mbp	使用 CSI 索引：samtools index -c input.bam
CRAM: reference not found	缺少参考 FASTA 或参考 FASTA 错误	设置 REF_PATH 环境变量或使用 -T ref.fa
重复标记不正确	跳过了 fixmate 步骤	运行完整流程：collate → fixmate → sort → markdup
flagstat 显示 0% 正确配对	双端 BAM 缺少 mate 信息	运行 samtools fixmate 以填充 mate 坐标
排序速度非常慢	每个线程的内存较低	增大 -m 4G；若内存受限则减小 -@
区域查询没有返回结果	BAM 未建立索引或坐标错误	运行 samtools index；使用从 1 开始的坐标：chr1:1000-2000
[E::hts_open_format] fail to open	文件路径错误或 BAM 已损坏	检查路径；使用 samtools quickcheck file.bam 测试

相关技能

• deeptools-ngs-analysis — samtools 下游的标准化 bigWig 覆盖度轨迹和 ChIP-seq 可视化
• pysam-genomic-files — 用于在自定义脚本中处理 BAM 的 Python API
• bedtools-genomic-intervals — 对 samtools 生成的 BAM/BED 文件进行基因组区间操作

参考资料

• samtools 文档 — 官方手册页和命令参考
• GitHub：samtools/samtools — 源码、发布版本和问题跟踪
• Danecek 等人（2021）《SAMtools 和 BCFtools 的十二年》— GigaScience 10(2)
• SAM 格式规范 — FLAG 位、CIGAR 字符串和可选标签