TCGA 数据分析实战 —— 富集分析

TCGA 数据分析实战 —— 富集分析

前言

通常，在识别完了差异基因之后，都会对差异基因进行功能富集，来获取差异基因参与的潜在生物学功能通路或生物学进程，有助于理解基因之间的作用关系以及发现基因在癌症发生发展过程中发挥的作用。

通路，通常是一些已知的功能相关的基因集合，而我们常说的基因集合，一般是忽略了基因之间互作关系的通路。

最常见的通路富集，是使用 GO 和 KEGG 数据库中预定义的生物学通路。

Gene Ontology (GO)

Gene Ontology（基因本体）定义了用于描述基因功能的类，以及这些类之间的结构关系，主要可以分为三类：

• Molecular Function(MF)：分子功能，基因产物的生物学活性，如催化或结合等
• Cellular Component(CC)：细胞组分，即基因产物发挥作用的地方，如内质网、高尔基体等
• Biological Process(BP)：由分子功能组成的一系列反应过程。

KEGG

KEGG 是系统分析基因功能和基因组信息的数据框，是一个整合了基因组、生物学通路、疾病、药物以及生物化学物质等信息的数据库。

KEGG 通路由一系列经手工绘制而成的通路图构成，每张通路图均包含分子之间相互作用和反应的网络，旨在将基因组中的基因与基因产物（主要是蛋白质）联系起来，记录了细胞中分子之间的相互作用网络以及具体生物所特有的变化形式。

这些通路主要分为 7 大类：

• 新陈代谢（Metabolism）
• 遗传信息处理（Genetic Information Processing）
• 环境信息处理（Environmental Information Processing）
• 细胞过程（Cellular Processes）
• 生物系统（Organismal Systems）
• 人类疾病（Human Diseases）
• 药物开发（drug development）

其他数据库

当然，除了我们最常用的 GO 和 KEGG，还有一些其他数据库定义的基因集，例如：

• Molecular Signatures Database (MSigDb)
• Reactome
• Disease Ontology (DO)
• Disease Gene Network (DisGeNET)

富集分析方法

富集分析方法主要可以分为四类：

过表达分析

通常是检验差异基因是否显著集中在预先定义的基因集。

使用累积超几何或 Fisher 精确检验
$$p=1-\sum _{i=0}^{k-1}\frac{\left(\genfrac{}{}{0ex}{}{M}{i}\right)\left(\genfrac{}{}{0ex}{}{N-M}{n-i}\right)}{\left(\genfrac{}{}{0ex}{}{N}{n}\right)}$$

式中，N 为背景基因的数量，M 为通路中的基因数，n 为兴趣基因的数量，k 为通路中兴趣基因的数量

显著性打分

对所有差异基因进行打分或排序，并评估基因集中的基因的富集分数

• GSEA：对于一个给定的已排序的差异基因列表 L，以及预定义的基因集 S。GSEA 算法的原理是，通过判断 S 中的基因 s 是随机分布还是主要集中在 L 的顶部或底部，来衡量该基因集合 S 对表型差异的贡献
• ssGSEA： 单样本 GSEA 分析
• GSVA

基于通路拓扑

上述两种方法将通路中的基因视为独立的，但是通路中基因之间具有紧密的互作关系，将这些信息考虑到富集分析中，比如，上游基因的改变对通路功能的影响比下游基因更大，所以，可以将基因的连接度以及调控类型信息以权重的方式加入富集分析当中。

• Pathway-Express
• NetGSA
• topologyGSA
• DEGraph
• PathNet 等

基于网络拓扑

而基于网络拓扑的富集方法，则是把基因在网络中的互作关系整合到富集分析当中

• EnrichNet
• NEA
• NOA

分析示例

我们使用上一节三种方法共同识别出的差异基因

common <- intersect(rownames(DEGs.limma), intersect(rownames(DEGs.DESeq2), rownames(DEGs.edgeR)))
gene_list <- DEGs.edgeR[common, 'gene_name']

GO

使用 TCGABiolinks 提供的 TCGAanalyze_EAcomplete 函数来进行 go 富集

system.time(
  ansEA <- TCGAanalyze_EAcomplete(
    TFname="gbm Vs lgg", gene_list
  )
)

TCGAvisualize_EAbarplot(
  tf = rownames(ansEA$ResBP),
  GOBPTab = ansEA$ResBP,
  GOCCTab = ansEA$ResCC,
  GOMFTab = ansEA$ResMF,
  PathTab = ansEA$ResPat,
  nRGTab = gene_list,
  nBar = 10,
  filename = NULL
)

或者使用 clusterProfiler 包进行富集分析，该包提供了两个函数

• enrichGO：过表达富集分析方法
• gseGO：GSEA 富集分析方法

enrichGO

可以先将基因转换为 entrez_gene_id，再进行富集分析。例如

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

# symbol to ID
gene.id <- bitr(
  gene_list, fromType = "SYMBOL",
  toType = "ENTREZID",
  OrgDb = org.Hs.eg.db
) 

这种方式会导致一些基因对应不上，即找不到对应的 ID。

