2012-01-29
修論提出したらやりたいこと
| 
「俺、この戦争が終わったら、○○するんだ......」リストの何番目かに載ってる懸案。思う様甘い物を食べに行きたいとか日本酒飲みに行きたいとか書店巡りしたいとか温泉浸かりたいとか自転車のタイヤに空気入れたいとか髪切りたいとかいろいろあるけども、それらはひとまず措いておく。
英語論文執筆用のコーパスを整備したい
わたしは英語ネイティブではないので、英語を書くときは基本「英借文」をすることになります。コーパスが云々というのは、要するにその英借文用のデータベースを作って、英借文のプロセスをシステム化しましょうという試み。（理解が間違ってたらすみません）
え、修論執筆前にやっておけよ、という感じですが、そこまで頭と手が回らなかった...
別の言い方をすると、Incremental PubMed/MEDLINE Expression Search: inMeXes
何を使うか→CasualConcとかどうか
MacOSXでのコーパス作成の強い味方。に思える。
CasualConc J
どうやって使えばいいかというと、論文PDFをテキストファイル化したものをCasualConcに食わせてあげればOK。
こんな感じで。
目的語の文脈中の頻度などの分析、keyword in contextという分野らしい。
今月か来月に参加する予定の温泉開発合宿であれこれ試してみたくなってきた。
まずは提出...そして発表...
PDFのテキスト化はどうするのっと
Xpdfに同梱されているpdftotextを使いましょう。
ただし、macportsで入るXpdfにはpdftotext入らなくなったらしい。
なので以下からバイナリでインストール。そうすると使える
Xpdf: Download
pdftotext hoge.pdf
のように打つとhoge.txtを生成してくれる。
pdftotext -enc UTF-8 hoge.pdf
のようにファイルエンコーディングを指定してあげたほうがいいかもしれない。
話はずれるけど
河本健先生、今年の生命科学夏の学校にお呼びできないか...
英語教師は理系に学ぼう：河本健 （編）(2007) 『ライフサイエンス論文作成のための英文法』羊土社
ううむ。
2011-12-31
つくばRのかくも長き不在、またはなぜつくばRは13ヶ月の時を経て再び開催されるに至ったかそしてありがとうさようなら2011年
そろそろid:myopommeさんのSweaveエントリ群への返答記事を書きたいし、今日は大晦日なので一年を振り返るエントリも書こうと思ったのだけれど、その前に書いておくべき事があったことを思い出したので。
ぼくたちとつくばRの419日戦争
去る2011年11月12日、フリーな統計解析・可視化環境であるRをテーマに、つくば周辺(半径40万キロくらい)をターゲットとした勉強会であるところのTsukuba.R #9が開催された。参加者17人、うちスピーカーが8名と、東京開催にしては小規模だがid:a_bickyさん@gentlementatsu君id:gepuro君が初発表してくれるなど、密度の濃い回であった。
宣伝してくれた方、Ustや会場で参加・発表してくれた方に厚く御礼申し上げたい。
いまさらながら、発表スライドやブログ記事へのリンクがまとめられている Tsukuba.R #9 のWikiページはこちらである。
また当日の様子は
- というわけでTsukuba.R#9@東京大学に参加してきた & 発表してきた - 糞ネット弁慶
- Tsukuba.R #9 終了 - 最近物忘れが激しいだれかさんの雑記帳
- Tsukuba.R #9で「Rデータフレーム自由自在」を発表してきました - あらびき日記
- Tsukuba.R #9 - 壊れた計算機
- 第９回 Tsukuba.R #TsukubaR #9 - Togetter
を参照いただきたい。
ところで今回のTsukuba.R、前回が2010年9月19日だった事と考え合わせると実に419日、月にして約13ヶ月もの時を経ての開催であった。
Tsukuba.R#8 おひらき！
Tsukuba.R終了。乙でした #tsukubar
> t1 <- as.Date("2010-9-19") > t2 <- as.Date("2011-11-12") > difftime(t2, t1, units="days") Time difference of 419 days
一体なぜかくも長くの間つくばRは開催されなかったのだろうか。
なぜ#8と#9の間がここまで空いたのか
Tsukuba.R #9は本来3月に開催されるはずであった。
が、震災の影響を鑑み一旦延期、その後も何度か開催を待つ声はあったのだが、
Tsukuba.Rを中心となって運営していた野郎どもの予定がなかなか合わず11月まで来てしまった。
というより、「そろそろやろうぜ」を有言実行し、カレンダーとにらめっこして場所抑えて告知する「主催者」が不在だったのだと思う。
こう書くと他人事のようなので書きなおすと、
勉強会の主催なりRに関するブログを書いたりするのに必要なある種の情熱のようなものが、まず自分(wakuteka)自身から枯渇していた。ずっと無気力だったとかそういうわけではないけれど、他に心を占めていたことがあったり院での研究が迷走しかけたりしてRに関して何か発信していくことに興味を失いかけていた。
そんなこんなで延ばし延ばしになっていたTsukuba.Rの開催を後押ししてくれたのは、id:dichikaやid:isseing333といったTokyo.Rで出会った人たちからのムチャ振り応援だった。一緒にごはんを食べている時に、「そういえばTsukuba.Rって次いつやるの」「ね、年内には......」「そっかそっかそれは良い、じゃあTokyo.Rで告知しよう」「えっ」という感じで見事に退路を断っていただき、
@yokkuns まだLT枠（5分）って空いてますか？一人LTをお願いしたい方がいるのですが。
2011-09-20 13:51:05 via Silver Bird to @yokkuns
@dichika いけますよ！
2011-09-20 13:51:29 via YoruFukurou to @dichika
@yokkuns じゃあ発表者 @wakuteka さんで枠だけとっておいてもらえますか？タイトル等は彼からおって連絡がいくと思いますので！
2011-09-20 13:53:19 via Silver Bird to @yokkuns
9月のTokyo.R#17で、LTをさせてもらったことでid:gepuro君が発表に名乗りを上げてくれたり、id:tsutatsutatsutaが乗ってくれてKashiwa.R#1につながったりもした。
興味を共有できる場は楽しい
そんなこんなでなんとなく再確認したのは、何かを面白がっている人の話を聞きにいくのは楽しいという事だった。
本を紐解けないほど学習のモチベーションが下がっていたとしても、勉強会に行きさえすれば（足を運ぶくらいの余裕さえあれば）眠っていた興味が起こされてくるものだなというのを実感した。
ひとりでなにかを学ぼうとしても、往々にしてその情熱は持続できないことがある。
それは単に意思が弱いとかそういう理由で片付けてしまうべきではなくて、自分の中の情熱の炎が弱まっていたら、近い興味を持っている、似たようなことを面白がっている人に会いに行けさえすればある種のやる気はまた湧いてくるものなんじゃないか。
ざっくりまとめると、みんな勉強会やろうぜ！（行こうぜ！） ということで。
興味と情熱を継続するための無理のない仕組み
最後にR Advent Calendarがとてもとても面白かった参加して良かったと書いて今日の記事を締めることにする。発起人@kos59125ありがとう。
興味と情熱を持続するために、人に会い勉強会に参加しに行くのは実はわりと骨だ。時間もお金もかかるしある程度の勇気がいる。言いたいこと聞きたいことをやりとりするだけのコミュニケーション力も要る。気分なりテンションなりが外向きにイケイケの時はいいけれど、内向きでいたいような時はハードルが高かったりするかもしれない。
要は、興味を共有できる「場」は同じ時間同じ場所に集まらなくても醸成できるのではないかって思うのだ。
そういう意味でこのR Advent Calendarは新鮮で刺激的で楽しかった。
12月に限らずなにか理由をつけて年に何回でもやりたいような気さえする。来年2012年、あと数分でやって来るけど、こういう感じでオンラインでも興味を共有できる「場」に加わりたいなぁ。
結局一年を振り返る記事になってしまった。来年もよろしくお願いいたします。
参考資料
http://lab.sakaue.info/wiki.cgi/JapanR2010?page=%CA%D9%B6%AF%B2%F1%C8%AF%C9%BD%C6%E2%CD%C6%B0%EC%CD%F7
http://pear.php.net/package/Calendar/
2011-12-24
Merry chRistmas!
## http://d.hatena.ne.jp/Zellij/20111205/p1 ## MerryChristmas.R library(grid) par(bg="blue4") plot.new() plot.window(ylim=c(-4,4),xlim=c(-3,3)) curve(sqrt(1-x^2)+(x^2)^(1/3), -1, 1, lwd=11,col="#FF0000",add=T) curve(-sqrt(1-x^2)+(x^2)^(1/3), -1, 1, lwd=11,col="#FF0000",add=T) axis(1,pos=0,col="grey",col.axis="#C0C0C0") axis(2,pos=0,las=2,col="#C0C0C0",col.axis="#C0C0C0") x <- stats::runif(50) y <- stats::runif(50) rot <- stats::runif(50, 0, 360) text <- c("\115\145\162\162\171\040\103\150\162\151\163\164\155\141\163\041") grid.text(text, x=x, y=y, rot=rot, gp=gpar(fontsize=10, col="#C0C0C0") )
静かな森に雪が降るよなアニメーションとかやりたかったけどまたいつかの機会に。
よく考えたらハートを描く必然性はなかった！！メリークリスマス！！！
参考
## Try it! ## install.packages("animation") library(animation) demo(Xmas2)
追記
id:tor_ozakiがワンライナーでgridパッケージ使わずに書いてくれたよ！ あと、こっちの方程式使ったほうがなんだかグッと来る（個人的に）
plot(0,typ="n",xli=c(-100,100),yli=c(-100,100));`for`(i,1:50,{t<-1:360;points(16*sin(t)^3,13*cos(t)-5*cos(2*t)-2*cos(3*t)-cos(4*t),co="red");text(runif(1,-100,100),runif(1,-100,100),lab="Merry Xmas",srt=runif(1,-90,90),family=c("Osaka"),col="grey")})
2011-12-06
文芸的な、あまりに文芸的な (R Advent Calendar 2011)
| 
一 「ドキュメント」らしいドキュメントの有るプログラム
僕は「ドキュメント」らしいドキュメントの無いRのプログラムを不完全なものとは思つてゐない。従つて「ドキュメント」らしいドキュメントの有るプログラムばかり書けとも言はない。第一僕のプログラムも大抵は「ドキュメント」を欠いてゐる。絵画は解説がなくとも成り立つ。それと丁度同じやうにプログラムはドキュメントがなくとも成り立つものである。（僕の「ドキュメント」と云ふ意味は単に「コメント」と云ふ意味ではない。）若し厳密に云ふとすれば、全然「コメント」のない所にさえプログラムは成り立つであらう。従つて僕は「ドキュメント」のない簡潔なプログラムにも勿論尊敬を表するものである。
しかし或解析の見通しを良くし、再利用性を高めるのは決してコードの短さではない。況や実行にかかる数秒程度の余分は評価の埒外にあるはずである。更に進んで考へれば、ソースコードの有無さへもかう云ふ問題には没交渉である。僕は上にも書いたように「ドキュメント」の無いソースコードを、――或は「ドキュメント」らしいドキュメントの無いソースコードを不完全なものとは思っていない。しかしこうでないRプログラミングも存在し得ると思ふのである。
二 如何にしてドキュメントを書くかに答ふ
次に僕は如何にしてRで「ドキュメント」らしいドキュメントのあるプログラムを書くかという問いに答へる責任を持つてゐる。試みにMASSパッケージのcatsというデータセットを用ゐて説明する。
昔々あるところに144匹の猫が居た。何でも薄暗いじめじめした所でニャーニャー泣いているうちに気づけば心の臓を抜き取られその重さを量られていたという。猫の犠牲に敬意を表し、いま心重量(Hwt)と体重(Bwt)の間に直線的な相関関係があるか、性別(Sex)で違いがあるかを見たいとする。
sample.R
library(lattice) library(xtable) data(cats, package = "MASS") lm1 <- lm(Hwt~Bwt*Sex, data = cats) lm1 print(xyplot(Hwt~Bwt|Sex, data = cats, type = c("p", "r")))
上のソースコードをRの統合開発環境であるRStudioで実行すると下のような画面が得られる。
めでたし。
ところが、データを解析して可視化すれば終わりというではない。
一週間、あるいは一ヶ月、一年後に同様の解析をしなければならない、ということはよくあるだろう。
8行程度のコードなら覚えていられるかもしれない。しかし実際の解析に用いるプログラムやデータは複数のファイルに複雑怪奇に入り組んでいることがある。時間が経つとどのデータファイルをどこに保存したか、どの関数のどのオプションを用いたか、そもそもどういう経緯でその解析を行ったのか分らなくなってしまうという事態を経験したのは僕だけではあるまい。
そこで、
- どのようなデータに対して
- どのような解析を行い
- どのような結果が得られたか
をこまめにまとめておくことが重要になってくる。例えば以下のように。
実物はここsample.pdf 
からダウンロードできるようにした。
驚くなかれ、このPDFは以下のコードから生成されたものである。
RのコードをLaTeXに埋め込んでPDFを生成するこのシステムは、WEBをもじってSweaveと呼ばれている。
sample.Rnw
\documentclass[a4paper]{article}
\usepackage[margin=20mm]{geometry}
\usepackage[noae]{Sweave}
\DefineVerbatimEnvironment{Sinput}{Verbatim}{
fontsize=\scriptsize,
fontfamily=courier,frame=single,framesep=2mm,
label={R Code},labelposition=topline} %numbers=left
\DefineVerbatimEnvironment{Soutput}{Verbatim}{
fontsize=\scriptsize,fontfamily=courier}
\DefineVerbatimEnvironment{Scode}{Verbatim}{fontshape=\small}
\begin{document}
<<echo=false,results=hide>>=
library(lattice)
library(xtable)
data(cats, package = "MASS")
@
\section{The Cats Data}
Consider the \texttt{cats} regression example from
Venables \& Ripley (2002)~\cite{venables2002modern}.
The data frame contains measurement of heart and body weight of \Sexpr{nrow(cats)}
cats (\Sexpr{sum(cats$Sex=="F")} female, \Sexpr{sum(cats$Sex=="M")} male).
A linear regression model of heart weight by sex and gender can be fitted in R
using the command
<<>>=
lm1 <- lm(Hwt~Bwt*Sex, data = cats)
lm1
@
Tests for siginificance of the coefficients are shown in Table~\ref{tab:coef},
a scatter plot including the regression lines is shown in Figure~\ref{fig:cats}.
\SweaveOpts{echo=false}
<<results=tex>>=
xtable(lm1, caption="Linear regression model for cats data.", label="tab:coef")
@
\begin{figure}[!h]
\begin{center}
<<fig=TRUE,width=12,height=6,echo=T>>=
print(xyplot(Hwt~Bwt|Sex, data=cats, type=c("p", "r")))
@
\caption{The cats data from package MASS.}
\label{fig:cats}
\end{center}
\end{figure}
This PDF file have been generated on \today
with \Sexpr{version$version.string} on a
\Sexpr{sub("_"," ",version$platform,)}.
\bibliographystyle{unsrt}
\bibliography{sample}
\end{document}
いや、Sweaveは知っている、ただLaTeXを含めたTeX環境の導入の壁が高いのだとの声があるかもしれない。しかしRStudioの登場によって状況は大きく変わった。だってRStudioには "Compile PDF" ボタンがあるのだ。
ポチッとな。
三 僕
最後に僕の繰り返したいのは僕も亦今後側目もふらずに「ドキュメント」らしいドキュメントのあるソースコードばかり作るつもりはないと云ふことである。僕等は誰も皆出来ることしかしない。僕の行なっている解析が常にかう云ふプログラムを作ることに適してゐるかどうかは疑問である。のみならずかう云ふプログラムを作ることは決して並み並みの仕事ではない。僕のRを好むのはRはあらゆるプログラミング言語の形式中、最もいいかげん包容力に富んでゐる為に何でもぶちこんでしまはれるからである。若しラムダ式の完成した(((括弧)の)(国))に生まれていれば、僕は或は院生よりもプログラマになってゐたかもしれない。僕はいろいろの言語たちに何度も色目を使つて来た。しかし今になつて考へて見ると、最も内心に愛してゐたのは統計解析兼可視化環境のGNU S――わがR言語だつた。
芥川龍之介 文芸的な、余りに文芸的な
Getting Started with Sweave: R, LaTeX, Eclipse, StatET, & TeXlipse
文芸的プログラミング - Wikipedia
Sweave マニュアル訳 - RjpWiki
CRAN Task View: Reproducible Research
R Advent Calendar 2011 : ATND
四 蜜柑
Sweaveコードの解説、実運用上のハマりどころ、学会発表一週間前にBlog更新にかまけていて大丈夫か、文体が保てていないぞ、Tsukuba.Rのレポートはどうした等に関しては、はい、がんばります。
蜜柑食べたい。
myopomme
2011/12/09 22:52
Sweaveは解析をしていると多様しがちですが、LaTexのminipageやlongtableに対応していないがね〜。結局、WEBはWEBで書きつつ、LaTexのコードはLaTexのコードで書きつつ、texファイルをrubyで自動整形したり...仕上げに時間がかかりません？
myopomme
2011/12/10 09:27
Rjpwikiに良いまとめがあります。ぜひこのブログをベースにLaTexにもしたしみたい。テフニシャンなろうとおもった時期もありし過去。今一度。
2011-11-26
第6回 Rユーザ会 (1) データ解析のためのGUIフロントエンド
| 
（鈴木 了太さん）
統計数理研究所 R研究集会（第6回 Rユーザ会）
R AnalyticFlowは、フローチャートを描くことで
R analytic flow
高度なデータ分析を実現するソフトウェアです。
フローとして構造化されたプロセスは容易に再現・共有が可能。
過去に蓄積されたデータ分析プロセスをフローに変換することで、
ナレッジの再利用もスムーズです。
個人・法人を問わず無償でご利用いただけます。
KNIME, Rapidminerからヒント
データ分析は基本的に反復的なプロセス、フローチャートでまとまってると便利
フローチャート的UIの利点：全体を把握しやすい・共有しやすい
対話的操作の利点：探索的（色々試行錯誤して最適なコードを組み上げるにはコンソールが便利）
↓
Ctrl+Enterで、実行履歴をフローチャート化できる
感想
コードはRのコンソール等でガリガリ書いて、ある程度固まりができたらフローチャート化していくのがよさそう。
その他
キャッシュ機能（途中のステップで.RData生成)
付属のコードエディタ、候補も出せる
英語版も、問い合わせの9割は英語
one more thing...（開発中の機能）
裏でRが動いていることをより意識させないUI
Sweaveの先へ(reporting、html出力など)
edit()likeにExcelで編集する機能も（tsv出力してExcel開いてる)
2011-11-20
FPKM値のパターンが類似する遺伝子群をプロットするには
| 
RでNGS解析 ~ RNA-Seqによる遺伝子発現量データを手軽に可視化する - ここからゲノムと計算機すか
発現パターンが類似する遺伝子群を取得するには
> mySimilar <- findSimilar(cuff, "PINK1", n = 7) > annotation(mySimilar)[,1:5] gene_id class_code nearest_ref_id gene_short_name locus 1 XLOC_000048 <NA> <NA> RER1 chr1:2323213-2344010 2 XLOC_000055 <NA> <NA> TPRG1L chr1:3541555-3546692 3 XLOC_000085 <NA> <NA> NMNAT1 chr1:10003485-10045555 4 XLOC_000147 <NA> <NA> NECAP2 chr1:16767166-16786582 5 XLOC_000172 <NA> <NA> PINK1 chr1:20959947-20978003 6 XLOC_000175 <NA> <NA> NBPF3 chr1:21766630-21811392 7 XLOC_001389 <NA> <NA> <NA> chr1:21749600-21754300
自分自身(PINK1)も含まれている。これは呼び出し時に除けたほうがいいのでは...
expressionPlot()関数に結果を投げてプロット
> p <- expressionPlot(mySimilar, showErrorbars = F) > p <- p + geom_line(data = fpkm(mySimilar), + aes(x = sample_name, y = fpkm, group = gene_id, colour = gene_id)) > p
現状、単にexpressionPlotを呼び出しただけでは各遺伝子が判別できないので、自分で色分け等してあげないといけない。もしggplotでの色分けの仕方を知らなければ、パッケージさえ入れれば可視化出来る、というお手軽感は薄れてしまう。
発現パターンが類似する遺伝子群を返してくれるfindSimilar()関数の関数定義を読む
| 
cummeRbundパッケージのfindSimilar()関数について深追いしてみた。findSimilar()関数は、ある遺伝子と発現パターンの類似する遺伝子群を出してくれるということでとても便利に見える。
けど、中でどんな処理をしているかよくわからないのは気持ち悪いので、BioconductorのサイトからcummeRbundのソースコードを落として読んでみた。
この記事はそのときのメモです。
ダウンロード＆展開
$ wget http://www.bioconductor.org/packages/2.9/bioc/src/contrib/cummeRbund_1.0.0.tar.gz $ tar zxvf cummeRbund_1.0.0.tar.gz $ cd cummeRbund/R $ export COLUMNS=80 #ブログでの表示用に $ ls AllClasses.R methods-CuffDist.R methods-CuffGeneSet.R AllGenerics.R methods-CuffFeature.R methods-CuffSet.R database-setup.R methods-CuffFeatureSet.R tools.R methods-CuffData.R methods-CuffGene.R
findSimilar()の関数定義 in methods-CuffSet.R
#Find similar genes .findSimilar<-function(object,x,n){ # x はgene_idかgene_short_nameか、適当に作ったFPKM値のベクトルであればよい # 単一の値しか与えられない if(is.character(x)){ myGene <- getGene(object,x) sig <- makeprobsvec(fpkmMatrix(myGene)[1,]) }else if(is.vector(x)){ sig <- makeprobsvec(x) } allGenes <- fpkmMatrix(object@genes) allGenes <- t(makeprobs(t(allGenes))) compare <- function(q){ JSdistVec(sig,q) # 距離(類似度)を計算 } myDist <- apply(allGenes,MARGIN=1,compare) mySimilarIds <- names(sort(myDist))[1:n] mySimilarGenes <- getGenes(object,mySimilarIds) } setMethod("findSimilar",signature(object="CuffSet"),.findSimilar)
やってること
- 1番目の引数でデータセットを、2番目の引数で目的の遺伝子名 or 発現データのベクトルを与える
- fpkm値の行列を、すべてのセルの和が1になるよう調整（計算のため）
- 距離(類似度)を計算(Jensen-Shannonダイバージェンスを利用)
- 類似度でソートして上から順に取ってくる
- 取ってきたリストのIDをデータセットに投げてGeneSetを取得
mySimilarIdsの部分ではIDのベクトルが類似度順になっているのに、最後のGetGenesに投げる所で折角の類似度ランキングが消えてしまっている。これはもったいないのでは？ さらに言えば、類似度だけじゃなくて有意性を評価するための指標もつけるべきだと思う。有意性が評価できるかどうかは知らない。どうなのだろう。
あと、vignetteには類似度の計算に Jensen-Shannon distance を使っていると書いてあったけど、上のリンクのように Jensen-Shannon divergence って表記の方が一般的な気がする。
fpkmMatrixを適用するとこうなる
> head(fpkmMatrix(cuff@genes)) hESC Fibroblasts iPS XLOC_000001 7.23836e-01 16.4011 54.06720 XLOC_000002 0.00000e+00 0.0000 0.00000 XLOC_000003 0.00000e+00 0.0000 0.00000 XLOC_000004 1.20000e+06 22616.4000 0.00000 XLOC_000005 1.13903e+03 41.1644 944.30800 XLOC_000006 0.00000e+00 0.0000 9.00455
findSimilar()内で呼び出している関数群の関数定義 in tools.R
JSdist<-function(mat){ res<-matrix(0,ncol=dim(mat)[2],nrow=dim(mat)[2]) colnames(res)<-colnames(mat) rownames(res)<-colnames(mat) for(i in 1:dim(mat)[2]){ # ここの処理、applyで書いたら速くなる？ for(j in i:dim(mat)[2]){ a<-mat[,i] b<-mat[,j] JSdiv<-shannon.entropy((a+b)/2)-(shannon.entropy(a)+shannon.entropy(b))*0.5 res[i,j] = sqrt(JSdiv) res[j,i] = sqrt(JSdiv) } } as.dist(res) } JSdistVec<-function(p,q){ JSdiv<-shannon.entropy((p+q)/2)-(shannon.entropy(p)+shannon.entropy(q))*0.5 JSdist<-sqrt(JSdiv) JSdist } makeprobsvec<-function(p){ phat<-p/sum(p) #FPKM行列の、全体の和を1にする phat[is.na(phat)] = 0 phat } shannon.entropy <- function(p) { if (min(p) < 0 || sum(p) <=0) return(Inf) p.norm<-p[p>0]/sum(p) -sum( log10(p.norm)*p.norm) } makeprobs<-function(a){ colSums<-apply(a,2,sum) b<-t(t(a)/colSums) b[is.na(b)] = 0 b }
コメント
手法を表す数式はわかったけれど、なぜこの手法で類似度を計算するのか、ヘの理解は持ち越し。
2011-11-18
多重検定で有意水準を調整する際の手法を調べる
| 
多重検定で有意水準を調整する際の手法を調べる ＆ Storey作のソフトQvalueをRのパッケージ化したqvalueパッケージ の使い所を把握
蛋白質 核酸 酵素 2009 Vol.54 No.10 1307-15 大規模データの解析における問題点 DNAマイクロアレイによる遺伝子発現量の測定を例として http://www.mbsj.jp/admins/ethics_and_edu/PNE/5_article.pdf # FDR, qvalueについて
Proc Natl Acad Sci U S A. 2003 Aug 5;100(16):9440-5. Epub 2003 Jul 25. Statistical significance for genomewide studies. http://www.ncbi.nlm.nih.gov/pubmed/12883005 Storey JD, Tibshirani R. # FDRの定義はここを参照
N Engl J Med. 2001 Feb 22;344(8):539-48. Gene-expression profiles in hereditary breast cancer. http://www.ncbi.nlm.nih.gov/pubmed/11207349 # Storey(2003)で使っているデータの出典
(_ _) zzz...
使い方解説記事はあまり頻繁には書けない
2011-11-17
RでNGS解析 ~ RNA-Seqによる遺伝子発現量データを手軽に可視化する
R | |
R(Bioconductor)のパッケージcummeRbundについてちょいと調べました。2011-10-31に ver. 1.0.0がリリースされたばかりのまだ新しいパッケージです。
cummeRbundを使うと、RNA-SeqのデータからCufflinks(Cuffdiff)を使って推定した複数サンプルの遺伝子発現量を、手軽に比較・解析することができます。パッケージ名はカマーバンドと読みます。
11/22追記:
RでNGSデータ（RNA-Seqによる発現データ）を可視化する - 牧場の朝に項目を追加してまとめ直しました。
注：cummeRbundはBioconductor2.9以上、R 2.14以上で動作します。
(original photo uploaded by rjrgmc28)
実行環境
> sessionInfo() R version 2.14.0 (2011-10-31) Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit) locale: [1] ja_JP.UTF-8/UTF-8/ja_JP.UTF-8/C/ja_JP.UTF-8/ja_JP.UTF-8 attached base packages: [1] grid stats graphics grDevices utils datasets methods [8] base other attached packages: [1] cummeRbund_1.0.0 ggplot2_0.8.9 proto_0.3-9.2 reshape_0.8.4 [5] plyr_1.6 RSQLite_0.10.0 DBI_0.2-5
cummeRbundインストール
> source("http://www.bioconductor.org/getBioC.R") > getBioC("cummeRbund") Using R version 2.14.0 (R-devel), biocinstall version 2.9.3. Installing Bioconductor version 2.9 packages: [1] "cummeRbund" Please wait... also installing the dependencies ‘itertools’, ‘iterators’, ‘DBI’, ‘plyr’, ‘proto’, ‘RColorBrewer’, ‘digest’, ‘colorspace’, ‘RSQLite’, ‘reshape’, ‘ggplot2’ ... > options(width = 80) # ブログの横幅に収めるために
パッケージの読み込み
> library(cummeRbund)
パッケージがインストールされたディレクトリを取得する
> packagePath <- system.file(package = "cummeRbund") # パッケージのディレクトリを取得 > packagePath [1] "/Library/Frameworks/R.framework/Versions/2.14/Resources/library/cummeRbund" > dir(packagePath, recursive = T) # パッケージディレクトリ以下の全ファイルを表示 [1] "data/sampleData.RData" "DESCRIPTION" [3] "doc/cuffData_schema.pdf" "doc/cummeRbund-manual.pdf" [5] "doc/cummeRbund-manual.R" "doc/cummeRbund-manual.Rnw" [7] "doc/index.html" "extdata/cds_exp.diff" [9] "extdata/cds.diff" "extdata/cds.fpkm_tracking" [11] "extdata/cuffData_schema.sql" "extdata/gene_exp.diff" [13] "extdata/genes.fpkm_tracking" "extdata/isoform_exp.diff" [15] "extdata/isoforms.fpkm_tracking" "extdata/promoters.diff" [17] "extdata/splicing.diff" "extdata/tss_group_exp.diff" [19] "extdata/tss_groups.fpkm_tracking" "help/aliases.rds" [21] "help/AnIndex" "help/cummeRbund.rdb" [23] "help/cummeRbund.rdx" "help/paths.rds" [25] "html/00Index.html" "html/R.css" [27] "INDEX" "Meta/data.rds" [29] "Meta/hsearch.rds" "Meta/links.rds" [31] "Meta/nsInfo.rds" "Meta/package.rds" [33] "Meta/Rd.rds" "Meta/vignette.rds" [35] "NAMESPACE" "NEWS" [37] "R/cummeRbund" "R/cummeRbund.rdb" [39] "R/cummeRbund.rdx"
$packagePath/dataと$packagePath/extdata以下にサンプルデータ、
$packagePath/doc以下にマニュアルがあるようです。
以下、doc/cummeRbund-manual.Rのソースを淡々と実行
Cuffdiffプログラムの出力ファイルを読み込む
> # extdata以下にあった *.diff ファイルを読み込む > extdataPath <- paste(packagePath, "/extdata", sep = "") > dir(extdataPath) [1] "cds_exp.diff" "cds.diff" [3] "cds.fpkm_tracking" "cuffData_schema.sql" [5] "gene_exp.diff" "genes.fpkm_tracking" [7] "isoform_exp.diff" "isoforms.fpkm_tracking" [9] "promoters.diff" "splicing.diff" [11] "tss_group_exp.diff" "tss_groups.fpkm_tracking" > cuff <- readCufflinks(extdataPath) > cuff CuffSet instance with: 3 samples # iPS, Fibroblasts, hESC 400 genes 1203 isoforms 575 TSS 545 CDS 960 promoters 1725 splicing 696 relCDS
これで発現量データをRで可視化する準備が整いました。
データセットから目的遺伝子のデータを抽出するには
例として、ミトコンドリア膜電位の異常を監視し異常な膜電位をもつミトコンドリアが細胞の中に蓄積するのを防ぐと言われているタンパク PINK1 をコードしている PINK1 (PTEN-induced kinase 1）のデータを抽出します。
単一データを抽出するにはgetGene()関数を用います。ここでは遺伝子のIDではなく、遺伝子名で指定してみます。
> myGene <- getGene(cuff, "PINK1") > myGene CuffGene instance for gene PINK1 Short name: PINK1 Slots: annotation fpkm diff isoforms CuffFeature instance of size 2 TSS CuffFeature instance of size 2 CDS CuffFeature instance of size 2
FPKM値はfpkm()で取得します。
> fpkm(myGene) gene_id sample_name fpkm conf_hi conf_lo quant_status 1 XLOC_000172 Fibroblasts 2919.340 4002.960 1835.730 OK 2 XLOC_000172 hESC 693.465 813.869 573.062 OK 3 XLOC_000172 iPS 1598.040 2282.380 913.710 OK
FPKMとはFragments Per Kilobase of transcript per Million mapped readsの略で、発現量の指標の一つです。詳しくは Cufflinks RNA-Seq analysis tools - FAQ および Mortazavi, Williams, et al. Nature Methods, 2008 を、日本語資料としてはKen Osakiさんやyag_aysさんの記事を御覧ください。
目的遺伝子のサンプル別FPKMを視覚的に比べるには
今回のデータセットにはFibroblasts, hESC, iPSの３サンプルについてのデータが格納されているので、これらを視覚的に比較してみます。
> expressionPlot(myGene)
> expressionBarplot(myGene)
> expressionBarplot(isoforms(myGene)) # isoformがある場合は並べてplotできる
データセットから複数遺伝子のデータを抽出するには
> sort(sampleIDs) # 解析対象の遺伝子ID一覧 [1] "XLOC_000069" "XLOC_000089" "XLOC_000105" "XLOC_000115" "XLOC_000132" [6] "XLOC_000151" "XLOC_000158" "XLOC_000170" "XLOC_001240" "XLOC_001262" [11] "XLOC_001263" "XLOC_001265" "XLOC_001297" "XLOC_001339" "XLOC_001348" [16] "XLOC_001359" "XLOC_001363" "XLOC_001369" "XLOC_001370" "XLOC_001411" > length(sampleIDs) [1] 20 # 20 genes
複数遺伝子のデータを抽出するにはgetGenes()関数を用います。ここでは、上記の遺伝子IDからデータを取得してみます。
> myGenes <- getGenes(cuff, sampleIDs) # sampleIDsに対応するデータをCuffSetから抽出 > class(myGenes) [1] "CuffGeneSet" # CuffSetから遺伝子毎に抽出したデータはCuffGeneSetオブジェクトとなる attr(,"package") [1] "cummeRbund" > > myGenes.fpkm <- fpkm(myGenes) # FPKMを取得 > dim(myGenes.fpkm) [1] 60 6 # 3 samples x 20 genes = 60行のデータ (6列) > head(myGenes.fpkm) gene_id sample_name fpkm conf_hi conf_lo quant_status 1 XLOC_000069 Fibroblasts 2.98269e-01 8.94808e-01 0.00 OK 2 XLOC_000089 Fibroblasts 1.16767e+04 1.68167e+04 6536.69 OK 3 XLOC_000105 Fibroblasts 2.55768e+03 3.65431e+03 1461.05 OK 4 XLOC_000115 Fibroblasts 0.00000e+00 0.00000e+00 0.00 OK 5 XLOC_000132 Fibroblasts 1.88151e+03 2.69236e+03 1070.66 OK 6 XLOC_000151 Fibroblasts 1.02775e+00 2.48121e+00 0.00 OK
目的の遺伝子セットをFPKMでクラスタリングしたヒートマップを見るには
> csHeatmap(myGenes, clustering = 'both') # clustering = 'row' なら行ごと
目的の遺伝子セットの遺伝子名を見るには
> featureNames(myGenes) tracking_id gene_short_name 1 XLOC_000069 ESPN 2 XLOC_000089 PGD 3 XLOC_000105 MFN2 4 XLOC_000115 PRAMEF1 5 XLOC_000132 EFHD2 6 XLOC_000151 PADI1 7 XLOC_000158 <NA> 8 XLOC_000170 FAM43B 9 XLOC_001240 UBE2J2 10 XLOC_001262 <NA> 11 XLOC_001263 C1orf86 12 XLOC_001265 <NA> 13 XLOC_001297 SLC2A7 14 XLOC_001339 <NA> 15 XLOC_001348 SPATA21 16 XLOC_001359 SDHB 17 XLOC_001363 TAS1R2 18 XLOC_001369 AKR7A3 19 XLOC_001370 AKR7A2 20 XLOC_001411 IL22RA1
IDじゃなくて遺伝子名でヒートマップ書きたいときはこれを使えばよさそう。
データセット全体の概要を掴むためのいろいろな図が簡単に描けます
ggplot2を描画ライブラリに使っている模様
(see also: http://had.co.nz/ggplot2/stat_density.html)
Density plot
> csDensity(genes(cuff))
Box plot
> csBoxplot(genes(cuff))
Scatter plot
> csScatter(genes(cuff),"hESC","Fibroblasts",smooth=T)
他にできること
findSimilar()で、ある遺伝子と同じようなFPKMのパターンを示す遺伝子を取ってくる、ってこともできます。
類似度の計算にはJensen-Shannon distanceを使っているそうな。
ヘルプのPDFを見るには以下のコードを実行、またはpackagePath/docから地道に辿ってみてください。
> vignette("cummeRbund-manual")
ソースコードは、Bioconductorのサイトからダウンロードできます。もしかしたら実際の解析にはパッケージの関数そのままでは物足りなくなる可能性があり、そんなときはソースコードを読みパッケージの関数定義を参考にして、自分にとって必要な関数を自作することになるでしょう。
cummeRbundパッケージのオブジェクト一覧
> ls("package:cummeRbund") [1] "addFeatures" "annotation" "CDS" [4] "csBoxplot" "csCluster" "csDensity" [7] "csHeatmap" "csScatter" "csVolcano" [10] "diffData" "distValues" "expressionBarplot" [13] "expressionPlot" "featureNames" "features" [16] "findSimilar" "fpkm" "fpkmMatrix" [19] "genes" "getGene" "getGenes" [22] "getLevels" "isoforms" "JSdist" [25] "JSdistVec" "makeprobs" "makeprobsvec" [28] "promoters" "readCufflinks" "relCDS" [31] "samples" "shannon.entropy" "splicing" [34] "TSS"
調べた経緯
↓
このpostをふぁぼる
↓
捕捉される
@bonohu 木曜にやります
2011-11-15 14:03:22 via YoruFukurou to @bonohu
というわけでやってみた。コメント・トラックバック歓迎です。
2011-07-18
PubMed THE FIRST
| 
original uploaded by hobby_blog
ふと思い立って、PubMedID1番の論文を調べてみました。
Biochem Med. 1975 Jun;13(2):117-26.
Formate assay in body fluids: application in methanol poisoning.
Makar AB, McMartin KE, Palese M, Tephly TR.
http://www.ncbi.nlm.nih.gov/pubmed/1
タイトル1語めのFormateってのはFormic acid、つまりギ酸のこと。
メタノール
ホルムアルデヒド
ギ酸
メタノールはヒトの体内でホルムアルデヒドを介してギ酸へと代謝されます。
それを利用して（？）メタノール中毒の人の体液中のギ酸を検査しました、みたいな内容。多分。
First PubMed Entry | Whizbang
↑上のブログ記事のようなネタがid:0とかにイースターエッグ的に仕込まれていることを密かに期待していたけど、そんなことはなかった。
公共DBにネタを仕込まれていても困るけど、明らかにネタと分かるものなら....あるいは4月1日限定なら...と思ってしまう今日この頃です。
画像は正義の味方PubMed仮面、Locusta migratoria var.birota expeditum。*1
*1：適当です
2011-07-04
「Rを教えるにあたってのクイック・ガイド」(日本語ユーザ編)
R | | 
イメージングのデータをよりによって，エクセルで，カット・アンド・ペーストを繰り返し，徹夜して解析している4年生を見て，こりゃいかんと思い，Rのスクリプトを書いて，「ほら，こんなに作業が簡単になるでしょ?」と，Rをオススメするはずだった．
「Rを教えるにあたってのクイック・ガイド」 - 勝虫日記
ところが，反応はイマイチ．こんな難しいことやってられません，と言いたげな顔になってしまったので，いきなり負荷をかけすぎたかと反省．やはり，教えるのってムズかしい．
に触発されて。
いきなり脱線しますがこの勝虫さんのエントリで紹介されていたA Quick Guide to Teaching R Programming to Computational Biology Students
オススメ本
Rによる医療統計学
元論文でも紹介されていたIntroductory Statistics with R (Statistics and Computing)の邦訳。
使用しているデータセットは医療分野のものだけれど、紹介されている手法とその解説は医療分野に限らず使えるものばかり。
Rグラフィックス ―Rで思いどおりのグラフを作図するために―
こちらはR Graphics, Second Edition (Chapman &Hall/CRC The R Series)の邦訳。
Traditional作図法とTrelis作図法の違いがよくわかる。原書のサポートページで、サンプルコードと実行結果が公開されてるのがありがたい。
またこちらは紹介されていなかったけど、
Rプログラミングマニュアル (新・数理工学ライブラリ 情報工学)
のように用途で関数を逆引きできる本が一冊手元にあると便利です。
役に立つサイト
seekR - 統計分析プログラミング言語 R のための検索エンジン
@hiratake55さんによるR専門の検索エンジン。どちらもGoogleカスタム検索を使っているけれど、日本語で検索した場合はRseek.orgよりもRに特化した情報が得やすいよ！
R Graphical Manual
遺伝研 小笠原 助教によるグラフィック用のコードと実行結果が閲覧できるデータベース。
R Graph Galleryがタグクラウドで用途別の検索インターフェースを誇るとしたら、こちらはほぼ全てのCRANパッケージに収録されるデモ用ソースコードのプロット結果を閲覧できるという網羅性の高さが武器だよ！
Rを学ぶ上で直面するよくある問題
実行したい処理をRでどう書けばいいのか分からない→だれかに質問する
- RjpWikiのQ&A(初級者コース)などで質問してみる
- ライフサイエンスQA（β）
を使ってみる
- http://twitter.com/#!/search?q=%23rstatsj にハッシュタグ付きでpostする
そもそもコマンドを打つのが面倒そう
Excel上でRを動かすRExcel - もうカツ丼でいいよな
ExcelでR自由自在
を参考に学部生さんが使い慣れた(?)Excelを魔改造して、Rに慣れてもらうことも...
とここまで書いて私も力尽きてしまいました。「もっといい方法知ってるよ！」とかありましたらぜひトラックバック欲しいな！
2011-05-13
DICOM形式のファイルを扱うC++ライブラリ
| 
RでMRIデータを可視化する - 牧場の朝
What C++ library can I use to read pixels from DICOM images? - Stack Overflow
結果、DCMTKというライブラリが良さそうなので、次回はこれを調べてみることに。
DCMTKはDICOM形式のビューワ・解析ツールとして有名なOsiriXの実装にも使われているようです。
macportを使うと簡単にインストールでき、最新版の3.6.0ではなくdcmtkdcmtk3.5.4_p2が入ります。(2011.5.13現在)
DCMTKのwebsite: dicom.offis.de - DICOM Software made by OFFIS - DCMTK - DICOM Toolkit
DCMTKのdocument: OFFIS DCMTK: Main Page
なおdsr2htmlなどのユーティリティを使うにはDCMDICTPATHにdicom.dicのpathを設定する必要がありました。
export DCMDICTPATH=/opt/local/lib/dicom.dic
2009-10-20
R苦手の会 #6が開催されました。
R | | 
今回からは可視化についてということで、
http://www1.doshisha.ac.jp/~mjin/R/05.html
を予定していましたが、
最初に作図に関する大まかな枠組みを掴もうということで、
http://cse.naro.affrc.go.jp/takezawa/r-tips/r/47.html
を手がかりに、plot()からはじまり高水準作図関数、低水準作図関数の辺りを扱いました。
2009-10-16
dokuwikiで、URL中のdoku.phpを消す方法
| 
http://www.example.com/dokuwiki.php?start のようなURLを、
http://www.example.com/start と表示させることができるようになる。(※URLはダミーです。)
サンプルはこんな感じ -> http://wakuteka.info/R/start
やりかた
conf/local.phpに以下の行を追加する。
$conf['userewrite'] = 1;
.htaccessの該当部分を、以下のようにコメントアウトする
(デフォルトでは、.htaccess.distという名前になっているかもしれません。その場合は、ファイル名を変更して下さい。)
## Enable this to restrict editing to logged in users only ## You should disable Indexes and MultiViews either here or in the ## global config. Symlinks maybe needed for URL rewriting. #Options -Indexes -MultiViews +FollowSymLinks ## make sure nobody gets the htaccess files <Files ~ "^[\._]ht"> Order allow,deny Deny from all Satisfy All </Files> # Uncomment these rules if you want to have nice URLs using # $conf['userewrite'] = 1 - not needed for rewrite mode 2 RewriteEngine on # Not all installations will require the following line. If you do, # change "/dokuwiki" to the path to your dokuwiki directory relative # to your document root. RewriteBase / # If you enable DokuWikis XML-RPC interface, you should consider to # restrict access to it over HTTPS only! Uncomment the following two # rules if your server setup allows HTTPS. RewriteCond %{HTTPS} !=on RewriteRule ^lib/exe/xmlrpc.php$ https://%{SERVER_NAME}%{REQUEST_URI} [L,R=301] RewriteRule ^_media/(.*) lib/exe/fetch.php?media=$1 [QSA,L] RewriteRule ^_detail/(.*) lib/exe/detail.php?media=$1 [QSA,L] RewriteRule ^_export/([^/]+)/(.*) doku.php?do=export_$1&id=$2 [QSA,L] RewriteRule ^$ doku.php [L] RewriteCond %{REQUEST_FILENAME} !-f RewriteCond %{REQUEST_FILENAME} !-d RewriteRule (.*) doku.php?id=$1 [QSA,L] RewriteRule ^index.php$ doku.php
追記
RewriteEngineをonにするとdoku.phpがURLに表示されなくなり、RewriteBaseをonにするとカレントディレクトリのdoku.phpを読み込んでくれる、という設定です。
2009-02-05
syou先生の高速化講座
| 
一時頃に待ち合わせて、卒研の目指すところと書いてるプログラムと、統計処理について忌憚のない意見を頂いた。
中心的な話題は、サンプル数が少ないけど次元数の多い対象への繰り返し処理をいかに高速化するか。
んで、以下の事を教えていただいた。
- 相関係数とかは十分速く書いてあるはずなのでRではもうそんなに速くならない
- expand.gridで組み合わせを生成できる！
- 繰り返し処理の関数applyの並列版がある(snowパッケージ)applypar
- 民族ごとの相関の平均をただ出すのはまずいかも→全部のデータでやればおk
- 遺伝子の特性を主成分分析して次元圧縮←これやってみたい
- id:rindai87先生がやってるらしい！次元の呪いについて詳しく聞きたいぜ！
rindai872009-02-04
体調戻った
| 
日曜にアキバに行った後から体調を崩しておりましたが、ようやく復活しました。
よかったよかった。
今日は各自の卒研進捗を発表し合う選択授業があって、
「疾患関連遺伝子知りたいのに正常人の発現データだけ見ても意味なくね？」
とか
「マイクロアレイ発現データの相関係数ってそもそもそうやって見るん？」
などという質問が飛び、自分の実験条件の前提がしっかりしてないことに気づかされるなどしました。
その他にも、資料を作っていくうちに「相関係数の2乗平均値でソートするべきでは？」とかアイデアが出てくる。
人に説明しようとする過程で、自分の理解不足が浮かび上がってくるので、こういう場は積極的に参加していきたいですね。
明日もsyou6162にコメントをもらいに行くよていなのので、データを出しておかないと。
