寫在前面的話
今天是2021年1月3號帘营，年底一番瞎忙活芬迄，停更了好長一段時間的簡書禀梳。人生中有個很重要伯樂導(dǎo)師告訴我算途，在這個事多的工作環(huán)境下扫外，沒人愿意看你太多的廢話筛谚。直奔重點才是最重要的刻获。后續(xù)將以簡潔和目標為主蝎毡，今日最后一話别垮，后共勉碳想。
對于基因突變（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體細胞的（snv and indel）移怯。采用的是貝葉斯體細胞基因分型模型脐帝。流程圖如下：
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)成的基因型）進行局部重新組裝來判斷SNV和InDel笆檀。 也就是說枷遂，當(dāng)Mutect2遇到顯示體細胞變異跡象的active region時痪伦，它會丟棄現(xiàn)有的映射信息笆制，并完全重新組裝該區(qū)域中的reads匆篓，以生成候選變異單倍型窗市。 像HaplotypeCaller一樣，Mutect2然后通過Pair-HMM算法將每個reads與每個單倍型對齊，以獲得似然矩陣。 最后交惯，它應(yīng)用貝葉斯體細胞似然模型來獲得體細胞變異與測序錯誤的對數(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樣本，得到的體細胞變異文件會產(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。默認情況下兴枯，出現(xiàn)在PoN中的變異位點不認為是tumor樣本的somatic變異也不會進行重組裝悠夯。如果tumor位點和PoN位點變異不是精確匹配夕膀，Mutect2也會重新組裝此區(qū)域，在變異文件中輸出此位點具壮，并在Filter列標記為“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標簽。
關(guān)于樣本間污染評估盔憨，分為3個步驟：
- 對腫瘤BAM運行GetPileupSummaries以總結(jié)tumor樣本在已知變異位點集上的reads支持情況徙菠。
- 如果存在配對樣本，會對正常樣本運行GetPileupSummaries以總結(jié)tumor樣本在已知變異位點集上的reads支持情況郁岩。
- 對已知變異位點集采用CalculateContamination來估計污染比例婿奔，segments.table文件的最后一列用來判斷突變位點是否為樣本污的位點。
需要指出幾點：在默認參數(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)如叼，這樣我們可以更準確地定量污染比例冰木。
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é)果準確度社证。這里用到的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ù)進行了解腺律，一般情況下有以下幾點奕短。
--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ù)需要進行設(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值來進行過濾的話，是不夠的周荐。若variant不可靠狭莱，則該項不為”PASS”或”.”。
Mutect的VCF中未提供終止位置概作，需要對存在多等位基因的情況進行分割后腋妙，根據(jù)參考基因的位置來獲得起始和終止位置。這里不做說明讯榕。
配對樣本和Tumor only樣本的INFO是不一樣的骤素。如下所示：
3.2 VCF FILTER項解釋
Mutect2一共有14個過濾標簽（vcf的FILTER列可能出現(xiàn)的tag），每個標簽對應(yīng)一個或者好幾個值瘩扼。vcf里每個點都有這14個值谆甜，值的意義在vcf的INFO列中列出，見下表的key集绰。這每個key或者FILTER對應(yīng)一個在運行Mutect時候設(shè)置的參數(shù)规辱。見下表的Argument。該描述參考佳期如夢你也是的簡書栽燕。個人感覺寫的讓人茅塞頓開罕袋！出現(xiàn)某個標簽標示該突變條目被過濾掉的原因改淑。為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：等位基因分數(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：方向偏差躏升，表示在一個或多個樣本中看到的方向偏差（以指定的偽影模式或補碼之一）。（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)檫_到閾值（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.似然模型中認為該點是體細胞變異的最小似然比钻弄，默認是5.3,若小于5.3則添加tlod標簽佃却。
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è)定值則過濾镣奋；還有第二層補充條件。
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：基因型的質(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模式補充之一（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.）
參考：
