GATK官方提供了一个SNP过滤的标准,howto-apply-hard-filters-to-a-call-set,如果你按照它的要求来过滤简化基因组中的SNP数据,也就是如下标准
QD > 2.0
FS > 60.0
MQ > 40.0
MQRankSum > -12.5
ReadPosRankSum > -8.0
SOR > 3.0
那么一顿操作之后,你会惊喜的发现,自己数据似乎都不见了。那么原因是什么呢?我们先来理解每个标准的含义
- QualByDepth(QD): 变异位点可信度除以未过滤的非参考read数
- FisherStrand (FS): Fisher精确检验评估当前变异是strand bias的可能性,这个值在0-60间
- RMSMappingQuality (MQ): 所有样本中比对质量的平方根
- MappingQualityRankSumTest (MQRankSum): 根据REF和ALT的read的比对质量来评估可信度
- ReadPosRankSumTest (ReadPosRankSum) : 通过变异在read的位置来评估变异可信度,通常在read的两端的错误率比较高
- StrandOddsRatio (SOR) : 综合评估strand bias的可能性
在解释原因之前,先让我们回顾下一个GBS数据比对后在IGV的情况
他们的比对位置并不随机,因此任何和strand bias有关的标准在过滤时,也就是FS > 60.0 SOR > 3.0
时会过滤掉90%的数据,因此过滤掉许多真实的变异。
因此,官方提供的标准,GBS数据只要用以下几个就行
QD > 2
MQ > 40.0
MQRankSum > -12.5
ReadPosRankSum > -8
当然具体标准,我建议用vcfR
导入VCF文件,通过柱状图分布来确定。
附上我的一批数据通过这些标准过滤的结果
> table(QD>2)
FALSE TRUE
1140 91348
> table(MQ > 40.0)
FALSE TRUE
8924 83565
> table(MQRankSum > -12.5)
TRUE
92465
> table(ReadPosRankSum > -8)
FALSE TRUE
103 92127
> table(FS >= 60.0)
FALSE TRUE
86515 5974
> table(SOR > 3)
FALSE TRUE
85705 6784