Dr.Bonoの生命科学データ解析
オミックスデータの解析を主軸とした共同研究がいくつか走り出したので、生命科学のビッグデータに真正面から向き合う必要が生じました。たまたま１ヶ月前くらいに坊農さんの新しい本が上梓されていて、その内容が今の自分にぴったりすぎて、嬉しくてこの記事を書いてます。読んでよかった。
自分はRなどを使って時系列データの解析をおこなうくらいのプログラムは書ける人間です。でも遺伝子発現・蛋白などのビッグデータの解析は全くの素人。もちろんトランスクリプション、NGS、データベースといった言葉は知っていたけれども、エンリッチメント解析、Biomart、VCFなんてのはこの本で知りました。本書は１章「生命科学データ解析の歴史」、２章「生命科学分野のの公共データベース」、３章「データの形式とその取り扱い方」、４章「基本データ解析」、５章「実用データ解析」という構成になっていて、生命科学データの基本的な予備知識と手法について触れられてます。どの章も面白かったけど、個人的には５章がドンピシャでした。いろんなことがわかった気持ちになりました。全体を通して内容は初心者にはとてもわかり易く、特に本文の脇に無数に散りばめられている注釈が素晴らしい。FASTAは欧米ではファスト・エーと呼ばれるそうです。
来年やってくる大学院生の必読書が増えました。3,240円はお買い得だと思う。坊農さんありがとう。
休眠論文がでました。
自分にとって冬眠関連の論文として1報目となる論文が発表されました。
Sunagawa, G.A., Takahashi M, 2016.
Hypometabolism during Daily Torpor in Mice is Dominated by Reduction in the Sensitivity of the Thermoregulatory System.
Sci. Rep., 2016, doi: 10.1038/srep37011
論文のリンク
理研のプレスリリース
今回の論文の生物学的な発見やそのために行った開発については日本語プレスリリースや論文を見ていただくといいかと思います。ここでは、そこに書ききれなかった、僕自身がどんな思いで研究をすすめたか２点ほど書き留めておこうと思います。
1． オスの休眠パターンがおもしろい
今までの多くの休眠論文はマウスのメスを実験対象としてきました。メスのほうが安定した低代謝を示すからです。一方で睡眠研究を行っていた自分はメスには性周期があり、3〜4日に1回の頻度で過活動（睡眠時間が減る）が生じることを知っていたたので、個体間のバラツキを減らしたいと考えて、迷わずにオスを用いました。その結果、1個体内では時間的にとても短い期間しか休眠は検出できなかったのですが、個体間のバラツキはメスよりも減ったのではないかと思います（→メスは現在測定中）。いずれ場所を変えて性差についての議論もしたいのですが、それよりも僕の興味をひきつけたのは、オスの体温パターンの変動です。図の最上段は近交系Aのメスの休眠パターンです。絶食の後半にずどんと体温が20度台に落ちている部分が休眠です。一方で同じ近交系Aのオスの休眠は様相が全く違います（2段目）。波打ちながら体温が落ちています。発振しているように見えます。さらに、他の近交系BやCでもオスは同様のパターンが見られます。なんとかして、この振動の意味を突き止められないか、と思っております。これは興味だけではなく、この振動の応じて生体内の代謝制御機構の「何か」が振動している可能性があり、能動的低代謝の制御のために有用な情報だと思われるからです。どんな数理モデルを用いればよいのか、アドバイス募集中です。
2． 解析をベイス統計モデルで行った
2015年4月に高橋研に入ラボしてから本格的にベイス統計モデルの勉強をはじめました。時系列では、いまリバイス中の別の論文で視覚機能の評価のために階層モデルを用いて、stanを使ったMCMCサンプリングでパラメータ推定を行ったのが最初なんですが、この休眠論文はそれに続いてベイズ統計モデルを用いた論文です。具体的にはFig 3で用いている休眠の判定で、代謝のベースラインの推定に2次トレンドモデル（2階差分のトレンド項を有する状態空間モデル）を用いています。また、哺乳類の体温制御モデルのパラメーター推定に、こちらは単回帰問題ではありますが、stanを用いたMCMCサンプリングを用いています。このモデリングおよび解析を行うために、そして英語の表記法なども含めて次の書籍を大いに参考にしました。わかりやすくユーモアに富んだ作品で、最後まで飽きずに読み通せました。本当に感謝です。また、McElreath氏の講義の動画もあります： http://xcelab.net/rm/statistical-rethinking/
Richard McElreath (著)
Statistical Rethinking: A Bayesian Course with Examples in R and Stan
(Chapman & Hall/CRC Texts in Statistical Science) (英語) ハードカバー – 2016/2/19
ベイス統計の「青本」になるかな
豊田秀樹先生の『はじめての統計データ分析 ベイス的＜ポストp値時代＞の統計学』が届いたので、パラパラ読みをしました。
久保先生の緑本に対して、勝手に「青本」と呼ばせていだきます。
McElreathの『Statistical Rethinking』を読み進めてますが、英語ということもあり、なかなかあっさり一気読みは難しく、いま7合目くらいです。そうこうしてるうちに、ベイス統計入門の日本語の良書が登場したので紹介いたします。この青本もRからStanを使って、MCMCサンプリングでベイズ的アプローチを勉強しましょう、という内容のようですが、完全に日本語ですし、統計学の前知識一切不要を謳っているとおり、入門書としてわかりやすそうです。ただし、RやStanの文法などはほとんど記載がないとのことです。まずは前書きのQ&Aを読むべしという指示通りに、巻末にあるQ&Aを読んでみましたが、Q1「有意性検定をなぜ使わないのですか？」やQ2「有意か否かだけで話が済んだほうが楽です」の答えは、そのまま印刷して廊下に貼っておきたい気分になりました。
まだ、全部は読んでいないですが、次の休みに一気読みします。
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)) |
