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 3893043GRanges オブジェクトの取得
トランスクリプトのデータの取得
トランスクリプトの名前、開始位置や終了位置などは上述のように二次元配置で取得できるが、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)