生信实战：教科书级手把手教你做GO富集分析！超级简单超级好看！

作者: xizzy | 来源:发表于2021-10-30 00:33 被阅读0次

2021都要过去了，你还在为做GO富集而烦恼吗？还在为垃圾公司做的GO效果差而苦恼吗？何不自己动手做一个呢（当然可以私信我给你代做）。下面我就手把手教你做sci级的GO富集分析+可视化吧

准备

做GO富集分析前，我们需要有我们的基因列表，通常我们的基因列表来自于我们差异分析得到的显著差异表达分析，有了这个列表后我们就可以做富集分析啦。不过在此之前你还需要做两件事情：

安装R和R包clusterProfiler、ggplot2，安装方法两个R包的方法:

if (!requireNamespace("BiocManager", quietly = TRUE))
              install.packages("BiocManager")
BiocManager::install("clusterProfiler")
install.packages("ggplot2")

当然你如果安装了conda也可以使用conda进行安装：

conda install r-ggplot2 -c bioconda -c conda-forge
conda install bioconductor-clusterProfiler -c bioconda -c conda-forge

2.对基因列表去重，方法也很简单，如果使用linux只需要sort -u gene.txt > unqi.gene.txt即可吧得到去重的列表。

分析

识别基因ID类别

以上步骤都完成后我们就可以进行分析啦，我们要知道主要进行富集的基因ID有很多种，比如ENSEMBL（ENXXX数字形式）、ENTREZID（纯数字）、SYMBOL（如TP53等），在富集的时候我们需要指定我们的输入类型，clusterProfiler的建议类型是ENTREZID，不过既然我们作为教科书级别的教程，我们就需要对每一种类型进行适配：R代码如下：

GetGeneType<-function(df,filename){
    if(str_detect(df[1,1],"EN")){
        return("ENSEMBL")
    }else if(is.integer( as.integer(df[1,1]))){
        return("ENTREZID")
    }else{
        return("SYMBOL")
    }
}

GO富集

需要进行富集非常简单，只需要使用enrichGO()函数即可，但是在此之前我们需要先把物种的Db指定好，如果是人的那需要安装org.Hs.eg.db，小鼠org.Mm.eg.db，猪org.Ss.eg.db，拟南芥：org.At.tair.db，都可以通过BiocManager::install()安装。以人为例：

suppressMessages(library(clusterProfiler))
suppressMessages(library(org.Hs.eg.db))
Db <- org.Hs.eg.db
df <- read.table(file="genelist.txt",sep='\t',header=F) #读取你的genelist文件
go <- enrichGO(df[,1], OrgDb = Db, ont='ALL',pAdjustMethod = 'BH',pvalueCutoff = 0.05,keyType = GetGeneType(df,genelist))

其中几个比较关键的参数就是ont指定富集那种类型，这里我们bp/cc/mf都做所以选择all，pvalueCutoff是说通过P值去筛选我们富集显著水平，数字为阈值，keyType就是输入的基因类型啦，这里载入我们上面写的函数GetGeneType()去自动识别。

查看结果和可视化

如果查看结果，需要一句代码即可：
go_res <- as.data.frame(go@result)
如果需要保存结果为xls，只需要：
write.table(go_res,file="go_result.xls",sep="\t",row.names=F)
那么如何得到我们在文章或者公司报告中常见的柱状图呢，如下：

image.png
虽然clusterProfiler也提供了一些可视化的方法，但是还是比较粗糙的，远达不到上面的效果。下面我手把手教你怎么实现：

library(ggplot2)
library(stringr)
#GO Term的名字太长，为了避免太长的名字影响到美观，我们设置控制名字长短的函数，通过指定max_char参数可以进行调整
shorten_names <- function(char_array, max_char=30){
  for(i in 1:length(char_array)){
      if(nchar(char_array[i]) >  max_char){
          char_array[i] <- paste0(substr(char_array[i],0,max_char),"...")
      }
  }
  return(char_array)
}
#根据P值选取TOP富集的基因，item控制选取的个数，如果不足指定的item时，就取全部
GetTop <- function(df,item){
    if(nrow(df) > item){
        return(df[1:item,])
    }else{
        return(df[1:nrow(df),])
    }
}
#将3个生物过程的TOP的GO term选出来
GetGoTopEnrich <- function(go_res,max_item=15){
    go_MF <- GetTop(go_res[go_res$ONTOLOGY=="MF",],max_item)
    go_CC <- GetTop(go_res[go$ONTOLOGY=="CC",][1:max_item,],max_item)
    go_BP <- GetTop(go_res[go_res$ONTOLOGY=="BP",][1:max_item,],max_item)
    go_enrich_df <- data.frame(ID=c(go_BP$ID, go_CC$ID, go_MF$ID),
                            Description=c(go_BP$Description, go_CC$Description, go_MF$Description),
                            GeneNumber=c(go_BP$Count, go_CC$Count, go_MF$Count),
                            GO_Ontology=factor(c(rep("biological process", nrow(go_BP)), rep("cellular component", nrow(go_CC)),rep("molecular function",nrow(go_MF))),levels=c("molecular function", "cellular component", "biological process")))
    return(go_enrich_df)
}
#得到用于画图的数据框用于画图
go_enrich_df$number <- factor(rev(1:nrow(go_enrich_df))) #sort the axis
#通过reorder重新排序，fill指定三个组成的颜色，~facet_wrap进行分页，并调整边距，旋转x轴标签
gobar <- ggplot(data = go_enrich_df, aes(x = reorder(number,GeneNumber), y =GeneNumber, fill = GO_Ontology)) + geom_bar(stat = 'identity', position = 'dodge')  + theme_bw() + theme(panel.grid=element_blank()) + 
            theme(axis.text.x = element_text(angle = 45, hjust = 1))  +  scale_x_discrete(labels=shorten_names(go_enrich_df$Description))+
            scale_fill_manual(values = c("#63B2EE", "#76DA91", "#F89588")) + facet_wrap(~GO_Ontology, scales = 'free_x', ncol = 3) + 
            theme(plot.margin=unit(c(1,1,2,4),'lines'))  +# distance to up right down left 
            xlab('GO terms') + ylab('Gene counts') 
