147-学习使用TxDb对象

Nov 17, 2019 4 min read cnposts

刘小泽写于19.11.17
教程来自Bioconductor：https://bioconductor.org/packages/release/bioc/vignettes/GenomicFeatures/inst/doc/GenomicFeatures.pdf【更新于19-10-29】

0 前言

需要一个包：GenomicFeatures ，方便从UCSC、BioMart数据库中提取、处理转录本相关的信息

suppressPackageStartupMessages(library('GenomicFeatures'))

这个包就是利用TxDb对象来存储转录本的metadata，包含了5’/3’非翻译区（UTRs）、蛋白编码区（CDSs）、外显子。TxDb对象的背后是一个SQLite数据库，包含了基因组坐标，前mRNA、外显子、蛋白编码序列的对应关系和相对应的基因ID

1 如何从TxDb对象获取数据

1-1 首先加载转录本数据

最简单的方式就是这样：

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

从中可以看到，这个数据中默认的基因ID是Entrez ID，因此如果自己使用的基因名不一致，还需要提前转换

1-2 可以根据想要的染色体进行过滤

查看目前包含的染色体

seqlevels(txdb)
# 你会发现除了常见的染色体，还有大量命名很奇怪的，比如”"chrUn_gl000233"“

如果只想要1号染色体的

seqlevels(txdb) <- "chr1"

又想复位怎么办？

想要原来全部的染色体，重新筛选，两种方法：

（简单粗暴）重新加载一遍那个包，把txdb覆盖掉
（有点技术）利用一个复位函数
```
seqlevels(txdb) <- seqlevels0(txdb)
```

1-3 提取想要的数据

TxDb和ChipDb、OrgDb都很像，都支持columns()、keytypes()、keys()、select() 等函数，对于提取子集最有帮助了

查看txdb数据库支持哪些ID转换

> columns(txdb)
 [1] "CDSCHROM"   "CDSEND"     "CDSID"      "CDSNAME"    "CDSSTART"  
 [6] "CDSSTRAND"  "EXONCHROM"  "EXONEND"    "EXONID"     "EXONNAME"  
[11] "EXONRANK"   "EXONSTART"  "EXONSTRAND" "GENEID"     "TXCHROM"   
[16] "TXEND"      "TXID"       "TXNAME"     "TXSTART"    "TXSTRAND"  
[21] "TXTYPE"

看看txdb允许查询的数据类型

> keytypes(txdb)
[1] "CDSID"    "CDSNAME"  "EXONID"   "EXONNAME" "GENEID"   "TXID"    
[7] "TXNAME"

用select挑选出来

需要提供 注释信息 + 查询ID + 转换后的ID类型 + 查询ID的类型

keys <- c("100033416", "100033417", "100033420")
> select(txdb, keys, columns="TXNAME", keytype="GENEID")
'select()' returned 1:1 mapping between keys and columns
     GENEID     TXNAME
1 100033416 uc001yxl.4
2 100033417 uc001yxo.3
3 100033420 uc001yxr.3

练习：还是上面👆的基因，找出它们所在的染色体和正负链

cols <- c("TXNAME", "TXSTRAND", "TXCHROM")

> select(txdb, keys=keys, columns=cols, keytype="GENEID")
'select()' returned 1:1 mapping between keys and columns
     GENEID     TXNAME TXCHROM TXSTRAND
1 100033416 uc001yxl.4   chr15        +
2 100033417 uc001yxo.3   chr15        +
3 100033420 uc001yxr.3   chr15        +

1-4 将提取的数据返回为GRanges对象

虽然用select提取数据很方便，但如果能保存成GRanges对象，后面提取信息会更加便利。比如提取外显子、转录本、编码区域的坐标信息

使用transcripts()、exons()、cds()函数，它们运行的结果就会返回GRanges对象

获取chr15的转录本的GRanges对象 => `transcripts()`

seqlevels(txdb) <- "chr15"
GR <- transcripts(txdb)
# chr15的转录本数量
> length(GR)
[1] 3337

> 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
  -------
  seqinfo: 1 sequence from hg19 genome

获取chr15的正负链信息 => `strand()`

tx_strand <- strand(GR)

> tx_strand
factor-Rle of length 3337 with 2 runs
  Lengths: 1732 1605
  Values :    +    -
Levels(3): + - *
# 正负链总数
> sum(runLength(tx_strand))
[1] 3337

利用transcripts函数进行筛选

利用的过滤条件存储在columns(txdb)中

比如取出3号染色体的正链转录本

# 先把txdb复位
seqlevels(txdb) <- seqlevels0(txdb)
# 然后用tx_chrom指定染色体编号，tx_strand指定链
GR_chr3 <- transcripts(txdb, filter=list(tx_chrom = "chr3", tx_strand = "+"))

> length(GR_chr3)
[1] 2179
# 看结果链的类型：全都是+链
> unique(strand(GR_chr3))
[1] +
Levels: + - *

利用promoters函数获得启动子区域

