📚 课程目录第2章参考基因组准备

参考基因组准备：地基打不牢，后面全白干

变异检测这活儿，说到底就是拿着你的测序数据去跟一个「标准答案」比对。这个标准答案，就是参考基因组。我入行那会儿吃过一次大亏——用了个旧版本的参考基因组，结果一堆变异位点位置对不上，后来全部重跑。嗯，从那以后，我再也不敢在基因组准备上偷懒了。

说白了，这一步就是给后续所有分析铺路。路铺不平，后面跑得再快也是白费油。今天咱们就把 GRCh38 的下载、索引建立、辅助文件整理，一条一条捋清楚。

2.1 选择哪个版本？GRCh38 还是 hg38？

其实这两个名字说的是同一个东西。GRCh38 是国际基因组参考联盟（GRC）给的名字，hg38 是 UCSC 的叫法。我个人习惯用 GRCh38.p13 或者 GRCh38.p14（补丁版本），因为 GATK 官方推荐的也是这个。

版本名称	组织来源	GATK 推荐度	常见文件名标记
GRCh38	GRC / NCBI	⭐⭐⭐⭐⭐	GRCh38.p13.genome.fa
hg38	UCSC	⭐⭐⭐	hg38.fa
hs38DH	BWA 作者	⭐⭐⭐	hs38DH.fa（含 decoy）

我的建议：直接去 GATK Resource Bundle 下载 GRCh38 的 FASTA。这个版本经过 GATK 团队验证，chromosome naming 统一，省去反复改名的麻烦。

2.2 下载参考基因组

常用的下载源有这几个：

NCBI 官方：ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/405/
Ensembl：ftp://ftp.ensembl.org/pub/release-X/fasta/homo_sapiens/dna/
GATK Resource Bundle（推荐）：gs://gatk-legacy-bundles/b37/ 或 gs://genomics-public-data/resources/broad/hg38/v0/

我一般用 wget 直接拉，速度快、稳定。举个例子：

# 从 GATK 官方 Google Cloud 下载（需要 gcloud auth）
gsutil cp gs://genomics-public-data/resources/broad/hg38/v0/Homo_sapiens_assembly38.fasta ./
gsutil cp gs://genomics-public-data/resources/broad/hg38/v0/Homo_sapiens_assembly38.fasta.fai ./

# 或者用 wget 从 NCBI 拉（慢一些但不用翻墙）
wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/405/GCF_000001405.39_GRCh38.p13/GCF_000001405.39_GRCh38.p13_genomic.fna.gz
gunzip GCF_000001405.39_GRCh38.p13_genomic.fna.gz
mv GCF_000001405.39_GRCh38.p13_genomic.fna GRCh38.p13.genome.fa

小提示：下载下来先做一遍 md5sum 校验。我曾经遇到过网络中断导致文件损坏，跑 bwa index 到一半报错，查了半小时才发现是文件本身有问题。

2.3 建立 bwa 索引

bwa index 是做比对（bwa mem）的前置条件。没有这一步，后面没法跑。命令很简单：

bwa index -a bwtsw GRCh38.p13.genome.fa

-a bwtsw 参数是专门为大型基因组（>10Mb）准备的。小基因组可以用 is，但人类基因组必须用 bwtsw。跑完之后你会看到这些文件：

GRCh38.p13.genome.fa.amb — 模糊碱基表（ambiguous base）
GRCh38.p13.genome.fa.ann — 注释信息
GRCh38.p13.genome.fa.bwt — 后缀数组（核心索引）
GRCh38.p13.genome.fa.pac — 打包的参考序列
GRCh38.p13.genome.fa.sa — 后缀数组索引

你想想看，这五个文件加起来大概 2-3GB，构建时间取决于机器内存和磁盘。我自己的服务器 64GB 内存，大约跑了 4 个小时。如果你的机器内存小于 32GB，建议先用 ulimit -s unlimited 解开栈限制，不然容易崩。

注意：bwa index 过程中不要中断进程。一旦中断，生成的 .bwt 和 .sa 文件可能就是坏的，得删掉重新跑。我之前在服务器上用 screen 跑，结果 session 超时断开，第二天回来看见报错，心态直接炸了。

2.4 samtools faidx：建立 FASTA 索引

这个索引文件主要是给 samtools、picard、GATK 这些工具用的。它们需要快速定位到基因组的某个区域，.fai 就是干这个的。

samtools faidx GRCh38.p13.genome.fa

生成的文件是 GRCh38.p13.genome.fa.fai。你可以用 less 看一下内容，格式大概是：

chr1	248956422	7	60	61
chr2	242193529	254136961	60	61
chr3	198295559	498433122	60	61
...

每一列分别是：染色体名、长度、在FASTA文件中的偏移量、单行碱基数、行间隔字符数。这个文件很小，但千万别丢了。很多变异检测流程（比如 GATK HaplotypeCaller）启动时会先读它。

2.5 创建 .dict 序列字典