ggsave("gobar.png",width = 12, height = 8, dpi = 300, units = "in", device='png') #保存图片

至此，GO富集大功告成！

总结

还是太麻烦？太高深?没关系，我已经为大家配置好了人和小鼠的，有其他的研究的也可以自行下载有关的包进行调整。将下面的代码保存为EasyEnrich.r，我们输入命令Rscript EasyEnrich.r hs genelist.txt即可完成，其中EasyEnrich.r改为你具体保存这个代码的路径，hs为你的物种，可以是人hs小鼠mm，而genelist.txt就是你的基因列表，记得先去重~
代码如下：

#Author: Xizzy txizzy#gmail.com
args <- commandArgs(T)
suppressMessages(library(clusterProfiler))
library(stringr)
suppressMessages(library(ggplot2))
ref <- args[1]
genelist <- args[2]

GetGeneType<-function(df,filename){
    if(str_detect(df[1,1],"EN")){
        return("ENSEMBL")
    }else if(is.integer( as.integer(df[1,1]))){
        return("ENTREZID")
    }else{
        return("SYMBOL")
    }
}

shorten_names <- function(char_array, max_char=30){
  for(i in 1:length(char_array)){
      if(nchar(char_array[i]) >  max_char){
          char_array[i] <- paste0(substr(char_array[i],0,max_char),"...")
      }
  }
  return(char_array)
}

GetTop <- function(df,item){
    if(nrow(df) > item){
        return(df[1:item,])
    }else{
        return(df[1:nrow(df),])
    }
}

GetGoTopEnrich <- function(go_res,max_item=15){
    go_MF <- GetTop(go_res[go_res$ONTOLOGY=="MF",],max_item)
    go_CC <- GetTop(go_res[go$ONTOLOGY=="CC",][1:max_item,],max_item)
    go_BP <- GetTop(go_res[go_res$ONTOLOGY=="BP",][1:max_item,],max_item)
    go_enrich_df <- data.frame(ID=c(go_BP$ID, go_CC$ID, go_MF$ID),
                            Description=c(go_BP$Description, go_CC$Description, go_MF$Description),
                            GeneNumber=c(go_BP$Count, go_CC$Count, go_MF$Count),
                            GO_Ontology=factor(c(rep("biological process", nrow(go_BP)), rep("cellular component", nrow(go_CC)),rep("molecular function",nrow(go_MF))),levels=c("molecular function", "cellular component", "biological process")))
    return(go_enrich_df)
}

df <- read.table(file=genelist,sep='\t',header=F)
if(ref=='hs'){
    suppressMessages(library(org.Hs.eg.db))
    Db <- org.Hs.eg.db
}else if(ref=='mm'){
    suppressMessages(library(org.Mm.eg.db))
    Db <- org.Mm.eg.db
}

go <- enrichGO(df[,1], OrgDb = Db, ont='ALL',pAdjustMethod = 'BH',pvalueCutoff = 0.05,keyType = GetGeneType(df,genelist))
go_res <- as.data.frame(go@result)
write.table(go_res,file="go_table.xls",row.names=F)
go_enrich_df <- GetGoTopEnrich(go_res)
go_enrich_df$number <- factor(rev(1:nrow(go_enrich_df))) #sort the axis
gobar <- ggplot(data = go_enrich_df, aes(x = reorder(number,GeneNumber), y =GeneNumber, fill = GO_Ontology)) + geom_bar(stat = 'identity', position = 'dodge')  + theme_bw() + theme(panel.grid=element_blank()) + 
            theme(axis.text.x = element_text(angle = 45, hjust = 1))  +  scale_x_discrete(labels=shorten_names(go_enrich_df$Description))+
            scale_fill_manual(values = c("#63B2EE", "#76DA91", "#F89588")) + facet_wrap(~GO_Ontology, scales = 'free_x', ncol = 3) + 
            theme(plot.margin=unit(c(1,1,2,4),'lines'))  +# distance to up right down left 
            xlab('GO terms') + ylab('Gene counts') 
ggsave("gobar.png",width = 12, height = 8, dpi = 300, units = "in", device='png')

最后：如果想了解更多和生信或者精品咖啡有关的内容欢迎关注我的微信公众号：生信咖啡，更多精彩等你发现！

网友评论

本文标题：生信实战：教科书级手把手教你做GO富集分析！超级简单超级好看！

本文链接：https://www.haomeiwen.com/subject/yzgwaltx.html

延伸阅读

深度阅读

您也可以注册成为美文阅读网的作者，发表您的原创作品、分享您的心情！