非模式生物本生烟草(Nicotiana benthamiana)

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)

结果一样。只是展示顺序不同