---

基因组生信分析教程

DNA 1. Germline Mutation Vs. Somatic Mutation 傻傻分不清楚

DNA 2. SCI 文章中基因组变异分析神器之 maftools

DNA 3. SCI 文章中基因组变异分析神器之 maftools

DNA 4. SCI 文章中基因组的突变信号（maftools）

DNA 5. 基因组变异文件VCF格式详解

DNA 6. 基因组变异之绘制精美瀑布图（ComplexHeatmap）

DNA 7. 基因组拷贝数变异分析及可视化 (GISTIC2.0)

DNA 8. 癌症的突变异质性及寻找新的癌症驱动基因(MutSigCV)

DNA 9. 揭秘肿瘤异质性与TMB, MSI之间的相关性

DNA 10. 识别癌症驱动基因 (OncodriveCLUST)

DNA 11. 识别肿瘤蛋白质三维结构上突变热点(HotSpot3D)

DNA 12. SCI 文章绘图之全基因组关联分析可视化(GWAS)

今天介绍肿瘤突变负荷 （TMB），这个指标在肿瘤治疗前的指导作用非常重要，无论是靶向还是免疫治疗，都是必须进行检测，以帮助患者获得最佳的治疗效果，那么这个值该怎么来计算呢？下面就给大家介绍一下。

简      介

肿瘤突变负荷(tumor mutation burden,TMB)是指肿瘤基因组去除胚系突变后的体细胞突变数量。定义为每百万碱基中被检测出的，体细胞基因编码错误、碱基替换、基因插入或缺失错误的总数。肿瘤突变负荷是一种定量的生物标志物，它反映了肿瘤细胞所携带的突变总数。具有高1MB的肿瘤细胞具有较高的新抗原水平，这被认为可以帮助免疫系统识别种瘤并刺激抗肿瘤T细胞的增殖和抗肿瘤反应。肿瘤细胞的突变会随时间推移而不断蓄积，这种现象在正常细胞中是看不到的。

目前市面上的检测方法是指患者靶向测序编码区每百万碱基(Mb)的体细胞突变数目，包括点突变和插入缺失。体细胞的非转义突变可表现为RNA和蛋白水平的变化，产生的新抗原(或新表位)、蛋白片段/段等被自身免疫系统识别为非自身抗原，激活T细胞，引起免疫反应。

免疫治疗(目前主要是免疫检查点抑制剂)在不加选择人群中的有效率相对较低，需要生物标记物(Biomarker)进行优势人群筛选。探寻最佳的疗效预测标志物是当下的研究热点。肿瘤突变负荷(tumor mutation burden,TMB)是一种生物标志物，可以帮助预测患者对免疫治疗的应答。

理论上，TMB 越高，能够被T细胞识别的新抗原产生也越多，接受免疫检查点抑制剂治疗的疗效可能就越好。即TMB越高，肿瘤免疫治疗的效果越好。

目前，有越来越多的证据支持肿瘤突变负荷(TMB)可能作为潜在的疗效预测标志物。多个瘤种的临床研究数据表明，在PD1未选择人群或PD-L1阳性人群，突变负荷越高，PD-1/PD-L1抑制剂治疗的疗效与TMB成正相关。TMB对于多种肿瘤的免疫治疗均有较好的预测价值。TMB预测疗效作用优于PD-L1表达。

PD-L1表达、TMB、多种驱动突变都是潜在的免疫治疗Biomarker。

T细胞活化通路中多个因素可能是免疫治疗潜在Biomarker。

上述Biomarker中仅有PD-L1表达和Mutational load有较高的临床可行性Presented By Jeffrey Weber at 2017 ASCO Annual Meeting

PD-L1表达检测应用广泛但目前存在较大问题TMB对免疫治疗的预测作用得到大型临床研究证实且具有临床可行性多种肺癌驱动突变对免疫治疗疗效预测与TMB高度相关。

图1.不同TMB阈值NSCLC患者治疗时间和生存率

TMB 的影响因素

TMB 可被很多因素影响，包括不同瘤种，内在遗传/分子因素，肿瘤微环境改变，外界致癌原，不同检测方法和不同界值都可能造成 TMB 缺乏异质性，而使之不能作为准确的反应预测因子（图 1）。

TMB 在不同瘤种，生物特征和患者特征都有不同。最高见于皮肤鳞癌，黑色素瘤，NSCLC，最低见于乳头状甲状腺癌；吸烟患者 TMB 常较高；不同种族和性别中 TMB 也不同。

MMR 通路基因突变、DNA 复制通路缺陷如 POLE/POLD 突变，以及 TP53，APOBEC，LRP1B，NF1，BRCA2 等基因突变和高 TMB 相关，而 EGFR，STK11，JAK1/2 和 B2M 等基因突变可能和更低的 TMB 相关。TMB 和拷贝数变异（CNV）也具有一定的相关性。