.dict 文件也是一个文本格式的索引，记录每条染色体的长度和排序。它主要由 Picard 和 GATK 使用。注意：有些工具（如 GATK4）对 .dict 的存在是硬性要求，没有它直接报错退出。

# 使用 Picard 工具（推荐）
java -jar picard.jar CreateSequenceDictionary \
    R=GRCh38.p13.genome.fa \
    O=GRCh38.p13.genome.dict

# 或者使用 GATK4 自带的命令（如果安装了）
gatk CreateSequenceDictionary \
    -R GRCh38.p13.genome.fa \
    -O GRCh38.p13.genome.dict

.dict 文件的内容长这样：

@HD	VN:1.5	SO:coordinate
@SQ	SN:chr1	LN:248956422	UR:GRCh38.p13.genome.fa	M5:6aef897c3d6ff0c78aff06ac189178dd
@SQ	SN:chr2	LN:242193529	UR:GRCh38.p13.genome.fa	M5:f98db6892b9d4a1b7a2c3c0a1b2c3d4e
...

关键信息：SN 是染色体名，LN 是长度，M5 是 MD5 校验值——这可以用来比对 FASTA 是否被改过。如果后续 GATK 报错说染色体名不匹配，多半就是 .dict 和 FASTA 里的 @SQ 对不上。我遇到过一次，后来养成习惯：每次更新基因组都重新生成 .dict，绝不复用旧文件。

2.6 资源文件获取（dbsnp、已知 indel 等）

变异检测通常需要一些已知变异数据库来做过滤和校准。这些文件虽然不是参考基因组本身，但跟参考基因组版本严格绑定。GRCh38 对应的资源文件推荐从以下渠道获取：

资源类型	文件名示例	下载来源	用途
dbsnp 151 / 155	Homo_sapiens_assembly38.dbsnp138.vcf.gz	GATK Resource Bundle	VQSR / 已知位点注释
已知插入缺失	Homo_sapiens_assembly38.known_indels.vcf.gz	GATK Resource Bundle	Indel 真实集合（BQSR）
Mills 和 1000G gold standard	Mills_and_1000G_gold_standard.indels.hg38.vcf.gz	GATK Resource Bundle	Indel 训练集
hapmap_3.3	hapmap_3.3.hg38.vcf.gz	GATK Resource Bundle	VQSR 训练集
1000G_omni2.5	1000G_omni2.5.hg38.vcf.gz	GATK Resource Bundle	VQSR 训练集

# 使用 gsutil 批量下载（需要先认证）
gsutil -m cp \
  "gs://genomics-public-data/resources/broad/hg38/v0/Homo_sapiens_assembly38.dbsnp138.vcf.gz" \
  "gs://genomics-public-data/resources/broad/hg38/v0/Homo_sapiens_assembly38.dbsnp138.vcf.gz.tbi" \
  "gs://genomics-public-data/resources/broad/hg38/v0/Homo_sapiens_assembly38.known_indels.vcf.gz" \
  "gs://genomics-public-data/resources/broad/hg38/v0/Homo_sapiens_assembly38.known_indels.vcf.gz.tbi" \
  ./

关键提醒：下载下来的 .vcf.gz 一定要配套下载 .tbi (tabix index) 文件。没有索引，GATK 很多模块会拒绝读取。我见过不止一个同事在跑 BQSR 的时候报错，最后发现是忘记下载 .tbi 了。

2.7 整理一下，检查清单

做完上面所有步骤后，你的参考基因组目录下应该至少有这七个文件：

GRCh38.p13.genome.fa — 主要 FASTA
GRCh38.p13.genome.fa.amb — bwa 模糊碱基表
GRCh38.p13.genome.fa.ann — bwa 注释
GRCh38.p13.genome.fa.bwt — bwa 后缀数组
GRCh38.p13.genome.fa.pac — bwa 压缩序列
GRCh38.p13.genome.fa.sa — bwa 后缀数组索引
GRCh38.p13.genome.fa.fai — samtools 索引
GRCh38.p13.genome.dict — Picard/GATK 序列字典

另外，资源文件放在单独的目录 resources/ 下，方便管理。我习惯这样组织目录结构：

project/
├── ref/
│   ├── GRCh38.p13.genome.fa
│   ├── GRCh38.p13.genome.fa.fai
│   ├── GRCh38.p13.genome.dict
│   └── (bwa index 五个文件)
├── resources/
│   ├── dbsnp138.vcf.gz
│   ├── dbsnp138.vcf.gz.tbi
│   ├── known_indels.vcf.gz
│   ├── known_indels.vcf.gz.tbi
│   └── (其他 VQSR 训练集)
├── data/
│   └── (原始测序数据)
└── output/
    └── (分析结果)

这样做的好处是：路径清晰，脚本里引用路径不会乱套，别人接手也能一眼看懂。嗯，这一块的基础工作虽然琐碎，但我可以负责任地说——前期做扎实了，后面每一个变异检测的步骤都会顺很多。不信你试试。