RNA 10. SCI 文章中基因表达富集之 KEGG 注释

已于 2022-03-16 07:46:58 修改
阅读量3.9k 收藏 32
点赞数
分类专栏： RNA数据分析 文章标签： r语言 数据挖掘 人工智能
于 2022-03-02 22:58:48 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/weixin_41368414/article/details/123242657
版权

全网最全 KEGG 注释结果绘图，直击 SCI 绘图标注，关注我，您最好的选择！

前言

1. KEGG 原理

KEGG(Kyoto Encyclopedia of Genes and Genomes)数据库是系统地分析基因功能、链接基因组信息和功能信息的数据库，包括代谢通路（pathway）数据库、分层分类数据库、基因数据库、基因组数据库等。KEGG的pathway数据库是应用最广泛的代谢通路公共数据库。

富集的含义：
这里pathway富集的含义与GO富集的含义相同，也是表示差异基因中注释到某个代谢通路的基因数目在所有差异基因中的比例显著大于背景基因中注释到某个代谢通路的基因数目在所有背景基因中的比例。因此，做pathway富集分析，也是涉及到前景基因和背景基因。前景基因就是你关注的要重点研究的基因集，背景基因就是所有的基因集。
**
富集显著性（P value）的计算：**
计算方法和公式与GO富集分析一样，也是利用超几何检验计算：

图片

其中，N为所有基因中具有Pathway注释的基因数目；n为N中差异表达基因的数目；M为所有基因中注释为某特定Pathway的基因数目；m为注释为某特定Pathway的差异表达基因数目。
计算得到的P value会进一步经过多重检验校正，得到corrected-pvalue（也就是Q value）。通常我们会以Q value≤0.05为阈值，满足此条件的pathway定义为在差异表达基因中显著富集的pathway。

2. 实例解析

1. 数据读取

数据的读取我们仍然使用的是 TCGA-COAD 的数据集，表达数据的读取以及临床信息分组的获得我们上期已经提过，我们使用的是edgeR 软件包计算出来的差异表达结果，提取上调基因 2832 的 ENSEMBL 号，

###########基因列表
DEG=read.table("DEG-resdata.xls",sep="\t",check.names=F,header = T)
geneList<-DEG[DEG$sig=="Up",]$Row.names
table(DEG$sig)


## 
## Down   Up 
## 1296 2832

读取样本分组信息，41个正常组织，478个癌症组织，如下：

######样本分组信息
group<-read.table("DEG-group.xls",sep="\t",check.names=F,header = T)
table(group$Group)


## 
##  NT  TP 
##  41 478

2. KEGG 注释结果

首先我们同样需要安装软件包并加载，这里面主程序就是 clusterProfiler 软件包，如下：

########
if(!require(clusterProfiler)){
  BiocManager::install("clusterProfiler")
}
if(!require(org.Hs.eg.db)){
  BiocManager::install("org.Hs.eg.db")
}
if(!require(DOSE)){
  BiocManager::install("DOSE")
}
if(!require(topGO)){
  BiocManager::install("topGO")
}
if(!require(pathview)){
  BiocManager::install("pathview")
}
if(!require(KEGG.db)){
  BiocManager::install("KEGG.db")
}
library(org.Hs.eg.db)
library(clusterProfiler)
library(DOSE)
library(topGO)
library(pathview)
library(KEGG.db)

根据数据比对我们找到1726 个KEGG数据中有的基因，同时我们获得到一个基因的三种命名方式，即"ENTREZID", “ENSEMBL”, ‘SYMBOL’，如下：

####KEGG 
eg <- bitr(geneList, 
           fromType="ENSEMBL", 
           toType=c("ENTREZID","ENSEMBL",'SYMBOL'),
           OrgDb="org.Hs.eg.db")

head(eg)

##           ENSEMBL ENTREZID  SYMBOL
## 1 ENSG00000062038     1001    CDH3
## 2 ENSG00000175832     2118    ETV4
## 3 ENSG00000167767   144501   KRT80
## 4 ENSG00000164283    11082    ESM1
## 5 ENSG00000120254    25902 MTHFD1L
## 6 ENSG00000129474    84962   AJUBA

