2. 环境搭建
环境搭建这事儿,说难不难,说简单吧,我第一次做的时候也卡了好几天。说白了,就是三件事:装好基础工具、装好 GATK、准备好参考基因组。顺序别搞反,不然依赖问题能让你怀疑人生。
这一节我会把每一步的关键点都撸一遍。你跟着敲命令就行,但注意版本和路径——群里经常有人因为版本对不上跑来问我。
2.1 基础软件安装
BWA、SAMtools、BCFtools 这三个是基因分型的基础三件套。我习惯用 conda 管理,干净又不怕冲突。当然,源码编译也行,但没必要折腾自己。
2.1.1 用 Conda 安装基础软件
conda create -n gvcf_env python=3.9
conda activate gvcf_env
conda install -c bioconda bwa=0.7.17 samtools=1.17 bcftools=1.17
嗯,这里要注意:BCFtools 1.17 以上版本改了部分参数行为,如果你后面用 sort 或 norm 发现报错,大概率是版本问题。我在项目里遇到过这种情况,后来统一锁了 1.17 才稳定下来。
2.1.2 验证安装
装完跑一下帮助命令,确认每个工具都能正常调出。
bwa
samtools --version
bcftools --version
正常输出类似这样(截取部分):
Program: bwa (alignment via Burrows-Wheeler transformation)
Version: 0.7.17-r1188
samtools 1.17
Using htslib 1.17
bcftools 1.17
Using htslib 1.17
如果提示 command not found,检查 conda 环境是否激活,或者路径有没有加。
which bwa 看了下路径才找到原因。建议每次新装工具后跑一下 which 确认。
2.2 GATK 安装
GATK(Genome Analysis Toolkit)是 GVCF 联合基因分型的核心。现在 GATK 4.x 已经全面采用 Java 运行,而且官方直接提供预编译 jar 包。你想想看,连编译都省了。
2.2.1 下载 GATK 包
我一般从官网或 GitHub Release 页下载。版本推荐 4.3.0.0 或 4.4.0.0,这两个我测试过最稳定。
# 直接下载 tar 包(以 4.4.0.0 为例)
wget https://github.com/broadinstitute/gatk/releases/download/4.4.0.0/gatk-4.4.0.0.tar.bz2
tar -xjf gatk-4.4.0.0.tar.bz2
cd gatk-4.4.0.0
其实你也可以用 conda install -c bioconda gatk4,但我觉得自己解压更清楚路径,后面写 pipeline 也好引用。
2.2.2 检查 Java 环境
java -version
# 应该看到类似:
# openjdk version "17.0.8" 2023-07-18 LTS
GATK4 需要 Java 11+,推荐 Java 17。如果你系统 Java 版本太老,就用 conda 装一个:
conda install -c conda-forge openjdk=17
2.2.3 验证 GATK 可运行
./gatk --list
能输出一大串工具列表就算成功。我第一次运行时报了 IllegalArgumentException,仔细一看是内存不够,默认只给了 4G。后来我加了个参数:
./gatk --java-options "-Xmx8G" --list
老项目里这个坑很常见,所以我习惯直接设 -Xmx8G 或更高,省得后面跑着跑着崩了。
./gatk 脚本启动,不要直接 java -jar,因为脚本内部会处理 classpath 和 jvm 参数。
2.3 参考基因组准备
参考基因组是比对和分型的基准。我用得最多的是 GRCh38(hg38)。注意要下载整套文件,不只是 FASTA,还要配套的索引。
2.3.1 下载 GRCh38 参考序列
# 从 UCSC 或 Ensembl 下载
mkdir -p ref/hg38
cd ref/hg38
wget ftp://ftp.ensembl.org/pub/release-110/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
gunzip Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
这个文件比较大,下载可能十几分钟。我建议用 nohup 或 screen 挂后台,否则网络波动断了又要重来。
2.3.2 构建 BWA 索引
bwa index -a bwtsw Homo_sapiens.GRCh38.dna.primary_assembly.fa
这一步会生成 .amb .ann .bwt .pac .sa 五个文件。耗时较长(约半小时),别中断。
2.3.3 构建 SAMtools FASTA 索引
samtools faidx Homo_sapiens.GRCh38.dna.primary_assembly.fa
生成 .fai 文件和 .gzi 文件(如果有 --gzip 参数)。另外,BCFtools 还需要 .dict 字典文件,用 Picard 生成(GATK 自带):
./gatk CreateSequenceDictionary -R Homo_sapiens.GRCh38.dna.primary_assembly.fa
2.3.4 整理参考基因组文件列表
最后的参考目录应该长这样:
| 文件 | 说明 |
|---|---|
Homo_sapiens.GRCh38.dna.primary_assembly.fa |
主要染色体 FASTA |
Homo_sapiens.GRCh38.dna.primary_assembly.fa.fai |
FASTA 索引(samtools 需要) |
Homo_sapiens.GRCh38.dna.primary_assembly.dict |
序列字典(GATK 需要) |
.amb, .ann, .bwt, .pac, .sa |
BWA 索引文件 |
ref.tar.gz 压缩包(Ensembl 或 NCBI 官方提供),解压即用,免去建索引的等待。但记得校验 MD5,我在某第三方站点下载过损坏的索引,跑了一天才发现比对结果全是乱码。
2.4 环境搭建整体流程
下面这张图把上面三个步骤的关系画清楚了。你可以当 checklist 用,每完成一个节点就心里打个勾。
流程图分三层:底层是基础工具(BWA、SAMtools、BCFtools),它们为 GATK 提供数据预处理和后处理能力;中间是 GATK 引擎,负责 HaplotypeCaller 和 GenotypeGVCFs 等核心算法;顶层是参考基因组,所有比对和变异检测都依赖它。
嗯,这一节就到这。你装好环境后,可以顺手跑一下 GATK 自带的测试数据(gatk-4.4.0.0/src/test/resources/ 目录下),看看能不能正常走完 HaplotypeCaller。遇到报错多是 Java 内存不够或路径没写对,排查思路我前面都写了。