拿到snp的rawdata后如何使用GATK进行筛选（GATK硬筛选文档翻译）

这篇文档旨在提供一些关于我们通常推荐的硬性过滤逻辑的见解，作为 VQSR（我们通常推荐用于过滤生殖系短变异的方法）的替代方法。希望这篇指南能够帮助您适应这些推荐，或者开发出适用于与我们常用数据集显著不同的数据集的新过滤方法。

需要注意的是，本文中描述的注释是针对在群体水平上调用的样本，用以过滤掉那些不可信赖的数据位点，即使使用 VQSR，过滤的也是变异位点，而不是单个样本的数据。

有关如何实际应用硬性过滤的说明，请参阅 JEXL 过滤表达式和 SelectVariants 及 VariantFiltration 工具文档。

1. 概述

硬性过滤是指为一个或多个注释选择特定的阈值，并剔除任何注释值超出该阈值的变异。注释是指描述每个变异的属性或统计信息，例如该变异位点周围的序列上下文如何、有多少读取覆盖了该位点、每个等位基因有多少读取覆盖、前向与反向方向上读取的比例等。

这种方法的局限性在于，它迫使我们单独考虑每个注释维度，最终可能因为一个注释表现不好而丢弃好的变异，或者为了保留好的变异而保留不好的变异。

相比之下，VQSR 更加强大，因为它使用机器学习算法从数据中学习出好变异（真阳性）和坏变异（假阳性）的注释特征。这使得我们能够根据不同维度上的聚类情况筛选变异，摆脱单一维度阈值的限制。

然而，这种方法需要大量变异和高质量的已知变异资源。对于使用小型基因面板或非模式生物的研究人员来说，这可能无法实现，因此只能依赖硬性过滤。

2. 可视化注释值的分布

一种帮助理解硬性过滤的方法是可视化特定管线调用出的真集（truth set）中的注释值分布。这些分布由管线方法学和序列数据的物理特性共同塑造；因此，对于某种数据生成技术和分析管线的组合，您可以根据真集的分布推导出过滤阈值。

可视化数据的来源

我们对一个完整基因组的三重体（样本 NA12878、NA12891 和 NA12892，之前已预处理）进行了变异调用，使用的是 HaplotypeCaller 的 GVCF 模式，生成了每个样本的 GVCF 文件。然后，我们使用 GenotypeGVCFs 对这些 GVCF 进行了联合调用，生成了三重体的未过滤 VCF 集合。最后，我们对三重体 VCF 运行了 VQSR，得到过滤后的集合。我们将通过的变异视为真变异，而过滤掉的变异视为假变异。注意，这个过程也可以通过使用 NA12878 的 Genome in a Bottle 真集作为资源对原始集合进行注释来实现（尽管结果不完全相同，但足够相似）。

接着，我们使用 VariantsToTable 提取了所有调用的注释值，并将其加载到 RStudio 进行绘图。

绘图方法和解释说明

我们为六个推荐的注释绘制了图，这些注释通常具有高度信息性：QD、FS、SOR、MQ、MQRankSum 和 ReadPosRankSum。同样的原则适用于 GATK 工具生成的大多数其他注释。

所有显示的图表都是使用 R 中的 ggplot2 库生成的密度图。x 轴是注释值，y 轴是密度值。密度图下的面积表示观察到注释值的概率。因此，图表下的总面积等于 1。如果您想知道在 0 到 1 之间观察到注释值的概率，则需要计算该区间下曲线下的面积。

简而言之，这些图显示了对于特定的变异集，注释值的分布是什么样的。需要注意的是，当我们在同一图表上比较两个或多个变异集时，需要记住它们包含的变异数量可能非常不同，因此在分布中某一部分的变异数量是不可比的；只有它们的比例是可比的。

3. QualByDepth (QD)

1. 什么是 QD？

QD 是一个注释，表示变异的置信度（从 QUAL 字段得出）除以非纯合参考（non-hom-ref）样本的未过滤深度。它的作用是对变异的质量进行归一化，避免由于深度覆盖过高而导致的置信度膨胀。

