TCGA 数据分析实战 —— GSVA、ssGSEA 和单基因富集分析

前言

前面，我们介绍过了差异基因的功能富集分析，今天，我们对这部分的内容作一些补充

主要介绍一下 GEVA、ssGSEA 和单基因的富集分析

GSVA

我们知道，GSEA 富集分析方法是针对两组样本来进行评估的，也就是说对基因列表的排列方式是根据基因与表型的相关度（例如，FC 值）来计算的，无法对单个样本使用

其富集分数（Enrichment Score，ES）的计算方式为

依次判断基因列表中的基因是否在基因集合中，如果在基因集合中，则 ES 加上该基因与表型的相关度，如果不是集合中的基因，则减去对应的值，最后可以计算出一个最大 ES。

Gene Set Variation Analysis（GSVA）与 GSEA 的原理类似，只是计算每个基因集合在每个样本中的 enrichment statistic（ES，或 GSVA score），其算法流程如下

不同于 GSEA 之处在于，对于不同的数据类型（只支持 log 表达值或原始的 read counts 值），假设了不同的累积密度函数（cumulative density function，CDF）

1. 芯片数据：正态分布密度函数
2. RNA-seq 数据：泊松分布密度函数

而且，GSVA 是为每个样本的每个基因计算对应的 CDF 值，然后根据该值对基因进行排序，这样，每个样本都有一个从大到小排序的基因列表

对于某一基因集合，计算其在每个样本中的 ES 值，也就是评估基因集合在基因列表中的富集情况。

例如，我们有一个排序后的样本

基因集合包含：B、E、H，我们可以绘制这样一张 K-S 分布图

x 轴为排序后的基因顺序，依照这一顺序，如果基因在集合内，则累积和会加上该基因的值（与基因的顺序有关，排名越靠前值越大），否则累积和不变。

将基因列表分为基因集内核基因集外两个集合，就可以绘制两个分布（红色和绿色曲线），分别计算两个分布之间的最大间距，以基因集内的分布值更大视为正间距（即红色曲线更高），两个最大间距之和即为该基因集的 ES 值

这样，就把基因水平的表达矩阵转换成了基因集水平的评分矩阵，可以使用差异表达基因识别算法，寻找显著差异的基因集，从而达到类似功能富集的作用

GSVA 分析

先获取基因表达矩阵，我们使用 TCGA 肺腺癌和肺鳞癌各 10 个样本的 read counts 数据

library(TCGAbiolinks)

get_prep <- function(cancer) {
  query <- GDCquery(
    project = cancer,
    data.category = "Transcriptome Profiling", 
    data.type = "Gene Expression Quantification", 
    workflow.type = "STAR - Counts",
    sample.type = c("Primary Tumor"),
  )
  # 选择 20 个样本
  query$results[[1]] <-  query$results[[1]][1:20,]
  GDCdownload(query)
  # 获取 read count
  exp.count <- GDCprepare(
    query,
    summarizedExperiment = TRUE,
  )
  return(exp.count)
}

luad.count <- get_prep("TCGA-LUAD")
lusc.count <- get_prep("TCGA-LUSC")

dataPrep_luad <- TCGAanalyze_Preprocessing(
  object = luad.count,
  cor.cut = 0.6,
  datatype = "unstranded"
)

dataPrep_lusc <- TCGAanalyze_Preprocessing(
  object = lusc.count,
  cor.cut = 0.6,
  datatype = "unstranded"
)
# 合并数据并使用 gcContent 方法进行标准化
dataNorm <- TCGAanalyze_Normalization(
  tabDF = cbind(dataPrep_luad, dataPrep_lusc),
  geneInfo = TCGAbiolinks::geneInfoHT,
  method = "gcContent"
)
# 分位数过滤
dataFilt <- TCGAanalyze_Filtering(
  tabDF = dataNorm,
  method = "quantile",
  qnt.cut =  0.25
)

