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

2021-07-31 本文已影响0人名本无名

前言

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

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

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

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

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

TMB

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

TCGA 的突变流程有 4 种，分别是：muse、varscan2、somaticsniper、mutect2

其中 muse 和 somaticsniper 只能计算点突变，无法识别 indel，所以不用考虑。我们就以 PAAD 的 mutect 为例来进行说明吧

library(TCGAbiolinks)
library(tidyverse)

GDCquery_Maf(
  tumor = "PAAD",
  save.csv = TRUE,
  directory = "~/Downloads/TCGA",
  pipelines = "mutect2"
)

默认会下载一个 .maf.gz 的文件，我们将其解析并保存为一个 csv 文件

有一个叫 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(
  "~/Downloads/TCGA/TCGA-PAAD/harmonized/Simple_Nucleotide_Variation/Masked_Somatic_Mutation/fea333b5-78e0-43c8-bf76-4c78dd3fac92/TCGA.PAAD.mutect.fea333b5-78e0-43c8-bf76-4c78dd3fac92.DR-10.0.somatic.maf.gz",
  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)
}

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

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

我们直接读取刚才下载 MAF 文件时保存的整理后的 CSV 文件

paad.csv <- get_TMB('~/Downloads/TCGA/TCGA.PAAD.mutect.fea333b5-78e0-43c8-bf76-4c78dd3fac92.DR-10.0.somatic.maf.csv')

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

> max(paad.csv$TMB)
[1] 443.4667

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

paad.tmb <- paad.csv %>% 
  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))

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

# 生存分析
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"
    )) %>%
  `names<-`(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)

差异分析

1 识别差异表达基因

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

get_count <- 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_count("TCGA-PAAD")

# 预处理
prepData <- TCGAanalyze_Preprocessing(
  object = paad.exp,
  cor.cut = 0.6,
  datatype = "raw_count"
)

# 标准化
dataNorm <- TCGAanalyze_Normalization(
  tabDF = prepData,
  geneInfo = TCGAbiolinks::geneInfo,
)
# 分位数过滤
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 个差异表达基因

> head(DEGs)
          logFC   logCPM       LR       PValue          FDR
ADCY1 -3.454446 4.629218 19.14697 1.210294e-05 0.0022529619
ADCY8 -7.009731 1.327014 18.77037 1.474400e-05 0.0024950866
AMY2B -3.243416 7.770663 18.08254 2.115327e-05 0.0031189554
APLP1 -3.541290 7.232156 21.57267 3.406721e-06 0.0009439648
ASCL1 -8.407121 4.683148 37.19467 1.069056e-09 0.0000142637
ASTN1 -3.441515 2.731279 15.56582 7.968209e-05 0.0067421913
> dim(DEGs)
[1] 209   5

2 通路富集

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

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

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

go 没有富集到通路，KEGG 富集到了 10 条通路

dotplot(kegg)
barplot(kegg)

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

再进行 GSEA 富集分析

gene_info <- DEGs %>%
  rownames_to_column(var = "SYMBOL") %>%
  filter(SYMBOL %in% gene_list) %>%
  inner_join(., gene.id[,1:2], by = "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.05 再看看结果，这里我就不再演示了，还是继续使用这个结果进行后续分析

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(rownames(DEGs), immune_gene$Symbol))
 [1] "C5"       "CALCA"    "CHGB"     "EDN3"     "IL13RA2"  "KLRC2"    "SSTR2"    "KL"       "SLC22A17"
[10] "ERAP2"

共交叠了 10 个基因，看看这些基因的表达差异情况

> DEGs[gene_list,]
             logFC    logCPM       LR       PValue          FDR
C5       -3.324999  4.020936 23.84367 1.044849e-06 0.0004149276
CALCA    -7.512828  3.813597 19.23037 1.158557e-05 0.0021839537
CHGB     -4.178971 10.264236 19.05373 1.270894e-05 0.0023080671
EDN3     -4.151774  7.076456 17.65441 2.649006e-05 0.0035480333
IL13RA2  -4.808872  2.562736 25.59899 4.202584e-07 0.0002342334
KLRC2    -4.662356  2.262040 18.94337 1.346564e-05 0.0023787247
SSTR2    -2.960156  3.227101 16.19731 5.707501e-05 0.0055727296
KL       -3.917833  4.283146 23.56430 1.208082e-06 0.0004387990
SLC22A17 -3.696153  7.487774 23.80980 1.063399e-06 0.0004149276
ERAP2     2.458962  5.496425 17.21783 3.332940e-05 0.0040042363

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

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