简单来说，QD 是对变异质量进行调整，以避免高覆盖率的数据人为增加变异的置信度。比直接使用 QUAL 或 DP 更好，因为它考虑了数据覆盖的变化情况。

2. `为什么要过滤 QD 低于 2 的变异？`

在硬性过滤中，常规推荐是过滤掉 QD 低于 2 的变异。让我们来看看原因：

QD 值的分布：

对于未经过滤的变异，QD 值可以在 0 到 40 之间分布。大多数变异的 QD 值分布在两个峰之间：大约在 QD = 12 和 QD = 32。这两个峰分别对应于在调用的样本中，异质性（杂合子，heterozygous）和纯合变异（homozygous-variant）的变异。这是因为纯合变异的样本提供的支持变异的读取数是杂合变异样本的两倍。

在分布的左边，可以看到 QD 在 0 到 5 之间的“肩部”区域。低 QD 值表明该变异位点可能存在质量问题。

过滤 QD 小于 2 的原因：

我们看到，VQSR（变异质量评分再校正）过滤掉的变异大多集中在 QD 低于 5 的区域。这意味着这些变异在数据中表现出较低的置信度或被认为是伪影。通常，建议使用 QD = 2 作为硬性过滤的阈值，以剔除质量较差的变异。虽然有些变异的 QD 值低于 2，但仍可能通过 VQSR，因此我们需要在硬性过滤中做出权衡。

3. 如何解读 QD 过滤中的异常情况？

虽然 QD 低于 2 的变异通常被认为质量较差，但有时我们会发现 VQSR 过滤掉了一些 QD 大于 30 的变异！这意味着尽管这些变异通过了硬性过滤的 QD 阈值，但在其他注释维度上表现不佳，被视为伪影。

同样，VQSR 也可能通过一些 QD 小于 2 的变异，说明这些变异可能在其他维度上表现良好。如果我们仅根据 QD 进行硬性过滤，可能会丢失这些变异。

总结：

- **QD 的作用**：用于归一化变异的质量，避免深度覆盖过高导致置信度膨胀。
- **过滤阈值**：通常推荐过滤掉 QD 低于 2 的变异，以移除质量较差的位点。
- **注意事项**：硬性过滤只能根据单一维度判断，可能会丢失一些在其他维度上表现良好的变异。

4. FisherStrand (FS)

1. 什么是 FS？

**FS（Fisher Strand）** 是一个 Phred 标度的概率值，用来评估是否在位点上存在链偏性（strand bias）。链偏性指的是备用等位基因在正向链（forward strand）和反向链（reverse strand）上出现的频率是否不同于参考等位基因。如果在正向和反向链上看到备用等位基因的频率明显不同，这可能表明该位点的数据有问题。

FS 通过 Fisher 精确检验（Fisher’s Exact Test）计算得出，这是一种用于处理小样本量数据的统计检验方法。

FS 的解读：

当位点没有或很少存在链偏性时，FS 值接近 0。
如果 FS 值较大，则表示该位点可能存在显著的链偏性，这通常是序列错误或测序技术问题的表现。

2. FS 值的分布

我们来看看 FS 值在未经过滤的变异中的分布：

FS 值的范围很广，有些变异的 FS 值接近 400。然而，大多数变异的 FS 值低于 10。
在未过滤的变异集中，几乎所有变异的 FS 值都低于 100。

3. `为什么推荐过滤掉 FS 大于 60 的变异？`

在硬性过滤中，我们通常建议过滤掉 FS 值大于 60 的变异。让我们来看看原因：

VQSR 过滤后的分布：

在通过 VQSR 的变异和未通过 VQSR 的变异中，我们可以看到，大多数未通过 VQSR 的变异的 FS 值都大于 55。换句话说，这些变异有显著的链偏性，是数据不可信的信号。
基于此，硬性过滤时建议过滤掉 FS 值大于 60 的变异，以去除具有显著链偏性的假阳性变异。

不过要注意，虽然将 FS 阈值设置为 60 能过滤掉许多假阳性变异，但仍会保留一些假阳性。如果我们将阈值设得更低（如低于 55），虽然能移除更多的假阳性，但也可能会失去一些真阳性。

