3. 输入数据准备：FASTQ、BAM、参考基因组
做GATK变异检测，第一步不是跑命令，是搞清楚你手里有什么数据。我见过不少人上来就对着BAM文件一顿操作，结果跑完才发现样本名不对、参考基因组版本搞混了。嗯，这节课我们就聊聊输入数据的门道。
说白了，GATK流水线吃三类东西：原始测序数据（FASTQ）、比对后的文件（BAM/CRAM）和参考基因组。这三样缺一不可，而且每样都得收拾干净。
一句话总结：花80%时间准备数据，20%时间跑分析。数据质量决定最终结果。
3.1 FASTQ文件：原始数据的“身份证”
FASTQ里存着每条read的序列和质量分数。一个标准的FASTQ记录长这样：
@EAS139:136:FC706VJ:2:2104:15343:197393 1:Y:18:ATCACG
CTCAAGGTTGTTGCAAGACGGGAGGTAGGTGTCTGGGTAG
+
IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII
第一行是read ID，后面跟着碱基序列，然后一行“+”分隔符，最后是质量分数（Phred+33编码）。我个人习惯拿到FASTQ先做几件事：
- 检查read数目：用
zcat file.fastq.gz | wc -l再除以4，确认样本测了多少对reads。 - 看质量分布：FastQC跑一轮，重点关注Q30比例、GC含量、过代表序列。
- 确认命名规范：样本名不能有空格、括号，项目组统一用“样本ID_批次_ Lane”格式。我在一个肿瘤项目里遇到过样本名带“/”的情况，结果GATK直接报错退出。
为什么要做这些？因为GATK的BQSR（碱基质量校正）严重依赖原始质量分数，如果FASTQ本身就有问题，后面全白搭。
3.2 BAM/CRAM文件：比对后的“地图”
BAM是SAM的二进制版本，CRAM是更压缩的格式。GATK官方推荐用CRAM，但我自己更习惯在分析前期用BAM，调试方便。
一个合格的BAM文件需要满足：
- 排序正确：GATK要求按坐标排序，而且必须带索引（.bai或.crai）。
- Read Group完整：@RG标签必须有ID、SM（样本名）、LB（文库名）、PL（测序平台）。缺少任何一个，GATK的HaplotypeCaller可能直接罢工。
- 标记过重复：用Picard MarkDuplicates或sambamba markdup处理过，不然变异检测结果会被PCR重复污染。
我曾经踩过的坑：有一次用GATK4 BaseRecalibrator，一直报“No reads found in interval”，查了两天发现BAM的@RG标签里SM字段写成了“Sample_1”，而参考基因组的字典里对应的样本名是“Sample_1_rerun”，不匹配。从那以后我建了个检查脚本，跑前先校验RG标签。
验证BAM完整性我常用这个命令：
samtools quickcheck sample.bam && echo "OK" || echo "Corrupt or incomplete"
samtools stats sample.bam | grep "reads mapped:"
3.3 参考基因组：变异的“字典”
参考基因组就是一把尺子，所有变异都是相对于它定义的。GATK需要三个东西：
- FASTA文件：比如hg38.fa，要求每个染色体一行，换行符统一（推荐用
picard NormalizeFasta处理）。 - 索引文件：
.fai由samtools faidx生成，.dict由picard CreateSequenceDictionary生成。 - 已知位点文件（可选但强烈推荐）：比如dbSNP、1000G的VCF，用于BaseRecalibration和VQSR。
不同项目用的参考版本可能不同。我常用的做法是把参考基因组统一放到一个公共目录，软链接到项目目录，避免到处拷贝浪费空间：
ln -s /ref/hg38/hg38.fa /project/ref/hg38.fa
ln -s /ref/hg38/hg38.fa.fai /project/ref/hg38.fa.fai
ln -s /ref/hg38/hg38.dict /project/ref/hg38.dict
一个建议：团队成员一定要统一参考基因组版本。我在一个跨国合作项目里吃过亏——中方用GRCh38，美方用hg38（其实是同一个东西但命名不同），最后合并VCF时染色体名字对不上，折腾了三天重跑。
3.4 本章知识结构图
下面这张图总结了输入数据准备的核心流程：从FASTQ到BAM再到参考基因组，每一步都有质控要点。
3.5 数据准备检查清单
为了不让项目卡在数据环节，我个人整理了一个自查表，每次启动新分析前过一遍：
| 检查项 | 具体内容 | 验证命令/工具 |
|---|---|---|
| FASTQ完整性 | reads数目匹配、无截断文件 | zcat | wc -l / FastQC |
| BAM排序 | 按坐标排序，有索引 | samtools flagstat / samtools idxstats |
| BAM RG标签 | ID、SM、LB、PL齐全 | samtools view -H | grep '@RG' |
| 重复标记 | 标记过的BAM有duplicate flag | samtools view -f 1024 | wc -l |
| 参考基因组索引 | .fai 和 .dict 文件存在且匹配 | samtools faidx / picard CreateSequenceDictionary |
| 已知位点VCF | 索引存在，染色体名称与参考一致 | tabix -p vcf known_sites.vcf.gz |
小技巧：把这些检查写成脚本，每天启动分析前跑一遍。check_inputs.sh 我从第一个项目用到现在，改了几十次，但核心逻辑没变过——宁可慢一点，不然后面出问题排查更浪费时间。
3.6 实战经验：数据目录结构
最后，说说目录组织。我习惯这样安排：
/project/
├── data/
│ ├── fastq/ # 原始FASTQ，按样本分文件夹
│ ├── bam/ # 比对后的BAM，含.bai索引
│ └── ref/ # 参考基因组软链接
├── scripts/
│ ├── check_inputs.sh # 数据检查脚本
│ └── run_gatk_pipeline.sh
└── results/
├── gvcf/
└── vcf/
为什么要强调目录结构？因为GATK的中间产物特别多，gVCF、BQSR表、VCF、过滤后的变异……如果乱放，出了错根本不知道哪个文件是哪个环节的。