GATK4插件Mutect2参数解析

作者按

最近在准备着换一个职业赛道,所以在做之前所有项目的回溯,遇到了最最基础的SNV+Indel的流程,给别人重新讲了一遍Mutect2的过滤规则和参数选择,发现这个,含金量比我之前写的SV和CNV高多了。贴出来给我考试攒人品啦。

Mutect2是基于somatic极大似然模型,通过寻找active region,迭代寻找可能突变的一款认可度很高的软件。

值得注意的是目前的Mutect2单独call变异属于BETA测试版,只有tumor-normal成对call是经过学术认可的。

但是目前各大医院、临检单位,为了压缩成本、减轻患者压力,普遍采取大panel像全外显子这种的采取成对call,小panel单独only tumor。通过SNP数据集,经验值或者其他一代二代验证方法等等过滤产生的假阳性。

那么过滤就需要细之又细,以下是Mutect2的所有参数。参考来源于Mutect官方文档之Mathematical notes on Mutect的chapter8。

概述

Mutect2一共有14个过滤标签(vcf的filter列可能出现的tag),每个标签对应一个或者好几个值。vcf里每个点都有这14个值,值的意义在vcf的info列列出,见下表的key。这每个key或者filter对应一个在运行Mutect时候设置的参数,见下表的Argument。

Mutect参数对应的过滤标签和INFO列的对应键值对

举个例子:

见下图的某位点,base_quality标签出现在了filter列,它在call的时候估算的参数为MBQ=0,而我们设置的min-median-base-quality参数为20,因为0<20,所以base_quality的标签出现在了tag里边。其余参数以此类推。

pos example

像这种与likelihood model相关性弱一点的参数很好理解。万一是t-lod标签呢?它代表什么意思呢?又是如何算出来的呢?这就涉及到它的数学模型了,我已经努力地在写一篇比较浅显易懂的介绍了,但是因为OneNote写公式实在太累,所以我决定手推,等下次更新的时候就可以看到一篇全篇公式的手写图片式科普23333

我一定是准备GMAT的数学太简单了才如此找虐的。。。

好了,言归正传,接下来依次介绍各个标签的含义:

参数含义

1. t_lod

  • tumor-lod is the minimum likelihood of an allele as determined by the somatic likelihoods model required to pass.
  • LOD threshold for calling tumor variant.
  • 似然模型中认为该点是体细胞变异的最小似然比,默认是5.3,若小于5.3则添加tlod标签。
  • Default value: 5.3
  • Example


    t-lod example

2.clustered_events

  • max-events-in-region is the maximum allowable number of called variants co-occurring in a single assembly region. If the number of called variants exceeds this they will all be filtered.

  • Variants coming from an assembly region with more than this many events are filtered.

  • 活动区域发生多次突变,且突变位点距离在3bp及以上。【为什么大于3bp呢?因为2个的话很有可能是可以合并的同一个。】

  • Default value: 2.

  • 因为这个曾经困扰了我很久,所以直接在GitHub扒了gatk源码。查看的时候关于这个参数的部分长这样


  • 然后又询问官方,得到的回复是建议认为是假性。


  • example



  • note
    如果是大panel,大部分位点都会带这个标签,根据经验,如果过滤,会过滤掉真阳性位点,我个人建议保留,或者至少验证一下。【真的不能粗暴过滤,我验证过的】

    以下是我当时调整参数--max-enents-in-region长度,得到clustered_events标签的个数。发现十几bp的时候,针对我的panel,带标签的点并不会减少多少。所以用官方的

3.duplicated_evidence

  • unique-alt-read-count is the minimum number of unique (start position, fragment length) pairs required to make a call. This count is a proxy for the number of unique molecules (as opposed to PCR duplicates) supporting an allele.
  • 设置单独支持该等位基因的最小分子数,默认值是0
  • Filter a variant if a site contains fewer than this many unique (i.e. deduplicated) reads supporting the alternate allele Default value: 0.
  • Example


4.multiallelic

  • max-alt-allele-count is the maximum allowable number of alt alleles at a site. By default only biallelic variants pass the filter.
  • 某一个点等位基因的最大数,默认情况下,只有双等位基因变异通过过滤器。
  • Default value: 1.
  • filter variants with too many alt alleles
  • Example


5.germline_risk

  • max-germline-posterior is the maximum posterior probability, as determined by the above germline probability model, that a variant is a germline event.
  • 该位点是germline event的最大后验概率, 根据模型计算P_GERMLINE值,大于设定值就会添加germline_risk的标签。
  • Maximum posterior probability that an allele is a germline variant.
  • Default value: 0.1.
  • Example


