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也能得到上述操作同样的效果。