0基础学习三代基因组测序组装（9）——GATK检测植物基因组SNP和INDEL变异

GATK的运行依赖于JAVA环境，目前（2023年3月6日）GATK更新到版本4.3.0，可以直接用conda下载。为了避免环境冲突，最好创建一个新环境专门用于GATK和相关变异检测工具运行。

1
2
3
4
5

conda create -n GATK
conda install -c bioconda gatk

conda install bwa-mem2
conda install samtools

bwa-mem2和samtools用于双端序列比对回基因组，需要单独下载这两个软件，后面再说。

安装完成之后可以通过gatk --help查看是否正常。

流程内容主要参考官网Best Practices Workflows的文章。

GATK工具的变异注释主要包括3个部分：数据预处理（Data Pre-processing）、变异检测（Variant Discovery）和变异优化（Callset Refinement）。

以Germline short variant discovery (SNPs + Indels) ，即胚系短变异的发现为例，官网对多个样本（群组数据，Cohort Data）的变异检测分析流程如下：

单个样本（Single-Sample Data）的变异检测标准流程如下：

分析流程基本类似，以我前面组装的三代植物基因组为例跑一下这两个流程，强调一下我用的是植物基因组，后续无法用Germline的注释数据资源对变异集进行功能注释！因此这里只跑到变异集过滤的步骤，拿到SNP和INDEL。

2.1 数据预处理

1. 构建参考基因组索引：组装的基因组作为reference参考基因组，首先需要对参考基因组建立索引，方便后续进行比对和对参考基因组进行查询。注意三个软件的索引文件不同，每个软件都必须建立索引。

1
2
3
4
5
6
7

#!/bin/bash
ref=$1

bwa-mem2 index $ref
samtools faidx $ref
referencename=`basename $ref | sed  "s/fasta/dict/" `		# .fasta文件后缀改为.dict
gatk CreateSequenceDictionary -R $ref -O $referencename 

2. fasta文件转化uBAM文件，标记adapter 序列：组装好的基因组格式是fasta格式，需要转换成uBAM格式（umapped的BAM文件），接着标记illumina二代测序的adapter序列。本质上都是调用了Picard工具，也可以直接用java写脚本。

1
2
3
4
5
6

#!/bin/bash
sampleName=$1

gatk FastqToSam -F1 raw_fastq/${sampleName}_1.fq.gz -F2 raw_fastq/${sampleName}_2.fq.gz -RG $sampleName -SM $sampleName -O ubam/${sampleName}.bam

gatk MarkIlluminaAdapters -I ubam/${sampleName}.bam -O markadapeters/${sampleName}.markadapeters.bam -M markadapeters/${sampleName}.metrics.txt

3. 标记后的序列转化成fastq格式，回比参考基因组，得到干净的BAM文件：第二和第三步来自官方的文章(How to) Map and clean up short read sequence data efficiently，通过联合SamToFastq、BWA - MEM和MergeBamAlignment三个程序，节省比对时间，绕过占用空间过大的SAM文件，MergeBamAlignment可以将已排序的SAM中的定义信息（我这里将SAM文件转换成BAM文件）与uBAM中的定义信息进行合并，得到干净的BAM并进行排序和构建索引。

1
2
3
4
5
6
7
8
9
10
11
12

#!/bin/bash
sampleName=$1
threads=50
ref=/public/home/wlxie/biosoft/GATK_file/gatk/ref/luobuma.fasta

gatk SamToFastq -I markadapeters/${sampleName}.markadapeters.bam -F  interleaved_fq/${sampleName}_1.interleaved.fq.gz -F2  interleaved_fq/${sampleName}_2.interleaved.fq.gz -CLIP_ATTR XT -CLIP_ACT 2

bwa-mem2 mem -M -t $threads $ref  interleaved_fq/${sampleName}_1.interleaved.fq.gz interleaved_fq/${sampleName}_2.interleaved.fq.gz | samtools view -Sb - > raw_bam/${sampleName}.bam

