1.加载所需R包和相关的向量

library(tidyverse)
library(clusterProfiler)
load(file = 'prepare.rdata')

2.关键差异基因列表获得，通过筛选后将结果传递到pull函数以id为筛选条件进行筛选

gene <- filter(de_result, 
               abs(log2FoldChange) > 1 & padj < 0.05) %>%
        pull(id)

3.TERM与基因对应关系信息

emapper <- read_delim("D:/share/R_data/data/rnaseq-apple/query_seqs.fa.emapper.annotations", 
                      delim = "\t", escape_double = FALSE, 
                      col_names = FALSE, comment = "#", trim_ws = TRUE) %>%
    #导入基因注释表
    dplyr::select(GID = X1, 
                  KO = X9, 
                  Pathway = X10)
    #对导入的基因注释表进行处理，只保留需要的三行
library(stringr)
pathway2gene <- dplyr::select(emapper, Pathway, GID) %>%
    #对emapper进行筛选,只保留GID和Pathway两列,Pathway这列一定要放前面
    separate_rows(Pathway, sep = ',', convert = F) %>%
    #对Pathway这一列进行拆分，拆分的分隔符是‘,’,convert = F表示是否改变数据类型,这里是字符型的所以不需要修改
    filter(str_detect(Pathway, 'ko')) %>%
    #str_detect()函数是stringr包中的函数，这是一个判断函数，即对结果是否进行判断。这一步是要筛选出Pathway这一列中的'ko'字符而过滤掉其他字符
    mutate(Pathway = str_remove(Pathway, 'ko'))
    #str_remove()函数剔除Pathway中所有ko字符

4.TERM与基因对应关系信息

library(magrittr)
get_path2name <- function(){
  keggpathid2name.df <- clusterProfiler:::kegg_list("pathway")
  keggpathid2name.df[,1] %<>% gsub("path:map", "", .)
  colnames(keggpathid2name.df) <- c("path_id","path_name")
  return(keggpathid2name.df)
}
pathway2name <- get_path2name()
#pathway2name <- mutate(pathway2name, path_id = str_remove(path_id, 'map'))

直接调用clusterProfiler中的包进行分析，会直接联网下载每一个pathway的名字

5.使用enricher()函数进行富集分析

library(clusterProfiler)
de_ekp <- enricher(gene,
                   TERM2GENE = pathway2gene,
                   TERM2NAME = pathway2name,
                   pvalueCutoff = 0.05,
                   qvalueCutoff = 0.05)
de_ekp_df <- as.data.frame(de_ekp)

enricher()函数是一个通用富集分析方法，可做pathway分析、GO富集分析、染色体富集分析。de_ekp是一个复杂数据，可用class()函数查看，使用as.data.frame()将其强制转换为dataframe。

6.所有基因的FC信息的获得与Pathway富集分析

geneList <- de_result$log2FoldChange
##先做gene_List列表，来自于de_result中的log2FoldChange
names(geneList) <- de_result$id
##给geneList中的vector取名字，名字来源于de_result中的id列
geneList <- sort(geneList, decreasing = T)
##对geneList由高到低进行排序
barplot(de_ekp, showCategory = 30)
##作Pathway富集分析柱状图
dotplot(de_ekp, showCategory = 10)
##作Pathway富集分析柱状图
cnetplot(de_ekp,
         FoldChange = geneList,
         #构建好的基因差异倍数  
         showCategory = 10,
         #展示前三个Pathway
         node_lable = "all",#category/gene/all/none
         #node_lable参数在老版本中应该是'T'/'F',新版本中可以进行指定,即显示标签的名称
         circular = TRUE,
         #将图片是否显示为圈图
         color.params = list(edge = TRUE))
##作Pathway富集分析圈图（推荐方式）    
emapplot(de_ekp, showCategory = 7, pie = 'count')    
##是圈图的简化版，省去了基因，只留下pathway,有关联的基因还是会通过线条连接

7.GO富集分析

##加载有关包
library(argparser, quietly=TRUE)
library(tidyverse)

##创建了一个参数解析器
p <- arg_parser("make OrgDB from emapper")

##向解析器中添加参数
p <- add_argument(p, "annotation", help="emapper annotation result", type="character")

##解析命令行参数，将参数值存储在argv变量中
argv <- list(annotation = "D:/bioinfotools/R/R/R-4.3.1/R-4.3.1/wkdir/rna-seq/query_seqs.fa.emapper.annotations")

##注释结果文件中读取数据并进行处理
emapper <- read_delim(argv$annotation, 
                      "\t", escape_double = FALSE, col_names = FALSE, 
                      comment = "#", trim_ws = TRUE) %>%
  dplyr::select(GID = X1, 
                Gene_Symbol = X6, 
                GO = X7, 
                KO = X9, 
                Pathway = X10, 
                OG = X21, 
                Gene_Name = X22)

##处理 emapper 数据框，提取基因信息并存储在gene_info变量中
gene_info <- dplyr::select(emapper,  GID, Gene_Name) %>%
  dplyr::filter(!is.na(Gene_Name))

##处理 emapper 数据框，提取基因与GO术语的关系，并存储在gene2go变量中
gene2go <- dplyr::select(emapper, GID, GO) %>%
  separate_rows(GO, sep = ',', convert = F) %>%
  filter(!is.na(GO)) %>%
  mutate(EVIDENCE = 'IEA') 

##使用 AnnotationForge 包的函数 makeOrgPackage 创建一个OrgDB包
AnnotationForge::makeOrgPackage(gene_info=gene_info,
                                go=gene2go,
                                maintainer='zhangsan <zhangsan@genek.tv>',
                                author='zhangsan',
                                outputDir="./",
                                tax_id=0000,
                                genus='M',
                                species='y',
                                goTable="go",
                                version="1.0")

##行使用 pkgbuild 包的函数 build 构建OrgDB包
pkgbuild::build('.//org.My.eg.db', dest_path = ".")

8.进行GO富集分析

安装org.My.eg.db_1.0软件包
install.packages('D:/bioinfotools/R/R/R-4.3.1/R-4.3.1/wkdir/rna-seq/org.My.eg.db_1.0.tar.gz',
                 repos = NULL,#从本地安装软件包
                 lib = 'D:/bioinfotools/R/R/R-4.3.1/R-4.3.1/R_lib')#安装文件夹

##加载软件包
library(org.My.eg.db, lib = 'D:/bioinfotools/R/R/R-4.3.1/R-4.3.1/R_lib')

##富集分析
#library(clusterProfiler)
de_ego <- enrichGO(gene = gene,
                  OrgDb = org.My.eg.db,
                  keyType = 'GID',
                  ont = 'ALL',
                  pvalueCutoff = 0.05,
                  qvalueCutoff = 0.05)

##数据转换(同pathway富集分析)
de_ego_df <- as.data.frame(de_ego)
head(de_ego_df)

9.作富集分析柱状图

library(ggplot2)
dotplot(de_ego, showCategory = 10, split = "ONTOLOGY") +
  facet_grid(ONTOLOGY~., scale="free")