GvizでCRPの周囲のゲノムブラウザの図を描いてみる #R

はじめに

GvisはRのパッケージで、ゲノムブラウザのような図を生成することができます。
この記事は基本的には下記の記事の焼き直しで新しい内容は特にありません。

UCSC Genome Browser をRのplotで再現する

ここでは、血液検査の結果で馴染み深い CRP（C反応性蛋白）の周囲の領域について、ゲノムブラウザのような画面を作成してみます。

ライブラリの読み込み

library(gviz)

染色体の番号と、リファレンスゲノムの指定。ここでは chr1:159,681,503-159,684,954 を表示します。

chr <- "chr1"
gen <- "hg19"
from <- 159681503
to   <- 159684954

イデオグラムを追加する

itrack <- IdeogramTrack(genome = gen, chromosome = chr)

plotTracks(list(itrack), from = from, to = to)

Gトラックを追加する

gtrack <- GenomeAxisTrack()

plotTracks(list(itrack, gtrack), from = from, to = to)

UCSCにアクセスして各種トラックを取得する

knownGenes <- UcscTrack(genome=gen, chromosome=chr, track="knownGene", from=from, to=to,
                        trackType="GeneRegionTrack", rstarts="exonStarts", rends="exonEnds",
                        gene="name", symbol="name", transcript="name", strand="strand", 
                        fill="#8282d2", name="UCSC Genes")

refGenes <- UcscTrack(genome=gen, chromosome=chr, track="xenoRefGene", from=from, to=to,
                      trackType="GeneRegionTrack", rstarts="exonStarts", rends="exonEnds", 
                      gene="name", symbol="name2", transcript="name", strand="strand", 
                      fill="#8282d2", stacking="dense", name="Other RefSeq")

ensGenes <- UcscTrack(genome=gen, chromosome=chr, track="ensGene", from=from, to=to,
                      trackType="GeneRegionTrack", rstarts="exonStarts", rends="exonEnds", 
                      gene="name", symbol="name2", transcript="name", strand="strand", 
                      fill="#960000", stacking="dense", name="Ensembl Genes")

plotTracks(list(idxTrack, axTrack, knownGenes,refGenes,ensGenes), from=from, to=to, showTitle=TRUE)

気になったのは、同じように書いたつもりなのに成功したり失敗したりする場合もあるみたいで、UCSCのサーバーにアクセスしているので、サーバの調子によってうまくいったり失敗するのかもしれません。

Biomartを使う

BioMartは接続が不安定なので、mirrorにはできるだけasiaを指定した方がいい

library(biomaRt)

bm <- useEnsembl(host = "grch37.ensembl.org", 
                 biomart = "ENSEMBL_MART_ENSEMBL", 
                 dataset = "hsapiens_gene_ensembl",
                 mirror="asia")

biomTrack <- BiomartGeneRegionTrack(genome = gen, chromosome = chr, 
                                    start = from, end = to,
                                    name = "ENSEMBL", biomart = bm)

plotTracks(list(idxTrack, axTrack, knownGenes,refGenes,ensGenes, biomTrack), from=from, to=to, showTitle=TRUE)