gatk MergeBamAlignment -R $ref -UNMAPPED  ubam/${sampleName}.bam -O align_bam/${sampleName}.bam -ALIGNED  raw_bam/${sampleName}.bam -MC true --CREATE_INDEX true

rm -rf markadapeters/${sampleName}.markadapeters.ba interleaved_fq/${sampleName}_1.interleaved.fq.gz  interleaved_fq/${sampleName}_2.interleaved.fq.gz raw_bam/${sampleName}.bam 	# 删除中间文件

4. 标记重复序列：做SNP分析前最重要的一点就是标记重复序列（mark duplicate），二代测序是在PCR扩增的基础上进行的，因此PCR扩增产生的多拷贝会结合到flowcell的不同位置，生成完全相同的测序cluster，最后得到重复序列。这一步就是标记这些重复序列（但是没有删除，对结果应该不影响），最后得到一个metrics文件（duplicate的统计信息）和一个bam文件（duplicate的详细信息，注意要创建索引）。

1
2
3
4
5
6

#!/bin/bash
sampleName=$1

gatk MarkDuplicates -I align_bam/${sampleName}.bam -O markdup/${sampleName}.markdup.bam -M markdup/${sampleName}.markdup_metrics.txt --CREATE_INDEX true

rm -rf  align_bam/${sampleName}.bam 	# 删除中间文件

最终生成的bam文件进行下一步变异检测，可以用samtools -view 查看bam文件的内容（也没啥好看的，感兴趣可以看看之前博客写的sam文件格式解读）。

• 注意：官网的数据预处理后续还有一步碱基质量重校正BQSR（Base Quality Scores Recallbrate），官方提供了两个工具BaseRecalibrator和ApplyBQSR。第一个工具计算需要校正的reads和特征值，输出校准表文件，需要注意的是，这个工具需要一个或者多个已知且可靠的物种变异位点数据库，比如人类就有千人基因组计划的各种变异位点数据库，1000G_omni2.5.hg38.vcf.gz、dbsnp_146.hg38.vcf.gz等等；第二个工具根据第一个工具生成校准表，重新调整原来BAM文件中的碱基质量值，重新输出到一个新的BAM文件中。

因为我这个植物没有已知的可靠变异位点数据，因此不用做最后这一步。

本部分程序需要运行10小时。

2.2 变异检测

因为我只有一个样本，所以可以按照单个样本（Single-Sample Data）的标准流程来做，也可以按照多样本（Cohort Data）的流程做。该部分主要用到的工具是HaplotypeCaller，两个流程只是对HaplotypeCaller工具产生的结果做了不同的处理，最后都是得到包含SNP和Indel的VCF文件。我也将分别跑两个流程来对比下差异。

首先要明白HaplotypeCaller这个工具具体做了什么，是怎么找出单碱基变异的：

• 1.定义活跃区域（Define active regions）：根据是否存在变异来确定需要操作的基因组的活跃区域。
• 2.通过组装活跃区域确定单倍型（Determine haplotypes by assembly of the active region）：对于每个活跃区域，构建一个类似De Bruijn图来重新组装活性区域，识别可能的单倍型，然后用Smith-Waterman算法将每个单倍型与参考基因组单倍型重新比对，发现潜在的变异位点。
• 3.确定read单倍型似然值（Determine likelihoods of the haplotypes given the read data）：对每个活跃区域使用PairHMM算法对每个read与每个单倍型进行两两比对，生成单倍型似然矩阵，将似然值边缘化，以获得每个潜在变异位点的等位基因的可能性。
• 4.指定样本基因型（Assign sample genotypes）：对每个潜在的变异位点使用贝叶斯算法，转化每个位点的基因型的似然值。然后将最可能的基因型指定为样本基因型。