对获得的基因进行KEGG注释，这里是人的结直肠癌的数据，故 organism = ‘hsa’， 阈值分别为 pvalueCutoff = 0.05 和 qvalueCutoff = 0.05， 注释的时候我们可以选择参数 keyType 中4种不同的选择，包括"kegg", ‘ncbi-geneid’, ‘ncib-proteinid’ and ‘uniprot’，下面我们看到利用基因列表注释到6个pathway通路，如下：

# Run KEGG enrichment analysis 
kegg <- enrichKEGG(eg$ENTREZID, 
                   organism = 'hsa',  
                   keyType = 'kegg', 
                   pvalueCutoff = 0.05,
                   pAdjustMethod = 'BH', 
                   minGSSize = 1,
                   maxGSSize = 500,
                   qvalueCutoff = 0.05,
                   use_internal_data = FALSE)
head(kegg)

##                ID                             Description GeneRatio
## hsa05034 hsa05034                              Alcoholism    38/384
## hsa05322 hsa05322            Systemic lupus erythematosus    31/384
## hsa04613 hsa04613 Neutrophil extracellular trap formation    33/384
## hsa04310 hsa04310                   Wnt signaling pathway    25/384
## hsa04657 hsa04657                 IL-17 signaling pathway    18/384
## hsa04060 hsa04060  Cytokine-cytokine receptor interaction    34/384
##           BgRatio       pvalue     p.adjust       qvalue
## hsa05034 187/8115 8.845476e-15 2.529806e-12 2.420867e-12
## hsa05322 136/8115 1.113336e-13 1.592070e-11 1.523512e-11
## hsa04613 190/8115 5.170512e-11 4.929222e-09 4.716958e-09
## hsa04310 167/8115 2.536698e-07 1.786110e-05 1.709196e-05
## hsa04657  94/8115 3.122570e-07 1.786110e-05 1.709196e-05
## hsa04060 295/8115 1.126797e-06 5.371067e-05 5.139777e-05
##                                                                                                                                                                                                         geneID
## hsa05034 2906/2904/7054/8367/8356/8364/85235/128312/8360/8366/440689/8348/8332/8335/8340/1813/2786/8331/8354/8346/8344/2792/8338/317772/8343/8336/810/3013/8359/653604/8351/8968/8970/8350/8368/6531/8342/3018
## hsa05322                                   2904/8367/8356/8364/85235/128312/8360/8366/440689/8348/8332/8335/8340/8331/8354/8346/8344/8338/317772/8343/8336/3013/8359/653604/8351/8968/8970/8350/8368/8342/3018
## hsa04613                          5582/8367/8356/8364/85235/128312/8360/8366/440689/8348/8332/8335/8340/366/8331/8354/8346/8344/8338/317772/8343/8336/2244/3013/8359/653604/8351/8968/8970/8350/8368/8342/3018
## hsa04310                                                            7472/4316/85409/7473/147111/54894/8840/7477/51176/85407/27121/8313/8549/8061/6424/5582/7479/7481/59352/22943/11211/89780/11197/7476/343637
## hsa04657                                                                                                        225689/4314/2921/4312/2919/4322/8061/2920/3576/6374/1437/27189/6372/4586/3934/6278/3605/112744
## hsa04060                  3624/51330/9518/3589/2921/2919/8744/51561/268/2920/655/3576/3552/6374/1437/5473/10913/1237/284340/3557/27189/6372/10344/5008/3625/4982/6373/338376/26525/4838/3623/56300/3605/112744
##          Count
## hsa05034    38
## hsa05322    31
## hsa04613    33
## hsa04310    25
## hsa04657    18
## hsa04060    34

3. 结果展示

针对 KEGG 富集结果我们这里给出了多种展示方式，根据自己的需求以及文章的设计，选择适合自己的即可。

绘图过程中我们需要安装并加载几个软件，如下：

library(stringr)
library(cowplot)
library(ggplot2)

