基于 gff 文件构建 TxDb 包

首先，TxDb 包是用 GenomicFeatures 包构建的，用于专门注释基因组中转录本、外显子、内含子等的包。

1. 如果是构建从 Ensembl 下载的参考基因组的 TxDb （transcript database）包，有两种构建方式，一种是直接利用 biomaRt 包构建，另一种是利用 .gff 注释文件构建，这里主要说的是第二种方式。

Introduction

# 方式一，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)