以上流程翻译自HaplotypeCaller – GATK (broadinstitute.org)，这个工具可用参数非常之多，下面跑的流程就展示一些常用的。

2.2.1 多样本的SNP和INDEL检测

1. 使用HaplotypeCaller的GVCF模式，找到每个样本SNP和INDEL变异。在GVCF模式下，每个样本的结果文件以gvcf（genomic vcf）格式文件呈现，实际上gvcf格式和vcf格式类似，gvcf记录所有位点的突变情况，并且提供这些位点是否是纯和的置信度，主要还是方便将所有样本的gvcf联合起来方便分析。

1
2
3
4
5
6

#!/bin/bash
sampleName=$1
threads=50
ref=/public/home/wlxie/biosoft/GATK_file/gatk/ref/luobuma.fasta

gatk HaplotypeCaller -R $ref --native-pair-hmm-threads ${threads} --emit-ref-confidence GVCF -I ../../pre_processing/markdup/${sampleName}.markdup.bam -O ${sampleName}.g.vcf.gz

这一步是整个SNP和IDEL检测中运行时间最长，需要算力最多的一步（主要是Pair-HMM算法花时间），Pair-HMM算法在本地调用的线程数是可以更改的，官方是给定默认值为4。GATK4.0版本开始放弃了多线程任务，这个参数可能是更新后遗漏的，因为我这里用的50线程和默认的4线程跑几乎没有区别，都是在28小时左右完成，此处参数存疑。

2. 建库合并所有样本的GVCF文件（单样本不用做这步），并将GVCF文件转化为VCF文件。对于多样本的GVCF合并，现在官方建议用GenomicsDBImport这个工具进行建库合并，速度会快很多。

1
2
3
4
5
6
7
8

#!/bin/bash
sampleName=$1
threads=50
ref=/public/home/wlxie/biosoft/GATK_file/gatk/ref/luobuma.fasta

# gatk GenomicsDBImport $(for file in `ls *.g.vcf.gz`; do echo "-V $file"; done)  --genomicsdb-workspace-path database -L chr01		# 单样本不用做这步，因为就一个GVCF。多样本注意-L参数是建库必须的,根据fasta参考基因组的染色体名称命名，拆分

gatk GenotypeGVCFs -R $ref -V ${sampleName}.g.vcf.gz -O raw_variants.vcf.gz		# 多样本使用参数-V gendb://database 也就是上一步建的数据库名，单样本直接用gvcf文件

可以看看最后生成的VCF文件：

文件解读放在最后过滤和拆分SNP和INDEL的时候再说，这里是初步得到变异信息，需要经过过滤和筛选。

本部分程序需要运行28小时。

2.2.2 单个样本的SNP和INDEL检测

使用HaplotypeCaller默认的single-sample模式，直接生成统计SNP和INDEL变异的VCF文件。

1
2
3
4
5
6

#!/bin/bash
sampleName=$1
threads=50
ref=/public/home/wlxie/biosoft/GATK_file/gatk/ref/baima.fasta

gatk HaplotypeCaller -R $ref --native-pair-hmm-threads ${threads} -I ../../pre_processing/markdup/${sampleName}.markdup.bam -O ${sampleName}.vcf.gz

同上一个步骤，此处参数--native-pair-hmm-threads对运算速度的提升存疑，最后同样是生成VCF文件，运行时间同样为28小时。

2.3 变异过滤（优化）

变异集过滤方法主要有两种：

• 1.软过滤：基于机器学习的方法，对原始vcf文件进行变异质量重矫正和过滤。比如有基于卷积神经网络CNN的CNNVariantTrain（有预训练的模型1D和2D），VariantRecalibrator、ApplyVQSR等可以用已知的人类变异数据集作为训练集，检测得到的SNP和INDEL的准确性（官方推荐用于人类变异过滤的方法，Variant Quality Score Recalibration，VQSR）。缺点显而易见，需要已知的真实变异数据集，除人类以外大多数生物都没有这方面的数据集。如果是研究人类基因组的话，可以从GATK官网资源处下载。
• 2.硬过滤：通过对6个指标的硬性阈值筛选质量合格的SNP和INDEL。