绘制气泡图，气泡图解读需要说明一下，GO富集程度通过Gene ratio、Pvalue和富集到此GO term上的基因个数来衡量。

  • 横坐标是Gene ratio，数值越大表示富集程度越大。Count 位于该GO term下的差异表达基因数

  • 纵坐标是富集程度较高的GO term（一般选取富集最显著的20条进行展示，不足20条则全部列出）。

  • Pvalue取值范围[0, 1]，以颜色表示，越红表示Pvalue越小，说明富集越明显。

  • 点的大小表示该term下差异基因的个数，点越大表示基因数越多。

dotplot(kegg,showCategory=10)

图片

由于KEGG 的描述字体重合，故我们利用scale_y_discrete设置，避免字体彼此重合，如下：

dotplot(kegg,showCategory=10)+
  scale_y_discrete(labels=function(x) stringr::str_wrap(x, width=60))

图片

绘制柱状图，如下：

barplot(kegg,showCategory=10)+
  scale_y_discrete(labels=function(x) stringr::str_wrap(x, width=60))


## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.

图片

绘制基因概念的网络图，Pathway 与差异基因关系网络图 Gene-Concept Network，对于基因和富集的Pathway 之间的对应关系进行展示说明：图中灰色的点代表基因，黄色的点代表富集到的Pathway ；如果一个基因位于一个Pathway 下，则将该基因与Pathway 连线；黄色节点的大小对应富集到的基因个数，top2富集到的Pathway ，如下：

cnetplot(kegg,showCategory=3,circular=TRUE,colorEdge=TRUE)

图片

热图可以看到哪些基因富集到哪些 Pathway，如下：

heatplot(kegg,showCategory=6)+
  scale_y_discrete(labels=function(x) stringr::str_wrap(x, width=60))


## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.

图片

软件安装并加载，如下：

if(!require(ggnewscale)){
  install.packages("ggnewscale")
}
library(enrichplot)
library(ggnewscale)

对于富集到的pathways之间的基因重叠关系进行展示，如果两个pathway的差异基因存在重叠，说明这两个节点存在overlap关系，在图中用线条连接起来，每个节点是一个富集到的pathway, 默认画top30个富集到的pathways, 节点大小对应该pathway下富集到的差异基因个数，节点的颜色对应p.adjust的值，从小到大，对应蓝色到红色。用法如下：

ed = enrichDO(eg$ENTREZID, pvalueCutoff=0.05)
met <- pairwise_termsim(ed)
emapplot(met,showCategory = 15)

图片

在pathway通路图上标记富集到的基因，代码如下 会给出一个url链接，示例：KEGG PATHWAY: Alcoholism - Homo sapiens (human)在浏览器中打开会看到如下所示的图片，如下：

browseKEGG(kegg, "hsa05034")

图片

KEGG的富集分析与GO差不多，同样的软件包同样的函数，唯一不同的是最后一个可以看到pathway通路上的基因，非常方便！

关注公众号，每日更新，扫码进群交流不停歇，马上就出视频版，关注我，您最佳的选择！

图片

桓峰基因

生物信息分析，SCI文章撰写及生物信息基础知识学习：R语言学习，perl基础编程，linux系统命令，Python遇见更好的你

46篇原创内容

公众号

References：

  1. Kanehisa M, Furumichi M, Tanabe M, Sato Y, Morishima K. KEGG: new perspectives on genomes, pathways, diseases and drugs. Nucleic Acids Res. 2017;45(D1):D353-D361. doi:10.1093/nar/gkw1092

  2. Yu G, Wang L, Han Y and He Q*. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS: A Journal of Integrative Biology. 2012, 16(5):284-287.

