📚 课程目录第4章预处理流程（一）

预处理流程（一）：BWA-MEM比对与SortSam排序、MarkDuplicates去重

预处理这一步，说白了就是把原始的 FASTQ 文件，一步步加工成 GATK 能直接吃进去的干净 BAM。很多朋友一上来就急着跑 HaplotypeCaller，结果报错一堆，排查半天发现源头在比对这一步就埋了雷。

我个人习惯把预处理拆分三个阶段：比对（BWA-MEM）→ 排序（SortSam）→ 去重（MarkDuplicates）。这三个阶段环环相扣，任何一个环节的参数没调好，后面都得加倍还债。

核心管线流程

原始 FASTQ → BWA-MEM 比对 → SAM/BAM → SortSam 排序 → MarkDuplicates 去重 → 分析就绪 BAM

1. BWA-MEM 比对：不只是跑个默认参数

BWA-MEM 是目前最主流的短读段比对工具，准确率和速度都很均衡。但我见过不少人直接 bwa mem ref.fa read1.fq read2.fq > out.sam 就完事了，嗯，这样 GATK 后面会很难受。

为什么？BWA 默认输出的 SAM 文件里，read group 信息是空的，而 GATK 的很多模块（比如 BaseRecalibrator）强制要求 ReadGroup。我之前有一次跑完整条流水线，结果在 BaseRecalibrator 报错，回头一看，好家伙，比对时的 RG 标签全丢了。

我的推荐命令

# 单样本 paired-end 示例
bwa mem -K 1000000 -M -t 8 \
    -R '@RG\tID:sample1\tLB:lib1\tPL:ILLUMINA\tSM:sample1' \
    reference.fasta \
    sample1_R1.fastq.gz \
    sample1_R2.fastq.gz \
    > sample1.raw.sam

参数说明：

-K 1000000：控制每次处理的碱基数，提高大内存场景下的稳定性
-M：将短拆分 read 标记为 次要（secondary） 而非 supplementary，避免 GATK 后续误判
-t 8：线程数，根据你的 CPU 核心调整，不是越多越好，我一般设 8-16
-R '@RG\tID:...'：ReadGroup 字符串，必须包含 ID、LB、PL、SM 四个字段

关于 -M 参数， 我多聊两句。GATK 的 HaplotypeCaller 对 supplementary alignment 的处理比较敏感。不加 -M 的话，那些跨过断裂点的 read 会被标成 supplementary，导致后续变异检测时被跳过或权重降低。说白了，加 -M 是让 GATK 从数据中借力，而不是被 BWA 的原始标注带沟里去。

BWA-MEM 与 GATK 的兼容性参数

参数	作用	我在项目中踩过的坑
`-K`	控制批处理碱基数	默认值太小会导致内存开销不稳定，尤其在 WGS 数据上
`-M`	标记短拆分 read 为次要	不加的话，GATK 的 ReadBackedPhasing 可能会漏检
`-R`	添加 ReadGroup 信息	漏掉这个字段，BaseRecalibrator 直接罢工
`-t`	线程数	线程数多到超出物理核心时，速度反而下降

曾经踩过的坑：BWA 的 -p 参数

我以前在混合数据时用过 -p 参数（单端模式），结果比对结果里出现了大量异常拆分，GATK 跑 VariantFiltration 时过不去。后来发现，-p 会导致 BWA 把双端当作单端处理，read pair 结构打乱。建议 paired-end 数据一定同时输入 R1 和 R2 文件，不要图省事。

2. SortSam 排序：为什么要排两次？

BWA 输出的 SAM 是按 read name 排序的（或者随机顺序）。但 GATK 所有的分析工具都要求 BAM 按 coordinate（基因组坐标）排序。所以排序这一步不是可选项，是强制项。

你可能会问：直接用 Picard 的 SortSam 还是用 sambamba/samtools 更快？我的建议是：如果在 GATK 生态内，就老老实实走 Picard。虽然 sambamba 跑排序确实快一截，但 GATK 团队对 Picard 输出的 BAM 兼容性验证最充分。

推荐的 SortSam 命令

java -jar picard.jar SortSam \
    I=sample1.raw.sam \
    O=sample1.sorted.bam \
    SORT_ORDER=coordinate \
    MAX_RECORDS_IN_RAM=5000000 \
    TMP_DIR=/path/to/tmp \
    VALIDATION_STRINGENCY=LENIENT