4. FS 与相关注释的区别

在评估链偏性时，除了 FS 之外，还有两个相关的注释：SB 和 SOR。

SB（Strand Bias）：提供每个等位基因在正向和反向链上支持的原始读取数。
SOR（Strand Odds Ratio）：使用另一种统计检验方法来评估链偏性。SOR 更适合处理高覆盖度的数据。

与 FS 不同，SB 是原始读取计数，而 FS 是基于 Fisher 精确检验的概率值。FS 会受到测序深度的影响，在某些情况下可能会过度惩罚那些位于外显子末端的变异，而 SOR 则可以更好地处理这种情况。

总结：

FS 的作用：用于检测链偏性，FS 值越大表示偏性越强，变异的可信度越低。
过滤阈值：通常建议过滤掉 FS 值大于 60 的变异。
注意事项：FS 是检测链偏性的一种常用方法，但在高覆盖度数据中，SOR 可能更为合适。

5. StrandOddsRatio (SOR)

1. 什么是 SOR？

SOR（Strand Odds Ratio） 是一种用来评估链偏性的新方法，使用类似对称优势比检验（symmetric odds ratio test）的方式进行估算。链偏性（strand bias）是指在正向和反向链上看到备用等位基因的比例差异较大，这可能表明数据质量存在问题。

SOR 之所以被引入，是因为 FS（Fisher Strand） 注释有时会过度惩罚那些发生在外显子末端的变异。在这些位置，常常只能在某个方向的链上覆盖到读取数，导致 FS 给出的结果显示出偏性。SOR 考虑了覆盖到两种等位基因的读取数的比例，从而更好地处理这种情况。

SOR 的特点：

SOR 值越高，表示链偏性越强。
SOR 是针对高覆盖度数据设计的，因此相比 FS，它在覆盖较高的情况下表现更好。

2. SOR 值的分布

让我们来看看未经过滤的变异集中 SOR 值的分布：

SOR 值的范围从 0 到大于 9，大多数变异的 SOR 值小于 3，几乎所有的变异 SOR 值都低于 9。
我们可以看到 SOR 值较大的变异数量较少，但存在一个长尾，SOR 值可以超过 9。

3. `为什么推荐过滤掉 SOR 大于 3 的变异？`

在硬性过滤中，我们通常推荐过滤掉 SOR 值大于 3 的变异。让我们来看看原因：

VQSR 过滤后的分布：

在通过 VQSR 的变异和未通过 VQSR 的变异中，我们看到大多数 SOR 值大于 3 的变异都未通过 VQSR 过滤。
尽管有些 SOR 值小于 3 的变异也未能通过 VQSR，但硬性过滤时使用 3 作为阈值，可以去除那些显示出显著偏性的长尾变异。

这个阈值的选择是基于观察到的注释值分布。设定 SOR = 3 的阈值可以移除表现出较强链偏性的变异，而不会过多地丢失潜在的好变异。

4. FS 与 SOR 的区别

尽管 FS 和 SOR 都是用来评估链偏性，但它们的适用场景有所不同：

FS：基于 Fisher 精确检验，适用于低覆盖度的数据，但在高覆盖度数据中可能会过度惩罚某些区域的变异。
SOR：使用优势比检验，专为高覆盖度数据设计，避免了 FS 在处理外显子末端变异时的不足。

因此，如果数据覆盖度较高，SOR 通常比 FS 更能准确反映变异的真实性。

总结：

SOR 的作用：用于评估链偏性，SOR 值越高，偏性越强，变异的可信度越低。
过滤阈值：通常建议过滤掉 SOR 值大于 3 的变异，以去除那些显示出显著偏性的假阳性。
注意事项：SOR 是针对高覆盖度数据设计的，通常比 FS 更适合处理这种数据。

6. RMSMappingQuality (MQ)

1. 什么是 MQ？

MQ（Root Mean Square Mapping Quality） 是一个注释，表示某个位点上所有读取的根均方（RMS）比对质量。与直接计算平均比对质量不同，MQ 通过计算每个读取比对质量的平方平均值再开方，以反映该位点的总体比对质量。这个值可以帮助我们评估该变异位点的整体比对准确性。

