GO和KEGG富集分析详细步骤

GO和KEGG富集分析


1. 将差异表达结果的基因名称转化为id

因为GO和KEGG分析需要用到id,所以这一步需要将基因名字转换为id。具体步骤如下:

  1. 新建空白文件夹,将差异分析得到的diff.xls复制粘贴到文件夹中

  2. 因为在这里只需要diff.xls中的基因名称和logFC两列,所以只复制这两列粘贴到新建的文本文件symbol.txt,如下图所示:
    在这里插入图片描述

  3. 新建R语言脚本文件symbol2id.R,代码如下:

    if (!requireNamespace("BiocManager", quietly = TRUE))
        install.packages("BiocManager")
    BiocManager::install("org.Hs.eg.db")
    
    
    setwd("C:\\Users\\Administrator\\Desktop\\cptac\\4_name2id")          #设置工作目录
    
    library("org.Hs.eg.db")          #引用包
    rt=read.table("symbol.txt",sep="\t",check.names=F,header=T)    #读取文件
    genes=as.vector(rt[,1])
    entrezIDs <- mget(genes, org.Hs.egSYMBOL2EG, ifnotfound=NA)    #找出基因对应的id
    entrezIDs <- as.character(entrezIDs)
    out=cbind(rt,entrezID=entrezIDs)
    write.table(out,file="id.txt",sep="\t",quote=F,row.names=F) 
    
  4. 设置好工作目录之后,打开R软件,运行上述代码即可。运行结束在文件夹中会有id.txt,打开后如下图所示:
    在这里插入图片描述

    可以看到后面已经有了id这一列了,至此本步骤结束。

2. GO富集分析

GO(gene ontology)是基因本体联合会(Gene Onotology Consortium)所建立的数据库,旨在建立一个适用于各种物种的、对基因和蛋白质功能进行限定和描述的、并能随着研究不断深入而更新的语言词汇标准。GO是多种生物本体语言中的一种,提供了三层结构的系统定义方式,用于描述基因产物的功能。在转录组项目中,GO功能分析一方面给出差异表达转录本的GO功能分类注释;另一方面给出差异表达转录本的GO功能显著性富集分析。

下面介绍GO分析的步骤:

  1. 将含有基因id的文本文件id.txt复制粘贴到新的文件夹中

  2. 新建R语言脚本,命名为GO.R,其代码如下:

    install.packages("colorspace")
    install.packages("stringi")
    install.packages("ggplot2")
    
    if (!requireNamespace("BiocManager", quietly = TRUE))
        install.packages("BiocManager")
    BiocManager::install("DOSE")
    
    if (!requireNamespace("BiocManager", quietly = TRUE))
        install.packages("BiocManager")
    BiocManager::install("clusterProfiler")
    
    if (!requireNamespace("BiocManager", quietly = TRUE))
        install.packages("BiocManager")
    BiocManager::install("enrichplot")
    
    library("clusterProfiler")
    library("org.Hs.eg.db")
    library("enrichplot")
    library("ggplot2")
    
    setwd("C:\\Users\\Administrator\\Desktop\\cptac\\5_GO分析")                  #设置工作目录
    rt=read.table("id.txt",sep="\t",header=T,check.names=F)           #读取id.txt文件
    rt=rt[is.na(rt[,"entrezID"])==F,]                                 #去除基因id为NA的基因
    gene=rt$entrezID
    
    #GO富集分析
    kk <- enrichGO(gene = gene,
                   OrgDb = org.Hs.eg.db, 
                   pvalueCutoff =0.05, 
                   qvalueCutoff = 0.05,
                   ont="all",
                   readable =T)
    write.table(kk,file="GO.txt",sep="\t",quote=F,row.names = F)                 #保存富集结果
    
    #柱状图
    pdf(file="barplot.pdf",width = 10,height = 8)
    barplot(kk, drop = TRUE, showCategory =10,split="ONTOLOGY") + facet_grid(ONTOLOGY~., scale='free')
    dev.off()
    
    #气泡图
    pdf(file="bubble.pdf",width = 10,height = 8)
    dotplot(kk,showCategory = 10,split="ONTOLOGY",orderBy = "GeneRatio") + facet_grid(ONTOLOGY~., scale='free')
    dev.off()
    

    这里GO分析用到的包为"clusterProfiler",画图用到的包为"enrichplot"。在代码中会设置p值和q值,设置的都是0.05,如果该条件下分析得到的可用基因较少,可将q设置为0,只看p值,但这样准确性也会降低一些。

  3. 打开R软件,运行上述代码,最终得到的结果如下图所示,下图按顺序分别是柱状图、气泡图以及GO分析结果。
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述

  4. 讲一下GO分析得到的文本文件,也就是上面三幅图中的最后一幅图,第一列是GO分析的分类,分别是BP,CC,MF;第二列是GO的id;第三列为对应的描述;第四列为基因背景的比例;第五列为p值,表示富集的显著性;第六列为p值得校正值;第七列为q值;第八列为基因id,也就是基因名称;最后一列就是富集在每个GO上的数目。对于柱状图和气泡图,会分为BP,CC,MF,每个类别颜色越红表示富集程度越高。

