基于 bioMart 构建绵羊（非常见物种） OrgDb 包/数据库

韩建刚（CAAS-UCD）

已于 2022-06-12 23:08:02 修改

阅读量3.3k

点赞数 3

分类专栏：注释文件&富集分析&ID转换文章标签：生物信息学

于 2022-06-12 13:48:35 首次发布

本文链接：https://blog.csdn.net/flashan_shensanceng/article/details/125223954

版权

注释文件&富集分析&ID转换专栏收录该内容

7 篇文章 1 订阅

订阅专栏

OrgDb （organism database）文件主要用于基因注释、ID转换、GO富集分析等，Bioconductor - BiocViews 仅提供部分物种正式发布的 OrgDb 包。此外还可通过 AnnotationHub 包检索一些未正式发布的 OrgDb 数据库，以绵羊为例，检索到的OrgDb 数据库是基于 Oar_v3.1 建立的，而我使用的参考基因组是Oar_rambouillet_v1.0，在下游分析或者ID转换很可能出现问题，因此建立基于Oar_rambouillet_v1.0 参考基因组的 OrgDb 包/数据库是非常有必要的

# 检索、下载、读取绵羊 orgdb 数据库
library(AnnotationHub)
library(AnnotationHubi)
ah=AnnotationHub()
query(ah,"org.Ovis_aries.eg.sqlite")
org.Oar.eg.db=ah[['AH96240']] 
saveDb(org.Oar.eg.db, file = "org.Oar.eg.db") # 保存
orgdb = loadDb(file = "org.Oar.eg.db") # 读取

以绵羊为例，有两种策略构建绵羊（非常见物种）的 OrgDb 包，一种是利用 eggnog-mapper 数据库比对参考基因组，将得到的注释信息用于建立 OrgDb 包；另一种方式是利用 bioMart 包检索相应参考基因组版本的注释信息，利用该注释信息构建 OrgDb 包。

第一种方式首先对参考基因组全部转录本的序列在 eggnog-mapper 数据库中进行比对，跨物种检索同源序列和同源基因，耗时较长，比较适合从 NCBI 下载的参考基因组进行后续分析。详细构建步骤可以参考下边几个文章。

1. 构建自己物种的orgDb - 简书

2. 生信 | 构建物种的OrgDb - 简书

3. AnnotationForge包构建非模式物种Orgdb包 - 简书

4. 应该是最好的eggnog-mapper功能注释教程 - 简书

# 下载eggnog-mapper软件，三种方式，推荐 conda 安装
# method 1, conda
# 先在 https://anaconda.org/ 搜索eggnog-mapper 安装方式
conda install -c bioconda eggnog-mapper

# method 2, git
# 非常不友好

# method 3, local
#

第二种方式利用的是 biomaRt 包的注释信息，biomaRt 包是基于 biomaRt 数据库和 Ensembl 数据库构建的，这种 OrgDb包构建方式适合从 Ensembl 下载的参考基因组，而且该方式比较灵活，而且信息量大，信息全面，可随意选择需要的信息用于构建OrgDb。

########################### biomaRt 包构建 org.Oaries.eg.db###################

biomaRt 可以在 bioconductor 网站下载，官网Accessing Ensembl annotation with biomaRt。
首先需要明确两个概念，database 和 dataset。database 是 bioMart 可用数据库，dataset 则是某个数据库下包含的所有数据集，我们访问的便是某个数据库中的某个数据集。使用基本流程就是：查看可用数据库，访问数据库，查看可用数据集，下载/访问数据集。

以绵羊为例：

一、绵羊参考基因组注释信息访问/下载

library(biomaRt)
# 1. 可用数据库
# display available Ensembl BioMart web services
listEnsembl()   #显示的是Ensembl最新版本，106
#          biomart                version
# 1         genes      Ensembl Genes 106
# 2 mouse_strains      Mouse strains 106
# 3          snps  Ensembl Variation 106
# 4    regulation Ensembl Regulation 106

# 2. 选择数据库
# select "genes" database
database = useEnsembl("genes")
database

# 3. 可用数据集
# display available dataset
dataset = listDatasets(database)

# 检索绵羊有关数据集（注释信息）
searchDatasets(mart = database, pattern = "Oar")  # 出现了两个版本，Oar_rambouillet_v1.0 和 Oar_v3.1
# sheep
# oarambouillet_gene_ensembl, Sheep genes (Oar_rambouillet_v1.0), Oar_rambouillet_v1.0
# oaries_gene_ensembl, Sheep (texel) genes (Oar_v3.1), Oar_v3.1

# 4. 访问/下载Oar_rambouillet_v1.0 数据集
# obtain sheep-biomart dataset
ensembl <- useDataset(dataset = "oarambouillet_gene_ensembl", mart = database)
# 或用useEnsembl
ensembl <- useEnsembl(biomart = "genes", dataset = "oarambouillet_gene_ensembl")
ensembl

