#这里需要对非模式物种制作ORG.DB包，如果是模式物种，“https://bioconductor.org/packages/release/BiocViews.html#___OrgDb”该网站有自带的成熟的包，自行下载使用就行。
#对上一个教程中得到的out.emapper.annotations文件，对表头修整下：
#windows上的R运行
library(dplyr)
library(stringr)
library(clusterProfiler)
library(AnnotationForge)
library(tidyr)
options(stringsAsFactors = F)#keepcharacter but not factor conversion
emapper <-read.delim("out.emapper.annotations")
emapper[emapper=="-"] <-NA#change "-" to "NA"
emapper <-emapper[-(49584:49586),]#remove the final 3 rows
gene_info <- dplyr::select(emapper,GID=query, Gene_Name=seed_ortholog)%>%
?dplyr::filter(!is.na(Gene_Name))
#gene_info表格
#提取GO信息
gene2go <- dplyr::select(emapper,GID=query,GO=GOs)%>%
?filter(!is.na(GO))%>%
?mutate(EVIDENCE='IEA')%>%
?separate_rows(GO, sep = ',', convert = F)
#gene2go表格，其实和实战8中，TBTOOLS做出来的是一样的。
#构建orgdb包
AnnotationForge::makeOrgPackage(gene_info=gene_info,
??????????????????????????????? go=gene2go,
??????????????????????????????? maintainer ='LJH',
??????????????????????????????? author = 'LJH',
??????????????????????????????? outputDir ="./",
??????????????????????????????? tax_id = 0000,
??????????????????????????????? genus = 'P',
??????????????????????????????? species ='tri',
??????????????????????????????? goTable ="go",
??????????????????????????????? version ="1.0")
#对新生成的org.Ptri.eg.db包中的DESCRIPTION，进行修改，Maintainer: LJH <abc@cba.com>，
#打包
pkgbuild::build('./org.Ptri.eg.db',dest_path = './')
#生成org.Ptri.eg.db_1.0.tar.gz将这个R包放到平时R包安装的路径中，D:\\R-4.3.1\\library，本地安装
install.packages('your_path',repos= NULL)
library(org.Ptri.eg.db)
#将实战5中的差异基因自行excel修改下基因名，使其与gene_info中的GID相对应。
#差异分析
DE <-read.delim("DE_genes_filter.txt")
ego <- enrichGO(gene = DE$GID,
??????????????? OrgDb = org.Ptri.eg.db,
??????????????? keyType = 'GID',
??????????????? ont = 'ALL',
??????????????? pvalueCutoff = 0.05,
??????????????? qvalueCutoff = 0.05)
#以下是自带的clusterprofiler的画图函数
dotplot(ego)
barplot(ego)
cnetplot(ego)
#这个富集文件要自己用EXCEL修改，我自己选了15条BP-4条CC-15条MF。GeneRatio自己做成百分比。#以下是用GGPLOT2画条形图了，各种函数，自己调节参数即可。
write.table(ego,file = "Ptri_GO_test",sep = '\t',quote = F)
ego2 <-read.delim("Ptri_GO_test")
library(ggplot2)
library(GOplot)
ggplot(ego2, aes(Description,-log10(p.adjust))) +
?geom_col(aes(fill = ONTOLOGY), width = 0.5, show.legend = FALSE) +
?scale_fill_manual(values = c('#D06660', '#5AAD36', '#6C85F5')) +
?facet_grid(ONTOLOGY~., scale = 'free_y', space = 'free_y') +
?theme(panel.grid = element_blank(), panel.background =element_rect(color = 'black', fill = 'transparent')) +
?scale_y_continuous(expand = expansion(mult = c(0, 0.1))) +
?coord_flip() +
?labs(x = '', y = '-log10(p.adjust)')
#GGPLOT2画气泡图
pp <- ggplot(ego2, aes(GeneRatio,Description))
pp + geom_point() +
?geom_point(aes(size = Count)) +
?geom_point(aes(size = Count, color = -1 * log10(qvalue))) +
?scale_colour_gradient(low = "green", high = "red") +
?labs(color = expression(-log[10](Qvalue)), size = "GeneNumber", x = "Rich Factor", y = "Pathway Name", title= "Top 30 of Pathway Enrichment") +
?theme_bw()
##########图片自己微调吧#######
#小林家的龙女仆
