基因突变检测之Mutect2
写在前面的话
今天是2021年1月3号，年底一番瞎忙活，停更了好长一段时间的简书。人生中有个很重要伯乐导师告诉我，在这个事多的工作环境下，没人愿意看你太多的废话。直奔重点才是最重要的。后续将以简洁和目标为主，今日最后一话，后共勉。
对于基因突变（snv or indel）检测，一般的思路是采用校正后的bam文件采用软件call变异。不同的软件其突变检测原理不一致。软件检测结果一般存储为VCF格式。不同的软件其VCF文件的infomation不一致。接下来解读一下行业内用的比较多的三款软件，分别为Mutect2, VarScan, VarDict2。今天主讲为Mutect2。创作不易，一个假期没啦！！！
对于Mutect2，讲述目录如下：
- Mutect2 突变检测原理
- Mutect2 过滤条件
- Mutect2 VCF文件格式解释
1. Mutect2 突变检测原理
1.1 Mutect2 更新情况
Mutect2是GATK4的模块，目前GATK4已经升级到4.1.9.0，不得不说，我的4.1.8.1版本也在几月前尚未焐热。这升级的速度让人有些头疼。但GATK官网推荐使用最新版本，github上展示更新如下：
4.1.9.0升级部分
那两个新Tools不是我们的重点，关于重点如下：
4.1.9.0 debug
其原理是通过对单倍型的局部组装call体细胞的（snv and indel）。采用的是贝叶斯体细胞基因分型模型。流程图如下：
image.png
1.2 Mutect2模块检测突变原理
Like HaplotypeCaller, Mutect2 calls SNVs and indels simultaneously via local de-novo assembly of haplotypes in an active region. That is, when Mutect2 encounters a region showing signs of somatic variation, it discards the existing mapping information and completely reassembles the reads in that region in order to generate candidate variant haplotypes. Like HaplotypeCaller, Mutect2 then aligns each read to each haplotype via the Pair-HMM algorithm to obtain a matrix of likelihoods. Finally, it applies a Bayesian somatic likelihoods model to obtain the log odds for alleles to be somatic variants versus sequencing errors.
与HaplotypeCaller一样，Mutect2通过在active region中对单倍型（决定同一性状的紧密连锁的基因构成的基因型）进行局部重新组装来判断SNV和InDel。 也就是说，当Mutect2遇到显示体细胞变异迹象的active region时，它会丢弃现有的映射信息，并完全重新组装该区域中的reads，以生成候选变异单倍型。 像HaplotypeCaller一样，Mutect2然后通过Pair-HMM算法将每个reads与每个单倍型对齐，以获得似然矩阵。 最后，它应用贝叶斯体细胞似然模型来获得体细胞变异与测序错误的对数比。
为加快分析速读，这里我根据染色体的bed文件拆分染色体，下面脚本以配对样本为例：
/software/GATK4/gatk-4.1.9.0/gatk --java-options "-XX:ParallelGCThreads=4 -Xmx20G -Djava.io.tmpdir=./tmp" Mutect2 \
-L chr1.bed --pcr-indel-model CONSERVATIVE \
-I sample_N.recal.bam -normal sample_N \
-I sample_T.recal.bam -tumor sample_T \
-R ucsc.hg19.fasta --native-pair-hmm-threads 8 \
--germline-resource Gatk_bundles/af-only-gnomad.raw.sites.hg19.vcf.gz \
--panel-of-normals PoN.vcf.gz \
--f1r2-tar-gz chr1.f1r2.tar.gz \
-O sample.chr1.raw.vcf
参数解释：
-L 分别为样本的Tumor 和Normal 的BAM文件，-tumor提供Tumor 样本的BAM文件对应的read group id（BAM头文件中的@RG->SM值），-normal提供正常样本的BAM的read group id。Mutect2可以据此排除Germline变异和个体特异性人工产物；如果tumor不能提供配对normal样本，得到的体细胞变异文件会产生更多假阳性。
-R：参考基因组序列文件。
--native-pair-hmm-threads ： 通过Pair-HMM算法将每个reads与每个单倍型对齐获得似然矩阵所用线程数。
--germline-resource：GATK数据库中hg19表明为Germline的突变位点。指定人群的germline信息来注释变异等位基因。
--panel-of-normals：PoN(正常样本callset)定义了一个预过滤变异位点集。PoN不仅呈现了通常的germline变异位点，也呈现了测序数据引入的噪音，如测序bias或比对bias。默认情况下，出现在PoN中的变异位点不认为是tumor样本的somatic变异也不会进行重组装。如果tumor位点和PoN位点变异不是精确匹配，Mutect2也会重新组装此区域，在变异文件中输出此位点，并在Filter列标记为“panel_of_normals”。
1.3 获得PoN文件
对于多少个样本可以构成PON，并没有确切的说法，但GATK论坛给的建议是最少用40个样本来构建，以normal01 为例。
/software/GATK4/gatk-4.1.9.0/gatk --java-options "-XX:ParallelGCThreads=4 -Xmx20G -Djava.io.tmpdir=./tmp" Mutect2 \
-R ucsc.hg19.fasta
-I normal01.bam \
-tumor normal01 \
-L bed.file \
-O normal01.vcf.gz
整理上一步生成的VCFs，运行命令。
/software/GATK4/gatk-4.1.9.0/gatk --java-options "-XX:ParallelGCThreads=4 -Xmx20G -Djava.io.tmpdir=./tmp" GenomicsDBImport \
-R ucsc.hg19.fasta
-L bed.file
--genomicsdb-workspace-path PoN.db
-V normal01.vcf.gz \
-V normal02.vcf.gz \
-V normal03.vcf.gz \
...
-V normaloN.vcf.gz \
-O PoN.vcf.gz
CreateSomaticPanelOfNormals形成PoN文件。
/software/GATK4/gatk-4.1.9.0/gatk --java-options "-XX:ParallelGCThreads=4 -Xmx20G -Djava.io.tmpdir=./tmp" CreateSomaticPanelOfNormals \
--germline-resource Gatk_bundles/af-only-gnomad.raw.sites.hg19.vcf.gz \
-R ucsc.hg19.fasta \
-V gendb://PoN.db \
-O PoN.vcf.gz
2. Mutect2 过滤条件
2.1 样本间污染评估
This step emits an estimate of the fraction of reads due to cross-sample contamination for each tumor sample and an estimate of the allelic copy number segmentation of each tumor sample. Unlike other contamination tools, CalculateContamination is designed to work well without a matched normal even in samples with significant copy number variation and makes no assumptions about the number of contaminating samples.
该步骤针对每个肿瘤样品发出由于样本间交叉污染引起的reads比例的估计，以及每个肿瘤样品的等位基因拷贝数分段的估计。 与其他污染工具不同，即使在拷贝数变化显着的样品中，CalculateContamination的设计也可以在没有匹配正常值的情况下很好地工作，并且不对污染样品的数量做出任何假设。在FilterMutectCalls后，会对疑似污染的突变条目给一个contamination标签。
关于样本间污染评估，分为3个步骤：
- 对肿瘤BAM运行GetPileupSummaries以总结tumor样本在已知变异位点集上的reads支持情况。
- 如果存在配对样本，会对正常样本运行GetPileupSummaries以总结tumor样本在已知变异位点集上的reads支持情况。
- 对已知变异位点集采用CalculateContamination来估计污染比例，segments.table文件的最后一列用来判断突变位点是否为样本污的位点。
需要指出几点：在默认参数中，此工具只考虑样本中纯和备用位点：等位基因频率范围在0.01-0.2（相关参数是--minimum-population-allele-frequency和--maximum-population-allele-frequency）,这样设计的理论基础是：如果某个纯和备用位点的人群频率较低，当发生样本交叉污染时，我们更容易观测到ref allele（或更常见allele）的出现，这样我们可以更准确地定量污染比例。
GetPileupSummaries默认参数
2.2 合并每一条Mutect中raw vcf
01 采用MergeVcfs合并所有染色体的raw.vcf
software/GATK4/gatk-4.1.9.0/gatk --java-options "-XX:ParallelGCThreads=4 -Xmx10G -Djava.io.tmpdir=./tmp" MergeVcfs \
-I sub_mutect2/sample.chr1.raw.vcf \
-I sub_mutect2/sample.chr2.raw.vcf \
-I sub_mutect2/sample.chr3.raw.vcf \
...
-I sub_mutect2/sample.chrX.raw.vcf \
-I sub_mutect2/sample.chrY.raw.vcf \
-O sample.MuTect2.raw.vcf
02 采用MergeMutectStats合并所有染色体的raw.vcf.stats。
/software/GATK4/gatk-4.1.9.0/gatk --java-options "-XX:ParallelGCThreads=4 -Xmx10G -Djava.io.tmpdir=./tmp" MergeMutectStats
--stats sub_mutect2/sample.chr1.raw.vcf.stats \
--stats sub_mutect2/sample.chr2.raw.vcf.stats \
--stats sub_mutect2/sample.chr3.raw.vcf.stats \
...
--stats sub_mutect2/sample.chrX.raw.vcf.stats \
--stats sub_mutect2/sample.chrY.raw.vcf.stats \
-O sample.MuTect2.raw.vcf.stats
03 测序偏好矫正
这一步主要是为了矫正测序产生的碱基偏好，对于FFPE样本来说这一步很重要。如果不是FFPE样本也可以做，并不会影响后续的结果准确度。这里用到的f1r2.tar.gz来源于第2步的输出。官网对该功能的表述如下：
This tool uses an optional F1R2 counts output of Mutect2 to learn the parameters of a model for orientation bias. It finds prior probabilities of single-stranded substitution errors prior to sequencing for each trinucleotide context. This is extremely important for FFPE tumor samples.
/software/GATK4/gatk-4.1.9.0/gatk --java-options "-XX:ParallelGCThreads=4 -Xmx10G -Djava.io.tmpdir=./tmp" LearnReadOrientationModel
-I sub_mutect2/chr1.f1r2.tar.gz \
-I sub_mutect2/chr2.f1r2.tar.gz \
-I sub_mutect2/chr3.f1r2.tar.gz \
...
-I sub_mutect2/chrX.f1r2.tar.gz \
-I sub_mutect2/chrY.f1r2.tar.gz \
-O sample.read-orientation-model.tar.gz
04 VCF过滤
首先，需要确定过滤什么。这个需要对该模块的参数进行了解，一般情况下有以下几点。
--contamination-table :Tables containing contamination information
--tumor-segmentation: sample.Tumor.segments.table
--min-allele-fraction: <Double>Minimum allele fraction required Default value: 0.0.
--min-median-base-quality:Minimum median base quality of alt reads Default value: 20.
--min-median-mapping-quality:Minimum median mapping quality of alt reads Default value: 30.
--min-median-read-position:Minimum median distance of variants from the end of reads Default value: 1. ###存在reads边缘的突变一般不可靠。
--min-reads-per-strand：Minimum alt reads required on both forward and reverse strands Default value: 0. ##为0的一般为链偏的突变类型了。
--normal-p-value-threshold： P value threshold for normal artifact filter Default value: 0.001.
--orientation-bias-artifact-priors：Sample.read-orientation-model.tar.gz
--unique-alt-read-count:Minimum unique (i.e. deduplicated) reads supporting the alternate allele Default value: 0
min-allele-fraction 和unique-alt-read-count一般要根据需要进行设定。
参考代码如下：
/software/GATK4/gatk-4.1.9.0/gatk --java-options "-XX:ParallelGCThreads=4 -Xmx10G -Djava.io.tmpdir=./tmp" FilterMutectCalls \
-R ucsc.hg19.fasta --min-median-mapping-quality 20 \
-V sample.MuTect2.raw.vcf \
-O sample.MuTect2.raw.filter.vcf \
--contamination-table contamina/sample.contamination.table \
--tumor-segmentation contamina/sample.Tumor.segments.table \
--orientation-bias-artifact-priors sample.read-orientation-model.tar.gz
3 VCF格式解释。
3.1 VCF表头解释
表头一共有10列，描述如下：
1) CHROM：表示染色体。
2) POS：表示突变位点的位置。
3) ID：表示突变的ID，比如在dbSNP中有该SNP的id，则会在此行给出；若没有，则用’.'表示其为一个novel variant。
4) REF：表示参考基因组上该位点的碱基。
5) ALT：与参考基因组序列比较发生突变的碱基。
6) QUAL：表示 Phred格式(Phred_scaled)的质量值，表示在该位点存在variant的可能性；该值越高，则variant的可能性越大；计算方法：Phred值 = -10 * log (1-p) p为variant存在的概率; 通过计算公式可以看出值为10的表示错误概率为0.1，该位点为variant概率为90%。
7) FILTER：表示过滤状态（filter statu），因为使用上一个QUAL值来进行过滤的话，是不够的。若variant不可靠，则该项不为”PASS”或”.”。
Mutect的VCF中未提供终止位置，需要对存在多等位基因的情况进行分割后，根据参考基因的位置来获得起始和终止位置。这里不做说明。
配对样本和Tumor only样本的INFO是不一样的。如下所示：
Tumor only的样本
有配对样本
3.2 VCF FILTER项解释
Mutect2一共有14个过滤标签（vcf的FILTER列可能出现的tag），每个标签对应一个或者好几个值。vcf里每个点都有这14个值，值的意义在vcf的INFO列中列出，见下表的key。这每个key或者FILTER对应一个在运行Mutect时候设置的参数。见下表的Argument。该描述参考佳期如梦你也是的简书。个人感觉写的让人茅塞顿开！出现某个标签标示该突变条目被过滤掉的原因。为PASS表示不满足所有的过滤条件。
Filter的解释
具体描述如下：
a) base_qual：alt的碱基质量值的中位数（alt median base quality）
clustered_events：在tumor中观察到的聚集事件（Clustered events observed in the tumor）
b) contamination：污染情况（contamination）
c) duplicate：（evidence for alt allele is overrepresented by apparent duplicates）
d) fragment：片段，表示片段长度的中位数（abs(ref - alt) median fragment length）
e) germline：胚系突变，有证据表明该位点是胚系的（Evidence indicates this site is germline, not somatic）
f) haplotype：单倍体，（Variant near filtered variant on same haplotype.）
g) low_allele_frac：等位基因分数是低于指定阈值的（Allele fraction is below specified threshold）
h) map_qual：比对质量值（ref - alt median mapping quality）
i) multiallelic：多等位基因，（Site filtered because too many alt alleles pass tumor LOD）。LOD：log odds ratio（对数差异比）。
j) n_ratio：N与alt的比值超过指定值（Ratio of N to alt exceeds specified ratio）N应该是normal的意思。
k) normal_artifact：在normal中的artifact（artifact_in_normal）
l) numt_chimera：NuMT variant有很多ALT reads是来自常染色体的（NuMT variant with too many ALT reads originally from autosome）。线粒体基因组插入（nuclear insertions of mitochondrial genome，NUMTs）
m) numt_novel：Alt深度低于常染色体NuMT的预期覆盖范围（Alt depth is below expected coverage of NuMT in autosome）
n) orientation：方向，表示通过方向偏差混合模型检测到的方向偏差（orientation bias detected by the orientation bias mixture model）
o) orientation_bias：方向偏差，表示在一个或多个样本中看到的方向偏差（以指定的伪影模式或补码之一）。（Orientation bias (in one of the specified artifact mode(s) or complement) seen in one or more samples.）
p) panel_of_normals：表示在normals中的黑名单位点（Blacklisted site in panel of normals）
q) position：表示alt variants与read 末端距离的中位数（median distance of alt variants from end of reads）
r) slippage：滑移，表示通过短的串联重复区域的收缩筛选出位点（Site filtered due to contraction of short tandem repeat region）
s) strand_bias：链偏移，表示仅来自于一条read方向的alt等位基因（Evidence for alt allele comes from one read direction only）
t) strict_strand：表示alt等位基因在两个方向均为显示（Evidence for alt allele is not represented in both directions）
u) weak_evidence：表示突变为达到阈值（Mutation does not meet likelihood threshold）
这里介绍几个主要的参数：
1：t_lod，tumor-lod is the minimum likelihood of an allele as determined by the somatic likelihoods model required to pass.似然模型中认为该点是体细胞变异的最小似然比，默认是5.3,若小于5.3则添加tlod标签。
image.png
2：clustered_events，Variants coming from an assembly region with more than this many events are filtered.活动区域发生多次突变，且突变位点距离在3bp及以上。【为什么大于3bp呢？因为2个的话很有可能是可以合并的同一个。】这个解释很赞的。
image.png
在igv中展示如下，在该突变位点周围存在其他类型的突变：
clustered_events igv示意图
3：multiallelic（多等位基因）
max-alt-allele-count is the maximum allowable number of alt alleles at a site. By default only biallelic variants pass the filter.
一般情况下multiallelic是假的突变类型，但别简单粗暴就给过滤掉了。有的时候其中一个突变类型可能频率很低，而另一个很高，这种有可能是真的。
4：strand_artifact：链偏好性的后验概率，根据计算的SA_POST_PROB，大于设定值则过滤；还有第二层补充条件。
3.3 VCF INFO项解释
a) AD：为样本中ref和alt的等位基因深度（Evidence for alt allele comes from one read direction only）
b) AF：在tumor中alt等位基因的频率（Allele fractions of alternate alleles in the tumor）
c) DP：通过筛选后，该位点的深度（Approximate read depth (reads with MQ=255 or with bad mates are filtered)）
d) F1R2：（Count of reads in F1R2 pair orientation supporting each allele）
F是比对到参考基因组正链，R是比对到参考基因组负链，1,2指的是read1和read2
e) F2R1：（Count of reads in F2R1 pair orientation supporting each allele）
f) FT：基因型水平筛选（Genotype-level filter）
g) GQ：基因型的质量值（Genotype Quality）表示在该位点该基因型存在的可能性；该值越高，则Genotype的可能性越大；计算方法：Phred值 = -10 * log (1-p) p为基因型存在的概率。
h) GT：基因型，两个数字中间用’/'分开，这两个数字表示双倍体的sample的基因型。0 表示样品中有ref的allele； 1 表示样品中variant的allele； 2表示有第二个variant的allele。因此： 0/0 表示sample中该位点为纯合的，和ref一致； 0/1 表示sample中该位点为杂合的，有ref和variant两个基因型； 1/1 表示sample中该位点为纯合的，和variant一致。
i) OBAM：该变异是否可以是给定的REF / ALT artifact模式之一（Whether the variant can be one of the given REF/ALT artifact modes）
j) OBAMRC：该变异是否可以是给定的REF / ALT artifact模式补充之一（Whether the variant can be one of the given REF/ALT artifact mode complements.）
k) OBF：方向偏差错误的altread比率（Fraction of alt reads indicating orientation bias error (taking into account artifact mode complement).）
l) OBP：给定的ref /alt artifact或者complement的方向偏差P值（Orientation bias p value for the given REF/ALT artifact or its complement.）
m) OBQ：给定的ref /alt 错误的方向偏差测量值（Measure (across entire bam file) of orientation bias for a given REF/ALT error.）
n) OBQRC：给定的ref /alt 补充的方向偏差测量值Measure (across entire bam file) of orientation bias for the complement of a given REF/ALT error.
o) PGT：单倍体信息（Physical phasing haplotype information, describing how the alternate alleles are phased in relation to one another）
p) PID：唯一ID（Physical phasing ID information, where each unique ID within a given sample (but not across samples) connects records within a phasing group）
q) PL：指定的三种基因型的可能性。这三种指定的基因型为(0/0,0/1,1/1)，这三种基因型的概率总和为1。（Normalized, Phred-scaled likelihoods for genotypes as defined in the VCF specification）
r) PS：（Phasing set (typically the position of the first variant in the set)）
s) SB：使用Fisher's Exact Test来检测链偏差情况（Per-sample component statistics which comprise the Fisher's Exact Test to detect strand bias.）
参考：