head(gene.id)
#      SYMBOL ENTREZID
# 1  EIF4A1P2   359792
# 2 HSP90AB3P     3327
# 3  GAPDHP65   389849
# 4    SETP14   389168
# 5  HNRNPKP2   389053
# 6  HNRNPKP4   644063

我选择使用 SYMBOL 进行富集分析，需要将 keyType 参数设置为 SYMBOL

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

富集结果

go[1:3, 1:5]
# ONTOLOGY         ID                                  Description GeneRatio   BgRatio
# GO:0099177       BP GO:0099177       regulation of trans-synaptic signaling  157/2147 490/18870
# GO:0050804       BP GO:0050804 modulation of chemical synaptic transmission  156/2147 489/18870
# GO:0050808       BP GO:0050808                         synapse organization  136/2147 483/18870

绘制富集分析结果的通路点图，点的大小表示通路中的基因数量，颜色表示的是显著性

dotplot(go)

条形图

barplot(go)

gseGO

GSEA 富集分析输入的基因列表需要排序，我们可以按照基因的 logFC 值对基因进行排序

gene_info <- DEGs.edgeR %>%
    filter(gene_name %in% gene_list) %>%
    # 必须降序
    arrange(desc(logFC))

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

GSEA 富集分析

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

富集结果

go2[1:3, 1:5]
# ONTOLOGY         ID                          Description setSize enrichmentScore
# GO:0045202       CC GO:0045202                              synapse     385       0.4872536
# GO:0007268       BP GO:0007268       chemical synaptic transmission     234       0.5365726
# GO:0098916       BP GO:0098916 anterograde trans-synaptic signaling     234       0.5365726

绘制第一条 GO item 的富集曲线

gseaplot(go2, geneSetID = "GO:0098916")

曲线表示富集分数的计算过程，根据基因排序，从左至由依次计算。从图中可以看出富集分数都是小于 0，表示基因富集在通路的下部。中间的竖线表示通路中的基因在列表中的位置

绘制第二种类型的富集曲线，添加了中间的基因与表型之间的相关矩阵热图，红色表示与第一个表型正相关，蓝色表示与第二个表型正相关

gseaplot2(go2, 1)

可以同时显示第 1-3 个富集分析结果的富集曲线

gseaplot2(go2, 1:3)

KEGG

KEGG 通路富集也有两种方法

• enrichKEGG
• gseKEGG

这两个函数只能用基因 ID，还是需要转换一下数据

gene_info <- DEGs.edgeR %>%
    filter(gene_name %in% gene_list) %>%
    inner_join(gene.id, by = c("gene_name" = "SYMBOL")) %>%
    # 必须降序
    arrange(desc(logFC))

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

enrichKEGG

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

富集结果

kegg[1:3, 1:5]
#                    category          subcategory       ID           Description GeneRatio
# hsa04727 Organismal Systems       Nervous system hsa04727     GABAergic synapse   41/1020
# hsa04724 Organismal Systems       Nervous system hsa04724 Glutamatergic synapse   47/1020
# hsa05032     Human Diseases Substance dependence hsa05032    Morphine addiction   41/1020

绘制点图

dotplot(kegg)

gseKEGG

kegg2 <- gseKEGG(
  geneList = geneList,
  organism     = 'hsa',
  minGSSize    = 120,
  pvalueCutoff = 0.05,
  verbose      = FALSE
)
kegg2[1, 1:5]
#                ID        Description setSize enrichmentScore      NES
# hsa01100 hsa01100 Metabolic pathways     189       0.2284664 1.859625
gseaplot(kegg2, geneSetID = "hsa01100")

gseaplot2(kegg2, geneSetID = "hsa01100")

结果可视化

通络网络结构

通过展示通路中基因之间的互作网络结构，可以看出基因在网络中的作用关系，以及潜在的生物学功能

ego <- setReadable(kegg, 'org.Hs.eg.db', 'ENTREZID')
p1 <- cnetplot(ego, showCategory = 2, foldChange = geneList)
# 设置分类的大小，可以是 pvalue 或 geneNum
p2 <-
  cnetplot(
    ego,
    showCategory = 2,
    categorySize = "geneNum",
    foldChange = geneList
  )
# 设置圆形布局
p3 <-
  cnetplot(
    ego,
    showCategory = 3,
    foldChange = geneList,
    circular = TRUE,
    colorEdge = TRUE
  )
# 合并图形
cowplot::plot_grid(
  p1, p2, p3, ncol = 3,
  labels = LETTERS[1:3],
  rel_widths = c(.8, .8, 1.2)
)

热图

使用热图的方式展示通路中基因的表达模式

heatplot(ego, foldChange=geneList)

Enrichment Map

Enrichment Map 是根据通路之间是否有基因交叠来确定通路间是否存在互作边

edo <- pairwise_termsim(kegg)
emapplot(edo, layout="kk") 

对通路关系网络进行聚类展示

emapplot_cluster(edo, node_scale=1.5, layout="kk") 

UpSet plot

使用 upsetplot 函数来绘制通路的 upset 图

upsetplot(edo)

也可以绘制 GSEA 分析的结果

upsetplot(go2)

山脊图

ridgeplot 函数可以绘制 GSEA 分析结果，可以很容易地看出上调和下调的通路

ridgeplot(go2)