log(exp.count[gene_list,]) %>% 
  as.data.frame() %>%
  rownames_to_column("Gene") %>%
  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()

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

exp.count <- cbind(dataHigh, dataLow)
gene.sur <- exp.count[c("CALCA", "ERAP2"),] %>%
  t() %>% as.data.frame() %>%
  mutate(
    CALCA = if_else(CALCA < median(CALCA), 0, 1),
    ERAP2 = if_else(ERAP2 < median(ERAP2), 0, 1)
  ) %>%
  rownames_to_column("Tumor_Sample_Barcode") %>%
  inner_join(data.surv, by="Tumor_Sample_Barcode")

km.CALCA <- survfit(Surv(time, status) ~ CALCA, data = gene.sur)
km.ERAP2 <- survfit(Surv(time, status) ~ ERAP2, data = gene.sur)

绘制生存曲线

ggsurvplot(
  km.CALCA, 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 值并不显著，可能由于样本量太少了，没有显著性，另外一个就不好了，结果不好。

也可以对这 10 个基因做一个 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

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

先获取样本的 FPKM 值，这个我们之前也是讲过的

query <- GDCquery(
  project = "TCGA-PAAD",
  data.category = "Transcriptome Profiling",
  data.type = "Gene Expression Quantification",
  workflow.type = "HTSeq - FPKM-UQ",
  barcode = colnames(paad.exp)
)
GDCdownload(query)
data <- GDCprepare(query)

exp <- assay(data) %>% as.data.frame() %>%
  rownames_to_column(var = "Ensembl_ID") %>% {
    # 如果一个 Ensembl ID 映射到多个 gene symbols，要将它删除
    unimap <- mapIds(
      org.Hs.eg.db, keys = .$Ensembl_ID, keytype = "ENSEMBL", 
      column = "SYMBOL", multiVals = "filter")
    filter(., Ensembl_ID %in% names(unimap))
  } %>%
  # 将 Ensembl ID 对应到基因 symbol
  mutate(Symbol = AnnotationDbi::select(
    org.Hs.eg.db, keys = .$Ensembl_ID, keytype = "ENSEMBL",
    columns = "SYMBOL")$SYMBOL, .before = Ensembl_ID) %>% 
  # 删除 Ensembl_ID 列和未配对到 symbol 的 NA 行
  dplyr::select(!Ensembl_ID) %>%
  filter(!is.na(Symbol)) %>%
  # 对基因进行分组取平均值
  group_by(Symbol) %>%
  summarise_all(mean) %>% 
  column_to_rownames(var = "Symbol") %>%
  `names<-`(substr(colnames(.), 1, 12))

转换为 TPM

exp.tpm <- exp(log(exp.fpkm) - log(sum(exp.fpkm)) + log(1e6))

2. TIMER

我们先进行 TIMER 分析

immune.timer <- deconvolute(
  exp.tpm, "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，再选择相应的癌型和基因，最后点击提交，就可以获取到最下方的分组箱线图

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

ERAP2 基因的结果如下

还可以对这 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 软件进行免疫浸润分析

write.table(exp.tpm, file = "~/Downloads/TCGA/PAAD_TPM.txt", sep = "\t")

source("~/Documents/WorkSpace/RStudio/Cibersort.R")

immune_cibersort <- CIBERSORT(
  sig_matrix = "~/Documents/WorkSpace/RStudio/LM22.txt", 
  mixture_file = "~/Downloads/TCGA/PAAD_TPM.txt",
  perm = 10, QN = F
)

可视化细胞占比

immune_cibersort[, 1:22] %>% as.data.frame() %>%
  rownames_to_column("sample") %>%
  pivot_longer(cols = -sample, names_to = "cell_type", 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[, 1:22] %>% as.data.frame() %>%
  mutate(label = if_else(rownames(immune_cibersort) %in% sample.high, 1, 0))
  
immune.cate <- immune.tmb %>% as.data.frame() %>%
  rownames_to_column("sample") %>%
  pivot_longer(cols = c(-sample, -label), names_to = "cell_type", values_to = "percent")

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

ggplot(immune.cate, 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_discrete(labels = c("Low", "High"), guide_legend(title = "TMB")) +
  theme(axis.text.x = element_text(angle=60, vjust=1, hjust=1))

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

完整的代码见：https://github.com/dxsbiocc/learn/blob/main/R/TCGA/TCGA_TMB_Immune.R

没怎么整理，将就着看看就行