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"
  )
  # 读取文件
  df <- read_csv(file, col_select = use_cols)
  data <- df %>% subset(!Variant_Classification %in% mut_type) %>%
    # 计算 VAF
    mutate(vaf = t_alt_count / t_depth) %>%
    group_by(Tumor_Sample_Barcode) %>%
    summarise(mut_num = n(), TMB = mut_num / 30, MaxVAF = max(vaf))
  return(data)
}
计算
paad.filtered <- get_TMB(maf)

我们对突变类型进行了过滤，删除一些无法反映到转录组上的突变。

外显子的总长度我使用的是 30M，我看取 38M、40M 的都有，不是很统一。如果我们以 TMB 中位值来区分高低 TMB 组的话，并不是选用一个固定的值来分，其实没什么影响。

在前面我们看到了一个样本的突变负荷特别高，我们来看看对应的 TMB 值为多少

max(paad.filtered$TMB)
# [1] 425.5333

现在，先删除该样本，并根据 TMB 中位值将样本分为高低两组

paad.tmb <- paad.filtered %>% 
  filter(TMB != max(TMB)) %>%
  mutate(
    label = if_else(TMB >= median(TMB), "High", "Low"), 
    .before = "MaxVAF"
  ) %>%
  mutate(Tumor_Sample_Barcode = substr(Tumor_Sample_Barcode, 1, 12))
max(paad.tmb$TMB)
# [1] 3.333333

然后对两类样本进行生存分析

# 生存分析
library(survival)
library(survminer)

data.surv <- clin %>% filter(!is.na(time)) %>%
  mutate(
    time = time / 30, 
    status = if_else(status == "Alive", 0, 1)
  ) %>% 
  inner_join(paad.tmb, by = "Tumor_Sample_Barcode") %>%
  dplyr::select(c("Tumor_Sample_Barcode", "time", "status", "label"))

km_fit <- survfit(Surv(time, status) ~ label, data = data.surv)

ggsurvplot(
  km_fit, data = data.surv,
  pval = TRUE, surv.median.line = "hv",
  legend.labs=c("TMB-High","TMB-Low"),
  legend.title="Group",
  title="Overall survival",
  xlab = "Time(month)",
  risk.table = TRUE
)

两组之间的生存率没什么差别，可能和胰腺导管腺癌本身的存活率就不高有关吧

这里还可以看看 TMB 的高低是否与其他临床特征相关，例如性别、年龄、TNM 分期和分级等指标。需要对前面的临床数据处理一下，这里我就不再举例了。

算了，还是举个例子吧，先把对应的数据挑出来，并整理一下列名

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"
    )) %>%
  rename_with(~ c("Tumor_Sample_Barcode", "time", "status", "T", "N", "M", "stage", "gender"))

处理一下临床信息，去掉 NA 或 不明确的样本

data.cate <- inner_join(paad.tmb, clin, by = "Tumor_Sample_Barcode") %>%
  filter_at(vars(T, N, M), all_vars(!endsWith(., "X") & !is.na(.))) %>%
  mutate(
    N = gsub("(N\\d).*", "\\1", N, perl = TRUE),
    stage = if_else(startsWith(stage, c("Stage III", "Stage IV")), "Stage III-IV", "Stage I-II")
  )

TMB 与性别之间的关系

ggboxplot(data.cate, x = "gender", y = "TMB",
          fill = "gender") +
  stat_compare_means(comparisons = list(
    c("female", "male")
  )) +
  stat_compare_means(label.y = 4.5) 

与 T 分期之间的关系

ggboxplot(data.cate, x = "T", y = "TMB",
          fill = "T") +
  stat_compare_means(comparisons = list(
    c("T2", "T3"), c("T2", "T4"), c("T3", "T4")
  )) +
  stat_compare_means(label.y = 4.5) 

差异分析

识别差异表达基因

根据 TMB 值将样本分组之后，就可以使用表达数据来识别两组之间的差异表达基因了。如何识别差异表达基因我们前面已经介绍过了，直接上代码

