Enrichment Analysis

本文包括过代表分析（ORA）和基因集富集分析（GSEA），均以最常用的 KEGG 为例。参考clusterProfiler说明文档。

过代表分析（ORA）

ORA 需输入一个基因 ID 向量。KEGG 支持的 ID 为 Entrez ID，如果使用蛋白质谱数据则需转换 ID。

library(clusterProfiler)
library(dplyr)

## choose organism ----
Organism <- 'Human' # 'Human' or 'Mouse'
ifelse(Organism == 'Human', library(org.Hs.eg.db), library(org.Mm.eg.db))

## prepare data ----
data(geneList, package="DOSE")
gene <- names(geneList)[abs(geneList) > 2]

## KEGG enrich ----
kk <- enrichKEGG(
  gene = gene,
  organism = ifelse(Organism == 'Human', 'hsa', 'mmu'), # human: hsa, mouse: mmu
  pvalueCutoff = 0.05,
  # qvalueCutoff = 2,
)
kk_result_p <- subset(kk@result, pvalue < 0.05)
View(kk_result_p)
# write.csv(kk_result_p, file = '_.csv')
# kk_result_p <- read.csv(file = '_.csv', row.names = 1)

# delete suffix ' - Mus musculus (house mouse)' when organism is mouse
if(Organism == 'Mouse') {
  kk_result_p$Description <- kk_result_p$Description |> stringr::str_sub(end = -30L)
}

## remove specific terms ----
## remove disease-related terms
tmp <- kk_result_p[stringr::str_detect(kk_result_p$category, 'disease', negate = T), ]
## remove specific terms
tmp <- tmp |> filter(Description != 'Alcoholism')

## bar plot ----
ggplot(tmp,
       aes(
         x = reorder(Descroption, -pvalue),
         y = -log10(pvalue),
         fill = FoldEnrichment
       )) +
  geom_bar(stat = "identity", width = 0.8) +
  scale_fill_gradient(low = "#95cefa",
                      high = "#234a6d",
                      name = 'FoldEnrichment') +
  coord_flip() + theme_bw() +
  theme(panel.grid = element_blank())

if(F) {
  ggsave(
    filename = '_.pdf',
    width = 200,
    height = 100,
    units = 'mm',
    dpi = 300
  )
}

使用setReadable()函数将富集结果的GeneID列转换为Gene Symbol.

kk <- setReadable(kk, OrgDb = org.Hs.eg.db, keyType="ENTREZID")
## The geneID column is translated to symbol
head(kk, 3)

基因集富集分析（GSEA）

GSEA 需要一个值为 log2FC，names 属性为 gene id 的向量，而且 log2FC 为降序排列。

library(clusterProfiler)
library(dplyr)
library(enrichplot)

# GSEA ----
## prepare gene list ----
data(geneList, package="DOSE")
head(geneList)

## choose organism ----
Organism <- 'Mouse' # 'Human' or 'Mouse'
ifelse(Organism == 'Human', library(org.Hs.eg.db), library(org.Mm.eg.db))

## GSEA KEGG ----
res_KEGG <- gseKEGG(
  geneList,
  organism = ifelse(Organism == 'Human', 'hsa', 'mmu'),
  pvalueCutoff = 0.05,
  pAdjustMethod = "BH"
)

GSEA_result <- subset(res_KEGG@result, pvalue < 0.05)
View(GSEA_result)
# readr::write_csv(GSEA_result, file = '_.csv')

## ridge plot ----
ridgeplot(
  res_KEGG,
  showCategory = 20,
  fill = "p.adjust",
  #填充色 "pvalue", "p.adjust", "qvalue"
  core_enrichment = TRUE,
  #是否只使用 core_enriched gene
  label_format = 30,
  orderBy = "NES",
  decreasing = FALSE
) +
  theme(axis.text.y = element_text(size = 8))

## line chart & hits plot ----
# select a specific term
n <- 1
gsea_description <- res_KEGG@result$Description[n]

if(Organism == 'Mouse') {
  gsea_description <- gsea_description |> stringr::str_sub(end = -30L)
}

gseaplot2(
  res_KEGG,
  geneSetID = n,
  title = gsea_description,
  subplots = 1:2, # choose subplots, including 1 line chart, 2 hits plot and 3 ?
  base_size = 10
)

