「简化基因组」如何过滤用GATK分析得到的SNP

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的情况

GBS比对情况

他们的比对位置并不随机，因此任何和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