Dr.Bonoの生命科学データ解析
オミックスデータの解析を主軸とした共同研究がいくつか走り出したので、生命科学のビッグデータに真正面から向き合う必要が生じました。たまたま１ヶ月前くらいに坊農さんの新しい本が上梓されていて、その内容が今の自分にぴったりすぎて、嬉しくてこの記事を書いてます。読んでよかった。
自分はRなどを使って時系列データの解析をおこなうくらいのプログラムは書ける人間です。でも遺伝子発現・蛋白などのビッグデータの解析は全くの素人。もちろんトランスクリプション、NGS、データベースといった言葉は知っていたけれども、エンリッチメント解析、Biomart、VCFなんてのはこの本で知りました。本書は１章「生命科学データ解析の歴史」、２章「生命科学分野のの公共データベース」、３章「データの形式とその取り扱い方」、４章「基本データ解析」、５章「実用データ解析」という構成になっていて、生命科学データの基本的な予備知識と手法について触れられてます。どの章も面白かったけど、個人的には５章がドンピシャでした。いろんなことがわかった気持ちになりました。全体を通して内容は初心者にはとてもわかり易く、特に本文の脇に無数に散りばめられている注釈が素晴らしい。FASTAは欧米ではファスト・エーと呼ばれるそうです。
来年やってくる大学院生の必読書が増えました。3,240円はお買い得だと思う。坊農さんありがとう。
Rにインストールしてるパッケージ
Rを3.3にあげたらパッケージを全部入れ直しになったので備忘録的にコードをかいておく。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 |
if ("h2o" %in% rownames(installed.packages())){ remove.packages("h2o") } pkgs<-c("ggplot2","MASS","glmmML","extrafont", "grid","reshape2","rstan","knitr", "scales","Cairo","pdfCluster","methods", "statmod","stats","graphics","RCurl", "jsonlite","tools","utils","reshape") for (pkg in pkgs){ if (!(pkg %in% rownames(installed.packages()))){ install.packages(pkg, dependencies = TRUE) } } install.packages("h2o", type="source", repos=(c("http://h2o-release.s3.amazonaws.com/h2o/rel-turing/3/R"))) library(extrafont) font_import() # run only once!! |
knitrをつかって教材作った
FASTER論文のセミナーを依頼されたのでknitrで講義資料を作ってみた → faster.tutorial.pdf
また、オリジナルのRnwファイルもアップしておきます → faster.Rnw
ということで、明日、生まれて初めて金沢に行ってきます！ さむいんだろうなぁ。。。
ggplot2のレジェンドをguidesで少しだけ改造
ggplot2の備忘録。
Rにて数万というサンプルを散布図にするとレジェンドを書き換えたくなることがある。
サンプルデータを用意する。
PNGに保存するときはコメントアウトを外してください。
|
1 2 3 |
.df<-function(mx, my, z) data.frame(x=rnorm(20000, mx), y=rnorm(20000, my), z) #png(file="sample%01d.png", width=1600, height=1200, res=200, pointsize=24) p<-ggplot(rbind(.df(-1,1, "Kobe"), .df(0,-1,"Osaka"), .df(1,1,"Tokyo"))) |
まずは普通に描画。
|
1 |
p+geom_point(aes(x=x, y=y, color=z)) |
点がでかすぎるので小さく。
|
1 |
p+geom_point(aes(x=x, y=y, color=z), size=0.2) |
色がこすぎるので透過に。
|
1 |
p+geom_point(aes(x=x, y=y, color=z), size=0.2, alpha=0.2) |
レジェンドが見えなくなっているのでguidesでデフォルトを上書きする。
|
1 2 3 |
p<-p+geom_point(aes(x=x, y=y, color=z), size=0.2, alpha=0.2) p+guides(color=guide_legend(override.aes=list(alpha=1, size=6))) #dev.off() |
詳しい使い方はSet guides for each scale.を参照してください。
Rとggplotでダブルプロットを作る
めちゃくちゃ久しぶりに備忘録。
Rにてダブルプロット。
ggplot2のgeom_rasterを使ってる。
こんなかんじ。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 |
library(ggplot2) .doublePlot<-function(z, start.datetime, one.day.size, max.count=NULL){ bin.size.in.secs<-(60*60*24)%/%one.day.size if (is.null(max.count)) max.count<-max(z) if (length(z)%%one.day.size!=0) z<-z[seq((length(z)%/%one.day.size)*one.day.size)] df1<-expand.grid(time=start.datetime+(seq(one.day.size)-1)*bin.size.in.secs, day=seq(length(z)%/%one.day.size)) df1$z<-z/max.count df2<-expand.grid(time=60*60*24+start.datetime+(seq(one.day.size)-1)*bin.size.in.secs, day=seq(length(z)%/%one.day.size)) df2$z<-c(z[seq(length(z)-one.day.size)+one.day.size], rep(NA, one.day.size))/max.count p<-ggplot(rbind(df1, df2)) p<-p+geom_raster(aes(x=time, y=day, fill=z)) p<-p+scale_x_datetime(breaks=start.datetime+60*60*6*(0:8), labels = scales::date_format("%H:%M"), expand=c(0,0)) p<-p+scale_y_reverse(breaks=c(seq(length(z)%/%one.day.size)), expand=c(0,0)) p<-p+scale_fill_gradient(low="white", high="black") } |
サンプルデータをかってにつくって、表示してみる。
|
1 2 3 4 5 6 7 |
z<-as.vector((rbind(matrix(rnorm(12*12*30, mean=1000, sd=1000), ncol=30), matrix(rnorm(12*12*30, mean=100, sd=500), ncol=30)))^2) z[sample(seq(12*24*30),3000)]<-z[sample(seq(12*24*30),3000)] start.datetime<-strptime("201311070600", format="%Y%m%d%H%M") one.day.size<-12*24 print(.doublePlot(z, start.datetime, one.day.size)) |
rglパッケージで3D散布図＆アニメーションをかく
久しぶりに備忘録。少し手を焼いたので。
Rで3D散布図。
rglパッケージは普通にインストール。
まずはデータ読み込み。ggplot2のmsleepデータセットを用いる。
|
1 2 |
library(rgl) data(msleep, package="ggplot2") |
plot3dコマンドで散布図作製。マウスでグリグリ動かすこともできる。
|
1 2 3 |
open3d(windowRect=c(0,0,512,512)) bg3d("white") plot3d(x=msleep$sleep_total-msleep$sleep_rem, y=msleep$sleep_rem, z=msleep$awake, col=msleep$conservation, xlab="NREM", ylab="REM", zlab="AWAKE") |
さらにグリグリうごかしている様子をアニメーションとして保存したい場合は、ImageMagickをインストールしておいて：
|
1 2 3 4 |
m<-par3d("userMatrix") Sys.setenv(PATH=paste("/opt/local/bin", Sys.getenv("PATH"), sep=":")) #change this depending on your env. movie3d(par3dinterp(userMatrix=list(m, rotate3d(m, pi/4, 1, 0, 0), rotate3d(m, pi/4, 0, 1, 0), rotate3d(m, pi/4, 0, 0, 1))), duration=6, movie="msleep", dir=".", convert="convert -delay 1x%d %s*.png %s.%s") |
NREMとREMとAWAKEの合計時間は24時かんなので見事に1つの平面にのっていますねｗ
open3dで作ったウィンドウを消すには：
|
1 |
rgl.close() |
TextmateにRのbundleいれた。MacBook Airにて。
S◯NYとの交渉がすっきりおわってからMBAをセットアップし始めた。生まれて初めてMacをさわる僕にとって何もかもが新しいわけだが、Windowsから移ってくる僕と同じような境遇の人達のためにも、自分のためにも、いくつかメモ。
最終的に何を入れたかはいつかアップするとして、最初に大きくつまずいたTextMateのbundleインストールについて備忘録。
TextMateのマニュアルの「5.7 Getting More Bundles」をみながら。
最初にsubversionのインストール。Xcodeをインストールしたためかすでに入っている。
~/.bash_profileに次の１行を追加：
|
1 |
export LC_CTYPE=en_US.UTF-8 |
ここが一番なやんだ。そもそも.bash_profileがどこにあるかわからなかった。そもそもそもそもこの段階では存在しないことも知らなかった。こちらをよみつつなんとか.bash_profileを作成。
一度bashを再起動して、次のコマンドを入力。
|
1 2 3 |
mkdir -p /Library/Application\ Support/TextMate/Bundles cd /Library/Application\ Support/TextMate/Bundles svn co http://svn.textmate.org/trunk/Bundles/R.tmbundle |
これで無事にTextMateでRスクリプトにも色がつくようになりました。
ggplot2でformatterを用いて軸ラベルを整形する
これも備忘録。いつも忘れてしまう。
まずはデフォルトのグラフ表示。
|
1 2 3 4 |
library(ggplot2) data(msleep) p<-ggplot(msleep, aes(x=sleep_total, fill=vore, group=vore))+geom_density(alpha=0.5) print(p) |
桁数を増やす。
|
1 2 |
fm1<-function(x) format(x, nsmall=4) print(p+scale_y_continuous(formatter=fm1)) |
Scientific表示。
|
1 2 |
fm2<-function(x) format(x, scientific=TRUE) print(p+scale_y_continuous(formatter=fm2)) |
組み込みのpercentもある。
|
1 |
<print(p+scale_y_continuous(formatter=percent)) |
ggplot2でviewportを用いて1枚の画像に複数のplotを配置する
毎回ggplotで複数のグラフを1枚の画像に収めるときに調べてる気がするので備忘録的にアップ。リゲス本を参照。
Source code:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 |
library(ggplot2) data(msleep) # remove rows including NAs in $vore msleep<-msleep[-which(is.na(msleep$vore)),] msleep$name<-factor(msleep$name) # make plots a<-qplot(data=msleep, x=name, y=sleep_total, fill=vore, geom="bar") +opts(axis.text.x=theme_text(angle=90, hjust=1, size =8, colour="grey50")) b<-qplot(data=msleep, x=bodywt, y=sleep_total, log="xy", color=vore) c<-qplot(data=msleep, x=brainwt, y=sleep_total, log="xy", color=vore) d<-qplot(data=msleep, x=brainwt, y=bodywt, log="xy", color=vore) # print plots grid.newpage() pushViewport(viewport(layout=grid.layout(2,3))) # 2x3 layout grid vplayout<-function(x,y) viewport(layout.pos.row=x, layout.pos.col=y) print(a, vp=vplayout(1,1:3)) print(b, vp=vplayout(2,1)) print(c, vp=vplayout(2,2)) print(d, vp=vplayout(2,3)) |
Cascading plots by viewport function under ggplot
Rの基礎とプログラミング技法
