一、基礎(chǔ)背景
根據(jù)研究對象和目的的不同，選擇后續(xù)進(jìn)行胚系突變還是體細(xì)突變的分析。
定義：
胚系突變是發(fā)生在生殖細(xì)胞（精子和卵子）中，由父親的精子或母親的卵子帶來的，是一出生就攜帶著的突變信息，通常全身細(xì)胞都攜帶，可能會遺傳給后代。
體系突變又叫體細(xì)胞突變，是指個體細(xì)胞中本來不攜帶基因突變，但在后天生長發(fā)育中由于各種內(nèi)外因素的影響造成變異，是癌組織樣本中特有的突變，通常身體內(nèi)只有部分細(xì)胞帶有，且不會遺傳給后代。
如何區(qū)分胚系突變和體系突變：
可以通過檢測結(jié)果、突變頻率、以往數(shù)據(jù)庫的報道情況、加做胚系驗證進(jìn)而確定。一般胚系突變的突變豐度理論值是50%（雜合）或者100%（純合），但在實際分析中，會把突變豐度50%左右（如40%-60%）認(rèn)為是胚系的雜合突變，突變豐度90%-100%認(rèn)為是胚系的純合突變；其他比例的突變是體細(xì)胞突變的可能性大。
但判斷體系、胚系突變最好的方法，還是用癌旁組織或血液作為配對樣本進(jìn)行胚系檢測，才能進(jìn)一步準(zhǔn)確區(qū)分。例如一個基因突變豐度為39.5%，從突變豐度去判斷不太準(zhǔn)確，需要進(jìn)行癌旁組織或血液胚系檢測來進(jìn)行確認(rèn)：如果胚系中未檢測到該突變，那么就是體系突變；如果胚系中檢測到該突變，那么就是胚系突變。
胚系、體系檢測樣本選擇：
胚系突變檢測一般采集血液、唾液或者口腔拭子等樣本，因為這些樣本中包含了生殖細(xì)胞的遺傳信息。
體系突變檢測一般采用腫瘤組織樣本進(jìn)行檢測，手術(shù)切除的腫瘤組織樣本最佳，包括新鮮組織、石蠟包埋組織、石蠟切片等，屬于大標(biāo)本；穿刺活檢或氣管鏡等取到的腫瘤組織次之，屬于小標(biāo)本。
惡性滲出液樣本，如腦脊液、胸水（胸腔積液）、腹水（腹腔積液）等也含有癌細(xì)胞，當(dāng)患者無法取得組織時也可使用此類樣本進(jìn)行體系檢測。
血液樣本中因含有腫瘤細(xì)胞釋放出的DNA，且方便快捷可反復(fù)取樣，對于很多晚期或年紀(jì)較大無法取得組織或其他類型樣本的患者，亦是優(yōu)選。
總結(jié)：
對于疾病研究，更為關(guān)注體細(xì)胞突變，下面我們的代碼以尋找體細(xì)胞突變?yōu)槔?，來展開后續(xù)分析。
二、代碼實戰(zhàn)
Step7: Mutect2突變檢測
由于該樣本只有疾病組沒有正常對照，所以我們需要借助一些數(shù)據(jù)庫來幫助我們過濾掉假陽性的突變，從不同維度（人群頻率、功能預(yù)測、臨床意義、癌癥關(guān)聯(lián)等）全面解讀一個變異，從而為后續(xù)的篩選和判斷提供綜合性的證據(jù)。。
· 人群頻率數(shù)據(jù)庫（gnomAD）：幫助Mutect2將常見的胚系多態(tài)性與真正的體細(xì)胞突變區(qū)分開。
· --panel-of-normalsPON：一個由多個正常樣本構(gòu)建的VCF文件，包含了在特定測序平臺、實驗流程中反復(fù)出現(xiàn)的技術(shù)性假陽性位點（比如某些難以比對的位置的系統(tǒng)錯誤）。PON文件可以從GATK官網(wǎng)下載Broad研究所提供的hg38 PON，或者如果有多個正常樣本，可以自己創(chuàng)建。
1.突變檢測
gatk Mutect2 \
-R ./0.ref_38/hg38.fasta \
-I ./6.gatk/EXON_LJJ_bqsr.bam -tumor "EXON_LJJ" \
-L ./0.ref_38/exon_data.bed \
--germline-resource ./tools/annovar/humandb/af-only-gnomad.hg38.vcf.gz \
--panel-of-normals /home/chencx/LM/project2305/exon/0.ref_38/somatic-hg38_1000g_pon.hg38.vcf.gz \
--af-of-alleles-not-in-resource 0.000001 \
-O ./8.Mutect2/EXON_LJJ_mutect2.vcf
2.過濾假陽性
FilterMutectCalls 工具基于一系列統(tǒng)計指標(biāo)（如鏈偏好性、讀段方位偏差、側(cè)翼序列復(fù)雜度、腫瘤深度和等位基因分?jǐn)?shù)等）對 Mutect2 的原始調(diào)用進(jìn)行過濾。
gatk FilterMutectCalls \
-R ./0.ref_38/hg38.fasta \
-V ./8.Mutect2/EXON_LJJ_mutect2.vcf \
-O ./8.Mutect2/EXON_LJJ_somatic.vcf
3.提取PASS突變
只保留可信度高的變異，過濾掉非 PASS 的變異
bcftools view -f PASS -o ./8.Mutect2/EXON_LJJ_filter.vcf ./8.Mutect2/EXON_LJJ_somatic.vcf
Step7: Annovar 注釋
僅僅知道基因組位置和堿基變化是不夠的，我們還需要知道這些突變有什么生物學(xué)意義，可以借助 ANNOVAR 或 snpEff等工具為每個變異添加詳細(xì)的信息（基因信息，功能影響，氨基酸突變，人群頻率等）
1.注釋
tools/annovar/table_annovar.pl ./8.Mutect2/EXON_LJJ_filter.vcf ./tools/annovar/humandb \
-buildver hg38 \
-out ./8.Mutect2/annotation/EXON_LJJ \
-remove \
-protocol refGene \
-operation g \
-nastring . \
-vcfinput
# 使用awk過濾gnomAD頻率 < 0.001的變異
awk -F'\t' '(NR==1 || $12 < 0.001 || $12 == ".")' \
./8.Mutect2/annotation/EXON_LJJ.hg38_multianno.txt > ./8.Mutect2/annotation/EXON_LJJ_fil0001.txt
2. vcf 注釋結(jié)果轉(zhuǎn) maf
grep -v '^Chr' ./8.Mutect2/annotation/EXON_FYY_fil0001.txt | cut -f 1-24 | awk -v T=EXON_FYY -v N=EXON_FYY_germline '{print $0"\t"T"\t"N}' > ./8.Mutect2/maf/EXON_FYY.annovar.vcf
# 創(chuàng)建正確的header（保持相同的列數(shù)）
head -1 ./8.Mutect2/annotation/EXON_LJJ_fil0001.txt | cut -f 1-24 | awk '{print $0 "\tTumor_Sample_Barcode\tMatched_Norm_Sample_Barcode"}' > ./8.Mutect2/maf/header.txt
3.合并所有vcf文件
cat ./8.Mutect2/maf/header.txt ./8.Mutect2/maf/*.annovar.vcf > ./8.Mutect2/maf/annovar_merge.vcf
ls -lh 8.Mutect2/maf
Step8: 結(jié)果可視化
這部分在R環(huán)境中進(jìn)行，對整體突變情況，以及關(guān)注基因的具體突變信息進(jìn)行可視化。
library(tidyverse)
rm(list = ls())
require(maftools)
options(stringsAsFactors = F)
## 1.讀入 annovar_merge.vcf 轉(zhuǎn)成 maf 格式變量 annovar.laml
annovar.laml <- annovarToMaf(annovar = "./8.Mutect2/maf/annovar_merge.vcf",
refBuild='hg38',
tsbCol = 'Tumor_Sample_Barcode',
table = 'refGene',
MAFobj = T)
write_csv(annovar.laml@data, file.path(outdir, 'annovar_1000.csv'))
table(annovar.laml@data$non_topmed_AF_popmax)
# EXON_FYY_germline EXON_JZY_germline EXON_LJJ_germline
# 2003 149 1762
## 2.繪制突變頻率匯總圖
pdf(file.path(outdir2, 'plotmafSummary_annovar.pdf'),
width = 5, height = 5)
plotmafSummary(maf = annovar.laml,
rmOutlier = TRUE, # 移除異常值
showBarcodes = T, # 顯示條形碼來表示每個樣本的突變頻率
textSize = 0.4, # 設(shè)置文本的大小
addStat = 'median', # 添加統(tǒng)計信息，這里選擇顯示中位數(shù)
dashboard = TRUE, # 顯示儀表盤式圖表
titvRaw = FALSE) # 不顯示硬件替代堿基轉(zhuǎn)換比率（Transition/Transversion）
dev.off()
## 3. 生成癌癥樣本的分析圖：可視化腫瘤樣本中的基因突變信息
pdf(file.path(outdir2, 'oncoplot_top30_annovar.pdf'),
width = 5, height = 5)
oncoplot(maf = annovar.laml,
top = 30, # 繪制前 30 個變異
fontSize=0.5,
sampleOrder = annovar.laml@clinical.data$Tumor_Sample_Barcode, # 樣本的排序順序
showTumorSampleBarcodes = T) # 顯示腫瘤樣本的名稱
dev.off()
該突變頻率匯總圖包括 6 個小圖，每個小圖解釋如下：
Variant Classification：這個小圖展示了不同變異分類的數(shù)量分布?？v軸顯示有 7 種變異分類。
Variant Type：這個小圖展示了不同變異類型的數(shù)量分布。變異類型可以是單核苷酸變異（SNV）、插入、缺失等。
SNV Class：這個小圖展示了單核苷酸變異（SNV）的亞類別數(shù)量分布。
Variants per sample (Median: 404)：這個小圖展示了 2 個 tumor 樣本中的變異數(shù)目分布。標(biāo)題中的 "Median: 404" 表示中位數(shù)變異數(shù)目為 404，這個值可以用來描述樣本間的變異負(fù)荷。
Variant Classification summary：這個小圖是變異分類的匯總圖。它將不同變異分類（例如，高風(fēng)險、低風(fēng)險等）的變異數(shù)目進(jìn)行了匯總并展示。
Top 10 mutated genes：這個小圖展示了變異頻率最高的前 10 個基因。
特定基因：
每個突變的信息都儲存在annovar.laml@data，可以具體查看。
pdf(file.path(outdir2, 'TNXB_annovar.pdf'),
width = 10, height = 3)
maftools::lollipopPlot(maf = annovar.laml,
gene = 'TNXB', # 指定要繪制的基因，這里是 "SP140"
AACol = 'AAChange.refGene', # 顯示的氨基酸改變信息的列，這里是 "AAChange.refGene"。
labelPos = 'all') # 指定標(biāo)簽顯示的位置，這里設(shè)置為 "all"，表示顯示所有標(biāo)簽。
dev.off()
有些基因在表格中存在，但是可視化的時候報錯：Error in maftools::lollipopPlot(maf = annovar.laml, gene = "EPHB4", AACol = "AAChange.refGene", : EPHB4 does not seem to have any mutations!
報錯原因是該基因所在的列的Variant_Type為NA，我們只需要把這一列補(bǔ)充完整即可繼續(xù)繪制：
annovar.laml@data <- annovar.laml@data %>% mutate(Variant_Type = replace_na(Variant_Type, "SNP"))
annovar.laml@data$Gene.refGene <- str_extract(annovar.laml@data$Gene.refGene, "^[^,]+")
annovar.laml@data$Hugo_Symbol <- str_extract(annovar.laml@data$Hugo_Symbol, "^[^,]+")
pdf(file.path(outdir2, 'EPHB4_annovar.pdf'),
width = 10, height = 3)
maftools::lollipopPlot(maf = annovar.laml,
gene = 'EPHB4',
AACol = 'AAChange.refGene',
labelPos = 'all')
dev.off()
該 lollipop 圖可視化了 EPHB4 基因的突變。基因 EPHB4 在所考慮的樣本中具有 66.67% 的體細(xì)胞突變率，并提供了該基因的 RefSeq 參考轉(zhuǎn)錄本編號 NM_004444。橫軸表示基因的各個氨基酸位置，縱軸表示基因在不同樣本中的突變信息。每個突變都用一個小的“棒棒糖”（lollipop）符號表示，其中糖棒的一端位于相應(yīng)的氨基酸位置，另一端指向具體的突變類型，具體信息可以再詳細(xì)了解。
參考帖子：
https://www.bohrium.com/notebooks/2112719645
https://mp.weixin.qq.com/s/k1GafYqt7tXzWysZqlCP_w
歡迎大家評論交流！
（每帖分享：人們總是會選擇自己相信的！）
