library(tidyverse)
library(clusterProfiler)
de_result <- read.table(
file = 'flower_sex/genes.counts.matrix.F_vs_M.DESeq2.DE_results',
header = TRUE) %>%
select(id,log2FoldChange,pvalue,padj)
这是de_result数据
library(readr)
emapper<- read_delim(
file = 'D:/geneomes/',
"\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 = X8)
emapper数据
### 候选基因向量(gene)
一个字符型向量,包含感兴趣的基因,如差异表达基因、WGCNA 得到的关键模块里的基因、PCA 中某个PC loadings 高的基因。
```{r}
gene <- filter(de_result,
abs(log2FoldChange) > 2 & padj < 0.01) %>%
pull(id)
```
### 基因差异向量(geneList)
一个有名向量,记录差异大小,名字为基因ID,值为 logFC,按 logFC 降序排列。
```{r}
geneList <- de_result$log2FoldChange
names(geneList) <- de_result$id
geneList <- sort(geneList, decreasing = T)
# 提取基因信息
gene_info <- dplyr::select(emapper, GID, Gene_Name) %>%
dplyr::filter(!is.na(Gene_Name))
# 提取GO信息
gene2go <- dplyr::select(emapper, GID, GO) %>%
separate_rows(GO, sep = ',', convert = F) %>%
filter(!is.na(GO)) %>%
mutate(EVIDENCE = 'IEA')
#过滤掉任何不以GO开头的GO
gene2go <- gene2go[grepl("^GO:", as.character(gene2go$GO)),]
# 构建 OrgDB
library(AnnotationForge)
AnnotationForge::makeOrgPackage(gene_info=gene_info,
go=gene2go,
maintainer='w<w@101214581.com>',
author='bank',
outputDir="./",
tax_id=0000,
genus='m',
species='y',
goTable="go",
version="1.0")
# 打包
pkgbuild::build('./org.My.eg.db', dest_path = "./enrich/")
```
安装自己构建的 OrgDb
```{r}
# 创建一个文件夹
dir.create('R_Library', recursive = T)
# 将包安装在该文件夹下
install.packages('enrich/org.My.eg.db_1.0.tar.gz',
repos = NULL, #从本地安装
lib = 'R_Library') # 安装文件夹
# 加载 OrgDB
library(org.My.eg.db, lib = 'R_Library')
### 富集分析
de_ego <- enrichGO(gene = gene,
OrgDb = org.My.eg.db,
keyType = 'GID',
ont = 'BP',
qvalueCutoff = 0.05,
pvalueCutoff = 0.05)
提示:
No gene can be mapped....
--> Expected input gene ID: evm.model.Chr6.25133,evm.model.Chr1.961.3.63ebc806,evm.model.Chr5.6899,evm.model.Chr6.1573,evm.model.Chr5.6420,evm.model.Chr3.2040
--> return NULL...