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

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)

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

构建绵羊（非常见物种）BSgenome参考基因组
1. Ensembl 下载绵羊参考基因组和注释文件,虽然这里用不到注释文件,但最好备份一下 2. .fa 格式参考基因组转为 .2bit 格式 # cd software directory wget ...
seer文献_文献精读|基于SEER数据库构建肝癌非转移Nomogram
文献精读|基于SEER数据库构建肝癌非转移Nomogram 作者:白介素2 摘要目的文章的目的比较明确,寻找肝癌肺转移相关的风险因素,并建立预后 Nomogram. 方法回顾性分析 SEER数据 ...
在微服务架构下基于 Prometheus 构建一体化监控平台的最佳实践
欢迎关注方志朋的博客,回复"666"获面试宝典随着 Prometheus 逐渐成为云原生时代的可观测事实标准,那么今天为大家带来在微服务架构下基于 Prometheus 构建一体 ...
付力力: 基于 ImpalaS 构建实时用户行为分析引擎
本文来自神策数据联合创始人&首席架构师付力力在 QCon 北京 2017 年全球软件开发者大会上的精彩分享,主题是"基于 ImpalaS 构建实时用户行为分析引擎". 付力 ...
技术干货｜阿里云基于Hudi构建Lakehouse实践探索
简介:阿里云高级技术专家王烨(萌豆)在Apache Hudi 与 Apache Pulsar 联合 Meetup 杭州站上的演讲整理稿件,本议题介绍了阿里云如何使用 Hudi 和 OSS 对象存储构建 ...
技术干货｜阿里云基于Hudi构建Lakehouse实践探索「内附干货PPT下载渠道」
简介: 阿里云高级技术专家王烨(萌豆)在Apache Hudi 与 Apache Pulsar 联合 Meetup 杭州站上的演讲整理稿件,本议题介绍了阿里云如何使用 Hudi 和 OSS 对象存储构 ...
基于yacto构建am5728 SDK
基于yacto构建am5728 SDK 构建SDK 1.1 介绍本页面提供了从源代码构建处理器SDK和各个组件的步骤. ProcessorSDK构建基于Arago项目,Arago项目为OpenEmb ...
Docker学习之六：基于Dockerfile构建镜像
镜像制作一般镜像的制作,通常需要修改镜像的配置文件,比如nginx的配置文件,可以通过以下的方式: 将配置文件做成存储卷,从宿主机编辑好之后,启动容器时应用程序加载配置文件的路径并和宿主机的目录建立 ...
Spark论文思想之-基于RDD构建的模型（Shark的来龙去脉）
3.1 介绍首先RDD提供以下功能: 跨集群的不可变存储(在Spark中,记录是指Java Object) 使用键对数据进行分区控制考虑分区的粗粒度运算符由于是内存计算,所以低延迟 3.2 在R ...

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

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

最新文章

热门文章