get_prep <- function(cancer) {
  query <- GDCquery(
    project = cancer,
    data.category = "Gene expression",
    data.type = "Gene expression quantification",
    platform = "Illumina HiSeq",
    file.type  = "results",
    sample.type = c("Primary Tumor"),
    legacy = TRUE
  )
  # 选择 20 个样本
  query$results[[1]] <-  query$results[[1]][1:20,]
  GDCdownload(query)
  # 获取 read count
  exp.count <- GDCprepare(
    query,
    summarizedExperiment = TRUE,
  )
  return(exp.count)
}
paad.exp <- get_prep("TCGA-PAAD")

prepData <- TCGAanalyze_Preprocessing(
    object = paad.exp,
    cor.cut = 0.6,
    datatype = "unstranded"
)
# 标准化
dataNorm <- TCGAanalyze_Normalization(
    tabDF = prepData,
    geneInfo = TCGAbiolinks::geneInfoHT,
)
# 分位数过滤
dataFilt <- TCGAanalyze_Filtering(
    tabDF = dataNorm,
    method = "quantile",
    qnt.cut =  0.25
)
colnames(dataFilt) <- substr(colnames(dataFilt), 1, 12)
sample.high <- subset(paad.tmb, label == "High")$Tumor_Sample_Barcode
sample.low <- subset(paad.tmb, label == "Low")$Tumor_Sample_Barcode

dataHigh <- dataFilt[,colnames(dataFilt) %in% sample.high]
dataLow <- dataFilt[,colnames(dataFilt) %in% sample.low]

DEGs <- TCGAanalyze_DEA(
    mat1 = dataLow,
    mat2 = dataHigh,
    metadata = FALSE,
    Cond1type = "Low",
    Cond2type = "High",
    fdr.cut = 0.01,
    logFC.cut = 1,
    method = "glmLRT"
)
# 差异水平
DEGsFiltLevel <- TCGAanalyze_LevelTab(
    DEGs, "TMB High", "TMB Low",
    dataHigh, dataLow
)

为了节省时间，我们只选了前 20 个样本来进行分析

使用 edgeR 包，控制 FDR 为 0.01 且 logFC 为 1，共识别出 209 个差异表达基因

DEGs[1:3, 1:5]
#                    logFC     logCPM       LR       PValue         FDR
# ENSG00000005073 4.507216  0.3712662 21.90683 2.862136e-06 0.005009341
# ENSG00000073282 3.688089  3.3028786 22.59479 2.000264e-06 0.004157298
# ENSG00000104827 5.072573 -0.6600451 24.51762 7.363358e-07 0.002040509
dim(DEGs)
# [1] 28  7

2 通路富集

通路富集分析我们在前面也已经介绍了，不再详细说明

library(clusterProfiler)
library(org.Hs.eg.db)
library(enrichplot)

gene.id <- bitr(
  DEGs$gene_name, fromType = "SYMBOL",
  toType = "ENTREZID",
  OrgDb = org.Hs.eg.db
) 
# go 
go <- enrichGO(
  gene = gene.id$ENTREZID,
  OrgDb = org.Hs.eg.db,
  ont = "ALL",
  pAdjustMethod = "BH",
  qvalueCutoff = 0.05,
  readable = T
)

dotplot(go)

这些通路是不是与 PAAD 有关呢？暂时还不知道，可能要查下文献或数据库

# kegg
kegg <- enrichKEGG(
  gene = gene.id$ENTREZID,
  organism = "hsa",
  pvalueCutoff = 0.05
)

没有富集到 KEGG 通路。

再进行 GSEA 富集分析

gene_info <- DEGs %>%
    filter(gene_name %in% gene.id$SYMBOL) %>%
    inner_join(., gene.id, by = c("gene_name" = "SYMBOL")) %>%
    arrange(desc(logFC))

geneList <- gene_info$logFC
names(geneList) <- gene_info$ENTREZID

go2 <- gseGO(
  geneList     = geneList,
  OrgDb        = org.Hs.eg.db,
  ont          = "ALL",
  minGSSize    = 100,
  maxGSSize    = 500,
  pvalueCutoff = 0.05,
  verbose      = FALSE
)

kegg2 <- gseKEGG(
  geneList = geneList,
  organism     = 'hsa',
  minGSSize    = 120,
  pvalueCutoff = 0.05,
  verbose      = FALSE
)