3. GO圈图绘制

话不多说,直接上步骤。

  1. 新建R语言脚本文件GOplot.R,脚本文件和GO分析得到的结果放在同一目录下,其代码如下:

    install.packages("digest")
    install.packages("GOplot")
    
    library(GOplot)
    setwd("C:\\Users\\Administrator\\Desktop\\cptac\\6_GO圈图绘制")              #设置工作目录
    
    ego=read.table("GO.txt", header = T,sep="\t",check.names=F)      #读取kegg富集结果文件
    go=data.frame(Category = "All",ID = ego$ID,Term = ego$Description, Genes = gsub("/", ", ", ego$geneID), adj_pval = ego$p.adjust)
    
    #读取基因的logFC文件
    id.fc <- read.table("id.txt", header = T,sep="\t",check.names=F)
    genelist <- data.frame(ID = id.fc$gene, logFC = id.fc$logFC)
    row.names(genelist)=genelist[,1]
    
    circ <- circle_dat(go, genelist)
    termNum = 5                                     #限定term数目
    geneNum = nrow(genelist)                        #限定蛋白数目
    
    chord <- chord_dat(circ, genelist[1:geneNum,], go$Term[1:termNum])
    pdf(file="circ.pdf",width = 11,height = 10.5)
    GOChord(chord, 
            space = 0.001,           #基因之间的间距
            gene.order = 'logFC',    #按照logFC值对基因排序
            gene.space = 0.25,       #基因名跟圆圈的相对距离
            gene.size = 4,           #基因名字体大小 
            border.size = 0.1,       #线条粗细
            process.label = 7.5)     #term字体大小
    dev.off()
    
    termCol <- c("#223D6C","#D20A13","#FFD121","#088247","#58CDD9","#7A142C","#5D90BA","#431A3D","#91612D","#6E568C","#E0367A","#D8D155","#64495D","#7CC767")
    pdf(file="cluster.pdf",width = 11.5,height = 9)
    GOCluster(circ.gsym, 
              go$Term[1:termNum], 
              lfc.space = 0.2,                   #倍数跟树间的空隙大小
              lfc.width = 1,                     #变化倍数的圆圈宽度
              term.col = termCol[1:termNum],     #自定义term的颜色
              term.space = 0.2,                  #倍数跟term间的空隙大小
              term.width = 1)                    #富集term的圆圈宽度
    dev.off()          
    
  2. 打开R软件运行上述代码即可。最终即可得到两个圈图,如下图所示:
    在这里插入图片描述

    左半圆圈为基因名字,从下到上按照logFC进行排序得,圆圈右半部分为GO的名称,基因与GO之间得连线表示这个基因存在于该GO上。
    在这里插入图片描述

    此图为聚类图,内部圆圈为基因或蛋白,颜色表示logFC的大小,内部的一个扇形表示一个基因,如果内部的一个扇形对应着外部的一个颜色的扇形,那么表示该基因只存在于这一个颜色对应的GO里面;如果内部一个扇形对应着外部三个扇形,那么表示内部的这个基因存在于三个GO里面。

4. KEGG富集分析

  1. 将差异分析得到的含有id的id.txt文件作为输入文件,新建文件夹,将id.txt拷贝到此文件夹下

  2. 新建R语言脚本文件,更改脚本文件的环境目录,代码如下:

    install.packages("colorspace")
    install.packages("stringi")
    install.packages("ggplot2")
    
    if (!requireNamespace("BiocManager", quietly = TRUE))
        install.packages("BiocManager")
    BiocManager::install("DOSE")
    
    if (!requireNamespace("BiocManager", quietly = TRUE))
        install.packages("BiocManager")
    BiocManager::install("clusterProfiler")
    
    if (!requireNamespace("BiocManager", quietly = TRUE))
        install.packages("BiocManager")
    BiocManager::install("enrichplot")
    
    
    library("clusterProfiler")
    library("org.Hs.eg.db")
    library("enrichplot")
    library("ggplot2")
    
    setwd("C:\\Users\\Administrator\\Desktop\\cptac\\7_KEGG分析")            #设置工作目录
    rt=read.table("id.txt",sep="\t",header=T,check.names=F)       #读取id.txt文件
    rt=rt[is.na(rt[,"entrezID"])==F,]                             #去除基因id为NA的基因
    gene=rt$entrezID
    
    #kegg富集分析
    kk <- enrichKEGG(gene = gene, organism = "hsa", pvalueCutoff =0.05, qvalueCutoff =2)      #富集分析
    write.table(kk,file="KEGGId.txt",sep="\t",quote=F,row.names = F)                          #保存富集结果
    
    #柱状图
    pdf(file="barplot.pdf",width = 10,height = 7)
    barplot(kk, drop = TRUE, showCategory = 30)
    dev.off()
    
    #气泡图
    pdf(file="bubble.pdf",width = 10,height = 7)
    dotplot(kk, showCategory = 30,orderBy = "GeneRatio")
    dev.off()
    
  3. 打开R软件运行上述代码,即可得到结果。
    在这里插入图片描述
    在这里插入图片描述

    KEGG因为数据库更新比较慢,而且分析时需要联网,因此富集到结果就会比较少。

  4. 运行完之后还会得到KEGGId.txt,里面的需要将里面的id转化为基因名字。因此新建perl脚本文件,代码太长,这里就不展示了。在该文件夹目录下打开powershell窗口,输入命令perl id2symbol.pl,运行完毕之后文件夹目录下就会产生新的含有基因名字的kegg文件,文件名为kegg.txt

  5. 至此,KEGG分析完毕

