AnnotationHub中没有本生烟草数据库，只好自己构建，再利用Y叔的clusterProfiler包做富集分析了。
首先参考非模式生物富集分析，进行数据库构建，这里我egg-NOG mapper上传的是本生烟草基因组蛋白序列进行注释，大概半小时就完成了。
注：parse_eggnog.py -O参数我给的是nta，即Nicotiana tabacum (common tobacco)，来过滤掉不靠谱的通路。
不知道自己想要物种的缩写可以去这查：https://www.genome.jp/kegg/catalog/org_list.html
构建完成后，R里面进行分析
library("clusterProfiler")
library("magrittr")
library("tidyverse")
library("RColorBrewer")
my_theme()
# > packageVersion("clusterProfiler")
# [1] ‘4.4.1’
# 读入根据推文步骤构建好的数据
KOannotation <- read.delim("KOannotation.tsv", stringsAsFactors=FALSE)
GOannotation <- read.delim("GOannotation.tsv", stringsAsFactors=FALSE)
GOinfo <- read.delim("go.tb", stringsAsFactors=FALSE)
GOannotation <- split(GOannotation, with(GOannotation, level))
# 读入差异富集基因，我之前做了聚类，所以加了一列簇的信息
kmeansRes <- read.table("degs_cl.txt")
# 由于注释得到的基因名字多了".1"，所以在原始的基因名字最后面加上".1"
kmeansRes$Gene_ID <- paste0(kmeansRes$Gene_ID,".1")
prefix <- 'kmeans8'
savepath <- "D:/your/pwd/enrichment"
for (i in kmeansRes$cl %>% unique) {
## BP
goBP <- enricher(gene = kmeansRes %>% filter(cl == i) %>% .$Gene_ID,
TERM2GENE=GOannotation[['BP']][c(2,1)],
TERM2NAME=GOinfo[1:2])
## check
write.table(as.data.frame(goBP),
paste0(prefix, '_cl', i, '_cp_BP.txt') %>% file.path(savepath, .),
quote = FALSE,
sep = "\t")
## MF
goMF <- enricher(gene = kmeansRes %>% filter(cl == i) %>% .$Gene_ID,
TERM2GENE=GOannotation[['MF']][c(2,1)],
TERM2NAME=GOinfo[1:2])
## check
write.table(as.data.frame(goMF),
paste0(prefix, '_cl', i, '_cp_MF.txt') %>% file.path(savepath, .),
quote = FALSE,
sep = "\t")
## CC
goCC <- enricher(gene = kmeansRes %>% filter(cl == i) %>% .$Gene_ID,
TERM2GENE=GOannotation[['CC']][c(2,1)],
TERM2NAME=GOinfo[1:2])
## check
write.table(as.data.frame(goCC),
paste0(prefix, '_cl', i, '_cp_CC.txt') %>% file.path(savepath, .),
quote = FALSE,
sep = "\t")
## KEGG
kegg <- enricher(gene = kmeansRes %>% filter(cl == i) %>% .$Gene_ID,
TERM2GENE=KOannotation[c(3,1)],
TERM2NAME=KOannotation[c(3,4)])
write.table(as.data.frame(kegg),
paste0(prefix, '_cl', i, '_cp_KEGG.txt') %>% file.path(savepath, .),
quote = FALSE,
sep = "\t")
}
kall <- lapply(kmeansRes$cl %>% unique, function(x) {
eachG <- kmeansRes %>% filter(cl == x) %>% .$Gene_ID
return(eachG)
}) %>%
set_names(kmeansRes$cl %>% unique %>% paste0('cl', .))
下面对所有的簇进行富集分析，这里是非常规的方法，语法和模式生物的不一样（比较见本段末尾）
# goBP_all
goBP_all <- compareCluster(geneCluster = kall,
fun = 'enricher',
TERM2GENE=GOannotation[['BP']][c(2,1)],
TERM2NAME=GOinfo[1:2])
dotplot(goBP_all)
ggsave("goBP_all.pdf", height = 9.5, width = 6)
# goMF_all
goMF_all <- compareCluster(geneCluster = kall,
fun = 'enricher',
TERM2GENE=GOannotation[['MF']][c(2,1)],
TERM2NAME=GOinfo[1:2])
dotplot(goMF_all)
ggsave("goMF_all.pdf", width = 6.5, height = 12)
# goCC_all
goCC_all <- compareCluster(geneCluster = kall,
fun = 'enricher',
TERM2GENE=GOannotation[['CC']][c(2,1)],
TERM2NAME=GOinfo[1:2])
dotplot(goCC_all)
ggsave("goCC_all.pdf", width = 6.5, height = 7)
# kegg all
kallKEGG <- compareCluster(geneCluster = kall,
fun = 'enricher',
TERM2GENE=KOannotation[c(3,1)],
TERM2NAME=KOannotation[c(3,4)],
pvalueCutoff = 0.05)
dotplot(kallKEGG)
ggsave("kallKEGG.pdf", height = 8.8, width = 7)
模式生物语法，以GO的MF（molecular functions）为例
kallGOMF <- compareCluster(geneCluster = kall,
fun = 'enrichGO',
OrgDb = sly.db,
keyType= 'SYMBOL',
ont = "MF",
pAdjustMethod = 'BH',
pvalueCutoff=0.01,
qvalueCutoff=0.1)
kallGOMFSim <- clusterProfiler::simplify(kallGOMF,
cutoff = 0.9,
by = 'p.adjust',
select_fun = min)
dotplot(kallGOMFSim, showCategory = 10)
完成，这里面还可以利用公式直接进行比较富集分析
# Gene_ID就是你的差异表达基因，cl是它所属的簇
kallKEGG <- compareCluster(Gene_ID ~ cl,
data = kmeansRes,
fun = 'enricher',
TERM2GENE=KOannotation[c(3,1)],
TERM2NAME=KOannotation[c(3,4)],
pvalueCutoff = 0.05)
dotplot(kallKEGG)
不用公式
用公式
结果一样。只是展示顺序不同
