TCGA 数据分析实战 —— TMB 与免疫浸润联合分析

TCGA 数据分析实战 —— TMB 与免疫浸润联合分析

前言

近年来，随着免疫检查点抑制剂的兴起，大大改变了传统的肿瘤治疗策略，尽管 PD-L1 和 dMMR 的检测都获得了 FDA 的批准，提高了免疫药物的响应和获益，但它们都有自身的不足。

各种检测方法判定的 PD-L1 水平不一致率较高，dMMR 在各种不同的癌种中的携带比例差异较大

而免疫治疗的效果主要是免疫细胞对癌细胞特异性抗原的识别，那么从理论上来说，那些携带基因突变越多的癌症患者，癌细胞产生的新抗原越多，被免疫细胞识别的可能性更高。也就是说，肿瘤组织的突变负荷（tumor mutational burden，TMB）越高，患者或许能从免疫治疗中获得更多的收益。

TMB 是指肿瘤细胞基因组中，基因的外显子编码区每兆碱基中发生置换和插入或缺失突变的总数。

接下来，我们就来看看，如何使用 TCGA 中的数据来分析 TMB 与免疫浸润之间的关系。

TMB

要计算样本的 TBM 值，必须先获取样本的突变数据，TCGA 的突变数据是 MAF 格式的，我们可以使用 TCGAbiolinks 下载对应癌型的突变数据。

TCGA 的突变流程有 4 种，分别是：muse、varscan2、somaticsniper、mutect2。其中 muse 和 somaticsniper 只能计算点突变，无法识别 indel。

我们就以 PAAD 为例来进行说明吧

library(TCGAbiolinks)
library(tidyverse)

query <- GDCquery(
    project = "TCGA-PAAD", 
    data.category = "Simple Nucleotide Variation",
    data.type = "Masked Somatic Mutation",
    access = "open"
)
GDCdownload(query)
maf <- GDCprepare(query)

maf 最后一列为该突变在哪种突变 calling 算法中识别出的，可以根据需要利用这一列的数据进行过滤。

有一个叫 maftools 的 R 包，可以用来处理 MAF 格式的数据。使用 read.maf 函数读取刚才获取的突变数据。

如果涉及到临床相关的分析，可以为 clinicalData 参数指定对应的临床信息。例如

library(maftools)

paad.clin <- GDCquery_clinic("TCGA-PAAD")
clin <- paad.clin %>% 
  dplyr::select(c(
    "submitter_id", "days_to_last_follow_up", "vital_status",
    "ajcc_pathologic_t", "ajcc_pathologic_n", "ajcc_pathologic_m",
    "ajcc_pathologic_stage", "gender"
    )) %>%
  `names<-`(c("Tumor_Sample_Barcode", "time", "status", "T", "N", "M", "stage", "gender"))

paad.maf <- read.maf(
    maf,
    clinicalData = clin,
    isTCGA = TRUE
)

如果还有 CNV 数据，也可以传递给 cnTable 参数。

有很多函数可以获取对应的信息，如

# 展示样本层面的信息
getSampleSummary(paad.maf)
# 展示基因层面的信息
getGeneSummary(paad.maf)
# 展示临床相关信息
getClinicalData(paad.maf)
# 获取 MAF 文件的所有列字段
getFields(paad.maf)
# 保存为文件
write.mafSummary(maf = paad.maf, basename = 'paad.maf')

1. 绘图

使用 plotmafSummary 来展示 MAF 文件的一些汇总信息，例如

plotmafSummary(maf = paad.maf, rmOutlier = TRUE, addStat = 'median', dashboard = TRUE, titvRaw = FALSE)

从图中可以看到，错义突变占据绝大部分的比例，C>T 的点突变类型是最多的，KRAS 和 TP53 两个基因的突变次数最多

绘制 oncoplots 图，只展示前 15 个基因

vc_cols <- RColorBrewer::brewer.pal(n = 8, name = 'Paired')
names(vc_cols) <- c(
  'Frame_Shift_Del',
  'Missense_Mutation',
  'Nonsense_Mutation',
  'Multi_Hit',
  'Frame_Shift_Ins',
  'In_Frame_Ins',
  'Splice_Site',
  'In_Frame_Del'
)
oncoplot(maf = paad.maf, colors = vc_cols, top = 15)

瀑布图展示了突变频率前 15 的基因在样本中的突变情况，上方的条形图展示的是样本的突变负荷，右侧的直方图展示了基因的每种突变类型的情况

图中，有一个样本的突变负荷特别高，在后续分析时应该考虑是否将其作为离群点去除

可以使用 somaticInteractions 函数检测互斥突变或同时突变的基因，该函数使用配对的 Fisher’s Exact 检验识别显著的基因对。

par(oma = c(3, 4, 5, 1))
somaticInteractions(maf = paad.maf, top = 25, pvalue = c(0.05, 0.1))

可以看到突变率前 10 的基因中，同时出现突变的基因对偏多

绘制基因的 Variant Allele Frequencies(VAF) 箱线图可以看出基因的克隆状态，在理想情况下，样本中克隆基因的平均等位基因频率约为 50%

paad.maf@data[["TumorVAF"]] <- paad.maf@data$t_alt_count / paad.maf@data$t_depth
plotVaf(maf = paad.maf, vafCol = 'TumorVAF')

2. 计算 TMB

使用 maftools 的 tcgaCompare 函数，可以将我们获取的 MAF 文件中的突变负荷与 TCGA MC3 项目中的 33 种癌型的突变负荷进行比较

maf.TMB <- tcgaCompare(
  maf = paad.maf, cohortName = 'PAAD-maf', 
  logscale = TRUE, capture_size = 50
)

可以看到，与 MC3 中胰腺导管腺癌的突变负荷差一些

下面，我们来计算每个样本的 TMB 值

get_TMB <- function(file) {
   
  # 需要用到的列
  use_cols <- c(
    "Hugo_Symbol", "Variant_Classification", "Tumor_Sample_Barcode", 
    "HGVSc", "t_depth", "t_alt_count"
  )
  # 删除这些突变类型
  mut_type <- c(
    "5'UTR", "3'UTR", "3'Flank", "5'Flank", "Intron", "IGR",
    "Silent", "RNA", "Splice_Region"