展开
全基因组(WGS) 思维导图
标签
  1. 首页
  2. 组学教程
  3. 全基因组(WGS)
  4. VCF质控

VCF质控

最后发布时间 : 2023-08-05 15:49:32 浏览量 :

学习资料

什么是质控？我们不妨先给它下个定义：质控的含义和目的是指通过一定的标准，最大可能地剔除假阳性的结果，并尽可能地保留最多的正确数据。有了这么一个定义之后，我们也就能够更加清晰地知道接下来该做些什么了。

在上次的文章里我已经说到在GATK HaplotypeCaller之后，首选的质控方案是GATK VQSR，它通过机器学习的方法利用多个不同的数据特征训练一个模型（高斯混合模型）对变异数据进行质控，然而不幸的是使用VQSR需要具备以下两个条件：

第一，需要一个精心准备的已知变异集，它将作为训练质控模型的真集。比如，对于我们人来说，就有Hapmap、OMNI，1000G和dbsnp等这些国际性项目的数据，这些可以作为高质量的已知变异集。GATK的bundle主要就是对这四个数据集做了精心的处理和选择，然后把它们作为VQSR时的真集位点。这里我强调一个地方：是真集的『位点』而不是真集的『数据』！还请大家多多注意。因为，VQSR并不是用这些变异集里的数据来训练的，而是用我们自己的变异数据。这个对于刚接WGS的同学来说特别容易搞混，不要因为VQSR中用了那四份变异数据，就以为是用它们的数据来训练模型。

实际上，这些已知变异集的意义是告诉我们群体中哪些位点存在着变异，如果在其他人的数据里能观察到落入这个集合中的变异位点，那么这些被已知集包括的变异就有很大的可能是正确的。也就是说，我们可以从数据中筛选出那些和真集『位点』相同的变异，把它们当作是真实的变异结果。接着，进行VQSR的时候，程序就可以用这个筛选出来的数据作为真集数据来训练，并构造模型了。

第二，要求新检测的结果中有足够多的变异，不然VQSR在进行模型训练的时候会因为可用的变异位点数目不足而无法进行。

由于条件1的限制，会导致很多非人的物种在完成变异检测之后没法使用GATK VQSR的方法进行质控。而由于条件2，也常常导致一些小panel甚至外显子测序，由于最后的变异位点不够，也无法使用VQSR。这个时候，我们就不得不选择硬过滤的方式来质控了。

那什么叫做硬过滤呢？所谓硬过滤其实就是通过人为设定一个或者若干个指标阈值（也可以叫数据特征值），然后把所有不满足阈值的变异位点采用一刀切掉的方法。

那么如何执行硬过滤？首先，需要我们确定该用哪些指标来评价变异的好坏。这个非常重要，选择对了事半功倍，选得不合理，过滤的结果有时还不如不过滤的。如果把这个问题放在从前，我们需要做比较多的尝试才能确定一些合适的指标，但现在就方便很多了，可以直接使用GATK VQSR所用的指标——毕竟这些指标都是经过精挑细选的。我想这应该不难理解，既然VQSR就是用这些指标来训练质控模型的，那么它们就可以在一定程度上描述每个变异的质量，我们用这些指标设置对应的阈值来进行硬过滤也将是合理的。VQSR使用的数据指标有6个（这些指标都在VCF文件的INFO域中，如果不是GATK得到的变异，可能会有所不同，但知道它们的含义之后也是可以自己计算的），分别是：

  • QualByDepth（QD）
  • FisherStrand (FS)
  • StrandOddsRatio (SOR)
  • RMSMappingQuality (MQ)
  • MappingQualityRankSumTest (MQRankSum)
  • ReadPosRankSumTest (ReadPosRankSum)

指标有了，那么阈值应该设置为多少？下面我想先给出一个硬过滤的例子，然后再逐个来对其进行分析，以便大家能够更好地理解变异质控的思路。值得注意的是不同的数据，有不同的情况，它的阈值有时是不同的。不过不用担心，当你掌握了如何做的思路之后完全有能力根据具体的情况举一反三。

执行硬过滤

