寫在前面的話
今天是2021年1月3號，年底一番瞎忙活，停更了好長一段時間的簡書。人生中有個很重要伯樂導(dǎo)師告訴我，在這個事多的工作環(huán)境下，沒人愿意看你太多的廢話。直奔重點才是最重要的。后續(xù)將以簡潔和目標(biāo)為主，今日最后一話，后共勉。
對于基因突變（snv or indel）檢測，一般的思路是采用校正后的bam文件采用軟件call變異。不同的軟件其突變檢測原理不一致。軟件檢測結(jié)果一般存儲為VCF格式。不同的軟件其VCF文件的infomation不一致。接下來解讀一下行業(yè)內(nèi)用的比較多的三款軟件，分別為Mutect2, VarScan, VarDict2。今天主講為Mutect2。創(chuàng)作不易，一個假期沒啦！??！
對于Mutect2，講述目錄如下：
- Mutect2 突變檢測原理
- Mutect2 過濾條件
- Mutect2 VCF文件格式解釋
1. Mutect2 突變檢測原理
1.1 Mutect2 更新情況
Mutect2是GATK4的模塊，目前GATK4已經(jīng)升級到4.1.9.0，不得不說，我的4.1.8.1版本也在幾月前尚未焐熱。這升級的速度讓人有些頭疼。但GATK官網(wǎng)推薦使用最新版本，github上展示更新如下：
那兩個新Tools不是我們的重點，關(guān)于重點如下：
其原理是通過對單倍型的局部組裝call體細(xì)胞的（snv and indel）。采用的是貝葉斯體細(xì)胞基因分型模型。流程圖如下：
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中對單倍型（決定同一性狀的緊密連鎖的基因構(gòu)成的基因型）進(jìn)行局部重新組裝來判斷SNV和InDel。 也就是說，當(dāng)Mutect2遇到顯示體細(xì)胞變異跡象的active region時，它會丟棄現(xiàn)有的映射信息，并完全重新組裝該區(qū)域中的reads，以生成候選變異單倍型。 像HaplotypeCaller一樣，Mutect2然后通過Pair-HMM算法將每個reads與每個單倍型對齊，以獲得似然矩陣。 最后，它應(yīng)用貝葉斯體細(xì)胞似然模型來獲得體細(xì)胞變異與測序錯誤的對數(shù)比。
為加快分析速讀，這里我根據(jù)染色體的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
參數(shù)解釋：
-L 分別為樣本的Tumor 和Normal 的BAM文件，-tumor提供Tumor 樣本的BAM文件對應(yīng)的read group id（BAM頭文件中的@RG->SM值），-normal提供正常樣本的BAM的read group id。Mutect2可以據(jù)此排除Germline變異和個體特異性人工產(chǎn)物；如果tumor不能提供配對normal樣本，得到的體細(xì)胞變異文件會產(chǎn)生更多假陽性。
-R：參考基因組序列文件。
--native-pair-hmm-threads ： 通過Pair-HMM算法將每個reads與每個單倍型對齊獲得似然矩陣所用線程數(shù)。
--germline-resource：GATK數(shù)據(jù)庫中hg19表明為Germline的突變位點。指定人群的germline信息來注釋變異等位基因。
--panel-of-normals：PoN(正常樣本callset)定義了一個預(yù)過濾變異位點集。PoN不僅呈現(xiàn)了通常的germline變異位點，也呈現(xiàn)了測序數(shù)據(jù)引入的噪音，如測序bias或比對bias。默認(rèn)情況下，出現(xiàn)在PoN中的變異位點不認(rèn)為是tumor樣本的somatic變異也不會進(jìn)行重組裝。如果tumor位點和PoN位點變異不是精確匹配，Mutect2也會重新組裝此區(qū)域，在變異文件中輸出此位點，并在Filter列標(biāo)記為“panel_of_normals”。
1.3 獲得PoN文件
對于多少個樣本可以構(gòu)成PON，并沒有確切的說法，但GATK論壇給的建議是最少用40個樣本來構(gòu)建，以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.
該步驟針對每個腫瘤樣品發(fā)出由于樣本間交叉污染引起的reads比例的估計，以及每個腫瘤樣品的等位基因拷貝數(shù)分段的估計。 與其他污染工具不同，即使在拷貝數(shù)變化顯著的樣品中，CalculateContamination的設(shè)計也可以在沒有匹配正常值的情況下很好地工作，并且不對污染樣品的數(shù)量做出任何假設(shè)。在FilterMutectCalls后，會對疑似污染的突變條目給一個contamination標(biāo)簽。
關(guān)于樣本間污染評估，分為3個步驟：
- 對腫瘤BAM運行GetPileupSummaries以總結(jié)tumor樣本在已知變異位點集上的reads支持情況。
- 如果存在配對樣本，會對正常樣本運行GetPileupSummaries以總結(jié)tumor樣本在已知變異位點集上的reads支持情況。
- 對已知變異位點集采用CalculateContamination來估計污染比例，segments.table文件的最后一列用來判斷突變位點是否為樣本污的位點。
需要指出幾點：在默認(rèn)參數(shù)中，此工具只考慮樣本中純和備用位點：等位基因頻率范圍在0.01-0.2（相關(guān)參數(shù)是--minimum-population-allele-frequency和--maximum-population-allele-frequency）,這樣設(shè)計的理論基礎(chǔ)是：如果某個純和備用位點的人群頻率較低，當(dāng)發(fā)生樣本交叉污染時，我們更容易觀測到ref allele（或更常見allele）的出現(xiàn)，這樣我們可以更準(zhǔn)確地定量污染比例。
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 測序偏好矯正
這一步主要是為了矯正測序產(chǎn)生的堿基偏好，對于FFPE樣本來說這一步很重要。如果不是FFPE樣本也可以做，并不會影響后續(xù)的結(jié)果準(zhǔn)確度。這里用到的f1r2.tar.gz來源于第2步的輸出。官網(wǎng)對該功能的表述如下：
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過濾
首先，需要確定過濾什么。這個需要對該模塊的參數(shù)進(jìn)行了解，一般情況下有以下幾點。
--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一般要根據(jù)需要進(jìn)行設(shè)定。
參考代碼如下：
/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：與參考基因組序列比較發(fā)生突變的堿基。
6) QUAL：表示 Phred格式(Phred_scaled)的質(zhì)量值，表示在該位點存在variant的可能性；該值越高，則variant的可能性越大；計算方法：Phred值 = -10 * log (1-p) p為variant存在的概率; 通過計算公式可以看出值為10的表示錯誤概率為0.1，該位點為variant概率為90%。
7) FILTER：表示過濾狀態(tài)（filter statu），因為使用上一個QUAL值來進(jìn)行過濾的話，是不夠的。若variant不可靠，則該項不為”PASS”或”.”。
Mutect的VCF中未提供終止位置，需要對存在多等位基因的情況進(jìn)行分割后，根據(jù)參考基因的位置來獲得起始和終止位置。這里不做說明。
配對樣本和Tumor only樣本的INFO是不一樣的。如下所示：
3.2 VCF FILTER項解釋
Mutect2一共有14個過濾標(biāo)簽（vcf的FILTER列可能出現(xiàn)的tag），每個標(biāo)簽對應(yīng)一個或者好幾個值。vcf里每個點都有這14個值，值的意義在vcf的INFO列中列出，見下表的key。這每個key或者FILTER對應(yīng)一個在運行Mutect時候設(shè)置的參數(shù)。見下表的Argument。該描述參考佳期如夢你也是的簡書。個人感覺寫的讓人茅塞頓開！出現(xiàn)某個標(biāo)簽標(biāo)示該突變條目被過濾掉的原因。為PASS表示不滿足所有的過濾條件。
具體描述如下：
a) base_qual：alt的堿基質(zhì)量值的中位數(shù)（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：片段，表示片段長度的中位數(shù)（abs(ref - alt) median fragment length）
e) germline：胚系突變，有證據(jù)表明該位點是胚系的（Evidence indicates this site is germline, not somatic）
f) haplotype：單倍體，（Variant near filtered variant on same haplotype.）
g) low_allele_frac：等位基因分?jǐn)?shù)是低于指定閾值的（Allele fraction is below specified threshold）
h) map_qual：比對質(zhì)量值（ref - alt median mapping quality）
i) multiallelic：多等位基因，（Site filtered because too many alt alleles pass tumor LOD）。LOD：log odds ratio（對數(shù)差異比）。
j) n_ratio：N與alt的比值超過指定值（Ratio of N to alt exceeds specified ratio）N應(yīng)該是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的預(yù)期覆蓋范圍（Alt depth is below expected coverage of NuMT in autosome）
n) orientation：方向，表示通過方向偏差混合模型檢測到的方向偏差（orientation bias detected by the orientation bias mixture model）
o) orientation_bias：方向偏差，表示在一個或多個樣本中看到的方向偏差（以指定的偽影模式或補(bǔ)碼之一）。（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 末端距離的中位數(shù)（median distance of alt variants from end of reads）
r) slippage：滑移，表示通過短的串聯(lián)重復(fù)區(qū)域的收縮篩選出位點（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：表示突變?yōu)檫_(dá)到閾值（Mutation does not meet likelihood threshold）
這里介紹幾個主要的參數(shù)：
1：t_lod，tumor-lod is the minimum likelihood of an allele as determined by the somatic likelihoods model required to pass.似然模型中認(rèn)為該點是體細(xì)胞變異的最小似然比，默認(rèn)是5.3,若小于5.3則添加tlod標(biāo)簽。
2：clustered_events，Variants coming from an assembly region with more than this many events are filtered.活動區(qū)域發(fā)生多次突變，且突變位點距離在3bp及以上。【為什么大于3bp呢？因為2個的話很有可能是可以合并的同一個?！窟@個解釋很贊的。
在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：鏈偏好性的后驗概率，根據(jù)計算的SA_POST_PROB，大于設(shè)定值則過濾；還有第二層補(bǔ)充條件。
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是比對到參考基因組負(fù)鏈，1,2指的是read1和read2
e) F2R1：（Count of reads in F2R1 pair orientation supporting each allele）
f) FT：基因型水平篩選（Genotype-level filter）
g) GQ：基因型的質(zhì)量值（Genotype Quality）表示在該位點該基因型存在的可能性；該值越高，則Genotype的可能性越大；計算方法：Phred值 = -10 * log (1-p) p為基因型存在的概率。
h) GT：基因型，兩個數(shù)字中間用’/'分開，這兩個數(shù)字表示雙倍體的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模式補(bǔ)充之一（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 補(bǔ)充的方向偏差測量值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.）
參考：