gene.counts <- as.data.frame(dataFilt) %>%
    tibble::rownames_to_column("gene_id") %>%
    inner_join(dplyr::select(TCGAbiolinks::get.GRCh.bioMart("hg38"), gene_id, gene_name) 
               %>% mutate(gene_id = str_remove(gene_id, "\\.\\d+"))) %>%
    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))
# 将数据拆分
luad.exp <- dplyr::select(gene.counts, substr(luad.count$barcode, 1, 12))
lusc.exp <- dplyr::select(gene.counts, substr(lusc.count$barcode, 1, 12))

我们使用 GSVA 包提供的 gsva 函数来将基因表达矩阵转换为基因集分数矩阵

library(GSVA)
library(GSEABase)

# 读取从 GSEA 官网下载的通路数据
c2gmt <- getGmt("~/Downloads/data/pathway/c2.cp.v7.2.symbols.gmt")
# 删选出常用的这三个数据库中的通路
gene.set <- c2gmt[grep("^KEGG|REACTOME|BIOCARTA", names(c2gmt)),]
# gsva 分析，read counts 使用泊松分布，通路至少包含 10 个基因
gs.exp <- gsva(as.matrix(gene.counts), gene.set, kcdf = "Poisson", min.sz = 10)

虽然 GSVAdata 包提供了通路数据 c2BroadSets 是基因 ID，但我们的基因表达数据的行是基因 Symbol，所以通路信息也必须是 Symbol 格式，要进行格式转换，比较麻烦

所以我们使用 GSEABase 包提供的 getGmt 函数来读取从 GSEA 官网下载的 C2 通路信息

得到结果如下，共包含 1642 条通路

然后，使用差异基因识别方法

差异分析

我们使用 limma 分析差异通路

DEA.gs <- TCGAanalyze_DEA(
  mat1 = gs.exp[, colnames(luad.exp)],
  mat2 = gs.exp[, colnames(lusc.exp)],
  metadata = FALSE,
  pipeline = "limma",
  Cond1type = "LUAD",
  Cond2type = "LUSC",
  fdr.cut = 0.05,
  logFC.cut = 0.5,
)

通过设置 FDR = 0.05，logFC = 0.5 共筛选出 13 条差异通路

查看通路的火山图

我们可以一起看下基因的火山图

DEA.gene <- TCGAanalyze_DEA(
  mat1 = luad.exp,
  mat2 = lusc.exp,
  metadata = FALSE,
  pipeline = "limma",
  Cond1type = "LUAD",
  Cond2type = "LUSC",
  fdr.cut = 0.01,
  logFC.cut = 1
)

总共识别出 3157 个差异表达基因

ssGSEA

single sample Gene Set Enrichment Analysis (ssGSEA) 是针对单个样本进行 GSEA 分析，其基因列表的排序方式和 ES 的计算方式都是依赖于样本中基因的表达值，而不再是依赖基因与表型的相关度

使用方式也很简单，只要在 gsva 函数中指定 method = "ssgsea"，例如

res.ssgsea <- gsva(as.matrix(gene.counts), gene.set, method = "ssgsea", kcdf = "Poisson", min.sz = 10)

也可以进行差异分析

DEA.ssgsea <- TCGAanalyze_DEA(
  mat1 = res.ssgsea[, colnames(luad.exp)],
  mat2 = res.ssgsea[, colnames(lusc.exp)],
  metadata = FALSE,
  pipeline = "limma",
  Cond1type = "LUAD",
  Cond2type = "LUSC",
  fdr.cut = 0.05,
  logFC.cut = 0.1,
)

或者绘制热图