首先是硬过滤的例子，这个过程我都用最新的GATK来完成。GATK 4.0中有一个专门的VariantFiltration模块（继承自GATK 3.x），它可以很方便地帮我们完成这个事情。不过，过滤的时候，需要分SNP和Indel这两个不同的变异类型来进行，它们有些阈值是不同的，需要区别对待。

# 使用SelectVariants，选出SNP
time /Tools/common/bin/gatk/4.0.1.2/gatk SelectVariants \
    -select-type SNP \
    -V ../output/E.coli/E_coli_K12.vcf.gz \
    -O ../output/E.coli/E_coli_K12.snp.vcf.gz

# 为SNP作硬过滤
time /Tools/common/bin/gatk/4.0.1.2/gatk VariantFiltration \
    -V ../output/E.coli/E_coli_K12.snp.vcf.gz \
    --filter-expression "QD < 2.0 || MQ < 40.0 || FS > 60.0 || SOR > 3.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" \
    --filter-name "Filter" \
    -O ../output/E.coli/E_coli_K12.snp.filter.vcf.gz

# 使用SelectVariants，选出Indel
time /Tools/common/bin/gatk/4.0.1.2/gatk SelectVariants \
    -select-type INDEL \
    -V ../output/E.coli/E_coli_K12.vcf.gz \
    -O ../output/E.coli/E_coli_K12.indel.vcf.gz

# 为Indel作过滤
time /Tools/common/bin/gatk/4.0.1.2/gatk VariantFiltration \
    -V ../output/E.coli/E_coli_K12.indel.vcf.gz \
    --filter-expression "QD < 2.0 || FS > 200.0 || SOR > 10.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" \
    --filter-name "Filter" \
    -O ../output/E.coli/E_coli_K12.indel.filter.vcf.gz

# 重新合并过滤后的SNP和Indel
time /Tools/common/bin/gatk/4.0.1.2/gatk MergeVcfs \
    -I ../output/E.coli/E_coli_K12.snp.filter.vcf.gz \
    -I ../output/E.coli/E_coli_K12.indel.filter.vcf.gz \
    -O ../output/E.coli/E_coli_K12.filter.vcf.gz

# 删除无用中间文件
rm -f ../output/E.coli/E_coli_K12.snp.vcf.gz* ../output/E.coli/E_coli_K12.snp.filter.vcf.gz* ../output/E.coli/E_coli_K12.indel.vcf.gz* ../output/E.coli/E_coli_K12.indel.filter.vcf.gz*

最后，只要符合了上面任意一个阈值的变异都会被设置为“Filter”，剩下的会被认为是正常的变异，并标记为“PASS”。流程的最后，我们需要把分开质控的SNP和Indel结果重新合并在一起，然后再把那些不必要的中间文件删除掉。

在具体的项目中，你如果需要使用硬过滤的策略，这个例子中的参数可以作为参考，特别是对于高深度数据而言。接下来我结合GATK所提供的资料与大家分享如何理解这些指标以及得出这些阈值的思路。

如何理解硬过滤的指标和阈值的计算

为了更好的理解，这里使用上一节生成的文件output/vcf/mother.sorted.markdup.vcf和下面多个样本生成的vcf文件(output/vcf/mutiple_sample.vcf)内容作为解释。完整的VCF文件

gatk HaplotypeCaller \
    -R output/index/ref.fasta \
    -I 2-germline/bams/mother.bam \
    -I 2-germline/bams/father.bam \
    -O output/vcf/mutiple_sample.vcf \
    -L 20:10,000,000-10,200,000

bcftools  view output/vcf/mutiple_sample.vcf  | less -S
bcftools  view  output/vcf/mother.sorted.markdup.vcf  | less -S

QualByDepth（QD）

QD是变异质量值（Quality）除以覆盖深度（Depth）得到的比值。这里的变异质量值就是VCF中QUAL的值——用来衡量变异的可靠程度，这里的覆盖深度是这个位点上所有含有变异碱基的样本的覆盖深度之和，通俗一点说，就是这个值可以通过累加每个含变异的样本（GT为非0/0的样本）的覆盖深度（VCF中每个样本里面的DP）而得到。举个例子：

