我用了10分钟，拆解了肿瘤免疫高分SCI的万用分析，简单又高级！

肿瘤突变负荷TMB的计算及其多组学分析

大家好，我是阿琛。随着免疫治疗以及个性化治疗方案的不断进展，肿瘤突变研究逐渐被大家所熟悉。其中，基于肿瘤突变以及产生的肿瘤新生抗原，被应用于各种临床试验和基础研究中。

Tumor mutational burden is the total number of mutations per coding area of a tumor genome.

根据定义，我们可以看到，肿瘤突变负荷（Tumor Mutational Burden， TMB）是指肿瘤中体细胞内编码蛋白的碱基突变数量，能够反映肿瘤细胞中所包含的突变数量，包括常见的替换，插入，缺失等各种类型。在多项临床研究中显示，TMB值与免疫治疗的效果存在显著的相关性。

下面，我们一起来看一下如何计算患者的TMB值。

1. TCGA突变数据的下载

对于TCGA数据库中数据的下载，除了直接从TCGA GDC（https://portal.gdc.cancer.gov/）外，一般常用的方法是借助R的TCGAbiolinks包来进行下载。

library(TCGAbiolinks) luad.varscan2.maf <- GDCquery_Maf("LUAD", pipelines = "varscan2") save(luad.varscan2.maf, file = "LUAD.maf.rda") #保存突变结果

这里，我们通过TCGAbiolinks包的GDCquery_Maf()函数来进行下载；其中pipelines参数可以是muse，varscan2，somaticsniper，mutect2中的一种，在此我们选择varscan2数据来进行下载。

2. 计算患者突变数量

TMB was calculated based on all somatic mutations, and then recalculated based on nonsynonymous mutations.

对于TMB值的计算，首先第一步需要获取可以引起氨基酸发生改变的突变数量。

library(maftools) rt <-read.maf(maf = "TCGA.LUAD.varscan.somatic.maf", isTCGA = T)

使用maftools包的read.maf()函数读取突变数据。

table(rt@data$Mutation_Status)

结果显示，所有突变类型都是somatic mutation。

接下来，我们进一步对每个样本的突变数量进行统计分析。

方法一

对于第一种方法，我们可以直接使用maftools包在读取数据过程中得到的统计结果。

variants_1 <-rt@variants.per.samplevariants_1[1:3,] dim(variants_1)

结果显示，其中包含了560名患者的突变数量总数。

方法二

当然，除了第一种方法外，我们也可以根据TMB的定义，自己计算患者的突变总数。

table(rt@data$Variant_Classification)

可以看到，在数据集中，包含了9种不同类型的突变形式。

library(tidyverse) variants_2 <- rt@data %>% group_by(Tumor_Sample_Barcode) %>% summarise(TCGA_sum = n()) head(variants_2)

通过group_by()和summarise()函数联用，从而统计每个患者对应的突变总数。

随后，我们将两个分析结果根据患者id进行合并，并比较一下两种方法得到的结果。

mut <- merge(variants_1, variants_2, by = "Tumor_Sample_Barcode") colnames(mut) <- c("Tumor_Sample_Barcode","variants_1","variants_2")

cor(as.numeric(mut$variants_1),as.numeric(mut$variants_2))

结果显示，两种方法计算得到的结果是完全一致的。

3. TMB值的计算

mut$TMB <- mut$variants_2%/%35

由于TMB值往往通过每MB突变数来进行衡量展示，因此，我们需要将计算得到的每个患者的突变综述除以总的MB数。目前，因为各家说法不一，对于人类总MB数存在多种理论，包括常见的35MB和38MB等。在此，我们使用35MB来进行计算。

根据统计得到的突变总数，最终得到每个患者对应的TMB值。

4. TCGA表达数据的下载

为了进一步计算患者某个基因的表达高低与TMB之间的关系，我们首先下载对应的基因表达集。

library(TCGAbiolinks) expquery <- GDCquery(project = "TCGA-LUAD", data.category = "Transcriptome Profiling", data.type = "Gene Expression Quantification", workflow.type = "HTSeq - FPKM") GDCdownload(expquery) expquery2 <- GDCprepare(expquery) save(expquery2,file = "luad.gdc.rda") #保存表达结果 #load("luad.gdc.rda") expMatrix <- TCGAanalyze_Preprocessing(expquery2)

通过TCGAbiolinks包的GDCquery()函数下载TCGA-LUAD患者对应的表达矩阵。

group_list <- ifelse(substr(colnames(expMatrix),14,14)=='0','tumor','normal') table(group_list)

结果显示，其中包含59例正常样品，535例肿瘤样品。

expMatrix_tumor <- expMatrix[,group_list=='tumor'] exp <- as.data.frame(expMatrix_tumor["ENSG00000141510",]) colnames(exp) <- c("ENSG00000141510")

提取表达矩阵中肿瘤组织的表达数据，并进一步提取基因TP53（ENSG00000141510）的表达。

exp$Type <- ifelse(exp$ENSG00000141510 > median(exp$ENSG00000141510), "TP53_high", "TP53_low") exp$Tumor_Sample_Barcode <- str_sub(rownames(exp), 1, 12) head(exp)

随后，根据TP53基因的表达中位值，将患者分成TP53_high组和TP53_low组。

5. 可视化展示

最后，根据表达分组，进行可视化展示，来形象地展示TP53基因表达与TMB值之间的关系。

TMB <- merge(mut, exp, by="Tumor_Sample_Barcode")

基于Tumor_Sample_Barcode中患者的id号，将突变数据和表达数据进行合并。

library(ggplot2) library(ggstatsplot) ggbetweenstats(TMB, x = "Type", y = "TMB", mean.ci = FALSE, type = "np", #默认为p，代表参数检验；np代表非参数检验 pairwise.comparisons = TRUE, pairwise.display = "s", p.adjust.method = "fdr", messages = FALSE)

使用ggstatsplot包的ggbetweenstats()函数绘制箱线图。关于ggstatsplot包的使用方法，大家可以参考之前的教程（小白绘图专题|ggstatsplot绘制多种发表级图形）。结果显示，在TCGA-LUAD患者中，TP53基因的高低表达对患者的TMB值无显著性影响。