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

USCS Genome Browserは、ゲノム中の遺伝子領域やCpG Island、配列保存度や一塩基多型部位などさまざまな情報（Track）を並べて可視化することができるサイトである。上の図はマウスゲノム(NCBI37/mm9)のX染色体65,921,878 - 65,980,988部分を表示している。この図を得るためには、ふつう上記のサイトで

• group: mammal, genome: Mouse, assembly: July 2007 (NCBI37/mm9) を選択
• テキストボックスに chrX:65921878 - 65980988 と入れてsubmit
• 表示したい個々のTrackを選択してrefresh

のような動作をGUI的に行う必要がある。

これらの一連の処理をRのスクリプトを書いて再現出来れば、ある基準によって抽出したゲノム中の領域について、知りたい個々の情報を自動的にplotすることができるのでとても便利な気がする。

...ということをGvizパッケージを用いるとできるようになるのでメモしておく。

##準備

インストール

source("http://bioconductor.org/biocLite.R")
biocLite("Gviz")

パッケージの読み込み

library(Gviz)

表示領域の指定

genome <- "mm9"
chr <- "chrX"
from <- 65921878
to <- 65980988

##データの読み込み

各Trackの読み込み

# UCSCにアクセスしているため、一つのTrackを読み込むのに30~40秒くらいかかる
knownGenes <- UcscTrack(genome=genome, 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=genome, 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=genome, 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")

cpgIslands <- UcscTrack(genome=genome, chromosome=chr, track="cpgIslandExt", from=from, to=to,
                         trackType="AnnotationTrack", start="chromStart", end="chromEnd", id="name",
                         shape="box", fill="#006400", name="CpG Islands")

gcContent <- UcscTrack(genome=genome, chromosome=chr, track="GC Percent", table="gc5Base",
                       from=from, to=to, trackType="DataTrack", start="start", end="end", data="score",
                       type="hist", window=-1, windowSize=1500, fill.histogram="black", col.histogram="black",
                       ylim=c(30, 70), name="GC Percent")

conservation <- UcscTrack(genome=genome, chromosome=chr, track="Conservation", table="phyloP30wayPlacental",
                          from=from, to=to, trackType="DataTrack", start="start", end="end", data="score",
                          type="hist", window="auto", col.histogram="darkblue", fill.histogram="darkblue", 
                          ylim=c(-3.7, 4), name="Conservation")

snpLocations <-  UcscTrack(genome=genome, chromosome=chr, track="snp128", from=from, to=to,
                            trackType="AnnotationTrack", start="chromStart", end="chromEnd", 
                            id="name", feature="func", strand="strand", shape="box", 
                            stacking="dense", fill="black", name="SNPs")

axTrack <- GenomeAxisTrack()

idxTrack <- IdeogramTrack(genome=genome, chromosome=chr)

各Trackの中身をデータフレーム形式で閲覧したいときはas()関数を使う。

result

> as(cpgIslands,"data.frame")
  X.seqnames  X.start    X.end X.width X.strand X.feature X.group     X.id X.density
1       chrX 65931358 65932587    1230        *   unknown       1 CpG: 117         1

##読み込んだデータの可視化
各Trackが読み込めたら、最後にplotTracks()でplotする。

plot

plotTracks(list(idxTrack, axTrack, knownGenes,refGenes,ensGenes, cpgIslands,
                gcContent, conservation, snpLocations), from=from, to=to, showTitle=TRUE)

なお、ここで読み込んだデータは、実はパッケージ内にデモデータとして組み込まれているので、
UCSCに時間をかけてアクセスするのが面倒な場合は以下のコードを実行すれば
同　じ　よ　う　な　効　果　が　得　ら　れ　る

忙しい人のためのGvizサンプルコード

library(Gviz)
data(ucscItems)
plotTracks(list(idxTrack, axTrack, knownGenes,refGenes,ensGenes, cpgIslands,
                gcContent, conservation, snpLocations),
           from=from, to=to, showTitle=TRUE)

##くらべてみよう
###UCSCで書いた図

###configure

• image width: 800px
• Base Position: dense
• Chromosome Band: pack
• GC Percent: full
• UCSC Genes: full
• RefSeq Genes: dense (Hide non-coding genes: check)
• Ensembl Genese: full
• CpG Islands: full
• Conservation: full (Data view scaling: auto-scaling to data view)
• SNPs (128): dense

###Gvizパッケージで書いた図

enjoy!