TCGA 数据分析实战 —— 富集分析
2021-07-17 本文已影响0人
名本无名
前言
通常，在识别完了差异基因之后，都会对差异基因进行功能富集，来获取差异基因参与的潜在生物学功能通路或生物学进程，有助于理解基因之间的作用关系以及发现基因在癌症发生发展过程中发挥的作用。
通路，通常是一些已知的功能相关的基因集合，而我们常说的基因集合，一般是忽略了基因之间互作关系的通路。
最常见的通路富集，是使用 GO 和 KEGG 数据库中预定义的生物学通路。
1. Gene Ontology (GO)
Gene Ontology（基因本体）定义了用于描述基因功能的类，以及这些类之间的结构关系，主要可以分为三类：
-
Molecular Function(MF)：分子功能，基因产物的生物学活性，如催化或结合等 -
Cellular Component(CC)：细胞组分，即基因产物发挥作用的地方，如内质网、高尔基体等 -
Biological Process(BP)：由分子功能组成的一系列反应过程。
2. KEGG
KEGG 是系统分析基因功能和基因组信息的数据框，是一个整合了基因组、生物学通路、疾病、药物以及生物化学物质等信息的数据库。
KEGG 通路由一系列经手工绘制而成的通路图构成，每张通路图均包含分子之间相互作用和反应的网络，旨在将基因组中的基因与基因产物（主要是蛋白质）联系起来，记录了细胞中分子之间的相互作用网络以及具体生物所特有的变化形式。
这些通路主要分为 7 大类：
- 新陈代谢（
Metabolism） - 遗传信息处理（
Genetic Information Processing） - 环境信息处理（
Environmental Information Processing） - 细胞过程（
Cellular Processes） - 生物系统（
Organismal Systems） - 人类疾病（
Human Diseases） - 药物开发（
Drug development）
3. 其他数据库
当然，除了我们最常用的 GO 和 KEGG，还有一些其他数据库定义的基因集，例如：
-
Molecular Signatures Database(MSigDb) Reactome-
Disease Ontology(DO) -
Disease Gene Network(DisGeNET)
富集分析方法
富集分析方法主要可以分为四类：
- 过表达分析：通常是检验差异基因是否显著集中在预先定义的基因集
- 累积超几何或
Fisher精确检验
式中，N 为背景基因的数量，M 为通路中的基因数，n 为兴趣基因的数量，k 为通路中兴趣基因的数量
- 显著性打分：对所有差异基因进行打分或排序，并评估基因集中的基因的富集分数
-
GSEA
对于一个给定的已排序的差异基因列表L，以及预定义的基因集S。GSEA算法的原理是，通过判断S中的基因s是随机分布还是主要集中在L的顶部或底部，来衡量该基因集合S对表型差异的贡献 -
ssGSEA
单样本GSEA分析 GSVA
- 基于通路拓扑：上述两种方法将通路中的基因视为独立的，但是通路中基因之间具有紧密的互作关系，将这些信息考虑到富集分析中，比如，上游基因的改变对通路功能的影响比下游基因更大，所以，可以将基因的连接度以及调控类型信息以权重的方式加入富集分析当中。
Pathway-ExpressNetGSAtopologyGSADEGraph-
PathNet等
- 基于网络拓扑：而基于网络拓扑的富集方法，则是把基因在网络中的互作关系整合到富集分析当中
EnrichNetNEANOA
分析示例
我们使用上一节三种方法共同识别出的 2749 个差异基因
> load("~/Downloads/gbm_lgg_deg.rda")
> gene_list <- rownames(DEGs.exp)
> DEGs.exp[1:6, 1:3]
TCGA-HT-A614-01A-11R-A29R-07 TCGA-HT-8104-01A-11R-2404-07 TCGA-TQ-A7RG-01A-11R-A33Z-07
A1BG 225 238 148
A2BP1 610 471 181
A2LD1 117 190 38
AATK 9071 3330 2785
ABAT 16024 17488 40945
ABCC3 11 192 9
1. 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 = "~/Downloads/go_enrichment.pdf"
)
或者使用 clusterProfiler 包进行富集分析，该包提供了两个函数
-
enrichGO：过表达富集分析方法 -
gseGO：GSEA富集分析方法
1.1 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
)
> head(gene.id)
SYMBOL ENTREZID
1 A1BG 1
4 AATK 9625
5 ABAT 18
6 ABCC3 8714
7 ABCC8 6833
8 ABCG4 64137
进行富集分析
go <- enrichGO(
gene = gene.id,
OrgDb = org.Hs.eg.db,
ont = "ALL",
pAdjustMethod = "BH",
qvalueCutoff = 0.05,
readable = T
)
> head(go)
ONTOLOGY ID Description GeneRatio BgRatio pvalue
GO:0050804 BP GO:0050804 modulation of chemical synaptic transmission 150/2256 454/18866 4.314421e-33
GO:0099177 BP GO:0099177 regulation of trans-synaptic signaling 150/2256 455/18866 5.706736e-33
GO:0050808 BP GO:0050808 synapse organization 136/2256 433/18866 1.356042e-27
GO:0042391 BP GO:0042391 regulation of membrane potential 134/2256 443/18866 1.907415e-25
GO:0050890 BP GO:0050890 cognition 97/2256 302/18866 8.988934e-21
GO:0099003 BP GO:0099003 vesicle-mediated transport in synapse 79/2256 220/18866 1.893468e-20
p.adjust qvalue
GO:0050804 1.716872e-29 1.281913e-29
GO:0099177 1.716872e-29 1.281913e-29
GO:0050808 2.719769e-24 2.030733e-24
GO:0042391 2.869229e-22 2.142329e-22
GO:0050890 1.081728e-17 8.076794e-18
GO:0099003 1.898833e-17 1.417776e-17
绘制富集分析结果的通路点图，点的大小表示通路中的基因数量，颜色表示的是显著性
dotplot(go)
条形图
barplot(go)
1.2 gseGO
GSEA 富集分析输入的基因列表需要排序，我们可以按照基因的 logFC 值对基因进行排序
gene_info <- DEGs.edgeR %>%
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) <- as.character(gene_info$ENTREZID)
GSEA 富集分析
go2 <- gseGO(
geneList = geneList,
OrgDb = org.Hs.eg.db,
ont = "ALL",
minGSSize = 100,
maxGSSize = 500,
pvalueCutoff = 0.05,
verbose = FALSE
)
> head(go2)
ONTOLOGY ID Description setSize enrichmentScore NES pvalue p.adjust
GO:0000228 CC GO:0000228 nuclear chromosome 188 -0.3746191 -3.036128 1e-10 5.903226e-10
GO:0000278 BP GO:0000278 mitotic cell cycle 198 -0.4102983 -3.349095 1e-10 5.903226e-10
GO:0001501 BP GO:0001501 skeletal system development 115 -0.4414610 -3.234870 1e-10 5.903226e-10
GO:0001525 BP GO:0001525 angiogenesis 117 -0.4349952 -3.197417 1e-10 5.903226e-10
GO:0001568 BP GO:0001568 blood vessel development 158 -0.4020008 -3.119252 1e-10 5.903226e-10
GO:0001775 BP GO:0001775 cell activation 263 -0.3259590 -2.814883 1e-10 5.903226e-10
qvalues rank leading_edge
GO:0000228 2.286361e-10 672 tags=53%, list=27%, signal=41%
GO:0000278 2.286361e-10 930 tags=70%, list=38%, signal=47%
GO:0001501 2.286361e-10 294 tags=36%, list=12%, signal=33%
GO:0001525 2.286361e-10 847 tags=67%, list=35%, signal=46%
GO:0001568 2.286361e-10 571 tags=48%, list=23%, signal=39%
GO:0001775 2.286361e-10 1220 tags=79%, list=50%, signal=44%
绘制第一条 GO item 的富集曲线
gseaplot(go2, geneSetID = "GO:0000228")
曲线表示富集分数的计算过程，根据基因排序，从左至右依次计算。从图中可以看出富集分数都是小于 0，表示基因富集在通路的下部。中间的竖线表示通路中的基因在列表中的位置
绘制第二种类型的富集曲线，添加了中间的基因与表型之间的相关矩阵热图，红色表示与第一个表型正相关，蓝色表示与第二个表型正相关
gseaplot2(go2, 1)
可以同时显示第 1-3 个富集分析结果的富集曲线
gseaplot2(go2, 1:3)
2. KEGG
KEGG 通路富集也有两种方法
enrichKEGGgseKEGG
2.1 enrichKEGG
kegg <- enrichKEGG(
gene = gene.id$ENTREZID,
organism = "hsa",
pvalueCutoff = 0.05
)
> head(kegg)
ID Description GeneRatio BgRatio pvalue p.adjust
hsa05033 hsa05033 Nicotine addiction 24/1087 40/8108 6.620414e-12 2.138394e-09
hsa04724 hsa04724 Glutamatergic synapse 43/1087 114/8108 4.747252e-11 7.666812e-09
hsa04727 hsa04727 GABAergic synapse 36/1087 89/8108 1.824548e-10 1.964430e-08
hsa04911 hsa04911 Insulin secretion 32/1087 86/8108 2.129231e-08 1.719354e-06
hsa04080 hsa04080 Neuroactive ligand-receptor interaction 81/1087 341/8108 8.656665e-08 4.646393e-06
hsa04721 hsa04721 Synaptic vesicle cycle 29/1087 78/8108 9.964211e-08 4.646393e-06
qvalue
hsa05033 1.561024e-09
hsa04724 5.596760e-09
hsa04727 1.434031e-08
hsa04911 1.255125e-06
hsa04080 3.391859e-06
hsa04721 3.391859e-06
dotplot(kegg)
2.2 gseKEGG
kegg2 <- gseKEGG(
geneList = geneList,
organism = 'hsa',
minGSSize = 120,
pvalueCutoff = 0.05,
verbose = FALSE
)
gseaplot(kegg2, geneSetID = "hsa05033")
gseaplot2(kegg2, geneSetID = "hsa05033")
3 结果可视化
3.1 通路网络结构
通过展示通路中基因之间的互作网络结构，可以看出基因在网络中的作用关系，以及潜在的生物学功能
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)
)
3.2 热图
使用热图的方式展示通路中基因的表达模式
heatplot(ego, foldChange=geneList)
3.3 Enrichment Map
Enrichment Map 是根据通路之间是否有基因交叠来确定通路间是否存在互作边
edo <- pairwise_termsim(kegg)
emapplot(edo, layout="kk")
对通路关系网络进行聚类展示
emapplot_cluster(edo, node_scale=1.5, layout="kk")
3.4 UpSet plot
使用 upsetplot 函数来绘制通路的 upset 图
upsetplot(edo)
也可以绘制 GSEA 分析的结果
upsetplot(go2)
3.5 山脊图
ridgeplot 函数可以绘制 GSEA 分析结果，可以很容易地看出上调和下调的通路
ridgeplot(go2)
先介绍这些吧，相关的内容如 GSVA、ssGSEA、单基因富集分析等后续再做介绍
