专栏名称: 生信技能树
生物信息学学习资料分析,常见数据格式及公共数据库资料分享。常见分析软件及流程,基因检测及癌症相关动态。
目录
相关文章推荐
深圳特区报  ·  深圳市政协七届五次会议定于2月24日召开 ·  19 小时前  
深圳大件事  ·  临时停运!深圳一景区发布通知 ·  2 天前  
深圳发布  ·  U20亚洲杯来啦!还有观赛指南→ ·  3 天前  
深圳特区报  ·  被举报教师唐某某,已被采取刑事强制措施! ·  3 天前  
51好读  ›  专栏  ›  生信技能树

如何让你的kegg注释结果图分门别类

生信技能树  · 公众号  ·  · 2024-05-09 09:13

正文

酒香也怕巷子深啊!好像大家不怎么知道我们 生信技能树 团队有一个 生物信息学入门课 ,详见; 生物信息学马拉松授课(买一得五) ,你值得拥有!

KEGG数据库是一个综合性的生物信息数据库,由日本京都大学生物信息学中心的Kanehisa实验室于1995年建立。它整合了基因组、化学和系统功能信息,旨在从分子水平上理解生物系统的高级功能和实用程序,特别是细胞、生物体和生态系统的功能。

KEGG数据库的分类情况相当广泛,它将信息分为三大类:系统信息、基因组信息和化学信息。进一步细分为16个主要的数据库。这些数据库包括:

  1. 系统信息 :涉及生物系统的高级功能和实用程序,如KEGG PATHWAY、KEGG BRITE、KEGG MODULE、KEGG DISEASE、KEGG DRUG和KEGG ENVIRON等。
  2. 基因组信息 :包括KEGG ORTHOLOGY、KEGG GENOME、KEGG GENES和KEGG SSDB等,这些数据库提供了基因组和基因相关的信息。
  3. 化学信息 :如KEGG COMPOUND、KEGG GLYCAN、KEGG REACTION等,这些数据库包含了关于化学物质、代谢物、生化反应等的信息。

虽然KEGG数据库看起来非常复杂,但是绝大部分让接触到KEGG数据库其实是里面的 PATHWAY数据库。

什么是KEGG PATHWAY数据库

这个是是一个手工绘制的代谢通路集合,它包含了分子间相互作用和反应网络的信息,覆盖了新陈代谢、环境信息加工、细胞过程、生物体系统和人类疾病等方面。

手工绘制的代谢通路集合