1     1429249 .       C       T       1044.77  .  .  GT:AD:DP:GQ:PL  0/1:48,15:63:99:311,0,1644  0/0:47,0:47:99:392,0,0  1/1:0,76:76:99:3010,228,0
  • 这个位点是1:1429249，VCF格式，但我把FILTER和INFO的信息省略了；
  • 它的变异质量值QUAL=1044.77；
  • 我们可以从中看到一共有三个样本，其中一个是杂合变异（GT=0/1），一个纯合的非变异（GT=0/0），最后一个是纯合的变异（GT=1/1）；
  • 每个样本的覆盖深度都在其各自的DP域上，分别是：63，47和76；
  • 按照定义，这个位点的QD值就应该等于质量值除以另外两个含有变异的样本的深度之和（排除中间GT=0/0这个不含变异的样本），也就是：
QD = 1044.77 / (63+76) = 7.516

变异检测质控和过滤（VQSR）

这是我们这个流程中最后的一步了。在获得了原始的变异检测结果之后，我们还需要做的就是质控和过滤。这一步或多或少都有着一些个性化的要求，我暂时就不做太多解释吧（一旦解释恐怕同样是一篇万字长文）。只用一句话来概括，VQSR是通过构建GMM模型对好和坏的变异进行区分，从而实现对变异的质控，具体的原理暂时不展开了。

## SNP Recalibrator
java -jar /path/to/GenomeAnalysisTK.jar \
   -T VariantRecalibrator \
   -R reference.fasta \
   -input sample_name.HC.vcf \
   -resource:hapmap,known=false,training=true,truth=true,prior=15.0 /path/to/gatk/bundle/hapmap_3.3.b37.vcf \ 
   -resource:omini,known=false,training=true,truth=false,prior=12.0 /path/to/gatk/bundle/1000G_omni2.5.b37.vcf \
   -resource:1000G,known=false,training=true,truth=false,prior=10.0 /path/to/gatk/bundle/1000G_phase1.snps.high_confidence.b37.vcf \ 
   -resource:dbsnp,known=true,training=false,truth=false,prior=6.0 /path/to/gatk/bundle/dbsnp_138.b37.vcf \ 
   -an QD -an MQ -an MQRankSum -an ReadPosRankSum -an FS -an SOR -an DP \ 
   -mode SNP \ 
   -recalFile sample_name.HC.snps.recal \
   -tranchesFile sample_name.HC.snps.tranches \ 
   -rscriptFile sample_name.HC.snps.plots.R

java -jar /path/to/GenomeAnalysisTK.jar -T ApplyRecalibration \
   -R  human_g1k_v37.fasta \
   -input sample_name.HC.vcf \ 
   --ts_filter_level 99.5 \ 
   -tranchesFile sample_name.HC.snps.tranches \ 
   -recalFile sample_name.HC.snps.recal \
   -mode SNP \
   -o sample_name.HC.snps.VQSR.vcf## Indel Recalibratorjava -jar /path/to/GenomeAnalysisTK.jar -T VariantRecalibrator \
   -R  human_g1k_v37.fasta \
   -input sample_name.HC.snps.VQSR.vcf \
   -resource:mills,known=true,training=true,truth=true,prior=12.0 /path/to/gatk/bundle/Mills_and_1000G_gold_standard.indels.b37.vcf \
   -an QD -an DP -an FS -an SOR -an ReadPosRankSum -an MQRankSum \
   -mode INDEL \
   -recalFile sample_name.HC.snps.indels.recal \
   -tranchesFile sample_name.HC.snps.indels.tranches \
   -rscriptFile sample_name.HC.snps.indels.plots.R

java -jar /path/to/GenomeAnalysisTK.jar -T ApplyRecalibration \ 
   -R human_g1k_v37.fasta\
   -input sample_name.HC.snps.VQSR.vcf \
   --ts_filter_level 99.0 \
   -tranchesFile sample_name.HC.snps.indels.tranches \
   -recalFile sample_name.HC.snps.indels.recal \
   -mode INDEL \
   -o sample_name.HC.snps.indels.VQSR.vcf

最后，sample_name.HC.snps.indels.VQSR.vcf 便是我们最终的变异检测结果。对于人类而言，一般来说，每个人最后检测到的变异数据大概在400万左右（包括SNP和Indel）。

GVCF 模式 GATK4 HaplotypeCaller