一、基础背景
根据研究对象和目的的不同，选择后续进行胚系突变还是体细突变的分析。
定义：
胚系突变是发生在生殖细胞（精子和卵子）中，由父亲的精子或母亲的卵子带来的，是一出生就携带着的突变信息，通常全身细胞都携带，可能会遗传给后代。
体系突变又叫体细胞突变，是指个体细胞中本来不携带基因突变，但在后天生长发育中由于各种内外因素的影响造成变异，是癌组织样本中特有的突变，通常身体内只有部分细胞带有，且不会遗传给后代。
如何区分胚系突变和体系突变：
可以通过检测结果、突变频率、以往数据库的报道情况、加做胚系验证进而确定。一般胚系突变的突变丰度理论值是50%（杂合）或者100%（纯合），但在实际分析中，会把突变丰度50%左右（如40%-60%）认为是胚系的杂合突变，突变丰度90%-100%认为是胚系的纯合突变；其他比例的突变是体细胞突变的可能性大。
但判断体系、胚系突变最好的方法，还是用癌旁组织或血液作为配对样本进行胚系检测，才能进一步准确区分。例如一个基因突变丰度为39.5%，从突变丰度去判断不太准确，需要进行癌旁组织或血液胚系检测来进行确认：如果胚系中未检测到该突变，那么就是体系突变；如果胚系中检测到该突变，那么就是胚系突变。
胚系、体系检测样本选择：
胚系突变检测一般采集血液、唾液或者口腔拭子等样本，因为这些样本中包含了生殖细胞的遗传信息。
体系突变检测一般采用肿瘤组织样本进行检测，手术切除的肿瘤组织样本最佳，包括新鲜组织、石蜡包埋组织、石蜡切片等，属于大标本；穿刺活检或气管镜等取到的肿瘤组织次之，属于小标本。
恶性渗出液样本，如脑脊液、胸水（胸腔积液）、腹水（腹腔积液）等也含有癌细胞，当患者无法取得组织时也可使用此类样本进行体系检测。
血液样本中因含有肿瘤细胞释放出的DNA，且方便快捷可反复取样，对于很多晚期或年纪较大无法取得组织或其他类型样本的患者，亦是优选。
总结：
对于疾病研究，更为关注体细胞突变，下面我们的代码以寻找体细胞突变为例，来展开后续分析。
二、代码实战
Step7: Mutect2突变检测
由于该样本只有疾病组没有正常对照，所以我们需要借助一些数据库来帮助我们过滤掉假阳性的突变，从不同维度（人群频率、功能预测、临床意义、癌症关联等）全面解读一个变异，从而为后续的筛选和判断提供综合性的证据。。
· 人群频率数据库（gnomAD）：帮助Mutect2将常见的胚系多态性与真正的体细胞突变区分开。
· --panel-of-normalsPON：一个由多个正常样本构建的VCF文件，包含了在特定测序平台、实验流程中反复出现的技术性假阳性位点（比如某些难以比对的位置的系统错误）。PON文件可以从GATK官网下载Broad研究所提供的hg38 PON，或者如果有多个正常样本，可以自己创建。
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 工具基于一系列统计指标（如链偏好性、读段方位偏差、侧翼序列复杂度、肿瘤深度和等位基因分数等）对 Mutect2 的原始调用进行过滤。
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 注释
仅仅知道基因组位置和碱基变化是不够的，我们还需要知道这些突变有什么生物学意义，可以借助 ANNOVAR 或 snpEff等工具为每个变异添加详细的信息（基因信息，功能影响，氨基酸突变，人群频率等）
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 注释结果转 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
# 创建正确的header（保持相同的列数）
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: 结果可视化
这部分在R环境中进行，对整体突变情况，以及关注基因的具体突变信息进行可视化。
library(tidyverse)
rm(list = ls())
require(maftools)
options(stringsAsFactors = F)
## 1.读入 annovar_merge.vcf 转成 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, # 设置文本的大小
addStat = 'median', # 添加统计信息，这里选择显示中位数
dashboard = TRUE, # 显示仪表盘式图表
titvRaw = FALSE) # 不显示硬件替代碱基转换比率（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：这个小图展示了不同变异分类的数量分布。纵轴显示有 7 种变异分类。
Variant Type：这个小图展示了不同变异类型的数量分布。变异类型可以是单核苷酸变异（SNV）、插入、缺失等。
SNV Class：这个小图展示了单核苷酸变异（SNV）的亚类别数量分布。
Variants per sample (Median: 404)：这个小图展示了 2 个 tumor 样本中的变异数目分布。标题中的 "Median: 404" 表示中位数变异数目为 404，这个值可以用来描述样本间的变异负荷。
Variant Classification summary：这个小图是变异分类的汇总图。它将不同变异分类（例如，高风险、低风险等）的变异数目进行了汇总并展示。
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') # 指定标签显示的位置，这里设置为 "all"，表示显示所有标签。
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，我们只需要把这一列补充完整即可继续绘制：
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% 的体细胞突变率，并提供了该基因的 RefSeq 参考转录本编号 NM_004444。横轴表示基因的各个氨基酸位置，纵轴表示基因在不同样本中的突变信息。每个突变都用一个小的“棒棒糖”（lollipop）符号表示，其中糖棒的一端位于相应的氨基酸位置，另一端指向具体的突变类型，具体信息可以再详细了解。
参考帖子：
https://www.bohrium.com/notebooks/2112719645
https://mp.weixin.qq.com/s/k1GafYqt7tXzWysZqlCP_w
欢迎大家评论交流！
（每帖分享：人们总是会选择自己相信的！）
