首先,TxDb 包是用 GenomicFeatures 包构建的,用于专门注释基因组中转录本、外显子、内含子等的包。
1. 如果是构建从 Ensembl 下载的参考基因组的 TxDb (transcript database)包,有两种构建方式,一种是直接利用 biomaRt 包构建,另一种是利用 .gff 注释文件构建,这里主要说的是第二种方式。
# 方式一,biomaRt
txdb <- makeTxDbFromBiomart(dataset="oarambouillet_gene_ensembl")
# 方式二,.gff 文件
# 1. 生成 metada,包含两列,name 和 value,必须文件
metadata <- data.frame(name="Resource URL", value="http://ftp.ensembl.org/pub/release-103/gtf/ovis_aries_rambouillet/Ovis_aries_rambouillet.Oar_rambouillet_v1.0.103.gtf.gz")
# 2. 利用 gff 文件生成 txdb 文件
txdb = makeTxDbFromGFF(file = "./Ovis_aries_rambouillet.Oar_rambouillet_v1.0.103.gtf", format="gtf", organism="sheep", taxonomyId=9940, dataSource="ensemblgenomes", metadata=metadata)
# 3. 保存和读取 txdb 文件
saveDb(txdb, file="TxDb.Oaries.Ensembl.Rambouilletv1.sqlite")
txdb <- loadDb("TxDb.Oaries.Ensembl.Rambouilletv1.sqlite")
# 4. 将 txdb文件 转成 txdb 包
# 需要注意命名规则,不能出现- _ 等符号
makeTxDbPackage(txdb, version="1.0", maintainer="Jiangang Han <jiangang_han@163.com>", author="Jiangang Han", destDir=".",license="Artistic-2.0",
pkgname="TxDb.Oaries.Ensembl.Rambouilletv1")
# 生成 TxDb.Oaries.Ensembl.Rambouilletv1 文件夹
# 5. 退出R,Linux 下 TxDb.Oaries.Ensembl.Rambouilletv1 同级目录运行
R CMD build ./TxDb.Oaries.Ensembl.Rambouilletv1 #生成TxDb.Oaries.Ensembl.Rambouilletv1.tar.gz包
R CMD INSTALL TxDb.Oaries.Ensembl.Rambouilletv1.tar.gz #安装
2. 检索 txdb 信息
library(TxDb.Oaries.Ensembl.Rambouilletv1)
library(GenomicFeatures)
txdb = TxDb.Oaries.Ensembl.Rambouilletv1
# 1. 染色体
# chromosomes that are active for the TxDb
seqlevels(txdb)
# only select chr1 or chr15
seqlevels(txdb) <- "chr1"
seqlevels(txdb) <- "chr15"
seqlevels(txdb)
# reset chromosome to original level
seqlevels(txdb) <- seqlevels0(txdb)
# 2. 四个经典访问命令
columns(), keytypes(), keys()
select(txdb, keys, columns, keytype)
# 3. 输出 GRange 对象
# return the coordinate information as a GRanges object.
# 转录本坐标,5列,seqnames,ranges,strand,tx_id,tx_name
GR <- transcripts(txdb)
GR[1:3]
## GRanges object with 3 ranges and 2 metadata columns:
## seqnames ranges strand | tx_id tx_name
## <Rle> <IRanges> <Rle> | <integer> <character>
## [1] chr15 20362688-20364420 + | 53552 uc001yte.1
## [2] chr15 20487997-20496811 + | 53553 uc001ytf.1
## [3] chr15 20723929-20727150 + | 53554 uc001ytj.3
# GRange 对象正负链
# tx_strand 和 GR 行数一致
tx_strand <- strand(GR)
# 进一步详细检索信息
GR <- transcripts(txdb, filter=list(tx_chrom = "chr15", tx_strand = "+"))
# 外显子坐标,四列
EX <- exons(txdb)
EX[1:4]
## GRanges object with 4 ranges and 1 metadata column:
## seqnames ranges strand | exon_id
## <Rle> <IRanges> <Rle> | <integer>
## [1] chr15 20362688-20362858 + | 192986
## [2] chr15 20362943-20363123 + | 192987
## [3] chr15 20364397-20364420 + | 192988
## [4] chr15 20487997-20488227 + | 192989
3. 计算启动子区域
感觉这个功能比较有意思,单独列了出来
# 将转录起始位点上游2000bp,下游400bp列为启动子区域
PR <- promoters(txdb, upstream=2000, downstream=400)
4. 分组展示基因,Working with Grouped Features
# 通过gene,exon 或 cds 分组展示转录本
transcriptsBy(x, by=c("gene", "exon", "cds"))
# 通过转录本分组展示外显子
exonsBy(x, by=c("tx", "gene"))
# 通过转录本或基因分组展示 cds 区域
cdsBy(x, by=c("tx", "gene"))
# 通过转录本分组展示内含子
intronsByTranscript(x)
# 通过转录本分组展示5' UTR
fiveUTRsByTranscript(x)
# 通过转录本分组展示3' UTR
threeUTRsByTranscript(x)
5. 获取序列信息
先读取预先构建的绵羊 BsGenome 包构建绵羊(非常见物种)BSgenome参考基因组_韩建刚(CAAS-UCD)的博客-CSDN博客
# 提取全部转录本序列信息
tx_seqs1 <- extractTranscriptSeqs(Oaries, TxDb.Oaries.Ensembl.Rambouilletv1, use.names=TRUE)
# 翻译成蛋白
suppressWarnings(translate(tx_seqs1))
cds_seqs <- extractTranscriptSeqs(Oaries,
cdsBy(txdb, by="tx", use.names=TRUE))
translate(cds_seqs)