奇怪，HaplotypeCaller中AD之和不等于DP

生信菜鸟团

发布于 2022-04-08 17:40:06

2.5K0

在基因组分析中，处理流程从上游测序数据到下游突变分析，中间的关键就是call突变。

GATK是最常使用的call突变工具包，其HaplotypeCaller和Mutect2命令分别用于call germline和somatic突变。我们相信GATK软件的可靠性以及其reassemble特性的优良表现，但是往往对其内部的复杂性有所忽视。

问题出现

多个样本走GATK Best Practice，通过HaplotypeCaller命令call出germline vcf文件。

看一下某突变在某样本中的详细信息。

$ bcftools view -r 22:16052126 -H -s SAMPLE1 temp.vcf.gz
22      16052126        .       T       C       223.71  .       AC=1;AF=0.017;AN=2;BaseQRankSum=2.2;DP=4056;ExcessHet=3.2494;FS=3.624;InbreedingCoeff=0.1267;MLEAC=8;MLEAF=0.019;MQ=53.36;MQRankSum=-2.091;QD=1.73;ReadPosRankSum=-0.719;SOR=0.947      GT:AD:DP:GQ:PGT:PID:PL:PS       0/1:48,7:67:99:.:.:109,0,1853:.

前几列分别是CHROM POS ID REF ALT QUAL FILTER，表示突变坐标等基础信息。

后面三列比较长的分别是INFO FORMAT SAMPLE1。INFO是多样本同时考虑计算出的突变质量信息，FORMAT表示样本列内容的格式，SAMPLE1表示该突变在SAMPLE1样本中的具体信息。FORMAT与样本列内容通过:分隔，逐个对应，比如，0/1的意思是GT=0/1。

我们会比较关注DP和AD信息。DP代表覆盖该位点的reads数，AD代表支持ref和alt的reads数，这些信息有助于评估该突变是否可信。

需要注意的是DP有两个，分别是INFO/DP和FORMAT/DP，在多样本vcf中，INFO/DP是FORMAT/DP的加和。

$ bcftools view -h temp.vcf.gz | grep -E "##FORMAT|##INFO" | grep -E "ID=AD|ID=DP"
##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Allelic depths for the ref and alt alleles in the order listed">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth (reads with MQ=255 or with bad mates are filtered)">
##INFO=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth; some reads may have been filtered">

在header中，可以看到不同tag的详细解释。

仔细看一下这个突变，INFO/DP=4056暂不用关注，因为是多样本的结果，AD=48,7表示支持T的reads有48个，支持A的有7个，FORMAT/DP=67表示该位点在筛选后有67个reads覆盖。

咦，等等，按理说覆盖位点的reads要么支持ref，要么支持alt，一共有67个reads覆盖，可是支持T和A的分别有48个和7个，加起来55，不相等啊。

看一看其他突变。

$ bcftools view -r 22:16122834 -H -s SAMPLE1 temp.vcf.gz
22      16122834        rs375374612     A       G       117347  .       AC=1;AF=0.259;AN=2;BaseQRankSum=1.69;DB;DP=26070;ExcessHet=91.9727;FS=1.288;InbreedingCoeff=-0.3437;MLEAC=117;MLEAF=0.261;MQ=40.49;MQRankSum=-1.067;QD=5.39;ReadPosRankSum=-0.129;SOR=0.783     GT:AD:DP:GQ:PGT:PID:PL:PS       0/1:779,412:1191:99:.:.:8034,0,17242:.

779+412=1191，是一致的。那为什么某些位点会出现AD之和不等于DP的情况呢。

官方解释

经过搜索，GATK Team的一篇文档给出了答案。Allele Depth (AD) is lower than expected[1]。

下面对其内容进行简要概括并予以评价。

问题的核心在于uninformative reads。reads在filter过后，会被区分成两类，uninformative read和informative read。

假设目前有两个read，和两个可能的allele，通过计算可以得到每个read支持每个allele的likelihood。

R	Likelihood of A	Likelihood of T
1	3.8708e-7	3.6711e-7
2	4.9992e-7	2.8425e-7

有意思，我们通常认为，read上是A，那这条read当然是支持A的，但是GATK不这么认为，他们选择根据haplotype对每个allele计算一个likelihood。具体细节请阅读How HaplotypeCaller works[2].

然后将likelihood转变成Phred-scaled likelihood。

