OrgDb才能进行富集分析
2021-06-09 本文已影响0人
MLD_TRNA
要进行GO或者KEGG富集分析，就需要知道每个基因对应什么样的GO/KEGG分类，OrgDb就是存储不同数据库基因ID之间对应关系，以及基因与GO等注释的对应关系的 R 软件包
如果自己研究的物种不在
http://bioconductor.org/packages/release/BiocViews.html#___OrgDb之列，很大可能就需要自己构建OrgDb，然后用clusterProfiler分析所以本次的内容都是基于非模式生物
分类
非模式生物要想找到自己的注释包，又分成两类：
- 一类是在AnnotationHub（https://bioconductor.org/packages/release/bioc/html/AnnotationHub.html）中存在的，例如棉铃虫
- 另一类是在AnnotationHub也不存在相应物种，就需要用AnnotationForge（ https://bioconductor.org/packages/release/bioc/html/AnnotationForge.html ）来自己构建
第一类：利用AnnotationHub得到org.db
下面以棉铃虫为例
首先下载并加载AnnotationHub
options("repos"= c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
options(BioC_mirror="http://mirrors.ustc.edu.cn/bioc/")
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("AnnotationHub", version = "3.8")
library(AnnotationHub)
然后加载现有的物种database
hub <- AnnotationHub() #这一步受限于网速，不成功的话多试几次
# 调用图形界面查看物种
display(hub)
# 或者根据物种拉丁文名称查找
query(hub,"helicoverpa")
# 'object[["AH66950"]]'
title
AH66950 | org.Helicoverpa_armigera.eg.sqlite
AH66951 | org.Heliothis_(Helicoverpa)_armigera.eg....
# 这里AH66950是我们需要的
# 然后下载这个sqlite数据库
ha.db <- hub[['AH66950']]
#查看前几个基因（Entrez命名）
head(keys(ha.db))
#查看包含的基因数
length(keys(ha.db))
#查看包含多少种ID
columns(ha.db)
#查看前几个基因的ID
select(ha.db, keys(ha.db)[1:3],
c("REFSEQ", "SYMBOL"), #想获取的ID
"ENTREZID")
#得到结果
ENTREZID REFSEQ SYMBOL
1 9977712 YP_004021052.1 COX1
2 9977713 YP_004021053.1 COX2
3 9977714 YP_004021054.1 ATP8
#保存到文件
saveDb(ha.db, "Harms-AH66950.sqlite")
#之后再使用直接加载进来
harms.db <- loadDb("Harms-AH66950.sqlite")
image
参考：https://mp.weixin.qq.com/s/lHKZtzpN2k9uPN7e6HjH3w
第二类：利用AnnotationForge得到org.db
通过上次的推送（功能注释整合流程），我们知道了怎样利用eggnog-mapper去人工构建注释、富集桥梁
上次选取的是一个病毒的小例子，这次可以用芝麻（Sesame）做演示
先下载蛋白序列
wget http://www.sesame-bioinfo.org/SesameFG/BLAST_search/G608_contig_2014-08-29.FgeneSH.pep.rar
# 解压后上传
因为统计了下有38406条序列，因此使用diamond来进行功能注释
# 前提是自己安装好eggnog-mapper并且下载好相应的数据库
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
image
关于结果解释：https://github.com/jhcepas/eggnog-mapper/wiki/Results-Interpretation
其中关于COG的介绍：倒数第二列
COG functional categories: COG functional category inferred from best matching OG 【会给出一个大写字母，每一个大写字母都有自己的解释：COG_explanation】
STEP1：自己构建的话，首先需要读入文件
egg_f <- "diamond.emapper.annotations"
egg <- read.csv(egg_f, 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
这一步不需要管代码什么意思，只需要知道它可以帮我们得到以上两个文件就好
options(stringsAsFactors = F)
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()
STEP6： 制作自己的Orgdb
# 查询物种的Taxonomy，例如要查sesame
# https://www.ncbi.nlm.nih.gov/taxonomy/?term=sesame
tax_id = "4182"
genus = "Sesamum"
species = "indicum"
makeOrgPackage(gene_info=gene_info,
go=gene2go,
ko=gene2ko,
pathway=gene2pathway,
version="0.0.1",
outputDir = ".",
tax_id=tax_id,
genus=genus,
species=species,
goTable="go")
sesame.orgdb <- str_c("org.", str_to_upper(str_sub(genus, 1, 1)) , species, ".eg.db", sep = "")
有了Orgdb就可以做富集分析了
# enrichGO最主要的目的就是将基因编号转换成GO号
ego <- enrichGO(gene = id.fc$ENTREZID,
#模式物种
#OrgDb = org.Mm.eg.db,
#非模式物种，例如芝麻
OrgDb = sesame.orgdb,
ont = "BP", #或MF或CC
pAdjustMethod = "BH",
#pvalueCutoff = 0.01,
qvalueCutoff = 0.01)
# 同理也能做enrichKEGG
添加于2019.4.8，鉴于很多小伙伴都在问如何使用自建数据库进行KEGG分析
########################################
# KEGG
########################################
if(F){
library(purrr)
library(tidyverse)
library(clusterProfiler)
########################################################################################
################################################
# 导入自己构建的 OrgDb
################################################
library(org.Harmigera.eg.db)
columns(org.Harmigera.eg.db)
########################################################################################
# 导入需要进行富集分析的基因列表，并转换为向量
#########################################################################################
gene_list <- DEGs[,1]
################################################
# 从 OrgDB 提取 Pathway 和基因的对应关系
################################################
pathway2gene <- AnnotationDbi::select(org.Harmigera.eg.db,
keys = keys(org.Harmigera.eg.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)
}
<article class="_2rhmJa">
最后，如果不想构建orgdb还想做富集分析
满足需求永远是第一生产力，如果你只想用一次，那么clusterProfiler的enricher可以学习一下
主要的两个参数需要注意：gene 是基因ID，TERM2GENE 是GO/KEGG terms与基因ID的对应，例如上面图片中的GO_terms、KEGG_KOs等eggnog-mapper结果，提取出来就好。对于没有对应terms的基因ID，那我们就把它们去掉。
参考：http://guangchuangyu.github.io/2015/05/use-clusterprofiler-as-an-universal-enrichment-analysis-tool/
