トランスクリプトのアノテーション Bioconductor

GenomicFeatures

GenomicFeatures はトランスクリプトのアノテーションを取得したり、管理したりするためのパッケージである。そのやりとりは TxDb オブジェクトを通して行う。

library("GenomicFeatures")

ヒトのトランスクリプトアノテーションを例に利用してみる。

library(TxDb.Hsapiens.UCSC.hg38.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene

染色体の種類を出力するとき seqlevel 関数を利用する。

seqlevels(txdb)[1:25]
##  [1] "chr1"  "chr2"  "chr3"  "chr4"  "chr5"  "chr6"  "chr7"  "chr8"  "chr9"
## [10] "chr10" "chr11" "chr12" "chr13" "chr14" "chr15" "chr16" "chr17" "chr18"
## [19] "chr19" "chr20" "chr21" "chr22" "chrX"  "chrY"  "chrM"

アノテーションデータの取得

ある遺伝子 ID を与えて、それに関連するアノテーションを取得することができる。具体的にどのようなデータが取得できるかは columns 関数で調べることができる。また、データを取得するときにあたえる検索条件に利用できる属性は keytypes で調べることができる。

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

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

以下、遺伝子の Ensembl ID を与えて、その遺伝子のトランスクリプトの開始位置と終了位置を取得する例。TxDb.Hsapiens.UCSC.hg38.knownGene は Ensembl ID を受け付けないので、org.Hs.eg.db パッケージを利用して、Ensembl ID を Entrez ID に変換しておく。

library(org.Hs.eg.db)

id <- select(org.Hs.eg.db,
             keys = c("ENSG00000004660", "ENSG00000007314"),
             keytype = "ENSEMBL",
             columns = c("ENSEMBL", "ENTREZID"))
id
##           ENSEMBL ENTREZID
## 1 ENSG00000004660    84254
## 2 ENSG00000007314     6329

res <- select(txdb,
              keys = id[, 2],
              keytype = "GENEID",
              columns = c("GENEID", "TXSTART", "TXEND", "TXTYPE"))
res
##   GENEID TXTYPE  TXSTART    TXEND
## 1  84254   <NA>  3860323  3890743
## 2  84254   <NA>  3860323  3893043
## 3  84254   <NA>  3865157  3893043
## 4   6329   <NA> 63938554 63972918

mat <- merge(id, res, by.x = 2, by.y = 1)
mat
##   ENTREZID         ENSEMBL TXTYPE  TXSTART    TXEND
## 1     6329 ENSG00000007314   <NA> 63938554 63972918
## 2    84254 ENSG00000004660   <NA>  3860323  3890743
## 3    84254 ENSG00000004660   <NA>  3860323  3893043
## 4    84254 ENSG00000004660   <NA>  3865157  3893043

GRanges オブジェクトの取得

トランスクリプトのデータの取得

トランスクリプトの名前、開始位置や終了位置などは上述のように二次元配置で取得できるが、Bioconductor を利用した解析の場合は、 GRanges クラスのオブジェクトとして取得したほうが便利である。

gr <- transcripts(txdb)
## GRanges object with 104178 ranges and 2 metadata columns:
##            seqnames               ranges strand   |     tx_id     tx_name
##               <Rle>            <IRanges>  <Rle>   | <integer> <character>
##        [1]     chr1       [11874, 14409]      +   |         1  uc001aaa.3
##        [2]     chr1       [11874, 14409]      +   |         2  uc001aab.3
##        [3]     chr1       [11874, 14409]      +   |         3  uc010nxq.1
##        [4]     chr1       [30366, 30503]      +   |         4  uc031tlb.1
##        [5]     chr1       [69091, 70008]      +   |         5  uc001aal.1
##        ...      ...                  ...    ... ...       ...         ...
##   [104174]     chrY [26315642, 26315875]      -   |     92692  uc033fht.1
##   [104175]     chrY [26578095, 26578328]      -   |     92693  uc033fhu.1
##   [104176]     chrY [57094007, 57094084]      -   |     92694  uc033fhw.1
##   [104177]     chrY [57094564, 57094605]      -   |     92695  uc033fhx.1
##   [104178]     chrY [57212178, 57214703]      -   |     92696  uc011ncc.1
##   -------
##   seqinfo: 455 sequences (1 circular) from hg38 genome

この例では、染色体 1 から性染色体まで合計 104,178 個のトランスクリプトが登録されていることがわかる。各トランスクリプトの名前、座標やストランドなども書かれている。

例えば、X 染色体のデータのみがほしい場合は以下のようにする。

gr[seqnames(gr) == "chrX"]
## GRanges object with 3772 ranges and 2 metadata columns:
##          seqnames                 ranges strand   |     tx_id     tx_name
##               <Rle>            <IRanges>  <Rle>   | <integer> <character>
##      [1]     chrX       [276324, 303355]      +   |     88361  uc011mgx.2
##      [2]     chrX       [281394, 303355]      +   |     88362  uc004cpc.3
##      [3]     chrX       [319145, 321319]      +   |     88363  uc011mgz.3
##      [4]     chrX       [319145, 321319]      +   |     88364  uc033dot.1
##      [5]     chrX       [624344, 646823]      +   |     88365  uc004cph.1
##      ...      ...                    ...    ... ...       ...         ...
##   [3768]     chrX [155489011, 155612961]      -   |     92128  uc004fnn.4
##   [3769]     chrX [155492535, 155612961]      -   |     92129  uc004fnp.4
##   [3770]     chrX [155907487, 155907564]      -   |     92130  uc033fck.1
##   [3771]     chrX [155908044, 155908109]      -   |     92131  uc033fcl.1
##   [3772]     chrX [156025658, 156028183]      -   |     92132  uc011nad.1
##   -------
##   seqinfo: 455 sequences (1 circular) from hg38 genome

同様な考えで、トランスクリプト名が uc004fnp.4 のデータのみを取得したい場合は以下のようにする。

 gr[elementMetadata(gr)$tx_name == "uc004fnp.4"]
## GRanges object with 1 range and 2 metadata columns:
##       seqnames                 ranges strand |     tx_id     tx_name
##          <Rle>            <IRanges>  <Rle>   | <integer> <character>
##   [1]     chrX [155492535, 155612961]      - |     92129  uc004fnp.4
##   -------
##   seqinfo: 455 sequences (1 circular) from hg38 genome

エクソンのデータの取得

トランスクリプトのデータの取得のように、エクソンの場合は exons 関数を用いる。

ex <- exons(txdb)
ex
## GRanges object with 327656 ranges and 1 metadata column:
##            seqnames               ranges strand   |   exon_id
##               <Rle>            <IRanges>  <Rle>   | <integer>
##        [1]     chr1       [11874, 12227]      +   |         1
##        [2]     chr1       [12595, 12721]      +   |         2
##        [3]     chr1       [12613, 12721]      +   |         3
##        [4]     chr1       [12646, 12697]      +   |         4
##        [5]     chr1       [13221, 14409]      +   |         5
##        ...      ...                  ...    ... ...       ...
##   [327652]     chrY [57094007, 57094084]      -   |    293786
##   [327653]     chrY [57094564, 57094605]      -   |    293787
##   [327654]     chrY [57212178, 57213357]      -   |    293788
##   [327655]     chrY [57213856, 57213964]      -   |    293789
##   [327656]     chrY [57214350, 57214703]      -   |    293790
##   -------
##   seqinfo: 455 sequences (1 circular) from hg38 genome

塩基配列データの取得

BSgenoem パッケージのデータと合わせて利用することで、トランスクリプトの塩基配列を取得することができる。ただし、TxDb と BSgenome のゲノムのバージョンが一致しているのを確認すること。次の例ではどちらも UCSC version 38 であることをパッケージの名前から確認できる。

library(TxDb.Hsapiens.UCSC.hg38.knownGene)
library(BSgenome.Hsapiens.UCSC.hg38)
txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
seq <- extractTranscriptSeqs(Hsapiens, txdb)
seq
##   A DNAStringSet instance of length 104178
##           width seq                                         names
##      [1]   1652 CTTGCCGTCAGCCTTTTCTT...AAGCACACTGTTGGTTTCTG uc001aaa.3
##      [2]   1595 CTTGCCGTCAGCCTTTTCTT...AAGCACACTGTTGGTTTCTG uc001aab.3
##      [3]   1488 CTTGCCGTCAGCCTTTTCTT...AAGCACACTGTTGGTTTCTG uc010nxq.1
##      [4]    138 GGATGCCCAGCTAGTTTGAA...ATTCAGAATTAAGCATTTTA uc031tlb.1
##      [5]    918 ATGGTGACTGAATTCATTTT...ATTCTAGTGTAAAGTTTTAG uc001aal.1
##      ...    ... ...
## [104174]     28 TGGGAAGTGAGGAGTGTCTCTGCCTGGC                uc033dok.1
## [104175]     30 CAGAATGGTAGATCCACTGACAGCTTGCAC              uc033dom.1
## [104176]    844 GCGGCACCCAGAGAGGGTCC...AGGAAACGAATAAATCTCTG uc033don.1
## [104177]   1420 TTGGGTATCCAGGCTCTGAG...AATAAATGTTCACCAAGGCA uc033doo.1
## [104178]     28 TGATTGGTGTGTTTACAAACCTTGAGCT                uc033dop.1

translate(seq[1:20])
##   A AAStringSet instance of length 20
##      width seq                                              names
##  [1]   550 LAVSLFFDLFFLFMCICCLLAQT...LTPRRLHPAQLEILY*KHTVGF uc001aaa.3
##  [2]   531 LAVSLFFDLFFLFMCICCLLAQT...LTPRRLHPAQLEILY*KHTVGF uc001aab.3
##  [3]   496 LAVSLFFDLFFLFMCICCLLAQT...DPETFASCTARDPLLKAHCWFL uc010nxq.1
##  [4]    46 GCPASLNFR*TTNNFVA*ICPKLSLGHTYAKKHYWLFI*DSELSIL   uc031tlb.1
##  [5]   306 MVTEFIFLGLSDSQELQTFLFML...KDMKTAIRQLRKWDAHSSVKF* uc001aal.1
##  ...   ... ...
## [16]  2163 VPGAVLLPSRRAAGTRVLSLLSQ...STFHLIFLYKHFSSE*KFFQTV uc031pjk.1
## [17]  1014 VLTRSWQWQVSGFRIKPRIRNFT...RLWRAKIAPLHPSLGNKSKTPS uc001abu.1
## [18]   163 QRLGLAGRSLRSEGKSLKTLMSK...MPEHQSRCEFQRGSLEIGLRPA uc001abv.1
## [19]   851 ADPCGVREGGTGSGLGSRAEGKV...AQPSSRYVWSFMLKNK**SLSR uc001abw.1
## [20]   706 GKV*RRLCPRGSCRCILRSATAR...PSRRMGPWLYFQGPPTLPSLCV uc031pjl.1
## There were 14 warnings (use warnings() to see them)