R	Phred-scaled likelihood of A	Phred-scaled likelihood of T
1	-6.4122	-6.4352
2	-6.3011	-6.5463

现在对read进行判断。计算方式很简单，如果最可能的allele与第二可能的allele相差大于0.2，就是informative，小于0.2，就是uninformative。

R	Difference	Tag
1	-6.4122-(-6.4352)=0.0230	uninformative
2	-6.3011-(-6.5463)=0.2452	informative

选择0.2作为cutoff，根据0.2=10^0.2=1.585，意味着最可能的allele比第二可能的allele，其真实的可能性大了约60%。

这里类似于引入了一个显著的概念，只有显著支持最可能的allele时，这个read才是informative的。

而这些和DP、AD有什么关系呢？

informative read会被AD和DP同时计数，而uninformative read只被AD计数，不被DP计数。GATK选择不在AD中包含uninformative read的原因是他们认为这些不够显著的read不可信。

所以AD之和与DP的差，是uninformative read的数量。

一些问题

该文档还提到一种情况。

2 151214 . G A 673.77 . AN=2;DP=10;FS=0.000;MLEAF=0.500;MQ=56.57;MQ0=0;NCC=0;SOR=0.693 GT:AD:DP:GQ:PL 0/1:0,0:10:38:702,0,38

这里更离谱，AD全是0，而DP是10，相当于覆盖的reads全是uninformative read，这样居然也能call出genotype，GATK的解释是uninformative read仍然被用于计算genotype，额，好吧，刚说不可信，你们就又拿来用了，我倾向于认为这是一种后备策略。

另外，文档刚开始贴出了INFO/DP，FORMAT/DP和FORMAT/AD的解释，说FORMAT/DP不包含uninformative read，并注明，The sum of all AD = DP，这里大概率写错了，毕竟文章的缘由就是AD之和与DP不匹配，而且给出的示例中，AD之和为0，和DP就不相等。评论区的Burak Alver也提到了该问题，然而作者没有回复。

影响

了解其成因后，我想看一下这种情况在我的数据上有多普遍。

使用bcftools获取AD之和不等于FORMAT/DP的突变，提取对应信息。

$ bcftools query -f '[%CHROM\_%POS\_%REF\_%ALT %SAMPLE %GT %DP %AD{0} %AD{1}\n]' \
    -i '(GT="alt")&(FORMAT/DP>=20)&(N_ALT=1)&((AD[:0]+AD[:1])!=FORMAT/DP)' temp.vcf.gz > temp_unequal
$ bcftools query -f '[%CHROM\_%POS\_%REF\_%ALT %SAMPLE %GT %DP %AD{0} %AD{1}\n]' \
    -i '(GT="alt")&(FORMAT/DP>=20)&(N_ALT=1)&((AD[:0]+AD[:1])=FORMAT/DP)' temp.vcf.gz > temp_equal
$ wc -l temp_unequal
278773 temp_unequal
$ wc -l temp_equal
281727 temp_equal

计算AD的和没找到好方法，所以只能逐个相加，如果写成sum(AD)就是计算所有sample的AD和了。逐个相加的话，必须同时设置N_ALT=1，不考虑multiallelic位点。

可以看到，不相等的突变还是挺多的。为了更直观的感受一下，读取做个分析。

import pandas as pd
data = pd.read_csv('temp_unequal',sep=' ',names = ['variant','sample','GT','DP','RD','AD'])
data['uninformative'] = data['DP'] - data['RD'] - data['AD']
(data['uninformative'].value_counts()/data.shape[0]).head(10)
---
1    0.331470
2    0.187870
3    0.114312
4    0.085080
5    0.063539
6    0.051024
7    0.040119
8    0.030627
9    0.022183
10   0.016121

前10逐渐下降，1占比1/3，最多，符合预期。

data.shape[0]
data[data['uninformative']/data['DP']>0.2].shape[0]
data[data['uninformative']/data['DP']>0.5].shape[0]
---
278773
42063
2461

占比DP方面，多数情况下，占比较低。

新的考虑：VAF如何计算

虽然大部分位点的uninformative read数都很少，但是也有一些位点，一半以上都是uninformative read。

这种情况，对分析有什么影响呢？

通常，突变下游分析会采用一些硬标准来过滤，比如VAF，之前，我的计算方式是VAF=AD[1]/DP。而在uninformative read占比很大的位点，使用VAF=AD[1]/DP，与VAF=AD[1]/sum(AD)的差距就会比较悬殊了。

