疑问:
当我想要通过VCF文件查看一个突变位点的变异reads支持数时,发现文件展示的DP (total depth) 和AD (allele depth)的数量对应不上是咋回事儿?比方说,ADs中记录的支持ref和alt的总reads数和DP中展示的总深度对应不起来,更奇怪的是,有时候一个位点的AD竟然是0,HaplotypeCaller为什么会把这样的变异call出来?有图有真相:
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
能看出来,在FORMAT这列(最后一列)显示的支持ref和alt的AD都是0,但是INFO列和FORMAT列中的DP都是10。因为INFO列的DP是未经过过滤的total depth,FORMAT列的DP是经过过滤的total depth,因此可以推断这个位置上没有任何一个read被GATK内置的过滤程序过滤掉。再进一步查看“bamout”文件,就会发现确实在这个位置上有10个reads的覆盖。所以为什么VCF报告的AD值都是0?
原因:uninformative reads
其实这个并不能算一个bug,问题出在uninformative reads上。
当该read给出的most likely allele的可能性并没有显著大于second most likely allele时,这个read就被称为uninformative reads。需要说明的是,该可能性转换为Phred格式(Phred scaled likelihoods)后需要大于0.2才会被软件认为是具有显著差异。换句话说,most likely allell可能性必须比second most likely allele的可能性大60%才行(因为0.2 = 10^0.2 = 1.585,约等于大60%)。
例如,我现在有两条reads和两个可能的等位基因,所有的reads都通过了HaplotypeCaller的过滤,这连个等位基因的可能性的值展示在下图中:
现在把上图中的值转换为Phred格式,即取log值:
现在,我们想看看read 1是否是imformative read,只需要看看most likely allele (A)和second most likely allele (T)的可能性差值,即-6.4122-(-6.4352)= 0.023,由于0.023少于0.2,read 1由此可判断为uninformative。同理可以判断出read 2 是informative。
小结:
所以,如果一个read 被判定为informative,它将记录在AD和DP中去。如果是uninformative,它只会记录在DP中,而不记录在AD中。AD实际上反应了究竟有多少reads支持了该位置基因型的判定,之所以没有把uninformative记录进去,就是因为这些read的提供的基因型判定证据不足。
参考:
https://software.broadinstitute.org/gatk/documentation/article?id=11096