5. KEGG圈图绘制

这里的圈图绘制和上面的GO圈图绘制步骤一样的。话不多说,直接上代码:

install.packages("digest")
install.packages("GOplot")

library(GOplot)
setwd("C:\\Users\\Administrator\\Desktop\\cptac\\8_KEGG圈图绘制")              #设置工作目录

ego=read.table("kegg.txt", header = T,sep="\t",check.names=F)       #读取kegg富集结果文件
go=data.frame(Category = "All",ID = ego$ID,Term = ego$Description, Genes = gsub("/", ", ", ego$geneID), adj_pval = ego$p.adjust)

#读取基因的logFC文件
id.fc <- read.table("id.txt", header = T,sep="\t",check.names=F)
genelist <- data.frame(ID = id.fc$gene, logFC = id.fc$logFC)
row.names(genelist)=genelist[,1]

circ <- circle_dat(go, genelist)
termNum = 2                                     #限定term数目
geneNum = nrow(genelist)                        #限定基因数目

chord <- chord_dat(circ, genelist[1:geneNum,], go$Term[1:termNum])
pdf(file="circ.pdf",width = 10,height = 9.6)
GOChord(chord, 
        space = 0.001,           #基因之间的间距
        gene.order = 'logFC',    #按照logFC值对基因排序
        gene.space = 0.25,       #基因名跟圆圈的相对距离
        gene.size = 4,           #基因名字体大小 
        border.size = 0.1,       #线条粗细
        process.label = 7.5)     #term字体大小
dev.off()

termCol <- c("#223D6C","#D20A13","#FFD121","#088247","#58CDD9","#7A142C","#5D90BA","#431A3D","#91612D","#6E568C","#E0367A","#D8D155","#64495D","#7CC767")
pdf(file="cluster.pdf",width = 10,height = 9.6)
GOCluster(circ.gsym, 
          go$Term[1:termNum], 
          lfc.space = 0.2,                   #倍数跟树间的空隙大小
          lfc.width = 1,                     #变化倍数的圆圈宽度
          term.col = termCol[1:termNum],     #自定义term的颜色
          term.space = 0.2,                  #倍数跟term间的空隙大小
          term.width = 1)                    #富集term的圆圈宽度
dev.off()

这里将代码的工作环境更改一下,然后将kegg分析所得到的kegg.txt和之前的id.txt复制到同一目录下,然后打开R软件运行代码即可。得到的圈图如下:
在这里插入图片描述
在这里插入图片描述

至此,KEGG圈图绘制结束。

  • 39
    点赞
  • 347
    收藏
    觉得还不错? 一键收藏
  • 11
    评论
GO富集分析是一种常用的生物信息学方法,用于确定一组基因在特定的生物学过程、细胞组分或分子功能中的富集情况。在R语言中,可以使用多个包来进行GO富集分析。 其中,clusterProfiler包是一个功能强大的包,主要用于进行GO和KEGG的功能富集分析,并提供了可视化功能。org.Hs.eg.db包用于转换不同数据库中的基因ID和符号之间的转换。enrichplot包提供了多种可视化方法来解释富集结果。而GOplot包则用于绘制功能富集的图形。\[1\] 在进行GO富集分析时,可以使用enrichGO函数来进行分析。该函数需要提供基因列表和参考基因组等参数。可以设置P值和q值的阈值,以及选择分析的层面(生物学过程、细胞组分或分子功能)。还可以选择是否将基因ID转换为基因名。\[2\] 在分析完GO富集后,可以根据显著性阈值筛选出显著富集的结果,并将结果保存为文件。可以构建数据框矩阵来存储GO富集的相关信息,包括分类、ID、术语、基因和校正的P值等。可以根据需要限定GO数目和基因数目。\[3\] 总之,使用R语言进行GO富集分析可以通过多个包来实现,包括clusterProfiler、org.Hs.eg.db、enrichplot和GOplot。可以根据具体需求设置参数并进行分析,最后可以筛选出显著富集的结果并进行可视化。 #### 引用[.reference_title] - *1* *2* *3* [R语言|GO富集分析](https://blog.csdn.net/weifanbio/article/details/124279953)[target="_blank" data-report-click={"spm":"1018.2226.3001.9630","extra":{"utm_source":"vip_chatgpt_common_search_pc_result","utm_medium":"distribute.pc_search_result.none-task-cask-2~all~insert_cask~default-1-null.142^v91^insert_down1,239^v3^insert_chatgpt"}} ] [.reference_item] [ .reference_list ]

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

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

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值