记录下硬过滤的6个指标，有些说明看不懂干脆都放英文了，参考自官网Hard-filtering germline short variants：

（1）QualByDepth (QD)：This is the variant confidence (from the QUAL field) divided by the unfiltered depth of non-hom-ref samples. 变异置信度，官方建议过滤该值小于2的变异。

（2）FisherStrand (FS)：This is the Phred-scaled probability that there is strand bias at the site.

（3）StrandOddsRatio (SOR)：This is another way to estimate strand bias using a test similar to the symmetric odds ratio test.

（4）RMSMappingQuality (MQ)：This is the root mean square mapping quality over all the reads at the site.比对reads质量的平方根。

（5）MappingQualityRankSumTest (MQRankSum)：This is the u-based z-approximation from the Rank Sum Test for mapping qualities.

（6）ReadPosRankSumTest (ReadPosRankSum)：This is the u-based z-approximation from the Rank Sum Test for site position within reads.

看不懂没关系，官方给出了6个硬过滤指标在SNP和INDEL中的阈值设置，详情可以看Hard-filtering germline short variants – GATK (broadinstitute.org)。以下例子均以官方的硬过滤指标为准，感兴趣可以去官网看各个参数的作用或者自己微调。

SelectVariants工具用于从vcf文件中提取SNP和INDEL信息，VariantFiltration工具用于硬过滤筛选：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

#!/bin/bash
sampleName=$1
threads=50
VARIANTS=/public/home/wlxie/biosoft/GATK_file/gatk/variants_discover/luobuma/raw_variants.vcf.gz

# SNP筛选、过滤和提取
gatk SelectVariants -select-type SNP -V $VARIANTS --restrict-alleles-to BIALLELIC  -O ${sampleName}_SNP.vcf.gz		# BIALLELIC 限制双等位基因，不考虑其他等位多态性
gatk VariantFiltration -V ${sampleName}_SNP.vcf.gz --filter-expression "QD < 2.0 || MQ < 40.0 || FS > 60.0 || SOR > 3.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" --filter-name "Filter" -O ${sampleName}_SNP.filter.vcf.gz
gatk SelectVariants  -V ${sampleName}_SNP.filter.vcf.gz --exclude-filtered true -O final.${sampleName}_SNP.vcf.gz		# 只显示通过过滤的变异（pass）
rm -rf ${sampleName}_SNP.vcf.gz*
rm -rf ${sampleName}_SNP.filter.vcf.gz*

# INDEL筛选、过滤和提取
gatk SelectVariants -select-type INDEL -V $VARIANTS --restrict-alleles-to BIALLELIC  -O ${sampleName}_INDEL.vcf.gz		# BIALLELIC 限制双等位基因，不考虑其他等位多态性
gatk VariantFiltration -V ${sampleName}_INDEL.vcf.gz --filter-expression "QD < 2.0 || FS > 200.0 || SOR > 10.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" --filter-name "Filter" -O ${sampleName}_INDEL.filter.vcf.gz
gatk SelectVariants  -V ${sampleName}_INDEL.filter.vcf.gz --exclude-filtered true -O final.${sampleName}_INDEL.vcf.gz		# 只显示通过过滤的变异（pass）
rm -rf ${sampleName}_INDEL.vcf.gz*
rm -rf ${sampleName}_INDEL.filter.vcf.gz*

得到的final.sampleName_SNP.vcf.gz和final.sampleName_INDEL.vcf.gz为最终的变异集结果文件。

因为结果文件时一个压缩过后的vcf文件，且vcf文件中前面带#部分的注释内容是用不到的，后面每一行代表一个变异位点信息，因此可以直直接统计行数来得到最终的SNP和INDEL的数量。