# 5. 还可以选择特定版本的参考基因组，比如我的是 103 版本
# select specific version
listEnsemblArchives()
listEnsembl(version = 103)
# 访问指定版本的数据集
ensembl.103 <- useEnsembl(biomart = 'genes', dataset = 'oarambouillet_gene_ensembl', version = 103)

二、脊椎动物参考基因组可以在 Ensembl 官网下载，biomaRt 还额外提供除脊椎动物外的，原生动物，后生动物，真菌，植物等多个物种的参考基因组序列

# 1. 查看可用参考基因组数据库
# display available Ensembl Genomes marts
listEnsemblGenomes()
# biomart                        version
# 1       protists_mart      Ensembl Protists Genes 53
# 2 protists_variations Ensembl Protists Variations 53
# 3          fungi_mart         Ensembl Fungi Genes 53
# 4    fungi_variations    Ensembl Fungi Variations 53
# 5        metazoa_mart       Ensembl Metazoa Genes 53
# 6  metazoa_variations  Ensembl Metazoa Variations 53
# 7         plants_mart        Ensembl Plants Genes 53
# 8   plants_variations   Ensembl Plants Variations 53

# 2. 选择数据库
database_protists <- useEnsemblGenomes(biomart = "protists_mart")

# 3. protists_mart 数据库下可用数据集
dataset_protists  = listDatasets(database_protists)

# 4. 下载数据集
ensembl_arabidopsis <- useEnsemblGenomes(biomart = "protists_mart", dataset = "aculicifacies_eg_gene")

三、检索注释信息

首先需要理解什么是属性 attribute，可以理解为“每一种对基因的描述信息可以看成一个属性，例如ensembl id, entrez id, refseq id, symbol 等”，所有的attribute 可以分为四类 page，需要注意的是在检索注释信息时，只能检索属于同一个 page 的多个attribute 的信息。

# list all attributes of ensembl
# 一共3074种，可以分为4类
attribute = listAttributes(ensembl)
table(attribute$page)

有两个命令可以获取相应的检索信息，getBM(attributes, filters, values, ensembl) 和 select(ensembl,keys, keytype, columns)。getBM() 是专属于 biomaRt 的命令，select() 是经典的读取 AnnotationDb 类文件的命令，两者本质上是一样的，构建绵羊的 OrgDb 用的是 select 命令。

ensembl <- useEnsembl(biomart = 'genes', dataset = 'oarambouillet_gene_ensembl', version = 103)

此时的 ensembl 文件已经成功建立的访问特定版本（103）绵羊参考基因组的所有注释信息的链接，不能用我们熟悉的数据框、矩阵、向量来理解 ensembl，可以把 ensembl想象成一个没有交互界面的网站，这个网站很简单，只能够通过命令进行访问和检索。那么如何查阅 ensembl 以及用 ensembl 检索注释信息？还是4个经典的命令，columns, keytypes 和 keys。个人理解：columns (ensembl)，所有与 oarambouillet_v1 有关的属性变量，例如ensembl_id，entrez_id，chromosome_name，exon_id，biotype 等等，变量的数量和内容与上边的 attributes 是一样的；keytypes (ensembl) 相当于对属性变量进行分组，比如不同的基因可以按照 biotype 分析proteion coding RNA，rRNA，miRNA等，可以按照 chromosome_name 分为 1号染色体，5号染色体，X染色体等，而有些变量是不能分组的，比如 ensembl_id 和 entrez_id，每一个 id 都是唯一的；keys 就是用于指定检索属于特定 keytype 的基因名字、染色体号等。

keytypes(ensembl) %>% as.data.frame()
columns(ensembl) %>% as.data.frame() # 与attribute一致

# 不是所有的keytype都可以在keys中使用
keys(ensembl, keytype="chromosome_name")
keys(ensembl, keytype="biotype")
keys(ensembl, keytype="chromosomal_region")

# keytype只能用一种
# 每个参数都必须赋值

知道用法之后，就可以利用 select() 命令检索需要的 attribute 信息了，这里我将所有的注释信息分为了三部分ensembl_bse，go 和 ensembl_transcripts_exons。

# 1. ensembl_base 
# ensembl_base 包含所有基因的基本信息
ensembl_base = select(ensembl,keys = keys(ensembl, keytype="chromosome_name"), 
                      keytype= "chromosome_name",
                      columns=c("ensembl_gene_id", "entrezgene_id", "external_gene_name", 
                                "description", "gene_biotype", "ensembl_gene_id_version",
                                "start_position", "end_position", "strand", "chromosome_name"
                      ))