桓峰基因
关注
  • 0
    点赞
  • 踩
  • 32
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
专栏目录
keggannot:KEGG注释Python库
05-01
打包以基于blastp对KEGG数据库注释蛋白质。 关于KEGG版本和支持的注意事项 该软件已经过2013年4月版KEGG数据库的测试。 如果不更改代码，则不能使用任何其他版本的KEGG。 该软件也按原样提供，很可能不会收到任何进一步的更新。 安装 python setup.py install 用法 usage: keggannot_genes2ko [-h] [-v] [-V] [-m] [kegg_dir] [blast_file] Annotate proteins with KEGG BLAST results positional arguments: kegg_dir KEGG data directory (default: None) blast_file Tabular blast results file (default: None)
获取KEGG数据库中的癌症Pathway的基因集
12-12
如何利用在线数据库检索各种癌症Pathway相关的已知的信息基因
ggkegg | KEGG数据可视化利器，就说你会用到吧...
最新发布
刘永鑫的博客——宏基因组公众号
02-09 83
前言一、ggkegg-利用图形语法分析KEGG信息并使其可视化~二、可视化学习圈子是干什么的？三、系统学习可视化四、猜你喜欢前言ggkegg-利用图形语法分析KEGG信息并使其可视化~今天有我们R语言可视化课程的同学，私信我有没有专门可视化KEGG信息数据的工具包，由于不是相关专业的，所以感觉查阅了资料，还别说，真知道了一个专门可视化KEGG信息的拓展工具包，这不！赶紧推荐给大家···「什么是KE...
KEGG（全基因组及代谢途径数据库）.doc
05-18
KEGG（全基因组及代谢途径数据库）.doc
KEGG圈图_R语言绘制SCI图的输入文件及代码.rar
07-30
KEGG圈图_R语言绘制SCI图的输入文件及代码.rar
kegg与go通路数据库介绍功能富集软件介绍.ppt
06-03
kegg与go通路数据库介绍功能富集软件介绍
linux转录组kegg注释,蛋白质组学KEGG注释详解
weixin_42367464的博客
05-16 3933
8种机械键盘轴体对比本人程序员，要买一个写代码的键盘，请问红轴和茶轴怎么选？一、KEGG简介KEGG is a database resource for understanding high-level functions and utilities of the biological system, such as the cell, the organism and the ecosyst...
KEGG注释（生信学习）
m0_70166478的博客
09-21 1625
学习
☞GO和KEGG富集倍数（Fold Enrichment）如何计算 enrich factor qvalue
生信小博士的博客
07-27 6625
比如说，某一类功能基因的总数是100个，其中差异基因有10个，那么只有实验只是影响了该功能基因的10%，但是如果该实验一共只鉴定到20个差异基因，那么这10个差异基因就占所有差异基因的50%，因此，需要同时考虑差异基因在所有功能分类中的分布，得到的结果才是准确的。其中，a为功能A中差异基因的数量，b为非差异基因的数量，c为所有功能分类中差异基因的总数量，d为非差异基因的总数量，n为识别到的所有差异基因的数目。今天我们来聊聊怎么来展示GO富集分析的结果，下面是一个GO富集分析的结果表。
KEGG功能注释工具 KofamKOALA 安装与使用
刘永鑫的博客——宏基因组公众号
04-06 9602
KEGG数据库，即京都基因和基因组百科全书（Kyoto Encyclopedia of Genes and Genomes），是系统分析基因功能、基因组信息的数据库。KofamKOALA...
富集分析原理和clusterProfiler包进行GO、KEGG富集分析详细说明
Tian問的博客
10-03 2万+
概念： 基因富集分析是指对于给定一组基因根据基因组注释信息（GO、KEGG）对基因进行聚类分析，即给定的基因是不是GO中的一个功能（或KEGG中的一个通路）。 基因的功能富集的目的说明给定的基因集对哪些功能的影响有针对性的，不是随机影响的。 原理： 基因富集分析是通过研究给定的基因集在功能节点上是否过出现来得到关注的基因集显著注释的功能节点。通常利用超几何分布等方法计算给定基因集在某个功能（或通路）上的P值，判断给定的基因集在功能（或通路）的基因数目超过了随机的期望，是一个小概率事件。 变
KEGG-GO注释分析代码.R
05-10
转录组分析差异表达基因注释时，用R语言的biocondcouter分析,根据代码，搜索需要分析的物种，就可出来结果
KEGG注释步骤（含KAAS网页运行，KofamKOALA本地安装及批量）
weixin_58269397的博客
04-01 886
再者，如果使用的是ensembl ID，则需要注意统一使用Transcript ID或Gene ID）。保留所有信息，包括Gene ID对应上的每一个KO Number，比对分数，E-value以及KO Number的详细信息等。：需要保持ID的一致性，否则注释不出结果（例如拿来做注释的往往是蛋白序列，一般会包括可变性剪切的标识符AT4G36920.1。需要填写邮箱信息，在网页中提交后需要再进入邮箱，点击邮件中的提交链接，才能开始计算。注释完毕后，从填写的邮箱中进入注释完毕后的网页，下载注释结果。
【R语言】——基因GO/KEGG富集分析!超级简单的保姆级教程！
weixin_54004950的博客
01-09 2万+
GO/KEGG功能富集分析中重要的是背景基因的选择，使用R clusterProfiler包对基因进行富集，需要导入目的基因（前景基因）相对应物种的参考基因组（背景基因），现阶段“bioconductor”已有十几种常见动物，如人类、小鼠等物种的OrgDb。”介绍如何使用DAVID在线分析工具对基因进行GO/KEGG功能富集分析和使用R ggplot包对获得的基因GO/KEGG功能富集结果进行可视化。：富集到该GO term/KEGG term中的基因数目/给定基因的总数目；
R绘图-KEGG功能注释组间差异分面条形图
qq_39859424的博客
04-27 954
一、 数据准备 KEGG功能注释的差异分析结果绘图。使用的是经过Pairwise Mann–Whitney U-tests检验的处理间class II功能丰度差异结果数据，包含class I作为分面信息进行条形图绘制。 library(ggplot2) library(ggsci) library(ggpubr) library(scales) # 1.1 读入数据 data = read.csv("KEGG.csv",header=TRUE,sep=",",stringsAsFactors = F
给定基因集gene list后，R语言GO、KEGG分析
m0_58549466的博客
03-25 6370
gene list到GO和KEGG
用R做GO功能注释和KEGG通路富集分析
热门推荐
Cassiel60的博客
04-19 2万+
KEGG -GO主要是使用R 中clusterProfiler包进行富集分析 以及使用pathview 包进行代谢途径整合和可视化。 进行在线分析的：https://pathways.embl.de/ 1.代码：https://github.com/Cassiel60/R/blob/master/enrichment.r 2.安装包：使用代码中的方法要多尝试几次，反正我第一次安装时就是安装...
0055-【生物数据库】-如何进行RNA差异基因KEGG注释分析-kobas在线分析
leadingsci的博客
06-29 1万+
1. 有参物种使用gene ID的方法 1. 差异基因文件准备 只需要用到两列 ENTREZ_GENE_ID logFC geneNames ENTREZ_GENE_ID normalAve tumorAve logFC pValue qValue CCL23 6368 95.05964624 5.566645819 -4.066608903 ...
基因富集分析
fengleqi的博客
06-15 7690
https://www.jianshu.com/p/a1f312b60569http://www.360doc.com/content/16/0711/10/26456292_574656968.shtmldiff_gene.entrez文件，是通过各种差异基因软件找出来的差异基因的entrez ID号列表，每一个ID号一行，几百个差异基因就几百行http://blog.sciencenet.cn...
R语言：GO富集和KEGG富集
08-23
R语言中可以使用许多不同的包来进行GO富集和KEGG富集分析。下面是一个使用clusterProfiler包进行富集分析的示例代码： 1. 安装和加载必要的包： ```R install.packages("clusterProfiler") library(clusterProfiler) ``` 2. 对基因列表进行GO富集分析： ```R # 假设有一个基因列表gene_list，包含感兴趣的基因 result_GO <- enrichGO(gene = gene_list, OrgDb = org.Hs.eg.db, # 数据库，例如Human keyType = "ENSEMBL", # 基因ID类型 ont = "BP", # 生物过程 pAdjustMethod = "BH", # 多重检验校正方法 pvalueCutoff = 0.05, # 显著性水平 qvalueCutoff = 0.2) # FDR阈值 # 查看富集结果 result_GO ``` 3. 对基因列表进行KEGG富集分析： ```R result_KEGG <- enrichKEGG(gene = gene_list, organism = "hsa", # 物种，例如Human pvalueCutoff = 0.05, # 显著性水平 qvalueCutoff = 0.2) # FDR阈值 # 查看富集结果 result_KEGG ``` 这些代码会返回富集分析的结果，包括富集的GO term和KEGG pathway，以及统计学上的显著性。你可以根据需要进一步处理和可视化这些结果。

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值