23-有参转录组实战9-差异基因GO富集分析

#这里需要对非模式物种制作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)#keep character 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 <[email protected]>，

#打包

pkgbuild::build('./org.Ptri.eg.db', dest_path = './')

#生成org.Ptri.eg.db_1.0.tar.gz将这个R包放到平时R包安装的路径中，本地安装

install.packages('your_path', repos = NULL)
library(org.Ptri.eg.db)

#将实战5中的差异基因自行excel修改下基因名，使其与gene_info中的GID相对应，我加尾缀.1吧，先看看效果如何，只做一列就行，加表头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 = "Gene Number", x = "Rich Factor", y = "Pathway Name", title = "Top 30 of Pathway Enrichment") +
  theme_bw()
##########图片自己微调吧#######

#小林家的龙女仆