[R]bioconductor之GenomicFeatures学
GenomicFeatures包主要提供了TxDb对象结构。基于此对象的操作，we can easily download the genomic locations of the transcripts, exons and cds of a given organism s from the UCSC Genome Bioinformatics and BioMart2 data resources.
官方教程：https://www.bioconductor.org/packages/release/bioc/html/GenomicFeatures.html
https://www.bioconductor.org/packages/release/bioc/vignettes/GenomicFeatures/inst/doc/GenomicFeatures.pdf
一、TxDb对象构建
常用的方法：加载对应R包
一般常用的物种基因注释信息都已经有人构建好TxDb对象，以R包的形式上传到Bioconductor里。我们使用时直接加载、赋值即可。
TxDb 包
BiocManager::install("TxDb.Hsapiens.UCSC.hg19.knownGene")
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
class(txdb)
txdb
txdb
如官网所说，由于TxDb对象基于SQL支持，因此可use the loadDb method to load the object directly from an appropriate .sqlite database file. 不过多学习了。
此外GenomicFeatures还支持自己手动从指定数据库中构建TxDb对象，例如从 UCSC提取构建的makeTxDbFromUCSC()函数；从BioMart提取构建的makeTxDbFromBiomart()函数等方法可以用作替补方案（没有对应TXDB的R包时）
二、子集筛选
1、指定染色体
seqlevel
head(seqlevels(txdb))
#只保留chr1
seqlevels(txdb) <- "chr1"
seqlevels(txdb)
#恢复初始设置
seqlevels(txdb) <- seqlevels0(txdb)
head(seqlevels(txdb))
seqlevels
2、key-column select模式
-
key指定某一列的任意范围；keytype为该列的列名；一般只能选取以下7列
keytypes(txdb)
#[1] "CDSID" "CDSNAME" "EXONID" "EXONNAME" "GENEID" "TXID" "TXNAME"
-
column为其它列的信息
#查看所有的column
columns(txdb)
columns(txdb)
- 示例操作
keys <- c("100033416", "100033417", "100033420")
select(txdb, keys = keys, keytype="GENEID", columns="TXNAME")
select
cols <- c("TXNAME", "TXSTRAND", "TXCHROM")
select(txdb, keys=keys, keytype="GENEID", columns=cols)
select2
3、提取transcripts, exons, cds, genes, promoters GRanges
- 主要是提取基因组里上述三种feature的起始坐标、strand作为Granges对象。
- Granges对象的相关操作可见之前的笔记。
- 分别对应
transcripts(),exons(),cds(),genes(),promoters()五个函数
GR.tr <- transcripts(txdb)
GR.tr
GR.tr
#配合filter参数，筛选指定信息
GR.ex.sub <- exons(txdb, filter=list(tx_chrom = "chr15", tx_strand = "-"))
GR.ex.sub
GR.ex.sub
关于 promoters启动子，由于其定义为transcription start site(TSS)位点上下游范围的一段序列，具体上下游的长度可自己定义，
promoters()函数默认设置为upstream=2000, downstream=200
head(promoters(txdb)) #2000+200,即2200的长度
promoters
4、Grouped Features & GRangesList
- 基因组的不同feature往往存在包含关系。比如某个基因可以有多个转录本；某个转录本又包含多个外显子。
- 可以通过对应的函数，得到存在如此包含关系的两种feature的情况，储存结果一般为
GRangesList -
transcriptsBy()，group by gene, exon or cds
##group transcripts by genes:
length(GR.tr);length(GR.ge)
#[1] 82960
#[1] 23056
GRList.tr2ge <- transcriptsBy(txdb, by = "gene")
class(GRList.tr2ge)
length(GRList.tr2ge)
#[1] 23459
head(GRList.tr2ge,2)
head(GRList.tr2ge,2)
-
exonsBy()，cdsBy()be grouped by transcript (tx) or gene.
length(GR.ex);length(GR.tr)
#[1] 289969
#[1] 82960
GRList.ex2tx <- exonsBy(txdb, by = "tx")
length(GRList.ex2tx)
#[1] 82960
此外GenomicFeatures包还提供了三种prespecified form function：
intronsByTranscript,fiveUTRsByTranscript,threeUTRsByTranscript，可通过名字直接了解包含关系，方便使用。
5、获取对应的序列信息
- 需要加载对应物种基因组的相关注释包，注意版本号要相同
- 利用之前我们使用的是TxDb.Hsapiens.UCSC.hg19.knownGene，那么就用对应的BSgenome.Hsapiens.UCSC.hg19
#BiocManager::install("BSgenome.Hsapiens.UCSC.hg19")
library(BSgenome.Hsapiens.UCSC.hg19)
BSgenome = BSgenome.Hsapiens.UCSC.hg19
-
getSeq()函数：需要提供对应物种的BSgenome，以及提取目标，一般是Granges对象（染色体、起始位点）
test <- head(GR.ex)
tmp <- getSeq(BSgenome,test)
tmp
class(tmp)
#[1] "DNAStringSet"
#attr(,"package")
#[1] "Biostrings"
tmp
- 这里又涉及到了Biostrings包的
DNAStringSet对象的学习，例如可以将核苷酸序列转换为氨基酸序列（如下）；以后有机会再深入了解下。
nchar(as.character(tmp[[1]]))
translate(tmp)
tx_seqs1 <- extractTranscriptSeqs(Hsapiens, TxDb.Hsapiens.UCSC.hg19.knownGene,
use.names=TRUE)
suppressWarnings(translate(tx_seqs1))
translate