if(F) {
  ggsave(
    filename = paste0('_', gsea_description, '.pdf'),
    width = 100,
    height = 80,
    units = 'mm',
    dpi = 300
  )
}

通用富集分析

即自定义基因集的富集分析。上述 ORA 和 GSEA 均使用了 KEGG 数据库中的基因集。如果不使用这些现成的基因集而使用自定义的基因集也可以进行富集分析，使用的函数分别是enricher()和GSEA()。以GSEA()为例。

首先需要一个值为 log2FC，names 属性为 gene id 的向量，而且 log2FC 为降序排列。此处仍以DOSE包自带的数据集作为示例。

data(geneList, package="DOSE")
head(geneList)

其次需要确定基因集。制作一个两列数据框my_gene_sets，一列为 term，即自定义的词条名，如CAF_activation，一列为 gene，即该词条内包含的基因，每行写一个基因。该数据框可以包含多种 term。

注意每个 term 一般不少 10 个基因，不多于 500 个基因。

CAF_activation <- c('ACTA2', 'FAP', 'VIM', 'PDGFRB', 'TAGLN', 'COL1A1', 'COL1A2', 'COL3A1', 'FN1', 'POSTN', 'THY1', 'MMP2', 'MMP9', 'TGFBI', 'CCN2')
ECF_remodeling <- c('COL1A1', 'COL1A2', 'COL3A1', 'COL5A1', 'COL6A1', 'FN1', 'POSTN', 'SPARC', 'TNC', 'VCAN', 'MMP2', 'MMP9', 'LOX', 'LOXL2', 'TGFBI')
tumor_stemness <- c('CD44', 'PROM1', 'ALDH1A1', 'LGR5', 'SOX2', 'NANOG', 'POU5F1', 'KLF4', 'MYC', 'BMI1', 'EPCAM')
immunosuppression <- c('TGFB1', 'IL6', 'IL10', 'CCL2', 'CCL5', 'CXCL8', 'CSF1', 'VEGFA', 'CD274', 'PDCD1LG2', 'ARG1', 'MRC1', 'STAT3', 'NFKB1', 'PTGS2')

# organize these vectors into a tibble
my_gene_sets <- tibble::lst(CAF_activation,
                            ECF_remodeling,
                            tumor_stemness,
                            immunosuppression) |> 
  tibble::enframe(name = 'term', value = 'symbol') |> 
  tidyr::unnest_longer(symbol)

# translate from gene symbol to entrez id
my_gene_ids <- bitr(my_gene_sets$gene, fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db")

my_gene_sets <- my_gene_sets |>
  inner_join(my_gene_ids, by = c("symbol" = "SYMBOL")) |> 
  dplyr::select(term, ENTREZID)
# note: make sure all gene symbols can be translated into entrez ids

# make sure all genes are covered in the gene list
my_gene_sets$symbol %in% names(geneList)

使用GSEA()进行富集。pvalueCutoff设为 1 即显示所有富集结果，设为 0.05 则仅显示统计学差异显著的富集结果。

# perform GSEA()
custom_gsea <- GSEA(
  geneList      = gene_list,
  TERM2GENE     = my_gene_sets,
  minGSSize     = 10,             # lower limit of genes contained in a gene set
  maxGSSize     = 500,            # upper limit
  pvalueCutoff  = 1,              # p value cutoff
  pAdjustMethod = "BH",
  verbose       = FALSE
)

custom_gsea_res <- as.data.frame(custom_gsea)
View(custom_gsea_res)
# readr::write_csv(custom_gsea_res, file = '~/result/_custom_GSEA.csv')

通用富集分析的绘图同前述类似。

可能遇到的问题

GSEA 分析的数据清洗阶段需要制作一个值为 log2FC，names 属性为 gene id 的向量。因为向量的 names 属性不能相同，这就要求 gene id 也是无重复的。

但是有些情况会遇到重复的基因，例如 TCGA 数据库以 ensembl id 作为基因的唯一识别符，转换成 gene symbol 或者 entrez id 会发生重复。即使是 ensembl id 去除版本号后自身也有重复。

遇到这些情况可以保留 log2FC 绝对值最大的基因，去除其它重复项，然后再做 GSEA。

# remove duplicate entrez ids
gene_diff <- gene_diff |> mutate(abs_lfc = abs(log2FoldChange)) |> 
  dplyr::group_by(ENTREZID) |> 
  slice_max(order_by = abs_lfc, n = 1, with_ties = FALSE)