MAX_RECORDS_IN_RAM：控制内存中缓存的记录数，根据你的可用 RAM 调整，避免大量磁盘 I/O
TMP_DIR：指定临时目录，保证磁盘空间充足，我一般设到一个单独的 SSD 盘
VALIDATION_STRINGENCY=LENIENT：对某些非致命格式错误睁一只眼闭一只眼，避免卡死在排序阶段

MAX_RECORDS_IN_RAM 给多少合适？ 我个人经验：如果你的服务器有 64GB RAM，设 5000000 比较保守。如果只有 16GB，我建议降到 2000000，否则 JVM 容易 OOM。记住，这个值不是越大越好——太大反而导致 GC 停顿变长。

3. MarkDuplicates 去重：分子条形码和 PCR 重复

MarkDuplicates 是预处理里的重头戏。PCR 扩增会产生大量重复 read，如果不标记，变异检测时会高估这些区域的覆盖度，导致假阳性 SNP 呼叫。

这里有一个关键点：REMOVE_DUPLICATES 参数。我建议 不要设为 true。GATK 团队官方推荐保留重复标记（用 0x400 flag 标记），而不是直接删除。为什么呢？因为某些高频变异检测（比如低频体细胞突变）需要原始覆盖信息，删除后数据就回不来了。

推荐的 MarkDuplicates 命令

java -jar picard.jar MarkDuplicates \
    I=sample1.sorted.bam \
    O=sample1.markdup.bam \
    M=sample1.markdup_metrics.txt \
    ASSUME_SORTED=true \
    REMOVE_DUPLICATES=false \
    OPTICAL_DUPLICATE_PIXEL_DISTANCE=2500 \
    TMP_DIR=/path/to/tmp

ASSUME_SORTED=true：既然已经用 SortSam 排好序，节省一次检查开销
OPTICAL_DUPLICATE_PIXEL_DISTANCE：光学重复判定距离，Illumina 数据默认 2500，实测够用
M=...：输出 metrics 文件，一定要保留，方便后续统计去重率

双端数据的去重逻辑：MarkDuplicates 默认使用 read pair endpoint 来判定重复，而不是单端。如果两个 read pair 两端比对位置完全一致，就标记为重复。这个逻辑在双端数据上很合理，但如果是单端数据，需要加 READ_NAME_REGEX=null 来关闭基于名字的匹配模式。

曾经踩过的坑：光学重复参数调太高

有一次我在 NovaSeq 数据上跑 MarkDuplicates，把 OPTICAL_DUPLICATE_PIXEL_DISTANCE 设成了 10000，结果去重率飙到 40%。后来一查文档，这个参数是针对不同测序平台的。NovaSeq 的 flow cell 密度更高，光学重复的距离更小，用 2500 就够了。设得过大，会把正常的 non-duplicate read 也误标记为光学重复。

性能优化技巧

MarkDuplicates 的处理速度，在 WGS 数据上是整个管线的瓶颈。我试过一些优化方法：

多线程模式：Picard 2.0+ 支持 --ADD_MATE_CIGAR 并行处理，但注意这需要额外内存
减少内存碎片：设置 -Xmx JVM 堆大小为物理内存的 75%，不要全占满
SSD 临时目录：把 TMP_DIR 放到 SSD 上，I/O 瓶颈明显改善

我的经验数据

用 32 核服务器跑一个 30x WGS 样本，整个过程（BWA + SortSam + MarkDuplicates）大约需要：

BWA-MEM（8 线程）：4-6 小时
SortSam：1-2 小时
MarkDuplicates：2-3 小时

如果你发现 MarkDuplicates 跑了一整天没结束，先检查 MAX_RECORDS_IN_RAM 是否设得太低，导致大量磁盘 spill。

4. 生成 BAM 索引

最后一步，千万别忘了建索引。GATK 所有工具都要求 BAM 文件附带 .bai 索引。我记得有一次需要强制重新比对一个区域，结果发现索引文件丢失，当场重新跑了一遍……

# 用 samtools 建索引
samtools index sample1.markdup.bam

索引文件虽然只多了几十 KB，但作用巨大。没有它，GATK 在读取特定区间时只能暴力遍历整个 BAM，速度慢到令人崩溃。