MQ 的特点：

MQ 值越高，表明在该位点上，读取的比对质量越好。
理想情况下，如果位点上的比对质量良好，MQ 值应该接近 60。

2. MQ 值的分布

我们可以看看未经过滤的变异的 MQ 值分布：

大多数变异的 MQ 值集中在大约 60 左右，显示出高质量的比对。
通常我们建议在硬性过滤时，过滤掉 MQ 值小于 40 的变异。这是因为低于 40 的 MQ 值表明在该位点上的比对质量较差。

3. `为什么推荐过滤掉 MQ 小于 40 的变异？`

MQ 反映了一个变异位点的比对质量，因此较低的 MQ 值通常表明该位点的数据不够可信。在硬性过滤中，推荐的阈值是过滤掉 MQ 值低于 40 的变异。

VQSR 过滤后的分布：

在通过 VQSR 的变异和未通过 VQSR 的变异中，几乎所有通过 VQSR 的变异的 MQ 值都在 60 附近，而未通过 VQSR 的变异通常有较低的 MQ 值。
通过使用 MQ = 40 作为硬性过滤的阈值，可以去除那些比对质量不佳的变异。

尽管有些研究者可能会选择更严格的阈值（例如过滤掉 MQ 小于 50 的变异），但我们推荐的阈值是较为宽松的 40，以保留更多可能是可信的变异。

4. 如何解释高 MQ 但被 VQSR 过滤的变异？

虽然大多数通过 VQSR 的变异的 MQ 值接近 60，但有些未通过 VQSR 的变异也有接近 60 的 MQ 值。这表明，即使比对质量良好，其他维度上的注释可能显示出这些变异存在问题。因此，仅依靠 MQ 进行硬性过滤无法捕捉所有伪影。

总结：

MQ 的作用：用于评估变异位点的比对质量，MQ 值越高，表示比对质量越好。
过滤阈值：通常建议过滤掉 MQ 值小于 40 的变异。
注意事项：虽然 MQ 能反映比对质量，但无法单独确定变异的真实性，因此应结合其他注释一起分析。

7. MappingQualityRankSumTest (MQRankSum)

1. 什么是 MQRankSum？

MQRankSum 是一个基于 u 值的 z 近似值，来自比对质量的秩和检验（Rank Sum Test）。它用于比较支持参考等位基因和支持备用等位基因的读取的比对质量。如果备用等位基因的读取比对质量明显低于参考等位基因，那么该变异可能不可靠。

MQRankSum 的解读：

正值：表示支持备用等位基因的读取的比对质量高于支持参考等位基因的读取。
负值：表示支持参考等位基因的读取的比对质量高于支持备用等位基因的读取。
接近 0 的值最好，表明支持参考和备用等位基因的读取在比对质量上没有明显差异。

2. MQRankSum 值的分布

让我们看看未经过滤的变异的 MQRankSum 值的分布：

MQRankSum 的值通常在 -10.5 到 6.5 之间分布。我们的推荐硬性过滤阈值是 -12.5。
在我们的数据集中，没有任何变异的 MQRankSum 小于 -10.5。这意味着硬性过滤使用 -12.5 作为阈值将不会移除任何变异。

3. `为什么选择 -12.5 作为 MQRankSum 的硬性过滤阈值？`

硬性过滤推荐的阈值是 -12.5，这个阈值相对宽松，意在保留大多数潜在的好变异。一般来说，如果您观察到 MQRankSum 值低于 -12.5，说明备用等位基因的比对质量比参考等位基因差得多，可能是错误的变异。因此，使用这个阈值来过滤掉这些可能的假阳性变异是合理的。

VQSR 过滤后的分布：

在通过 VQSR 的变异和未通过 VQSR 的变异中，负值较低（例如小于 -2.5）的变异往往未通过 VQSR 过滤，而值接近 0 的变异通常会通过 VQSR。
这表明，尽管选择 -12.5 作为阈值不会过滤掉我们数据集中的任何变异，但将阈值设置得更严格（如 -2.5）可能会移除更多不可信的变异。

