一句话：我们经常会面对非模式物种的GO或者KEGG富集与注释。
1、载入我们所需要的包
if("clusterProfiler" %in% rownames(installed.packages()) == FALSE) {source("http://bioconductor.org/biocLite.R");biocLite("clusterProfiler")}
biocLite("colorspace")
suppressMessages(library(clusterProfiler))
ls("package:clusterProfiler")
library(AnnotationHub)
library(biomaRt)
library(enrichGO)
2、收索目标物种
hub <- AnnotationHub() unique(hub$dataprovider) #查看数据注释来源 query(hub, "Apis cerana")
3、抓取目标OrgDb
Apis_cerana.OrgDb <- hub[["AH62635"]] columns(Apis_cerana.OrgDb) Apis_cerana.OrgDb
4、了解这个Apis_cerana.OrgDb对象
4.1、以symbol形式展示
head(keys(Apis_cerana.OrgDb,keytype = "SYMBOL"),30)
4.2、可以查看gene类型
keytypes(Apis_cerana.OrgDb)
4.3、默认基因是以ENTREZID，例如查看前90个基因。
keys(Apis_cerana.OrgDb)[1:50]
4.4、格式转换。例如将上面前50个基因转换为SYMBOL格式对照
bitr(keys(Apis_cerana.OrgDb)[1:50], \'ENTREZID\', "SYMBOL", Apis_cerana.OrgDb)
4.5、多个格式转化，例如将上面的50个基因转换为多个形式
bitr(keys(Apis_cerana.OrgDb), \'ENTREZID\', c("SYMBOL","REFSEQ", "GO", "ONTOLOGY"), Apis_cerana.OrgDb)
5、测试
5.1、取前100个基因集做测试
sample_genes <- keys(Apis_cerana.OrgDb,keytype = "ENTREZID")[1:100]
5.2、做GO富集分析
enrich.go = enrichGO( gene = sample_genes,
OrgDb = Apis_cerana.OrgDb,
keyType = \'SYMBOL\',
ont= "ALL",
pAdjustMethod = \'BH\',
pvalueCutoff = 0.05,
qvalueCutoff = 0.1,
readable = T)
head(enrich.go)dim(hainan_enrich.go)
dim(hainan_enrich.go[hainan_enrich.go$ONTOLOGY==\'BP\',])
dim(hainan_enrich.go[hainan_enrich.go$ONTOLOGY==\'CC\',])
dim(hainan_enrich.go[hainan_enrich.go$ONTOLOGY==\'MF\',])
画图
barplot(hainan_enrich.go,showCategory=20,drop=T) dotplot(hainan_enrich.go,showCategory=50) write.table(as.data.frame(hainan_enrich.go), file="hainan_enrich_go.xls")
5.3、画GO分类图
5.3.1、重新取样，然后做GO富集分析
rm(list=ls()) #重头开始
sample_genes <- keys(Apis_cerana.OrgDb,keytype = "ENTREZID")[1:100] #采样100个基因
ego_ALL = enrichGO( gene =sample_genes,
OrgDb = Apis_cerana.OrgDb,
keyType = \'ENTREZID\',
ont= "ALL",
pAdjustMethod = \'BH\',
pvalueCutoff =1,
qvalueCutoff =1,
readable = T)
5.3.2、对富集结果进行分类
ego_result_BP <- ego_ALL[ego_ALL$ONTOLOGY==\'BP\',] dim(ego_result_BP) ego_result_CC <- ego_ALL[ego_ALL$ONTOLOGY==\'CC\',] dim(ego_result_CC) ego_result_MF <- ego_ALL[ego_ALL$ONTOLOGY==\'MF\',] dim(ego_result_MF)
构建数据框
display_number = c(48, 26, 39)
go_enrich_df <- data.frame(ID = c(ego_result_BP$ID, ego_result_CC$ID, ego_result_MF$ID),
Description = c(ego_result_BP$Description, ego_result_CC$Description, ego_result_MF$Description),
GeneNumber = c(ego_result_BP$Count, ego_result_CC$Count, ego_result_MF$Count),
type = factor(c(rep("biological process", display_number[1]), rep("cellular component", display_number[2]), rep("molecular function", display_number[3])), levels=c("molecular function", "cellular component", "biological process")))
缩短名称
## numbers as data on x axis
go_enrich_df$number <- factor(rev(1:nrow(go_enrich_df)))
#shorten the names of GO terms
shorten_names <- function(x, n_word=4, n_char=40){
if (length(strsplit(x, " ")[[1]]) > n_word || (nchar(x) > 40))
{
if (nchar(x) > 40) x <- substr(x, 1, 40)
x <- paste(paste(strsplit(x, " ")[[1]][1:min(length(strsplit(x," ")[[1]]), n_word)], collapse=" "), "...", sep="")
return(x) }
else { return(x) }
}
设置画图轴
labels=(sapply(levels(go_enrich_df$Description)[as.numeric(go_enrich_df$Description)], shorten_names)) names(labels) = rev(1:nrow(go_enrich_df))
画图
library(ggplot2)
CPCOLS <- c("#8DA1CB", "#FD8D62", "#66C3A5")
library(ggplot2)
CPCOLS <- c("#8DA1CB", "#FD8D62", "#66C3A5")
p <- ggplot(data=go_enrich_df, aes(x=number, y=GeneNumber, fill=type)) +
geom_bar(stat="identity", width=0.8) +
coord_flip() +
scale_fill_manual(values = CPCOLS) +
theme_bw() +
scale_x_discrete(labels=labels) +
xlab("GO term") +
theme(axis.text=element_text(face = "bold", color="gray50")) +
labs(title = "ALL Enriched GO Terms")
p
pdf("go_enrichment_of_sample_targets.pdf")
dev.off()
svg("go_enrichment_of_sample_targets.svg")
p
dev.off()
参考文献
https://blog.csdn.net/sinat_30623997/article/details/79250940https://blog.csdn.net/sinat_30623997/article/details/79250940