TMB 的检测

TMB 的检测受样本质量、检测方法和分析方法等多种因素影响，临床应用前应充分了解 TMB 检测的条件。

样本：肿瘤纯度要求 ≥ 20%，需要有正常对照为检测提供胚系变异信息。

方法：全外显子测序（WES）是 TMB 检测的「金标准」。但是 WES 价格昂贵，检测时间长，需要新鲜标本，因而应用受限。靶向测序 panel 已经成为 WES 的有效替代，为准确性考虑，其覆盖范围应 ≥ 1.0 Mb，测序深度 ≥ 500×。

分析：TMB 的中位值和分布范围在不同癌种中有所不同，因此，在各个癌种中分别确定界值十分重要。应使用相同的筛选策略，选择排序在 20% 以上的病例定义为 TMB-H，而前瞻性的临床疗效才是确定 TMB 界值的最佳标准。不同靶向测序 panel 的 TMB 不能通用。

TMB Cut-off值怎么定?

按照在受检者瘤种中的排列顺序将TMB分为高、中、低。从高到低排序，位于0%-25%区间的TMB为高，26%-75%为中，76%-100%为低。即按照人群分布(3分位法/4分位法)。人群中TMB排序最高的前33%或前25%被划定为高TMB因为各公司/研究单位的TMB高、中、低的评估以内部数据库为准。所以，不同公司/研究单位的结果可能没有可比性。

TMB是否也存在优势人群呢?

研究显示，鳞癌比腺癌具有更高的突变负荷，吸烟者具有更高的突变负荷，可能提示鳞癌和吸烟者是优势人群。

TMB检测存在题注意事项:

第一，TMB有别于PD-L1表达，PD-L1表达低的病人不一定TMB低，也可能存在TMB高的情况，在未来研究中，TMB和PD，L1同时检测是非常必要的。

第二，目前涉及TMB检测的研究均为回顾性研究，亟需前瞻性的临床研究来进一步落实TMB的真正作用。

第三，相比于PD-L1的检测，TMB检测有自身不足的地方，TMB检测需要的组织量相对较多，并且需要用二代测序做检测耗时比较长，耗费也比较高。

肿瘤TMB的计算原理和数学模型

美国FDA目前批准了3个检测肿瘤TMB（ tumor mutation burden）的方法：Omisc Core WES、MSK-IMPACT、FMI。如果按实验流程及操作复杂度来讲，这可能是有史以来最复杂的病理诊断方法。流程一般如下：手术取样、提取DNA、建测序文库、靶向捕获、NGS测序、数据比对、变异检测、变异注释、结果解读。复杂方法下的简单目的：数一数肿瘤特有的变异有几个，多数情况下平均每100万个位点还不超过10个变异。

TMB定义是：每Mb区域里非同义的体细胞变异个数。

图 1. somatic:检测到的非同义突变的体细胞变异个数。L:有效覆盖区域。

有效覆盖区域L容易知道，如果覆盖度足够好，可以用bed文件区间长度代替，那么计算TMB主要的工作就是怎么计算somatic。在肿瘤基因检测中计算somatic处于最核心的位置，靶向治疗、免疫治疗等都围绕着somatic展开。

Somatic定义为个体生长过程中产生的变异。目标就是要找在测序中肿瘤组织测到而正常组织没有测到的变异。因为测序有错误，如果测序足够深，肿瘤及正常组织大多数位点都会有变异。那么通过判断正常组织有无变异将不适合。我们可以模型化为肿瘤组织与正常组织某一变异有相同的频率为种系变异germline，否则为somatic。方法是对肿瘤组织及正常组织同时外显子靶向高通量测序，对每个位点进行比较。通过测序深度、频率等参数比较判断是否somatic。业内有很多计算软件可以用来计算somatic。下表是业内使用最多的开源软件，现在选取业界使用最多的VarScan2 的算法来介绍somatic及其TMB的计算原理。

表1 热门软件采用的数学模型

Software	Mathematical model
Variant caller	Statistical approach
VarScan2	Fisher's exact test
Strelka	Bayesian algorithm
Mutect2	Bayesian algorithm

以下是VarScan2输出VCF格式文件位点的结果

SPV=8.8337E-2是用Fisher’s Exact Test计算的P-value值，用来判断是否somatic。那这值是怎么计算得到的，源文献没有具体说明，本人在这里给出计算过程。Fisher’s Exact Test接受2×2行列表作为计算对象，从上述结果提取数据可以得到下表。

表2. 肿瘤样本与正常样本覆盖深度2×2行列表。RD: 支持参考序列的碱基数；AD: 支持变异的碱基数。

	tumor	normal
RD	47	99
AD	74	111

