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


![hgt_genome_408_262e50.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/12353/6f6670f8-b331-f7f1-af9b-9d484f2325a9.png)

[USCS Genome Browser](http://genome.ucsc.edu/cgi-bin/hgGateway)は、ゲノム中の遺伝子領域や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パッケージ](http://www.bioconductor.org/packages/2.11/bioc/html/Gviz.html)を用いるとできるようになるのでメモしておく。

##準備

```r:インストール
source("http://bioconductor.org/biocLite.R")
biocLite("Gviz")
```

```r:パッケージの読み込み
library(Gviz)
```

```r:表示領域の指定
genome <- "mm9"
chr <- "chrX"
from <- 65921878
to <- 65980988
```

##データの読み込み

```r:各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()`関数を使う。

```r: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する。

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

![ucscItems.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/12353/8f3a38c2-3558-9f81-3e21-4f604456224c.png)

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

```r:忙しい人のためのGvizサンプルコード
library(Gviz)
data(ucscItems)
plotTracks(list(idxTrack, axTrack, knownGenes,refGenes,ensGenes, cpgIslands,
                gcContent, conservation, snpLocations),
           from=from, to=to, showTitle=TRUE)
```

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


![hgt_genome_408_262e50.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/12353/6f6670f8-b331-f7f1-af9b-9d484f2325a9.png)


###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パッケージで書いた図
![ucscItems.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/12353/8f3a38c2-3558-9f81-3e21-4f604456224c.png)

enjoy!