annotation_col <- data.frame(sample = rep(1:2, each = 20))
rownames(annotation_col) <- c(colnames(luad.exp), colnames(lusc.exp))
pheatmap::pheatmap(
    res.ssgsea[rownames(DEA.ssgsea), rownames(annotation_col)],
    show_colnames = F,
    # 不展示行名
    cluster_rows = F,
    # 不对行聚类
    cluster_cols = F,
    # 不对列聚类
    annotation_col = annotation_col,
    # 加注释
    cellwidth = 5,
    cellheight = 5,
    # 设置单元格的宽度和高度
    fontsize = 5
)

单基因富集分析

单基因富集分析并不是说拿单个基因来进行富集分析，单个基因怎么能进行富集分析呢？一个基因根本没法进行统计检验。

其实，这里说的单基因并不是拿单个基因来富集，而是基于单个基因来进行富集分析，这个“基于”，就是以单个基因为基础，向外扩展，抓取与其相关的基因，然后用这些相关的基因来进行功能富集

所以，要理解这个单基因富集分析的意思，这样一说就已经很明了了。针对单个基因我们可以做什么？

主要有两种做法：

1. 定性分组：我们可以根据给定基因的表达值对样本进行分组，然后识别在两组样本之间差异表达的基因，最后用这些差异表达基因来进行功能富集

2. 定量相关：通过计算其他基因与目标基因表达之间的相关性，将具有显著相关的基因作为一个集合，也可以进行富集分析

定性分组

我们以 CCDC134 基因为例，以该基因表达值的中位值来对样本进行分组

gene <- "CCDC134"
gene.exp <- gene.counts[gene,] %>%
    t() %>% as.data.frame() %>%
    mutate(label = ifelse(!!sym(gene) < median(!!sym(gene)), 0, 1))

group.low <- gene.counts[,gene.exp$label == 0]
group.high <- gene.counts[,gene.exp$label == 1]

识别两组样本之间的差异表达基因

DEGs <- TCGAanalyze_DEA(
  mat1 = group.low,
  mat2 = group.high,
  metadata = FALSE,
  pipeline = "limma",
  Cond1type = "CCDC134_Low",
  Cond2type = "CCDC134_High",
  fdr.cut = 0.01,
  logFC.cut = 1,
)

共识别出 710 个差异表达基因

富集分析

对差异基因进行富集分析

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

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

GSEA 富集

gene_info <- DEGs %>%
  rownames_to_column(var = "SYMBOL") %>%
  inner_join(., gene.id[,1:2], by = "SYMBOL") %>%
  # 必须降序
  arrange(desc(logFC))

# 构造输入数据格式
geneList <- gene_info$logFC
names(geneList) <- as.character(gene_info$ENTREZID)

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

定量相关

我们计算 CCDC134 基因与其他基因的相关性，筛选出相关系数的绝对值大于 0.6 的基因，共筛选出 37 个基因。

cor.genes <- cor(t(gene.counts), as.numeric(gene.exp$CCDC134)) %>%
    as.data.frame() %>%
    tibble::rownames_to_column("gene") %>%
    rename_with(.cols = 2, ~"corr") %>%
    arrange(-corr) %>%
    dplyr::filter(gene != "CCDC134")

high.cor <- cor.genes %>%
    dplyr::filter(abs(corr) > 0.6)

富集分析

使用高度相关的基因进行富集分析

go <- enrichGO(
    gene = high.cor$gene,
    OrgDb = org.Hs.eg.db,
    keyType = "SYMBOL",
    ont = "ALL",
    pAdjustMethod = "BH",
    qvalueCutoff = 0.05,
    readable = T
)
dotplot(go)

根据相关性的大小排序，然后进行 GSEA 分析

# 构造输入数据格式
geneList <- cor.genes$corr
names(geneList) <- as.character(cor.genes$gene)

go2 <- gseGO(
    geneList     = geneList,
    OrgDb        = org.Hs.eg.db,
    keyType      = "SYMBOL",
    ont          = "ALL",
    minGSSize    = 10,
    maxGSSize    = 500,
    pvalueCutoff = 0.1,
    verbose      = FALSE
)
dotplot(go2)