4. 如何调整 MQRankSum 的阈值？

每个数据集的注释分布可能不同，因此在使用硬性过滤时，建议根据您特定的数据集调整阈值。虽然 -12.5 是推荐的宽松阈值，但如果在您的数据集中发现很多 MQRankSum 小于 -2.5 的变异不可信，您可以考虑将阈值调整为更严格的 -2.5。

总结：

MQRankSum 的作用：比较支持参考等位基因和备用等位基因的读取的比对质量，判断它们之间是否有显著差异。
过滤阈值：通常建议使用 -12.5 作为硬性过滤阈值，但可以根据数据集进行调整。
注意事项：该注释主要用于检测备用等位基因和参考等位基因在比对质量上的显著差异，因此结合其他注释进行判断更为可靠。

这就是关于 MappingQualityRankSumTest (MQRankSum) 的详细介绍。如果没有问题，我们可以继续讲解最后一个注释 ReadPosRankSumTest (ReadPosRankSum)。
很好！接下来我们讲解最后一个注释 ReadPosRankSumTest (ReadPosRankSum)。

8. ReadPosRankSumTest (ReadPosRankSum)

1. 什么是 ReadPosRankSum？

ReadPosRankSum 是一个基于 u 值的 z 近似值，来自位点位置的秩和检验（Rank Sum Test）。它用于比较读取中的参考等位基因和备用等位基因在读取中的位置差异。如果备用等位基因的读取位置明显集中在读取的末端，而参考等位基因的读取分布均匀，则该变异可能是不可靠的。

ReadPosRankSum 的解读：

负值：表示备用等位基因的读取更常出现在读取的末端位置。因为在读取的末端区域，测序错误更为常见，所以这通常是不好的信号。
正值：表示参考等位基因的读取更多出现在读取的末端位置。
接近 0 的值是最佳情况，表明参考等位基因和备用等位基因在读取中的位置分布相似，数据质量较为可信。

2. ReadPosRankSum 值的分布

我们来看一下未经过滤的变异集中 ReadPosRankSum 值的分布：

未经过滤的变异的 ReadPosRankSum 值大多数在 -4 到 4 之间。
硬性过滤的推荐阈值是过滤掉 ReadPosRankSum 值小于 -8 的变异。然而在我们的数据集中，没有任何变异的 ReadPosRankSum 值小于 -8。

3. `为什么推荐过滤掉 ReadPosRankSum 小于 -8 的变异？`

负值：负的 ReadPosRankSum 值表明备用等位基因的读取位置通常在读取的末端，这在测序中容易引入错误。因此，负的 ReadPosRankSum 值较低的变异（如小于 -8）可能是不可信的变异。
阈值设定：虽然在我们的数据集中没有值小于 -8 的变异，但在不同数据集上，使用 -8 作为阈值可以过滤掉那些备用等位基因读取集中在读取末端的伪变异。

VQSR 过滤后的分布：

当我们查看通过 VQSR 的变异和未通过 VQSR 的变异时，可以看到，VQSR 过滤掉的变异通常具有较低的 ReadPosRankSum 值（小于 -1.0）或较高的值（大于 3.5）。
硬性过滤使用 -8 作为阈值，可以保留大部分可信的变异，同时移除那些位置分布明显异常的变异。

4. 如何调整 ReadPosRankSum 的阈值？

正如其他注释一样，您可以根据数据集的具体情况调整 ReadPosRankSum 的阈值。例如，如果在您的数据集中有许多 ReadPosRankSum 值较高或较低的变异表现不佳，您可以将阈值调整得更严格，过滤掉更多不可信的变异。

总结：

ReadPosRankSum 的作用：用于比较参考和备用等位基因在读取中的位置分布，负值通常表示备用等位基因更多出现在读取末端，可能是伪变异。
过滤阈值：通常建议使用 -8 作为过滤阈值，但您可以根据数据集的具体情况调整这一阈值。
注意事项：ReadPosRankSum 注释帮助识别备用等位基因是否集中在读取的末端，这是一个测序错误可能集中的区域，因此结合其他注释一起使用可以提高变异的鉴别能力。