第2章：数据预处理流程（Fastq → BAM）
这一章，咱们聊聊最基础但也最关键的环节——从原始测序数据到比对后的BAM文件。很多人急着跑GATK，结果在数据预处理上栽了跟头。我见过太多案例，最后发现问题出在源头。
说白了，这一步你做得有多扎实，直接决定了后面变异检测的准确性。
核心原则：数据预处理是GATK最佳实践的第一道关。它的目标很明确——把混乱的原始数据整理成GATK能理解的标准格式。这个过程容不得半点马虎。
2.1 Fastq质检——拒接坏数据
拿到测序公司给的数据，我习惯先用FastQC扫一遍。你看似可以跳过这步，但万一遇到接头污染或者质量波动，后面全白干。
我的经验是：不要只看总数据量。
避坑指南：我曾经遇到过一个项目，数据总量看着漂亮，结果跑完GATK才发现绝大部分序列都在前50bp有严重的碱基质量衰减。原因就是没做FastQC。后来重做，白烧了三周时间。
运行FastQC很简单：
fastqc sample_R1.fastq.gz sample_R2.fastq.gz -o ./qc_report/
重点看这几项指标：
| 指标 | 含义 | 合格标准 |
|---|---|---|
| Per base sequence quality | 每个位点的平均质量 | Q30占比 > 80% |
| Per sequence GC content | GC含量分布 | 与实际基因组GC%偏差 < 10% |
| Adapter Content | 接头污染比例 | < 5% |
| Sequence Duplication Levels | 重复序列比例 | PCR重复常偏高，但异常高需注意 |
小技巧：如果发现接头污染，别慌。用Trimmomatic或cutadapt切掉就行。我个人觉得cutadapt更轻量，但Trimmomatic支持滑动窗口裁切，适合低质量碱基集中在前段的情况。
# 裁切低质量和接头
java -jar trimmomatic-0.39.jar PE \
sample_R1.fastq.gz sample_R2.fastq.gz \
sample_R1_trimmed.fastq.gz sample_R1_unpaired.fastq.gz \
sample_R2_trimmed.fastq.gz sample_R2_unpaired.fastq.gz \
ILLUMINACLIP:adapters.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36
2.2 BWA-MEM比对——把reads送回基因组
这一步说白了就是把测序片段“贴”回参考基因组。目前主流是BWA-MEM算法，它的速度和准确度都够用。
我的习惯是先建索引：
# 为参考基因组建立索引（含bwa索引和fasta索引）
bwa index -a bwtsw hg38.fa
samtools faidx hg38.fa
然后执行比对：
bwa mem -t 8 -M -R "@RG\tID:sample1\tLB:lib1\tSM:sample1\tPL:illumina" \
hg38.fa sample_R1_trimmed.fastq.gz sample_R2_trimmed.fastq.gz \
> sample.sam
必须注意：Read Group（RG）标签不能省略。@RG中的ID、LB、SM、PL四个字段必须填写完整。GATK后续依赖这些信息进行样本分离和bias校正。你不填，GATK会报错，甚至拒绝对BAM进行处理。
还有个参数容易忽略：我用 -M 把多比对reads标记为“secondary”，这样MarkDuplicates才能正确处理它们。你想想看，如果不标记，下游去重就会漏掉情况。
2.3 排序与去重——整理战场
比对出来的SAM文件是乱的——按reads的顺序而不是按基因组坐标。GATK要求BAM是排序好的，所以必须做SortSam。
# SAM → 排序BAM
java -jar picard.jar SortSam \
I=sample.sam \
O=sample_sorted.bam \
SORT_ORDER=coordinate \
CREATE_INDEX=true
排序完成，接下来去重。为什么要去重？
说白了，PCR扩增会产生很多完全一模一样的reads，这些不是真实的基因组片段，而是技术噪声。如果不去重，变异等位基因频率会被抬高，假阳性就来了。
# 标记并去除PCR重复
java -jar picard.jar MarkDuplicates \
I=sample_sorted.bam \
O=sample_dedup.bam \
M=marked_dup_metrics.txt \
REMOVE_DUPLICATES=true \
CREATE_INDEX=true
避坑指南：我曾经做过一个全外显子组项目，样本质量本身不太好，PCR重复率飙到30%。我当时没注意metrics里的数字，直接跑GATK，结果一堆假阳性。后来回头一看，好家伙，30%的reads都是无用数据。所以，不看metrics就去重等于白做。
去重后的metrics文件里，重点看：
- LIBRARY：文库的统计结果
- READ_PAIRS_EXAMINED：被检查的reads对总数
- PERCENT_DUPLICATION：重复率，WGS一般<5%，WES可能到10-20%
- ESTIMATED_LIBRARY_SIZE：有效文库大小
2.4 最终QC——确认数据可用性
排序去重之后的BAM还不能直接给GATK，得先确认质量。我习惯用samtools flagstat看看基本统计量：
samtools flagstat sample_dedup.bam
重点看：
- total reads：初始reads数
- mapped reads：比对上的比例，一般>95%才合格
- properly paired：正确比对的比例，>90%算正常
- duplicates：被标记为重复的比例
我个人会把这三个关键比例放在一起看：
| 指标 | 优秀 | 及格 | 需要重做 |
|---|---|---|---|
| 比对率 (mapped%) | > 98% | > 95% | < 90% |
| 正确配对率 (properly paired%) | > 95% | > 90% | < 80% |
| 去重率 (不重复reads%) | < 5% (WGS) | < 15% (WES) | > 20% |
提示：如果发现比对率偏低，检查一下是不是参考基因组版本选错了。比如人类样本用hg19去比对hg38，匹配率会惨不忍睹。我遇到过好几个人犯这个错误。
最终，生成一个干净的标准BAM文件。你可以直接用这个BAM开始GATK的BaseRecalibrator步骤了。
嗯，到这里，数据预处理这块就梳理完了。记住一句话——你花在数据预处理上的每一分钟，都会在后面的分析中十倍回报给你。