# 填充NA："entrezgene_id"或"external_gene_name"只有一个NA，则用另外一个填充
# 删除NA："entrezgene_id"或"external_gene_name"均为NA，删除行
ensembl_base$external_gene_name[ensembl_base$external_gene_name==""] = NA
for (i in 1:nrow(ensembl_base))
  if (is.na(ensembl_base[i,2]))  ensembl_base[i,2]=ensembl_base[i,3]
for (i in 1:nrow(ensembl_base))
  if (is.na(ensembl_base[i,3]))  ensembl_base[i,3]=ensembl_base[i,2]
ensembl_base = ensembl_base[complete.cases(ensembl_base$external_gene_name), ]
dim(ensembl_base)
colnames(ensembl_base)[1]="GID"
write.table(ensembl_base, file = "ensembl_base.txt", sep = "\t", 

# 2.
# 转录本和外显子
# transcript and exons
ensembl_ts_ex = select(ensembl,keys = keys(ensembl, keytype="chromosome_name"), 
                       keytype= "chromosome_name",
                       columns=c("ensembl_gene_id", "transcript_count", "ensembl_exon_id", 
                                 "ensembl_transcript_id", "transcript_biotype",
                                 "transcript_start", "transcript_end", "transcription_start_site"
                          ))
colnames(ensembl_ts_ex)[1]="GID"
write.table(ensembl_ts_ex, file = "ensembl_ts_ex.txt", sep = "\t", 
            row.names = FALSE, col.names = TRUE, quote = FALSE)

# 3.
# go 信息
ensembl_go = select(ensembl,keys = keys(ensembl, keytype="chromosome_name"),
                    keytype= "chromosome_name", 
                    columns=c("ensembl_gene_id", "go_id", "name_1006"))
ensembl_go$go_id[ensembl_go$go_id==""] = NA
ensembl_go = ensembl_go[complete.cases(ensembl_go$go_id), ]
colnames(ensembl_go) = c("GID",'GO',"EVIDENCE")
write.table(ensembl_go, file = "ensembl_go.txt", sep = "\t", 
            row.names = FALSE, col.names = TRUE, quote = FALSE)

四、构建 OrgDb

# 1. 读取对应文件（可略过）
gene_info = read.csv("ensembl_base.txt", sep="\t", header = TRUE)
transcript_exon = read.csv("ensembl_ts_ex.txt", sep="\t", header = TRUE)
go = read.csv("ensembl_go.txt", sep="\t", header = TRUE)

# 2. 构建 OrgDb 包
# 第一个参数是指定用于生成 OrgDb 的多个 dataframe，每个 dataframe 第一列列名必须是GID
# 每个文件的第一列 GID 不需要完全一致
# 所有的行都不能有行名
# 不同的 dataframe 不能包含相同的列名
makeOrgPackage(gene_info = gene_info, transcript_exon = transcript_exon,
               go = go,
               version="0.1",
               maintainer="Jiangang Han <jiangang_han@163.com>",
               author="Jiangang Han <jiangang_han@163.com>",
               outputDir = ".",
               tax_id="9940",
               genus="Ovis",
               species="aries",
               goTable="go")

# 3. 命令运行过程信息
Populating genes table:
genes table filled
Populating protein_coding table:
protein_coding table filled
Populating go table:
go table filled
table metadata filled

'select()' returned many:1 mapping between keys and columns
Dropping GO IDs that are too new for the current GO.db
Populating go table:
go table filled
Populating go_bp table:
go_bp table filled
Populating go_cc table:
go_cc table filled
Populating go_mf table:
go_mf table filled
'select()' returned many:1 mapping between keys and columns
Populating go_bp_all table:
go_bp_all table filled
Populating go_cc_all table:
go_cc_all table filled
Populating go_mf_all table:
go_mf_all table filled
Populating go_all table:
go_all table filled
Creating package in ./org.Oaries.eg.db 
Now deleting temporary database file
[1] "./org.Oaries.eg.db"

# 4. 成功生成"./org.Oaries.eg.db"文件夹

# 5. 退出Ｒ，org.Oaries.eg.db　同级目录运行
R CMD build ./org.Oaries.eg.db  #生成org.O.eg.db包
R CMD INSTALL BSgenome.Oaries.Ensembl.rambouilletv1.tar.gz　＃安装

四、其他（非必要信息）

　　可以根据　biotype　提取指定类型的基因，例如　protein coding　或　lncRNA。需要主义的是生成　OrgDb 包的时候，不同的 dataframe 不能包含相同的列名，因此只能选择下方某个文件或者对某些文件进行合并。

# protein coding
ensembl_base_coding = subset(ensembl_base, subset = gene_biotype=="protein_coding")
ensembl_base_coding$external_gene_name[ensembl_base_coding$external_gene_name==""] = NA
for (i in 1:nrow(ensembl_base_coding))
  if (is.na(ensembl_base_coding[i,2]))  ensembl_base_coding[i,2]=ensembl_base_coding[i,3]
for (i in 1:nrow(ensembl_base_coding))
  if (is.na(ensembl_base_coding[i,3]))  ensembl_base_coding[i,3]=ensembl_base_coding[i,2]
ensembl_base_coding = ensembl_base_coding[complete.cases(ensembl_base_coding$external_gene_name), ]
colnames(ensembl_base_coding)[1]="GID"
write.table(ensembl_base_coding, file = "ensembl_base_coding.txt", sep = "\t", 
            row.names = FALSE, col.names = TRUE, quote = FALSE)


# lncRNA
ensembl_base_lnc = subset(ensembl_base, subset = gene_biotype=="lncRNA")
ensembl_base_lnc = ensembl_base_lnc[,c(1,5:10)]
colnames(ensembl_base_lnc)[1]="GID"
write.table(ensembl_base_lnc, file = "ensembl_base_lnc.txt", sep = "\t", 
            row.names = FALSE, col.names = TRUE, quote = FALSE)
# miRNA
ensembl_base_mi = subset(ensembl_base, subset = gene_biotype=="miRNA")
ensembl_base_mi$external_gene_name[ensembl_base_mi$external_gene_name==""] = NA
ensembl_base_mi=ensembl_base_mi[complete.cases(ensembl_base_mi$external_gene_name), ]
for (i in 1:nrow(ensembl_base_mi))
  if (is.na(ensembl_base_mi[i,2]))  ensembl_base_mi[i,2]=ensembl_base_mi[i,3]
colnames(ensembl_base_mi)[1]="GID"
write.table(ensembl_base_mi, file = "ensembl_base_mi.txt", sep = "\t", 
            row.names = FALSE, col.names = TRUE, quote = FALSE)
# pseudogene
ensembl_base_pseudo = subset(ensembl_base, subset = gene_biotype=="processed_pseudogene")
ensembl_base_pseudo = ensembl_base_pseudo[,c(1,5:10)]
colnames(ensembl_base_pseudo)[1]="GID"
write.table(ensembl_base_pseudo, file = "ensembl_base_pseudo.txt", sep = "\t", 
            row.names = FALSE, col.names = TRUE, quote = FALSE)
# rRNA
ensembl_base_rRNA = subset(ensembl_base, subset = gene_biotype=="rRNA")
ensembl_base_rRNA$external_gene_name[ensembl_base_rRNA$external_gene_name==""] = NA
ensembl_base_rRNA = ensembl_base_rRNA[complete.cases(ensembl_base_rRNA$external_gene_name), ]
ensembl_base_rRNA = ensembl_base_rRNA[,-2]
colnames(ensembl_base_rRNA)[1]="GID"
write.table(ensembl_base_rRNA, file = "ensembl_base_rRNA.txt", sep = "\t", 
            row.names = FALSE, col.names = TRUE, quote = FALSE)
# Small nucleolar RNA
ensembl_base_sno = subset(ensembl_base, subset = gene_biotype=="snoRNA")
ensembl_base_sno$external_gene_name[ensembl_base_sno$external_gene_name==""] = NA
ensembl_base_sno = ensembl_base_sno[complete.cases(ensembl_base_sno$external_gene_name), ]
ensembl_base_sno = ensembl_base_sno[,-2]
colnames(ensembl_base_sno)[1]="GID"
write.table(ensembl_base_sno, file = "ensembl_base_sno.txt", sep = "\t", 
            row.names = FALSE, col.names = TRUE, quote = FALSE)
# spliceosomal RNA
ensembl_base_sn = subset(ensembl_base, subset = gene_biotype=="snRNA")
ensembl_base_sn$external_gene_name[ensembl_base_sn$external_gene_name==""] = NA
ensembl_base_sn = ensembl_base_sn[complete.cases(ensembl_base_sn$external_gene_name), ]
ensembl_base_sn = ensembl_base_sn[,-2]
colnames(ensembl_base_sn)[1]="GID"
write.table(ensembl_base_sn, file = "ensembl_base_sn.txt", sep = "\t", 
            row.names = FALSE, col.names = TRUE, quote = FALSE)


protein_coding = read.csv("ensembl_base_coding.txt", sep="\t", header = TRUE)
lncRNA = read.csv("ensembl_base_lnc.txt", sep="\t", header = TRUE)
miRNA = read.csv("ensembl_base_mi.txt", sep="\t", header = TRUE)
pseudogene = read.csv("ensembl_base_pseudo.txt", sep="\t", header = TRUE)
rRNA = read.csv("ensembl_base_rRNA.txt", sep="\t", header = TRUE)
snoRNA = read.csv("ensembl_base_sno.txt", sep="\t", header = TRUE)
snRNA = read.csv("ensembl_base_sn.txt", sep="\t", header = TRUE)