PATHWAY数据库中的通路被分为7类:

  • 细胞过程(Cellular Processes)
  • 环境信息处理(Environmental Information Processing)
  • 遗传信息处理(Genetic Information Processing)
  • 人类疾病(Human Diseases)
  • 新陈代谢(Metabolism)
  • 生物体系统(Organismal Systems)
  • 药物相关( Drug Development

对于人类(Homo sapiens)这个物种,KEGG数据库提供了物种特异性的通路,这些通路被详细地分类和注释。例如,在KEGG PATHWAY数据库中,可以通过特定的物种名称缩写(如hsa代表人类)来检索人类的特定通路,如hsa00010代表人类的糖酵解/糖异生通路。

少量基因做KEGG PATHWAY数据库注释

如果只有少量基因,其实每个基因去查询一下,就可以看到它的KEGG PATHWAY数据库注释信息啦, 每个基因都是独立的查询即可。比如大名鼎鼎的TP53,TP53基因大家都有所耳闻,而且也大概都知道它是  tumor protein p53的简称,其实它还有很多别名,比如BCC7;LFS1;P53;TRP53; :

  • https://www.genecards.org/cgi-bin/carddisp.pl?gene=TP53

可以看到其实是可以做很多数据库注释,但是绝大部分人只知道KEGG PATHWAY数据库。

很多数据库注释

如果仅仅是想做KEGG PATHWAY数据库注释,也是可以使用r代码的,借助这个 library(org.Hs.eg.db) 包即可:

rm(list=ls())
library(org.Hs.eg.db)
eg2symbol=toTable(org.Hs.egSYMBOL)
eg2name=toTable(org.Hs.egGENENAME)
eg2alias=toTable(org.Hs.egALIAS2EG)
eg2alis_list=lapply(split(eg2alias,eg2alias$gene_id),function(x){paste0(x[,2],collapse = ";")})
GeneList=mappedLkeys(org.Hs.egSYMBOL)

GeneList[1]
if( GeneList[1] %in% eg2symbol$symbol ){
  symbols=GeneList
  geneIds=eg2symbol[match(symbols,eg2symbol$symbol),'gene_id']
}else{
  geneIds=GeneList
  symbols=eg2symbol[match(geneIds,eg2symbol$gene_id),'symbol']
}
geneNames=eg2name[match(geneIds,eg2name$gene_id),'gene_name']
geneAlias=sapply(geneIds,function(x){ifelse(is.null(eg2alis_list[[x]]),"no_alias",eg2alis_list[[x]])})

借助于 org.Hs.eg.db 包,我们已经拿到了全部的人类全部基因的全名和别名,就是如下所示的4个变量:

  • geneIds
  • symbols
  • geneNames
  • geneAlias

同理,任意基因都可以去查询在 org.Hs.eg.db 包里面的有的信息:

> ls("package:org.Hs.eg.db")
 [1] "org.Hs.eg"                "org.Hs.eg_dbconn"         "org.Hs.eg_dbfile"        
 [4] "org.Hs.eg_dbInfo"         "org.Hs.eg_dbschema"       "org.Hs.eg.db"            
 [7] "org.Hs.egACCNUM"          "org.Hs.egACCNUM2EG"       "org.Hs.egALIAS2EG"       
[10] "org.Hs.egCHR"             "org.Hs.egCHRLENGTHS"      "org.Hs.egCHRLOC"         
[13] "org.Hs.egCHRLOCEND"       "org.Hs.egENSEMBL"         "org.Hs.egENSEMBL2EG"     
[16] "org.Hs.egENSEMBLPROT"     "org.Hs.egENSEMBLPROT2EG"  "org.Hs.egENSEMBLTRANS"   
[19] "org.Hs.egENSEMBLTRANS2EG" "org.Hs.egENZYME"          "org.Hs.egENZYME2EG"      
[22] "org.Hs.egGENENAME"        "org.Hs.egGENETYPE"         "org.Hs.egGO"             
[25] "org.Hs.egGO2ALLEGS"       "org.Hs.egGO2EG"           "org.Hs.egMAP"            
[28] "org.Hs.egMAP2EG"          "org.Hs.egMAPCOUNTS"       "org.Hs.egOMIM"           
[31] "org.Hs.egOMIM2EG"         "org.Hs.egORGANISM"        "org.Hs.egPATH"           
[34] "org.Hs.egPATH2EG"         "org.Hs.egPFAM"            "org.Hs.egPMID"           
[37] "org.Hs.egPMID2EG"         "org.Hs.egPROSITE"         "org.Hs.egREFSEQ"         
[40] "org.Hs.egREFSEQ2EG"       "org.Hs.egSYMBOL"          "org.Hs.egSYMBOL2EG"      
[43] "org.Hs.egUCSCKG"          "org.Hs.egUNIPROT"     

单个基因或者多个基因都可以通过代码,一个个的去注释,但是很浪费时间,因为如果基因很少量,大概率上大家 就直接使用网页工具了。

多个基因做注释需要统计学(超几何分布检验)

如果是少量的几个基因,确实是可以一个个去网页工具里面查询获取具体的信息,但是如果是几十个或者上百个基因,比如我们做完差异分析后的基因列表,就需要借助统计学(超几何分布检验),在KEGG数据库注释中使用统计学方法,尤其是超几何分布检验,是为了识别差异基因集中富集的生物学通路。以下是该过程的一般步骤:

  1. 差异基因集的获取 :首先,你需要从一个实验条件(如疾病状态)与另一个实验条件(如正常状态)的比较中获得一组差异表达的基因。

  2. 基因功能注释 :将这些差异基因与KEGG数据库中的通路进行匹配。这通常涉及到将基因标识符(如Ensembl ID、GenBank Accession号码等)转换为KEGG Orthology (KO)编号。

  3. 超几何分布检验 :利用超几何分布检验来确定某个通路中的基因在差异基因集中是否显著富集。超几何分布检验的基本原理是评估从特定大小的群体(背景基因集)中随机抽取一定数量的样本(前景基因集,即差异基因集)时,获得特定特征(如属于某个KEGG通路)的基因数目的概率。

    检验的统计公式通常为:

    其中:

  • ( N ) 是背景基因集中所有基因的总数。
  • ( M ) 是背景基因集中具有某种特定功能的基因数。
  • ( n ) 是前景基因集中的基因数。
  • ( k ) 是前景基因集中具有该特定功能的基因数。
  • ( p ) 是在背景基因集中抽取到具有该功能基因的概率,( p = \frac{M}{N} )。
  • 多重检验校正 :由于一次分析中可能涉及多个通路的检验,因此需要使用如Bonferroni校正或False Discovery Rate (FDR)等方法来校正P值,以减少假阳性结果。

  • 结果解释 :最后,根据校正后的P值确定哪些通路是显著富集的,并据此解释差异基因集中的生物学功能。

  • 可视化 :使用图表或富集图(如气泡图、条形图)来直观展示富集分析的结果。

  • 这个过程涉及到对统计学原理的深刻理解以及对KEGG数据库结构的熟悉。在实际应用中,许多生物信息学工具和软件包(如R包clusterProfiler)可以自动化这些步骤,简化富集分析的过程。

    当然了,看不懂上面的统计学文字描述很正常, 你可以试试看朴素的理解它,当然可以。超几何分布检验可以用一个更朴素的方式来理解,通过一个简单的例子来说明。

    假设你有一个装满球的箱子,其中一些球是红色的,其余的是蓝色的。红色球和蓝色球的比例代表了整个群体(也就是背景)中特定特征(比如红色)的比例。现在,你闭上眼睛从箱子里随机取出几个球,组成一个较小的样本组(也就是前景)。你感兴趣的是,取出的样本组中红色球的比例是否显著高于整个箱子中红色球的比例。

    超几何分布检验的朴素理解

    1. 随机抽取 :你从箱子里随机抽取球(不把抽出来的球放回去),看能抽到多少红球。

    2. 期望比例 :整个箱子中红球和蓝球的比例给你一个期望值,也就是在随机抽取的情况下,你期望抽到的红球的比例。

    3. 实际观察 :你实际抽出来的样本中红球的比例与你的期望值进行比较。

    4. 检验统计量 :超几何分布检验会计算在随机情况下,抽取的样本中红球比例至少和你实际观察到的一样高(或更高)的概率。

    5. 显著性 :如果你观察到的红球比例远高于期望值,并且这种情况发生的概率非常低,那么你可能得出结论:抽取的样本中红球是“显著富集”的。

    应用到KEGG富集分析

    在KEGG富集分析中,这个例子可以这样对应:

    • 箱子和球 :箱子代表你研究的生物体的所有基因(背景基因集),球代表这些基因。红色球代表参与某个特定生物学过程或通路的基因。

    • 抽取样本 :你从箱子中抽取的球代表你通过实验发现的差异基因集(前景基因集)。

    • 期望与观察 :你期望前景基因集中参与特定通路的基因比例与背景基因集中的比例相似。但如果你观察到前景基因集中该通路的基因显著更多,这可能意味着该通路在你的实验条件下起到了特别的作用。

    • 超几何分布检验 :通过这个检验,你可以统计上判断前景基因集中某个通路的富集是随机发生的,还是很可能反映了真实的生物学现象。

    当然了,如果你尝试去理解这个故事本身,也很难朴素的理解超几何分布检验,另外一个比较好的方式就是通过图像来展示,比如下面的韦恩图 ;

    image-20240509085216820

    其中的deg就是你需要注释的基因列表,通常是几十个甚至上百个,然后FRG就是kegg等数据库里面的一个pathway啦,上面的韦恩图(Venn diagram)展现了两个基因列表的交集的数量,就可以通过超几何分布检验判定这个交集数量是否有统计学显著性。

    我们通常是不需要一个个去手工整理每个pathway的基因列表,然后一个个去跟我们的需要注释的基因列表去做超几何分布检验的, 因为在代码里面都是批量而且自动化的!如下所示的代码就是一个很好的例子:

    library(clusterProfiler)
    library(ggplot2)
    data(geneList, package='DOSE')
    de 1:100]
    yy 0.05)
    head(yy@result[,1:3])
    dotplot(yy)

    如果是初次接触r语言,可能会对上面的代码比较困惑。其实就是使用了 clusterProfiler 包来进行KEGG富集分析,并且利用 ggplot2 包来生成可视化图形。下面是对每一行代码的详细解释:

    1. library(clusterProfiler) 这行代码加载了 clusterProfiler 包,这是一个用于生物数据富集分析的R包,可以进行基因本体(GO)和KEGG通路富集分析。
    2. library(ggplot2) 加载 ggplot2 包,这是一个流行的R包,用于创建高质量的统计图形。
    3. data(geneList, package='DOSE') 使用 data() 函数从名为 DOSE 的R包中加载名为 geneList 的数据集。 DOSE 包是一个包含多个用于富集分析的数据集的包。
    4. de 这行代码从 geneList 数据集中提取前100个基因的名称,并将它们存储在名为 de 的变量中。 names(geneList) 获取 geneList 中所有基因的名称, [1:100] 选取前100个。






    请到「今天看啥」查看全文