自己构建物种包OrgDb，然后用clusterProfiler富
2019-09-25 本文已影响0人
生信小白2018
利用clusterProfiler进行富集分析时，发现网上没有物种包OrgDb，利用网上教程构建了一个。主要参考了刘小泽的简书，特别感谢一下！！！
emapper首先得到注释文件，前提是自己安装好eggnog-mapper并且下载好相应的数据库
emapper.py --cpu 20 -i filter.pep.fa_16 --output filter.pep.fa_16.out -d virNOG -m diamond
or
emapper.py -m diamond \
-i sesame.fa \
-o diamond \
--cpu 19
得到注释文件后进行处理，只保留表头query_name这一行的注释信息，去掉头尾的# 等信息
sed -i '/^# /d' diamond.emapper.annotations
sed -i 's/#//' diamond.emapper.annotations
开始一定要加上options(stringsAsFactors = F)，否则所有的字符在数据框中都会被R默认设置为factor！！！
library(clusterProfiler)
library(dplyr)
library(stringr)
options(stringsAsFactors = F)
#hub <- AnnotationHub::AnnotationHub()
#query(hub,"Theobroma cacao")
#STEP1：自己构建的话，首先
需要读入文件
setwd("G:/xxx")
egg_f <- "xx_annotations1"
egg <- read.table(egg_f, header=TRUE, sep = "\t")
#egg_f <- "test111"
#egg <- read.table(egg_f, header=TRUE, sep = "\t")
#gene_info_new<-read.table("id.txt" , header=TRUE, sep = "\t")
egg[egg==""]<-NA #这个代码来自花花的指导(将空行变成NA，方便下面的去除)
#STEP2: 从文件中挑出基因query_name与eggnog注释信息
gene_info <- egg %>%
dplyr::select(GID = query_name, GENENAME = eggNOG_annot) %>% na.omit()
#STEP3-1：挑出query_name与GO注释信息
gterms <- egg %>%
dplyr::select(query_name, GO_terms) %>% na.omit()
#STEP3-2：我们想得到query_name与GO号的对应信息
# 先构建一个空的数据框(弄好大体的架构，表示其中要有GID =》query_name，GO =》GO号， EVIDENCE =》默认IDA)
# 关于IEA：就是一个标准，除了这个标准以外还有许多。IEA就是表示我们的注释是自动注释，无需人工检查http://wiki.geneontology.org/index.php/Inferred_from_Electronic_Annotation_(IEA)
# 两种情况下需要用IEA：1. manually constructed mappings between external classification systems and GO terms； 2.automatic transfer of annotation to orthologous gene products.
gene2go <- data.frame(GID = character(),
GO = character(),
EVIDENCE = character())
# 然后向其中填充：注意到有的query_name对应多个GO，因此我们以GO号为标准，每一行只能有一个GO号，但query_name和Evidence可以重复
for (row in 1:nrow(gterms)) {
gene_terms <- str_split(gterms[row,"GO_terms"], ",", simplify = FALSE)[[1]]
gene_id <- gterms[row, "query_name"][[1]]
tmp <- data_frame(GID = rep(gene_id, length(gene_terms)),
GO = gene_terms,
EVIDENCE = rep("IEA", length(gene_terms)))
gene2go <- rbind(gene2go, tmp)
}
#STEP4-1: 挑出query_name与KEGG注释信息
gene2ko <- egg %>%
dplyr::select(GID = query_name, KO = KEGG_KOs) %>%
na.omit()
#STEP4-2: 得到pathway2name, ko2pathway
# 需要下载 json文件(这是是经常更新的)
# https://www.genome.jp/kegg-bin/get_htext?ko00001
# 代码来自：http://www.genek.tv/course/225/task/4861/show
if(F){
# 需要下载 json文件(这是是经常更新的)
# https://www.genome.jp/kegg-bin/get_htext?ko00001
# 代码来自：http://www.genek.tv/course/225/task/4861/show
library(jsonlite)
library(purrr)
library(RCurl)
update_kegg <- function(json = "ko00001.json") {
pathway2name <- tibble(Pathway = character(), Name = character())
ko2pathway <- tibble(Ko = character(), Pathway = character())
kegg <- fromJSON(json)
for (a in seq_along(kegg[["children"]][["children"]])) {
A <- kegg[["children"]][["name"]][[a]]
for (b in seq_along(kegg[["children"]][["children"]][[a]][["children"]])) {
B <- kegg[["children"]][["children"]][[a]][["name"]][[b]]
for (c in seq_along(kegg[["children"]][["children"]][[a]][["children"]][[b]][["children"]])) {
pathway_info <- kegg[["children"]][["children"]][[a]][["children"]][[b]][["name"]][[c]]
pathway_id <- str_match(pathway_info, "ko[0-9]{5}")[1]
pathway_name <- str_replace(pathway_info, " \\[PATH:ko[0-9]{5}\\]", "") %>% str_replace("[0-9]{5} ", "")
pathway2name <- rbind(pathway2name, tibble(Pathway = pathway_id, Name = pathway_name))
kos_info <- kegg[["children"]][["children"]][[a]][["children"]][[b]][["children"]][[c]][["name"]]
kos <- str_match(kos_info, "K[0-9]*")[,1]
ko2pathway <- rbind(ko2pathway, tibble(Ko = kos, Pathway = rep(pathway_id, length(kos))))
}
}
}
save(pathway2name, ko2pathway, file = "kegg_info.RData")
}
update_kegg(json = "ko00001.json")
}
#STEP5: 利用GO将gene与pathway联系起来，然后挑出query_name与pathway注释信息
load(file = "kegg_info.RData")
gene2pathway <- gene2ko %>% left_join(ko2pathway, by = "KO") %>%
dplyr::select(GID, Pathway) %>%
na.omit()
library(AnnotationForge)
#STEP6： 制作自己的Orgdb
# 查询物种的Taxonomy，例如要查sesame
# https://www.ncbi.nlm.nih.gov/taxonomy/?term= sesame
tax_id = "4182"
genus = "Sesamum"
species = "indicum"
#gene2go <- unique(gene2go)
#gene2go<-gene2go[!duplicated(gene2go),]
#gene2ko<-gene2ko[!duplicated(gene2ko),]
#gene2pathway<-gene2pathway[!duplicated(gene2pathway),]
makeOrgPackage(gene_info=gene_info,
go=gene2go,
ko=gene2ko,
maintainer = "xxx <xxx@163.com>",
author = "",
pathway=gene2pathway,
version="0.0.1",
outputDir = ".",
tax_id=tax_id,
genus=genus,
species=species,
goTable="go")
ricenew.orgdb <- str_c("org.", str_to_upper(str_sub(genus, 1, 1)) , species, ".eg.db", sep = "")
options(stringsAsFactors = F) 的作用
如果说一个data.frame中的元素是factor，你想转化成numeric，你会怎么做？比如d[1,1]是factor
正确答案是 先as.character(x) 再as.numeric(x)
哈哈，我刚发现如果直接as.numeric,就不是以前的数字了，坑爹啊。
原来as.data.frame()有一个参数stringsAsFactors
如果stringAsFactor=F
就不会把字符转换为factor 这样以来，原来看起来是数字变成了character，原来是character的还是character
KEGG富集分析
setwd("G:/xxx")
library(purrr)
library(tidyverse)
library(clusterProfiler)
################################################
# 导入自己构建的 OrgDb
################################################
install.packages("org.xxx.db", repos=NULL, type="sources")
library(org.xxx.db)
columns(org.xxx.db)
# 导入需要进行富集分析的基因列表，并转换为向量
#########################################################################################
DD<-"DEGs"
DEGs<- read.table(DD, header=TRUE, sep = "\t")
gene_list <- DEGs[,1]
################################################
# 从 OrgDB 提取 Pathway 和基因的对应关系
################################################
pathway2gene <- AnnotationDbi::select(org.xxx.db,
keys = keys(org.xxx.db),
columns = c("Pathway","KO")) %>%
na.omit() %>%
dplyr::select(Pathway, GID)
################################################
# 导入 Pathway 与名称对应关系
################################################
load("kegg_info.RData")
#KEGG pathway 富集
ekp <- enricher(gene_list,
TERM2GENE = pathway2gene,
TERM2NAME = pathway2name,
pvalueCutoff = 1,
qvalueCutoff = 1,
pAdjustMethod = "BH",
minGSSize = 1)
ekp_results <- as.data.frame(ekp)
barplot(ekp, showCategory=20,color="pvalue",
font.size=10)
dotplot(ekp)
emapplot(ekp)
GO富集分析
#########################################################################################
GO 分析
#########################################################################################
ego <- enrichGO(gene = gene_list, #差异基因 vector
keyType = "GID", #差异基因的 ID 类型，需要是 OrgDb 支持的
OrgDb = org.xxx.db, #对应的OrgDb
ont = "CC", #GO 分类名称，CC BP MF
pvalueCutoff = 1, #Pvalue 阈值 （pvalue=1指输出所有结果，pvalue=0.05指输出符合要求的结果）
qvalueCutoff = 1, #qvalue 阈值 pAdjustMethod = "BH", #Pvalue 矫正方法
readable = FALSE) #TRUE 则展示SYMBOL，FALSE 则展示原来的ID（选false是因为不是所有gene都有symbol的)
ego_results<-as.data.frame(ego) ###生成的ego文件转换成data.frame格式即可。
write.table(ego_results, file = "ego_results.txt", quote = F) ###让保存的字符串不用“”引起来
pdf(file = "ego_barplot.pdf") ##打开一个PDF文件
barplot(ego, showCategory=20, x = "GeneRatio") ##把图画到这个PDF文件里
dev.off() ##关闭PDF
dotplot(ego)
emapplot(ego)
参考
https://www.jianshu.com/p/9c9e97167377
https://www.jianshu.com/p/5d5394e0774f