6.artifact_in_normal

  • normal-artifact-lod is the maximum acceptable likelihood of an allele in the normal by the somatic likelihoods model. This is different from the normal likelihood that goes into the germline model, which makes a diploid assumption. Here we compute the normal likelihood as if it were a tumor in order to detect artifacts.
  • 当tumor和control成对call的时候,会对control组的normal样本单独设置对数比阈值,该阈值越高,过滤标准越严格,因为认为normal全部是假阳性,所以会设置较低的LOD值。
  • LOD threshold for calling normal artifacts
  • Default value: 0.0.
  • example
    我没有tumo-normal对,所以无图可用,可怜兮兮~~~

7.strand_artifact

  • max-strand-artifact-probability is the posterior probability of a strand artifact, as determined by the model described above, required to apply the strand artifact filter.
  • 链偏好性的后验概率,根据计算的SA_POST_PROB,大于设定值则过滤;还有第二层补充条件;
  • This is necessary but not sufficient – we also require the estimated max a posteriori allele fraction to be less than min-strand-artifact-allele-fraction.The second condition prevents filtering real variants that also have significant strand bias, i.e. a true variant that also has some artifactual reads.
  • 如果链偏好性的最大后验概率比SA_MAP_AF(MAP estimates of allele fraction given
  • 变异频率的最大后验概率)值小,会保留,以防将真阳性位点加上strand_artifact标签。
  • Filter a variant if the probability of strand artifact exceeds this number
  • Default value: 0.99.
  • Example


8.base_quality

  • min-median-base-quality is the minimum median base quality of bases supporting a SNV.
  • 支持SNV的最小碱基质量。( median base quality of bases 碱基质量衡量标准,在fastqc软件中该值为25)
  • filter variants for which alt reads' median base quality is too low.
  • Default value: 20.
  • Example


9.mapping_quality

  • min-median-mapping-quality is the minimum median mapping quality of reads supporting an allele.
  • 支持SNV的最小mapping质量。
  • Filter variants for which alt reads' median mapping quality is too low.
  • Default value:30.
  • Example


10.fragment_length

  • max-median-fragment-length-difference is the maximum difference between the median fragment lengths reads supporting alt and reference alleles. Note that fragment length is based on where paired reads are mapped,not the actual physical fragment length.
  • 是支持alt和ref片段长度之间的最大差异。片段长度是基于mapping成对读取时的长度,而不是实际物理片段长度。过滤掉ref和alt比对片段差异巨大的点。
  • Filter variants for which alt reads' median fragment length is very different from the median for ref reads.
  • Default value: 10000.
  • example
    这个也是tumor-normal成对call才会出现的,我没有例子可展示,可怜兮兮X2

11.read_position

  • min-median-read-position is the minimum median length of bases supporting an allele from the closest end of the read. Indels positions are measured by the end farthest from the end of the read.
  • 位点到read末尾的最近读取端的最小中值长度。DENELS的位置是由读数末尾最远的一端测量的。
  • filter variants for which the median position of alt alleles within reads is too near the end of reads.
  • Default value: 5.
  • Example


12.panel_of_normals

  • One of the two unadjustable filters. the panel of normals filter removes all alleles at a site belonging to the panel of normals, which is a vcf of blacklisted artifact sites. It can be disabled by not passing a panel of normals to Mutect2.
  • 该位点也在PON中存在
  • Example


13.contamination

  • If FilterMutectCalls is passed a contamination-table from CalculateContamination it will filter alleles with allele fraction less than the whole-bam contamination in the table.
  • 过滤样品污染
  • --contamination-table:File
  • example:
    单个样品就不需要做啦

14.str_contraction

  • One of the two unadjustable filters, an STR contraction filter which removes variants that are the deletion of a single repeat unit of an STR when this repeat unit contains more than one base.

  • 当该重复大于一个base,过滤短串联重复区,STR( Variant is a short tandem repeat )。

  • RPA=Number of times tandem repeat unit is repeated, for each allele (including reference)

  • RU=Tandem repeat unit (bases)

  • Example


写在最后

很多参数,比如t-lod是模型最终的检验T值,artifact_in_normal是normal的后验概率相关,如果不了解模型,可能无法理解其假阳性的中间推导过程。建议各位同学仔细读一读Mathematical Notes on Mutect(David Benjamin� and Takuto Sato†Broad Institute, 75 Ames Street, Cambridge, MA 02142
(Dated: September 26, 2018)。

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 216,324评论 6 498
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 92,356评论 3 392
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 162,328评论 0 353
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 58,147评论 1 292
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 67,160评论 6 388
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 51,115评论 1 296
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 40,025评论 3 417
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 38,867评论 0 274
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 45,307评论 1 310
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,528评论 2 332
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 39,688评论 1 348
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 35,409评论 5 343
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 41,001评论 3 325
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,657评论 0 22
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,811评论 1 268
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 47,685评论 2 368
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 44,573评论 2 353