一条都没富集到，这时，可以将 FDR 阈值设置为 0.1 再看看结果，这里我就不再演示了，还是继续使用这个结果进行后续分析

3 免疫相关基因

首先，我们从 IMMPORT 网站 https://www.immport.org/shared/home 中下载免疫相关基因列表

点击 Resources

往下翻，找到 Gene Lists

然后点击下载

下载后的文件后缀为 json，但是格式还是 txt 格式的，所以我将后缀改成了 txt 格式，再读取基因列表，

immune_gene <- read_delim("~/Downloads/GeneList.txt", delim = "\t")

找出差异基因中的免疫相关基因

(gene_list <- intersect(DEGs$gene_name, immune_gene$Symbol))
# [1] "CGB3"   "SLURP1" "IL20RB"

看看这些基因的表达差异情况

dplyr::filter(DEGs, gene_name %in% gene_list)
#                    logFC     logCPM       LR       PValue         FDR gene_name      gene_type
# ENSG00000104827 5.072573 -0.6600451 24.51762 7.363358e-07 0.002040509      CGB3 protein_coding
# ENSG00000126233 5.041713 -0.3448266 22.27137 2.367095e-06 0.004630317    SLURP1 protein_coding
# ENSG00000174564 3.151819  3.0399704 21.07785 4.409955e-06 0.006015362    IL20RB protein_coding

两个下调，一个上调。这里需要说明的是，在识别差异基因时，两个分组的前后顺序表明了差异的方向，这里的上调，是 TMB High 组相对于 Low 组上调了，记住上下调是谁相对于谁的。是参数中 Cond2type 相对于 Cond1type 而言。

我们来看看这些基因的表达情况

genes <- dplyr::filter(DEGs, gene_name %in% gene_list) %>%
    dplyr::select(gene_name)

dataFilt[rownames(genes),] %>%
    log1p %>% as.data.frame() %>%
    mutate(Gene = genes$gene_name) %>%
    pivot_longer(cols = -Gene, names_to = "sample", values_to = "expression") %>%
    mutate(TMB = if_else(sample %in% sample.low, "Low", "High")) %>%
    ggplot(aes(Gene, expression, fill = TMB)) +
    geom_boxplot()

我们挑两个基因来看看它们与生存之间有没有关系，就拿 CGB3 和 IL20RB 两个基因吧，logFC 值最大和最小

gene.sur <- dataFilt[rownames(genes),] %>%
    log1p %>% t() %>% as.data.frame() %>%
    rename_with(~genes$gene_name) %>%
    mutate(
        CGB3 = if_else(CGB3 < median(CGB3), 0, 1),
        IL20RB = if_else(IL20RB < median(IL20RB), 0, 1)
    ) %>%
    rownames_to_column("Tumor_Sample_Barcode") %>%
    inner_join(data.surv, by="Tumor_Sample_Barcode")
km.CGB3 <- survfit(Surv(time, status) ~ CGB3, data = gene.sur)
km.IL20RB <- survfit(Surv(time, status) ~ IL20RB, data = gene.sur)

绘制生存曲线

ggsurvplot(
  km.CGB3, data = gene.sur,
  pval = TRUE, surv.median.line = "hv",
  legend.labs=c("High expression", "Low expression"),
  legend.title="",
  title="Overall survival",
  xlab = "Time(month)",
  risk.table = TRUE
)

虽然从图像上看，生存曲线是分开的，但是 P 值并不显著，可能由于样本量太少了，没有显著性，也是类似的。

也可以对这几个基因做一个 cox 回归分析，看下生存会不会好一些。

免疫浸润

我们使用 immunedeconv 包来进行免疫浸润分析，该包提供了很多免疫浸润算法的支持，包括

• quantiseq
• timer
• cibersort
• cibersort_abs
• mcp_counter
• xcell
• epic

该包需要从 GitHub 上安装

install.packages("remotes")
remotes::install_github("icbi-lab/immunedeconv")

library(immunedeconv)

因为需要输入的是 TPM 标准化且未 log 转换的表达值，我们需要重新整理一下数据

1. 获取基因的 TPM

本来想要通过基因的长度来标准化，但是整了一下觉得太麻烦了。想想还是用 TPM 来转换比较方便。