使用单样本vcf看一下两者的区别。

$ bcftools view -i 'GT="alt" & FORMAT/AD[:1]/FORMAT/DP>0.2' -H one_sample.vcf.gz | wc -l
61825
$ bcftools view -i 'GT="alt" & FORMAT/AD[:1]/sum(AD)>0.2' -H one_sample.vcf.gz | wc -l
68941

按照GATK的说法，uninformative read是不可信的，那，理论上，VAF计算为AD[1]/sum(DP)更可靠。

然而，在multiallelic位点会出问题。

$ bcftools norm -m-both -f ${hg19} one_sample.vcf.gz -Oz -o one_sample_split.vcf.gz # split multiallelic sites into biallelic records
$ tabix one_sample_split.vcf.gz
$ bcftools view -r 22:16349288 -H one_sample.vcf.gz
22      16349288        rs142693008     GAC     G,GACAC 14642.6 .       AC=0,1;AF=0.096,0.127;AN=2;BaseQRankSum=0.207;DB;DP=14134;ExcessHet=65.5628;FS=1.621;InbreedingCoeff=-0.3018;MLEAC=42,60;MLEAF=0.094,0.134;MQ=54.75;MQRankSum=-2.245;QD=2.07;ReadPosRankSum=0.282;SOR=0.898     GT:AD:DP:GQ:PGT:PID:PL:PS       0/2:73,11,22:124:99:.:.:422,334,3231,0,2219,2806:.
$ bcftools view -r 22:16349288 -H one_sample_split.vcf.gz
22      16349288        rs142693008     GAC     G       14642.6 .       AC=0;AF=0.096;AN=2;BaseQRankSum=0.207;DB;DP=14134;ExcessHet=65.5628;FS=1.621;InbreedingCoeff=-0.3018;MLEAC=42;MLEAF=0.094;MQ=54.75;MQRankSum=-2.245;QD=2.07;ReadPosRankSum=0.282;SOR=0.898      GT:AD:DP:GQ:PGT:PID:PL:PS       0/0:73,11:124:99:.:.:422,334,3231:.
22      16349288        rs142693008     G       GAC     14642.6 .       AC=1;AF=0.127;AN=2;BaseQRankSum=0.207;DB;DP=14134;ExcessHet=65.5628;FS=1.621;InbreedingCoeff=-0.3018;MLEAC=60;MLEAF=0.134;MQ=54.75;MQRankSum=-2.245;QD=2.07;ReadPosRankSum=0.282;SOR=0.898      GT:AD:DP:GQ:PGT:PID:PL:PS       0/1:73,22:124:99:.:.:422,0,2806:.

原本，这个位点sum(AD)=73+11+22=106，而split之后，sum(AD)就丢失了支持另外一个allele的informative reads数，如果采用VAF=AD[1]/sum(AD)，其值就与真实情况相差甚远了。

解决这个问题，最好的办法是通过AD重新计算DP，然后延用之前的VAF=AD[1]/DP。

在计算方面，bcftools通过plugin对此提供了支持，fill-tags for DP[3]。

$ bcftools +fill-tags one_sample.vcf.gz -Oz -o recalculated_DP.vcf.gz -- -t 'DP:1=int(sum(FORMAT/AD))' # bcftools=1.15测试
$ bcftools query -f '[%CHROM\_%POS\_%REF\_%ALT %SAMPLE %GT %DP %AD{0} %AD{1}\n]' \
    -i '(GT="alt")&(FORMAT/DP>=20)&(N_ALT=1)&((AD[:0]+AD[:1])!=FORMAT/DP)' recalculated_DP.vcf.gz | wc -l
10539
$ bcftools query -f '[%CHROM\_%POS\_%REF\_%ALT %SAMPLE %GT %INFO/DP %AD{0} %AD{1}\n]' \
    -i '(GT="alt")&(INFO/DP>=20)&(N_ALT=1)&((AD[:0]+AD[:1])!=INFO/DP)' recalculated_DP.vcf.gz | wc -l
0

美中不足的是，这里更新的是INFO/DP，对于单样本vcf而言，更新INFO/DP足够了，然而，考虑到多样本vcf，更新FORMAT/DP才是正确的做法。

经过测试，目前并不能指定FORMAT/DP，因此，我在作者此功能的commit下进行了评论[4]，希望作者能够进行更进一步的改进。