48.《Bioinformatics Data Skills》之
GenomicFeatures包
通过R包存储与提取数据是一个很好的选择，有利于数据版本管理与代码可重复性。GenomicFeatures就是这样的一个R包，此包嵌入了一系列的转录组注释数据（如图1），提供了一系列的统一方便的数据提取操作，下面让我们了解一下。
图1 GenomicFeatures包含的数据
当我们安装完GenomicFeatures包之后，一系列的依赖数据包就已经存在，这些包的命名方式统一为：TxDb.<物种>.<来源>.<版本号>.<名称>，可以直接加载。例如加载来自UCSC的人类基因组hg19的已知基因数据：
> library(TxDb.Hsapiens.UCSC.hg19.knownGene)
> txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
> txdb
TxDb object:
# Db type: TxDb
# Supporting package: GenomicFeatures
# Data source: UCSC
# Genome: hg19
# Organism: Homo sapiens
# Taxonomy ID: 9606
# UCSC Table: knownGene
# Resource URL: http://genome.ucsc.edu/
# Type of Gene ID: Entrez Gene ID
# Full dataset: yes
# miRBase build ID: GRCh37
# transcript_nrow: 82960
# exon_nrow: 289969
# cds_nrow: 237533
# Db created by: GenomicFeatures package from Bioconductor
# Creation time: 2015-10-07 18:11:28 +0000 (Wed, 07 Oct 2015)
# GenomicFeatures version at creation time: 1.21.30
# RSQLite version at creation time: 1.0.0
# DBSCHEMAVERSION: 1.1
GenomicFeatures包的数据统一采用TxDb（即transcriptDb）数据类型，通过一系列地函数对转录组，蛋白编码区域，外显子区域，启动子区域等（图2）进行统一的提取。
图2 基因
提取特征
通过transcripts(), exons(), cds(), genes(), promoters() （采用help(transcripts)查看）方法提取不同区域数据，例如提取所有的基因区域：
> hm_genes <- genes(txdb)
> length(hm_genes)
[1] 23056
> head(hm_genes)
GRanges object with 6 ranges and 1 metadata column:
seqnames ranges strand | gene_id
<Rle> <IRanges> <Rle> | <character>
1 chr19 58858172-58874214 - | 1
10 chr8 18248755-18258723 + | 10
100 chr20 43248163-43280376 - | 100
1000 chr18 25530930-25757445 - | 1000
10000 chr1 243651535-244006886 - | 10000
100008586 chrX 49217763-49233491 + | 100008586
-------
seqinfo: 93 sequences (1 circular) from hg19 genome
已知的人类基因数目有2万多个。
相比于直接提取区域，更常见的操作是按照某个特征进行分组提取，比如提取所有转录组（或基因）的外显子区域：
> hm_exons_by_tx <- exonsBy(txdb, by = "tx")
> hm_exons_by_gn <- exonsBy(txdb, by = "gene")
> length(hm_exons_by_tx)
[1] 82960
> length(hm_exons_by_gn)
[1] 23459
注：by="tx"代表转录组，根据需要可以替换为"exon"和"cds".
与直接提取区域的函数类似，这些函数家族包括transcriptsBy(), exonsBy()，cdsBy()等（通过help(exonsBy)函数查看）。
限定查询区域
我们可以将查询范围限定在某个区域，例如限定在1号染色体上查询，可以采取如下操作：
> seqlevels(txdb) <- "chr1"
提取特征的结果就只包含1号染色体区域了：
> chr1_exons <- exonsBy(txdb, by = "tx")
> all(unlist(seqnames(chr1_exons)) == "chr1")
[1] TRUE
之后可以如下命令恢复所有区域的查询：
> txdb <- restoreSeqlevels(txdb)
查询重叠区域
查询与某个区域重叠的特征同样有一个家族的函数，包括transcriptsByOverlaps(), exonsByOverlaps(), cdsByOverlaps()等（通过help(transcriptsByOverlaps)了解）。举例说明，假设8号染色体[123000000, 123500000]区间为特征区域，想知道这段区域上都有哪些转录组，可以采取如下操作：
> region <- GRanges("chr8", IRange(123000000, 123500000))
> region_expand <- GRanges("chr8", IRanges(123000000, 123500000)) + 10e3
> transcriptsByOverlaps(txdb, region_expand)
GRanges object with 1 range and 2 metadata columns:
seqnames ranges strand | tx_id tx_name
<Rle> <IRanges> <Rle> | <integer> <character>
[1] chr8 122966845-123139423 - | 33971 uc003ypj.3
-------
seqinfo: 93 sequences (1 circular) from hg19 genome
这里我们将特征区域的上下游都扩大了1Kb，事实上这个家族的函数都包括了一个参数叫做maxgap，默认为-1L（若两个范围相信相邻则gap为0），通过设定这个参数为1000也能得到上述操作同样的效果。
