今天,只分享一个精选内容。全文3116字,阅读8分钟。
----/ BEGIN /----
变异的质控,是我们在得到变异数据之后,接下来最重要的一个步骤。通常我们都是使用GATK VQSR模块来完成这个事情,关于VQSR的基本原理我在这篇文章中有写,但暂时不算详细。下面是大家经常都会用到的VQSR基本命令(以GATK4为例):
## 首先是SNP mode
time $gatk VariantRecalibrator \
-R $reference/Homo_sapiens_assembly38.fasta \
-V $outdir/poplation/${outname}.HC.vcf.gz \
-resource:hapmap,known=false,training=true,truth=true,prior=15.0 $GATK_bundle/hapmap_3.3.hg38.vcf \
-resource:omini,known=false,training=true,truth=false,prior=12.0 $GATK_bundle/1000G_omni2.5.hg38.vcf \
-resource:1000G,known=false,training=true,truth=false,prior=10.0 $GATK_bundle/1000G_phase1.snps.high_confidence.hg38.vcf \
-resource:dbsnp,known=true,training=false,truth=false,prior=2.0 $GATK_bundle/dbsnp_146.hg38.vcf \
-an DP -an QD -an FS -an SOR -an ReadPosRankSum -an MQRankSum \
-mode SNP \
-tranche 100.0 -tranche 99.9 -tranche 99.0 -tranche 95.0 -tranche 90.0 \
-rscriptFile $outdir/poplation/${outname}.HC.snps.plots.R \
--tranches-file $outdir/poplation/${outname}.HC.snps.tranches \
-O $outdir/poplation/${outname}.HC.snps.recal && \
time $gatk ApplyVQSR \
-R $reference/Homo_sapiens_assembly38.fasta \
-V $outdir/poplation/${outname}.HC.vcf.gz \
--ts_filter_level 99.0 \
--tranches-file $outdir/poplation/${outname}.HC.snps.tranches \
-recalFile $outdir/poplation/${outname}.HC.snps.recal \
-mode SNP \
-O $outdir/poplation/${outname}.HC.snps.VQSR.vcf.gz && echo "** SNPs VQSR done **"
## 然后是Indel mode
time $gatk VariantRecalibrator \
-R $reference/Homo_sapiens_assembly38.fasta \
-input $outdir/poplation/${outname}.HC.snps.VQSR.vcf.gz \
-resource:mills,known=true,training=true,truth=true,prior=12.0 $GATK_bundle/Mills_and_1000G_gold_standard.indels.hg38.vcf \
-resource:dbsnp,known=true,training=false,truth=false,prior=2.0 $GATK_bundle/dbsnp_146.hg38.vcf \
-an DP -an QD -an FS -an SOR -an ReadPosRankSum -an MQRankSum \
-mode INDEL \
--max-gaussians 6 \
-rscriptFile $outdir/poplation/${outname}.HC.snps.indels.plots.R \
--tranches-file $outdir/poplation/${outname}.HC.snps.indels.tranches \
-O $outdir/poplation/${outname}.HC.snps.indels.recal && \
time $gatk ApplyVQSR \
-R $reference/Homo_sapiens_assembly38.fasta \
-input $outdir/poplation/${outname}.HC.snps.VQSR.vcf.gz \
--ts_filter_level 99.0 \
--tranches-file $outdir/poplation/${outname}.HC.snps.indels.tranches \
-recalFile $outdir/poplation/${outname}.HC.snps.indels.recal \
-mode INDEL \
-O $outdir/poplation/${outname}.HC.VQSR.vcf.gz && echo "** SNPs and Indels VQSR (${outname}.HC.VQSR.vcf.gz finish) done **"
在WGS系列文章中大家也可以到类似的程序操作命令。但是大多数初学者可能并不完全理解GATK VQSR中训练集参数(-resource)的内在含义, 。前段时间在知识星球里有几位小伙伴反复问到了这个问题,我一一作了回答,最后进行了总结。虽然这种参数是GATK VQSR模块所特有,但如果能理解好这些参数以及它们背后包含的意义,应该有助于我们更好地理解基因组的变异质控,更恰当地使用GATK。另外,对于生信算法开发者来说,还可以从这样的策略中或多或少地得到一些启发。
所以,这篇文章里,就把我在知识星球中有关这一块的内容整理出来分享给大家。
SNP和Indel的特征是不同的,评价它们的好坏需要使用不同的标准 它们各自的评估模型就需要分开训练和计算 不需要自己去把VCF的SNP和Indel分出来
SNP
对于SNP来说,VQSR的训练集数据目前主要有四个:HapMap,Omni,1000G和dbSNP,参数一般如下:
-resource:hapmap,known=false,training=true,truth=true,prior=15.0 hapmap_3.3.hg38.vcf \
-resource:omni,known=false,training=true,truth=false,prior=12.0 1000G_omni2.5.hg38.vcf \
-resource:1000G,known=false,training=true,truth=false,prior=10.0 1000G_phase1.snps.high_confidence.hg38.vcf \
-resource:dbsnp,known=true,training=false,truth=false,prior=2.0 dbsnp_146.hg38.vcf \
-an DP -an QD -an FS -an SOR -an ReadPosRankSum -an MQRankSum \
-mode SNP \
这个参数的格式,从左到右,按照逗号(,)隔开,分别是:
训练集名字,这个名字是可以随便改动的,但是为了便于交流,一般还是默认按照数据集的名字来设置( );
known:该数据是否作为已知变异数据集,用于对变异数据的标注;
training:该数据是否作为模型训练的数据集,用于训练VQSR模型;
truth:该数据是否作为验证模型训练的真集数据,这个数据同时还是VQSR训练bad model时自动进行参数选择的重要数据;
prior:该数据集在VQSR模型训练中的权重,或者叫Prior likelihood( )。
VQSR训练集的数据可以用-resource参数继续往下添加,有多少就加多少,对于人类WGS/WES数据来说,目前用的主要还是上面的这4个,
但这些差异的实质和所代表的意义什么呢?在使用的时候到底是基于什么原则设定了这些参数?
毫无疑问的是这些差异的实质一定源自于不同的数据集本身 只要能够弄清楚这几点,那么参数的设置也会不言自明
第一个是HapMap, 由于这个数据集包含了大量家系数据,并且有非常严格的质控和严密的实验验证,因此它的准确性是目前公认最高的
所以VQSR进行质控模型训练的时候,会将其作为一个很重要的训练集(training=true)。它的权重也会被设置得很高,比如在WGS数据分析中常常设为prior=15——这里的Prior是Prior likelihood的Phred-scale,我们如果把15转换为likelihood,那么就是0.96838。此外,由于它的高准确性,通常还将作为模型验证的一个真集数据(truth=true)。
第二个是Omni 它也是一个高可信的变异结果 通常情况下也可以把它作为验证结果的真集数据 你没太大的“洁癖”把它设置为true都是没问题的(这也是GATK最佳实践的一般做法)。
第三个是1000G
第四个是dbSNP。 这是一个绝对不可以作为训练集位点的数据 dbSNP收集的数据,实际都是研究者们发表了相关文章提交上来的变异,这些变异很多是没做过严格验证的,很多甚至还是假的,在没被反复验证之前,是不可信的
Indel
对于Indel来说,VQSR模型的训练参数只有两个:Mill和dbSNP。
-resource:mills,known=true,training=true,truth=true,prior=12.0 Mills_and_1000G_gold_standard.indels.hg38.vcf \
-resource:dbsnp,known=true,training=false,truth=false,prior=2.0 dbsnp_146.hg38.vcf \
-an DP -an QD -an FS -an SOR -an ReadPosRankSum -an MQRankSum \
-mode INDEL \
第一个是Mills
第二个还是dbSNP,WGS系列的这篇文章里,给出的流程就是这种情况。
一些需要额外注意的地方
GATK VQSR在执行的时候要基于全基因组的所有变异数据,而不能拆分不同的染色体分别取执行,否则会导致各个染色体所训练的质控模型不一致。另外,除了要区分SNP和Indel模式之外,GATK VQSR分为两个步骤进行:VariantRecalibrator 和 ApplyVQSR,两者缺一不可。VariantRecalibrator用来进行模型计算,获得数据的情况,ApplyVQSR则是根据我们设定的ts_filter_level参数,最终过滤得到数据,这个参数基于我们对模型真集数据的灵敏度和特异性来确定,一般会设置为99.0%(比如上面例子)。但事实上,这个比例,并不是说它是最合适的,设定为99或者99.5本身并没有什么理论证明,更多还是一种约定俗成,或者大部分情况下,看到设置为这个值,结果看起来是合理的,而且是在正常人样本数据得到的认识。所以,有时候还是需要具体问题具体分析,多看VQSR得到的tranche图,特别是对于非人物种,如果你觉得95.0%也很合适,那么也不一定非得是99%。
----/ END /----
※ ※ ※
你还可以读
这是知识星球:『达尔文星球』(原名:解螺旋技术交流圈),是一个我与读者朋友们的私人朋友圈。我有9年前沿而完整的生物信息学、NGS领域的工作经历,在该领域发有多篇Nature级别的科学文章,我也希望借助这个知识星球把自己的一些微薄经验分享给更多对组学感兴趣的伙伴们。
这是知识星球上第一个真正与基因组学和生物信息学强相关的圈子。我希望能够借此营造一个高质量的组学知识圈和人脉圈,通过提问、彼此分享、交流经验、心得等,彼此更好地学习生信知识,提升基因组数据分析和解读的能力。
在这里你可以结识到全国优秀的基因组学和生物信息学专家,同时可以分享你的经验、见解和思考,有问题也可以向我提问和圈里的星友们提问。