Help us understand the problem. What is going on with this article?

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

More than 1 year has passed since last update.

hgt_genome_408_262e50.png

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)

ucscItems.png

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

忙しい人のための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

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

enjoy!

lifull
日本最大級の不動産・住宅情報サイト「LIFULL HOME'S」を始め、人々の生活に寄り添う様々な情報サービス事業を展開しています。
https://lifull.com/
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
No comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
ユーザーは見つかりませんでした