KEGG富集.txt

#########################加载所需软件包####################################
#BiocManager::install('clusterProfiler')                    #安装一次就行 

#每次都加载以下包
library(tidyverse)                                                        #                   
library(clusterProfiler)                                                  #
library(dplyr)                                                            #
###########################################################################
#读取eggnog-mapper注释结果

eggnogmapper<- read_delim(
  file = 'Iyun.emapper.annotations',
  "\t",
  escape_double =FALSE,
  col_names=FALSE,
  comment = "#",trim_ws = TRUE)%>%
  dplyr::select(GID=X1,
                Gene_Symbol=X9,
                GO=X10,
                KO=X12,
                Pathway=X13,
                OG=X7,
                Gene_Name=X21)




#################构建orgdb库


# 提取基因信息
gene_info <- dplyr::select(eggnogmapper,  GID, Gene_Name) %>%
  dplyr::filter(!is.na(Gene_Name))

# 提取GO信息
gene2go <- dplyr::select(eggnogmapper, GID, GO) %>%
  separate_rows(GO, sep = ',', convert = F) %>%
  filter(!is.na(GO)) %>%
  mutate(EVIDENCE = 'IEA')
gene2go <- subset(gene2go,grepl("^.*GO.*$",GO))


# 构建 OrgDB,参数十分玄学，请按照以下格式
library(AnnotationForge)
AnnotationForge::makeOrgPackage(gene_info=gene_info,
                                go=gene2go,
                                maintainer='zhangsan <zhangsan@genek.tv>',
                                author='zhangsan',
                                outputDir="/home/lx_sky6/yyj/R/KEGG/KEGG/myresearch",
                                tax_id=0000,
                                genus='I',                     #只改动这里
                                species='yun',               #和这里，进行命名
                                goTable="go",
                                version="1.0")

pkgbuild::build('/home/lx_sky6/yyj/R/KEGG/KEGG/myresearch/org.Iyun.eg.db', dest_path = "/home/lx_sky6/yyj/R/KEGG/KEGG/myresearch")       #这两个路径自己创建





#安装生成的压缩包

install.packages('/home/lx_sky6/yyj/R/KEGG/KEGG/myresearch/org.Iyun.eg.db_1.0.tar.gz',
                 repos = NULL,
                 lib='/home/lx_sky6/yyj/R/KEGG/KEGG/R_Library')











#--------------------------------------分割线，以上orgdb包的构建————————————————————————————————————



#加载org.Wtaiwan.eg.db
library(org.Iyun.eg.db,lib='/home/lx_sky6/yyj/R/KEGG/KEGG/R_Library')

# 导入基因ID
library(readr)
shui_gene <- read_table2("~/yyj/genome/selection/3.1/9.3.1.id")

#degene
de_gene <- pull(shui_gene, GID)



########### ############ ################   run__GO富集

de_ego <- enrichGO(gene = de_gene,
                   OrgDb = org.Iyun.eg.db,
                   keyType = 'GID',
                   ont = 'ALL',
                   qvalueCutoff = 0.05,
                   pvalueCutoff = 0.05)

#以上得到的de_ego结果可视化放到de_ego_df
de_ego_df <- as.data.frame(de_ego)

#取前多少个？？？head(de_ego_df)


######################################可视化

#绘制点状图（内容同上）
enrichplot::dotplot(de_ego,showCategory = 20)      #数字设置展示点的个数











########### ############ ##############      run__KEGG富集 

#根据eggnogmapper列表获取pathway2gene
pathway2gene <- dplyr::select(eggnogmapper, Pathway, GID) %>%
  separate_rows(Pathway, sep = ',', convert = F) %>%
  filter(str_detect(Pathway, 'ko')) %>%
  mutate(Pathway = str_remove(Pathway, 'ko'))


#通过以下代码（不用做改变）从线上获取生成pathway2name
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()




#kegg富集分析
library(clusterProfiler)
de_ekp <- enricher(de_gene,
                   TERM2GENE = pathway2gene,
                   TERM2NAME = pathway2name,
                   pvalueCutoff = 0.05,
                   qvalueCutoff = 0.05)
#可视化
de_ekp_df <- as.data.frame(de_ekp)

#head(de_ekp_df)


#点状图
enrichplot::dotplot(de_ekp, showCategory = 20)







#导出富集的结果文件
write.table(de_ego_df,file = '/home/lx_sky6/yyj/genome/selection/3.1/3.1_GO.csv',sep = ',',quote = FALSE)

write.table(de_ekp_df,file = 'home/lx_sky6/yyj/genome/selection/3.1/3.1_Kegg.csv',sep = ',',quote = FALSE)