1

zcat final.luobuma_SNP.vcf.gz | grep -v -P "^#" -c

用2.2.1多样本流程的SNP和INDEL结果文件为final.luobuma_SNP.vcf.gz和final.luobuma_INDEL.vcf.gz；用2.2.2单个样本流程的SNP和INDEL结果为final.luobuma_sm_SNP.vcf.gz和final.luobuma_sm_INDEL.vcf.gz。可以看到两者在统计SNP和INDEL数量上的差距非常小，说明这两个流程对单样本来说都是可以用的。

以final.luobuma_SNP.vcf.gz文件进行解读，通过less -S命令一行一条信息查看文件内容：

每列信息如下：

• 1.CHROM：染色体信息
• 2.POS：变异所在参考基因组的位置
• 3.ID：变异的ID，如果有参考变异集，会给出id，否则为.表示新发现的变异
• 4.REF：变异在参考基因组上的信息，必须为ATCGN五个之一
• 5.ALT：突变之后的情况，类型同上，.表示缺失
• 6.QUAL：突变后的质量值，质量值越高越可靠，通常只用pass的数据
• 7.FILTER：是否通过过滤
• 8.INFO：每个位点的详细信息（包括硬过滤的指标，详细可以到header里找）
• 9.FORMAT：格式
• 10.样本名（实际是前面格式的具体值）

主要解释一下第九列和第十列，就是上图中红色框框起来的部分，两列值是用冒号分隔且一一对应的，需要注意的值是GT：

• GT：0表示和参考序列一致（REF allele），1表示和样本序列一致（ALT allele），双等位基因只有0和1，0/1和0|1表示杂合，1/1和1|1表示纯和。“|”和“/”区别是前者是phased genotype，就是知道REF/ALT allele是来自于父本还是母本，在这里对我这个植物基因组没有什么意义，全都统计进杂合和纯和SNP个数就行。
• AD：REF和ALT allele的覆盖度，在二倍体是是用逗号分割的两个值表示，前一个代表参考基因组的基因型，后者代表样本基因型。
• DP：样本中该位点覆盖度，AD两个数字的和。

3.1 杂合率统计

分别统计SNP和INDEL文件中杂合单碱基变异的个数：

1
2
3
4

zcat final.luobuma_SNP.vcf.gz | grep -c "0/1"
zcat final.luobuma_SNP.vcf.gz | grep -c "0|1"
zcat final.luobuma_INDEL.vcf.gz | grep -c "0/1"
zcat final.luobuma_INDEL.vcf.gz | grep -c "0|1"

杂合率 = (杂合SNP数 + 杂合INDEL数) / 基因组大小

此处计算的杂合率可以和前面做的基因组Survey做个比较，说明基因组Survey的可靠性。

3.2 单碱基准确度计算

我们知道测序过程中不可避免地存在错误，三代测序数据单碱基变异的来源，包括真实的单碱基变异和测序错误导致的单碱基变异。当测序错误导致的单碱基变异存在于参考基因组上时，利用二代测序数据进行单碱基变异的检测时，会将其识别为纯合单碱基变异。因此，可以将三代数据组装的最终版基因组作为参考基因组，利用二代数据将纯合子单碱基变异率作为组装结果的错误率，即:

组装结果的准确率 = 1 - 纯合子单碱基变异率

1
2
3
4

zcat final.luobuma_SNP.vcf.gz | grep -c "1/1"
zcat final.luobuma_SNP.vcf.gz | grep -c "1|1"
zcat final.luobuma_INDEL.vcf.gz | grep -c "1/1"
zcat final.luobuma_INDEL.vcf.gz | grep -c "1|1"

这种单碱基准确度的计算结果也可以作为基因组组装质量的评估指标之一，即序列一致性评估——利用高质量的二代测序数据来评估三代测序数据组装结果在单碱基水平上的准确性。