Fisher’s Exact Test接受2×2行列表作为输入，并按如下方式计算P-value值：
1、计算行列表中每一行，每一列的总和以及观察总数。
2、给定行和列总和，如果原假设为真，则使用超几何概率函数计算条件概率，以观察行列表中的准确结果。条件概率为

图2.条件概率计算公式

其中R1和R2是行总和，C1和C2是列总和，N是行列表中观测的总数，nij是表中第i行和第j列的值。
3、查找与行和列之和一致的所有可能的非负整数矩阵。对于每个矩阵，使用P公式计算相关的条件概率。
4、根据感兴趣的替代假设，使用这些值来计算检验的p值。
(1)对于双边测试，对于观察到的行列表，将所有小于或等于Pcutoff的条件概率求和。
(2)对于左侧测试，将（1,1）象限频率小于或等于n11的所有矩阵的条件概率求和。
(3)对于右侧测试，在观察到的行列表中求和（1,1）单元频率大于或等于n11的所有矩阵的条件概率。

VarScan2 somaitc计算的是单侧测试，如果按下表排列数据是左侧测试。

表2. 肿瘤样本与正常样本覆盖深度2×2行列表。

	tumor	normal
RD	47	99
AD	74	111

那怎么计算P-value呢？根据P-value的定义：在原假设为真的前提下，出现该样本或比该样本更极端的结果的概率之和。上表就是所有tumor Ref depth小于等于47的概率之和。

(P根据公式1计算)

P-value=0.0883，这P的意思是假设tumor Alt与normal Alt有同样的频率，但这的假设可能性为0.0883。按通常0.05或0.01阈值来cutoff的话，应该接受这一假设。肿瘤样本中测到的位点是种系变异而非somatic。

再看一个测序位点

表3.肿瘤样本与正常样本覆盖深度2×2行列表。

	tumor	normal
RD	144	230
AD	39	0

P-value=1.3226e-15，这P值就足够小，认为是种系变异就不合适了，可以判断为肿瘤特有的变异somatic。

按照上述方法对肿瘤样本与正常样本所有同时覆盖的位点都计算一个P-value值，当P-value值小于设定值时判定位肿瘤体细胞变异。因影响NGS测序结果的因素多，reads覆盖的概率分布并不跟理想模型一致。如果按照常见0.05值来判定somatic，结果会比预期多。可以通过更低的P-value值或其他条件过滤掉可能的假阳性位点。按本人计算经验，P-value设0.001跟mutect2设的默认参数差不多。如果somatic是10个，覆盖到的区域2M，那么TMB=10/2=5。

TMB的计算实操

1. 软件安装

在安装这个软件maftools时，需要先安装BioManager,然后在安装maftools，如下：

if (!require("BiocManager")) {install.packages("BiocManager")
}
if (!require(maftools)) {BiocManager::install("maftools")
}
library(maftools)

2. 数据读取

maftools工具需要读入两个文件，如下：

1.MAF文件-可以是gz压缩。必需的；

2.与MAF中每个样本/肿瘤样本条码相关的可选但推荐的临床数据；

maf文件格式

MAF文件包含许多字段，从染色体名称到cosmic注释。然而，maftools中的大多数分析使用以下列如下：1.Hugo_Symbol；

2.Chromosome；

3.Start_Position；

4.End_Position；

5.Reference_Allele；

6.Tumor_Seq_Allele2；

7.Variant_Classification；

8.Variant_Type；

9.Tumor_Sample_Barcode.

同时读取maf文件和临床信息文件，看下结果，如下：

# path to TCGA LAML MAF file
laml.maf = system.file("extdata", "tcga_laml.maf.gz", package = "maftools")
# clinical information containing survival information and histology. This is
# optional
laml.clin = system.file("extdata", "tcga_laml_annot.tsv", package = "maftools")
laml = read.maf(maf = laml.maf, clinicalData = laml.clin)
## -Reading
## -Validating
## -Silent variants: 475
## -Summarizing
## -Processing clinical data
## -Finished in 1.940s elapsed (0.290s cpu)
print(laml@data[1:5, ])
##    Hugo_Symbol Entrez_Gene_Id           Center NCBI_Build Chromosome
## 1:      ABCA10          10349 genome.wustl.edu         37         17
## 2:       ABCA4             24 genome.wustl.edu         37          1
## 3:      ABCB11           8647 genome.wustl.edu         37          2
## 4:       ABCC3           8714 genome.wustl.edu         37         17
## 5:       ABCF1             23 genome.wustl.edu         37          6
##    Start_Position End_Position Strand Variant_Classification Variant_Type
## 1:       67170917     67170917      +            Splice_Site          SNP
## 2:       94490594     94490594      +      Missense_Mutation          SNP
## 3:      169780250    169780250      +      Missense_Mutation          SNP
## 4:       48760974     48760974      +      Missense_Mutation          SNP
## 5:       30554429     30554429      +      Missense_Mutation          SNP
##    Reference_Allele Tumor_Seq_Allele1 Tumor_Seq_Allele2 Tumor_Sample_Barcode
## 1:                T                 T                 C         TCGA-AB-2988
## 2:                C                 C                 T         TCGA-AB-2869
## 3:                G                 G                 A         TCGA-AB-3009
## 4:                C                 C                 T         TCGA-AB-2887
## 5:                G                 G                 A         TCGA-AB-2920
##    Protein_Change i_TumorVAF_WU i_transcript_name
## 1:        p.K960R      45.66000       NM_080282.3
## 2:       p.R1517H      38.12000       NM_000350.2
## 3:       p.A1283V      46.97218       NM_003742.2
## 4:       p.P1271S      56.41000       NM_003786.1
## 5:        p.G658S      40.95000    NM_001025091.1