# 比如要获得chr3的信息
seqlevels(txdb) <- "chr3"
PR <- promoters(txdb, upstream=2000, downstream=400)
> PR
GRanges object with 4328 ranges and 2 metadata columns:
             seqnames              ranges strand |     tx_id     tx_name
                <Rle>           <IRanges>  <Rle> | <integer> <character>
  uc003bot.3     chr3       236279-238678      + |     13060  uc003bot.3
  uc003bou.3     chr3       236279-238678      + |     13061  uc003bou.3
  uc003bov.2     chr3       237326-239725      + |     13062  uc003bov.2
  uc003bow.2     chr3       237326-239725      + |     13063  uc003bow.2
  uc011asi.2     chr3       359366-361765      + |     13064  uc011asi.2
         ...      ...                 ...    ... .       ...         ...
  uc003fyp.3     chr3 197686487-197688886      - |     17383  uc003fyp.3
  uc003fyq.4     chr3 197676397-197678796      - |     17384  uc003fyq.4
  uc021xjx.1     chr3 197773630-197776029      - |     17385  uc021xjx.1
  uc003fyx.3     chr3 197807143-197809542      - |     17386  uc003fyx.3
  uc010iat.2     chr3 197807143-197809542      - |     17387  uc010iat.2
  -------
  seqinfo: 1 sequence from hg19 genome

1-5 整合分组（grouped features）

有时关心的不止是某个基因/转录本的坐标，还有：某个基因包含多少转录本；某个转录本有多少外显子

就可以利用transcriptsBy()、exonsBy()或cdsBy()

转录本按照基因分组

transcriptsBy()函数返回的是GRangesList对象。

GRList <- transcriptsBy(txdb, by = "gene")
length(GRList)
[1] 1271

> names(GRList)[10:13]
[1] "100128025" "100128164" "100128378" "100128640"

> GRList[11:12]
GRangesList object of length 2:
$`100128164`
GRanges object with 2 ranges and 2 metadata columns:
      seqnames              ranges strand |     tx_id     tx_name
         <Rle>           <IRanges>  <Rle> | <integer> <character>
  [1]     chr3 169661772-169684522      - |     17029  uc011bpp.2
  [2]     chr3 169672015-169684522      - |     17030  uc003fgf.3
  -------
  seqinfo: 1 sequence from hg19 genome

$`100128378`
GRanges object with 1 range and 2 metadata columns:
      seqnames            ranges strand |     tx_id     tx_name
         <Rle>         <IRanges>  <Rle> | <integer> <character>
  [1]     chr3 52096110-52099128      - |     16010  uc010hmb.2
  -------
  seqinfo: 1 sequence from hg19 genome

除了按照基因分组，还可以按照exon/cds分组

除了对转录本进行分组，还可以使用exonsBy / cdsBy，不过它们只能按照转录本或基因分组tx/gene

获取每个转录本的全部外显子 => `exonsBy`

> GRList <- exonsBy(txdb, by = "tx")
> length(GRList)
[1] 4328
> names(GRList)[10:13]
[1] "13069" "13070" "13071" "13072"

> GRList[[12]]
GRanges object with 25 ranges and 3 metadata columns:
       seqnames          ranges strand |   exon_id   exon_name exon_rank
          <Rle>       <IRanges>  <Rle> | <integer> <character> <integer>
   [1]     chr3 2140550-2140662      + |     47793        <NA>         1
   [2]     chr3 2142242-2142323      + |     47795        <NA>         2
   [3]     chr3 2380862-2380917      + |     47797        <NA>         3
   [4]     chr3 2613100-2613242      + |     47799        <NA>         4
   [5]     chr3 2777899-2778025      + |     47800        <NA>         5
   ...      ...             ...    ... .       ...         ...       ...

147-学习使用TxDb对象

0 前言

1 如何从TxDb对象获取数据

1-1 首先加载转录本数据

1-2 可以根据想要的染色体进行过滤

查看目前包含的染色体

如果只想要1号染色体的

又想复位怎么办？

1-3 提取想要的数据

查看txdb数据库支持哪些ID转换

看看txdb允许查询的数据类型

用select挑选出来

练习：还是上面👆的基因，找出它们所在的染色体和正负链

1-4 将提取的数据返回为GRanges对象

获取chr15的转录本的GRanges对象 => `transcripts()`

获取chr15的正负链信息 => `strand()`

利用transcripts函数进行筛选

利用promoters函数获得启动子区域

1-5 整合分组（grouped features）

转录本按照基因分组

获取每个转录本的全部外显子 => `exonsBy`

Yunze Liu

Bioinformatics Sharer

Related

147-学习使用TxDb对象

0 前言

1 如何从TxDb对象获取数据

1-1 首先加载转录本数据

1-2 可以根据想要的染色体进行过滤

查看目前包含的染色体

如果只想要1号染色体的

又想复位怎么办？

1-3 提取想要的数据

查看txdb数据库支持哪些ID转换

看看txdb允许查询的数据类型

用select挑选出来

练习：还是上面👆的基因，找出它们所在的染色体和正负链

1-4 将提取的数据返回为GRanges对象

获取chr15的转录本的GRanges对象 => transcripts()

获取chr15的正负链信息 => strand()

利用transcripts函数进行筛选

利用promoters函数获得启动子区域

1-5 整合分组（grouped features）

转录本按照基因分组

获取每个转录本的全部外显子 => exonsBy

Yunze Liu

Bioinformatics Sharer

Related

获取chr15的转录本的GRanges对象 => `transcripts()`

获取chr15的正负链信息 => `strand()`

获取每个转录本的全部外显子 => `exonsBy`