例如，直接在 TCGAanalyze_Preprocessing 函数中指定要获取的数据类型

padd.tpm <- TCGAanalyze_Preprocessing(
    object = paad.exp,
    cor.cut = 0.6,
    datatype = "tpm_unstrand"
)

exp <- as.data.frame(padd.tpm) %>%
    tibble::rownames_to_column("gene_id") %>%
    inner_join(dplyr::select(TCGAbiolinks::get.GRCh.bioMart("hg38"), gene_id, gene_name)) %>%
    dplyr::select(-gene_id) %>%
    group_by(gene_name) %>%
    summarise_all(mean) %>%
    tibble::column_to_rownames(var = "gene_name") %>%
    rename_with(function(x) substr(x, 1, 12))

2. TIMER

我们先进行 TIMER 分析

library(immunedeconv)

immune.timer <- deconvolute(
  exp, "timer", 
  indications = rep("PAAD", 20)
)

以堆积条形图来展示样本中免疫细胞的占比

immune.timer %>%
  pivot_longer(cols = -cell_type, names_to = "sample", values_to = "percent") %>%
  ggplot(aes(sample, percent, fill = cell_type)) +
  geom_bar(stat='identity') +
  coord_flip()

TIMER 官网的 SCNA 模块可以根据给定基因，来比较不同肿瘤免疫浸润水平下的体细胞拷贝数变异（CNV）

首先，进入官网: https://cistrome.shinyapps.io/timer/

然后，选择 SCNA，再选择相应的癌型和基因，最后点击提交，就可以获取到最下方的分组箱线图

我们前面找到了几个免疫相关的差异表达基因，我们可以分别下载这几个基因的 CNV 与免疫细胞的分析结果

IL20RB 基因的结果如下

还可以对这 6 种免疫细胞的比例与生存时间做一个 cox 回归模型，如

immune.os <- immune.timer %>%
  column_to_rownames("cell_type") %>% 
  t() %>% as.data.frame() %>%
  rownames_to_column("Tumor_Sample_Barcode") %>%
  inner_join(data.surv, by = "Tumor_Sample_Barcode")
  
cox <- coxph(Surv(time, status) ~ `B cell` + `T cell CD4+` + `T cell CD8+` + 
               Neutrophil + Macrophage + `Myeloid dendritic cell`, 
             data = immune.os)
cox_fit <- survfit(cox)

但是看了一下没啥用，也可以当做一种尝试吧

3. CIBERSORT

再使用 CIBERSORT 软件进行免疫浸润分析

set_cibersort_binary("~/Documents/WorkSpace/RStudio/Cibersort.R")
set_cibersort_mat("~/Documents/WorkSpace/RStudio/LM22.txt")

immune_cibersort <- deconvolute(exp, "cibersort")

可视化细胞占比

immune_cibersort %>%
    pivot_longer(cols = -cell_type, names_to = "sample", values_to = "percent") %>%
    ggplot(aes(sample, percent, fill = cell_type)) +
    geom_bar(stat='identity') +
    theme(axis.text.x = element_text(angle=60, vjust=1, hjust=1)) +
    guides(fill = guide_legend(ncol = 1, byrow = TRUE))

我们可以比较 TMB 高低两组样本之间，免疫细胞占比是否存在差异。我们使用 wilcox 秩和检验

immune.tmb <- immune_cibersort %>%
    pivot_longer(cols = -cell_type, names_to = "sample", values_to = "percent") %>%
    mutate(label = if_else(sample %in% sample.high, 1, 0))


label_high <- group_by(immune.tmb, cell_type) %>%
    summarise(y = max(percent))

ggplot(immune.tmb, aes(cell_type, percent, fill = factor(label))) +
    geom_boxplot(
        position = position_dodge2(preserve = "single", width = 0.5)) +
    stat_compare_means(
        aes(group = label, label = paste0("p = ", ..p.format..)), label.y = label_high$y) +
    scale_fill_manual(name = "TMB", labels = c('0' = "Low", '1' = "High"), values = 2:3) +
    theme(axis.text.x = element_text(angle=60, vjust=1, hjust=1))

虽然我们从图上看有些细胞的占比之间差别比较明显，但是 p 值并不显著，可能还是由于样本量太少的原因吧