maftools 计算 TMB

仅计算TMB

我们只需要调取maftools软件包里面的tmb()函数即可完成计算，如下：

tmb1 = tmb(maf = laml)  ##取log值

head(tmb1)
##    Tumor_Sample_Barcode total total_perMB total_perMB_log
## 1:         TCGA-AB-2903     0        0.00            -Inf
## 2:         TCGA-AB-2823     1        0.02        -1.69897
## 3:         TCGA-AB-2834     1        0.02        -1.69897
## 4:         TCGA-AB-2840     1        0.02        -1.69897
## 5:         TCGA-AB-2842     1        0.02        -1.69897
## 6:         TCGA-AB-2848     1        0.02        -1.69897

第三列就是TMB值了,可以用四分位数看看tmb值的分布情况：

quantile(tmb1$total_perMB)
##   0%  25%  50%  75% 100%
## 0.00 0.10 0.18 0.24 0.68

不取TMB log值，则设为F

tmb2 = tmb(maf = laml, logScale = F)  ##不去log值

head(tmb2)
##    Tumor_Sample_Barcode total total_perMB total_perMB_log
## 1:         TCGA-AB-2903     0        0.00            -Inf
## 2:         TCGA-AB-2823     1        0.02        -1.69897
## 3:         TCGA-AB-2834     1        0.02        -1.69897
## 4:         TCGA-AB-2840     1        0.02        -1.69897
## 5:         TCGA-AB-2842     1        0.02        -1.69897
## 6:         TCGA-AB-2848     1        0.02        -1.69897
quantile(tmb2$total_perMB)
##   0%  25%  50%  75% 100%
## 0.00 0.10 0.18 0.24 0.68

与TCGA数据库比较TMB

tcgaCompare()使用突变负荷比较33个TCGA队列的突变负荷，但是这里需要指出我们测序使用的panel大小是多少Mb,比如这里例子中是 capture_size = 50 也就是50mb大小的区域。

laml.mutload = tcgaCompare(maf = laml, cohortName = "Example-LAML", logscale = TRUE,capture_size = 50)

该函数得到的结果可以获取TMB，如下：

head(laml.mutload$mutation_burden_perSample)
##    Tumor_Sample_Barcode total       cohort total_perMB
## 1:         TCGA-AB-3009    34 Example-LAML        0.68
## 2:         TCGA-AB-2807    25 Example-LAML        0.50
## 3:         TCGA-AB-2959    23 Example-LAML        0.46
## 4:         TCGA-AB-3002    21 Example-LAML        0.42
## 5:         TCGA-AB-2849    20 Example-LAML        0.40
## 6:         TCGA-AB-2923    20 Example-LAML        0.40

利用maftools软件包完全可以轻松实现TMB的计算，后续我们在分享 MSI的计算方法，这样我们就可以充分利用两种指标来评估患者的PD-L1的治疗！！

References:

1. Spigel, David R. et al. “Total mutation burden (TMB) in lung cancer (LC) and relationship with response to PD-1/PD-L1 targeted therapies.” Journal of Clinical Oncology 34 (2016): 9017-9017.

2. Rizvi NA, Hellmann MD, Snyder A, et al. Cancer immunology. Mutational landscape determines sensitivity to PD-1 blockade in non-small cell lung cancer. Science. 2015;348(6230):124-128. doi:10.1126/science.aaa1348

3. Lawrence MS, Stojanov P, Polak P, et al. Mutational heterogeneity in cancer and the search for new cancer-associated genes. Nature. 2013;499(7457):214-218. doi:10.1038/nature12213

桓峰基因，铸造成功的您！

未来桓峰基因公众号将不间断的推出单细胞系列生信分析教程，

敬请期待！！

有想进生信交流群的老师可以扫最后一个二维码加微信，备注“单位+姓名+目的”，有些想发广告的就免打扰吧，还得费力气把你踢出去！