igraphでクリーク探索
igraphはRとPythonのbindingがあるのにチュートリアルがほとんどなさげなので、ネットワーク分析の本があると手っ取り早く知れて良い。
クリーク探索
こんなネットワークからクリークを探してみる

g <-matrix(c(0,1,1,0,1,0,0,1,0,1,0,0,0,1,1,1,0,1,1,1,0,0,
0,1,0,1,1,1,1,0,1,1,0,1,0,0,0,1,1,1,0,0,0,1,
0,1,0,0,0),nrow=7,ncol=7,byrow=TRUE)
a <- graph.adjacency(g)
maximal.cliques(a)
[[1]]
[1] 0 1 2
[[2]]
[1] 1 6
[[3]]
[1] 0 2 4
[[4]]
[1] 2 3 4 5
[[5]]
[1] 3 6
追記 100818
Pythonの場合
a = [[0,1,1,0,1,0,0],
[1,0,1,0,0,0,1],
[1,1,0,1,1,1,0],
[0,0,1,0,1,1,1],
[1,0,1,1,0,1,0],
[0,0,1,1,1,0,0],
[0,1,0,1,0,0,0]]
g = GraphBase.Adjacency(a,mode=ADJ_UNDIRECTED)
g.maximal_cliques() # [(0, 1, 2), (1, 6), (0, 2, 4), (2, 3, 4, 5), (3, 6)]
g.largest_cliques() # [(2, 3, 4, 5)]
化合物の類似性ネットワークのために閾値とedgeの数をプロットする。
昨日の化合物の類似性ネットワークのやつは何がメリットなのかなぁと考えてみたけど、一つにはedgeの数を自由に決められるあたりなんじゃないかと。あんまり低い閾値だとなんでもかんでもつながって何見てるかわからんし、逆に厳しい閾値にしちゃうとシングルトンばかりになっちゃうし。
というわけで、プロット欠かせてそれ見て閾値決める必要はあるんだろうなと。最終的にはeddgesの数とか、ノードあたりに生えているedgeの数の平均とかを統計的に処理して自動で閾値を決めるのが楽なんだろうけど、、、(一般的な方法あるのかな?)。X-meansに倣って、BICとかAIC使ってみようかなぁ。
で、プロット描いた。最初に書いたコードがナイーブすぎて遅くていらっときたので書き直した。まだかなり遅いけど、耐えられるのでこれでよしとした。
import pybel
mols = list(pybel.readfile("sdf", "pc_sample.sdf"))
fps = [x.calcfp() for x in mols]
similarity_matrix = []
for i in range(len(fps)):
for j in range(i+1,len(fps)):
similarity = fps[i] | fps[j]
similarity_matrix.append(similarity)
threshold = 0.0
print "threshold\tnodes"
while threshold < 100:
t = threshold / 100.0
new_matrix = [e for e in similarity_matrix if e > t]
nodes = len(new_matrix)
print "%2.1f\t%d" % (threshold, nodes)
similarity_matrix = new_matrix
threshold += 0.1
あとはRで
library(ggplot2)
nodes <- read.table("xxxx",header=T)
qplot(threshold,edgess, data=nodes, log="y")
dev.copy(device=png,filename="xxxxxx")
dev.off()

node数が1158なのでedgeの数が10000から100あたりの間でちょうどいい閾値を探したいんだけど、これだけじゃわからん。やっぱシングルトンの数とかノードあたりに生えているエッジ数の平均とか数えないとだめかな。
PypeR
the Journal of Statistical SoftwareにPypeRというpythonからRを使うパッケージの論文が出てた。
Most Bioconductor tools are presented as R packages. Therefore, R can be a good complement to Python in bioinformatics computation.
メモリの使用量が低く抑えられる感じ?RPyやRPy2とかと比較している
でも、なんかpythonの中でRの構文埋め込んでるみたいな感じだなぁ。PypeRってのはそういうことなのかな。
Bioinformatics Programming Using PythonMitchell L. Model
Oreilly & Associates Inc / ¥ 5,584 ()
在庫あり。
Building Bioinformatics Solutions: With Perl, R and MysqlConrad Bessant,Ian Shadforth,Darren Oakly
Oxford Univ Pr (Txt) / ¥ 4,163 ()
在庫あり。
R Programming for Bioinformatics (Chapman & Hall/Crc Computer Science & Data Analysis)Robert Gentleman
Crc Pr I Llc / ¥ 6,532 ()
通常2~5週間以内に発送
Rでfoldlとfoldr
Rにおける高階関数はCommon Higher-Order Functions in Functional Programming Languagesにあるように一通りそろっている
Reduce uses a binary function to successively combine the elements of a given vector and a possibly given initial value. Filter extracts the elements of a vector for which a predicate (logical) function gives true. Find and Position give the first or last such element and its position in the vector, respectively. Map applies a function to the corresponding elements of given vectors. Negate creates the negation of a given function.
foldl,foldrも簡単
foldl <- function(f,x,i){Reduce(f,x,i)}
foldr <- function(f,x,i){Reduce(f,x,i,right=TRUE)}
実行してみる
> foldl("-",list(1,2,3),1) # ((1-1)-2)-3
[1] -5
> foldr("-",list(1,2,3),1) # 1-(2-(3-1))
[1] 1
あとRの高階関数は大文字で始める規則にでもなってんのかな?
RのNAMESPACEファイル
Rでパッケージ作るときにパッケージ変数として扱いたいものを
pit_dir <- file.path(Sys.getenv("HOME"),".pit")
pit_config <- file.path(pit_dir,"pit.yaml")
pit_profile_file <- NULL
と先頭に書いておきたいのだけどNAMESPACEファイルを作成すると
> library(pit)
要求されたパッケージ yaml をロード中です
> pit.switch("test")
以下にエラー pit.switch("test") :
'pit_profile_file' に対するロックされたバインディングの値は変更できません
とでる。ちなみにpit.switchはこれ
pit.switch <- function(name=NULL){
if(is.null(name)) name <- "default"
pit_profile_file <<- file.path(pit_dir, paste(name, ".yaml", sep=""))
config <- pit.config()
ret <- config$profile
config$profile <- name
write(as.yaml(config),pit_config)
if(name != ret) write(paste("Profile switch to",name,"from",ret), stderr())
return(ret)
}
環境をとってきてunlockBindingすればいいんだろうなぁと思うんだけどやり方が分からなくて、
pit_dir <- file.path(Sys.getenv("HOME"),".pit")
pit_config <- file.path(pit_dir,"pit.yaml")
#pit_profile_file <- NULL
assign("pit_profile_file",NULL,pos =".GlobalEnv")
とやって、直接大域環境に変数をセットしてるんだけど、これって本当はどうするのがいいんだろうか?
追記10.07.08
環境オブジェクトを作るという方法もあるそうなので、こっちでやってみた。
pit_dir <- file.path(Sys.getenv("HOME"),".pit")
pit_config <- file.path(pit_dir,"pit.yaml")
conf<-new.env()
conf$pit_profile_file <- NULL
こっちのほうが精神衛生上よろしいかも。それにしてもRの環境ってなんだろか?Lispみたいな感じで考えとけばいいのかな?あとでソースコードを探ってみよう。
参考
Shiga.R
Shiga.Rお疲れ様でした。内容盛りだくさんで面白かったのと、色々な人に出会えて有意義でした。というかこの三日間は充実してた。
plotter.Rはpit使わなくても、usernameとpassword直接入力すれば動きます。twitpicのAPIを使っているので、画像をポストしつつtwitterに投稿するようになっている(デモの時にはtwitterにはポストされなかったけど)。
帰りの新幹線でRいじってたので、irisをMDSとspeで次元縮約したのをポストしてた。
require("RCurl")
require("pit")
profile=pit.get("twitter.com",require=list(username="username",password="password"))
pweet <- function(status){
tempfile <- tempfile(pattern = "file", tmpdir = tempdir())
dev.copy(device=png,filename=tempfile)
dev.off()
twitpicreq <- "http://twitpic.com/api/uploadAndPost"
params = list(media=fileUpload(filename=tempfile),username=profile$username,password=profile$password,message=status)
twitpicresponse <- postForm(twitpicreq,.params = params)
}
あと、R版pit
R CMD INSTALL pit_0.8.tar.gz
でインストールできます。ただ、pit.setするときにエディタが起動するように鳴っているのだけど、terminalでRを起動しているときにはいいのだけど、それ以外の時(RのコマンドインターフェースとかESS)ではエディタが立ち上がらずおかしなことになるので、そこのとこをちょっとどうにかしないといけない。
あと、NAMESPACEというファイルを作ったら<<-で変数の代入が出来なくて、面倒になってファイルを削除してしまったのだけど、今日話を聞いてたらなんとなくやり方がわかったのであとでなおす。
Population Stochastic Modelling
違うことをしていたらPSM:PSM: Non-Linear Mixed-Effects modelling using Stochastic Differential Equationsという面白そうなパッケージを見つけたのでサンプルで遊んでみた。

集団から個人をサンプリングしてきてシミュレーションできるという理解でいいのかな?
Rのパッケージの internalな関数の文書化
R CMD checkすると、ドキュメントの存在しない関数は警告がでるのだけど、パッケージ内部で使うための関数も引っかかってしまう。
if a package pkg contains user-level objects which are for “internal” use only, it should provide a file pkg-internal.Rd which documents all such objects, and clearly states that these are not meant to be called by the user.
gridのRdを読めと書いてあったので読んでみる(R-2.11.0/src/library/grid/man/grid-internal.Rd)。
\name{grid-internal}
\alias{grid.legend}
\alias{grid.multipanel}
\alias{grid.panel}
\alias{grid.strip}
\alias{is.unit}
\alias{layoutRegion}
\alias{viewport.layout}
\alias{viewport.transform}
\alias{layout.heights}
\alias{layout.widths}
\alias{layout.torture}
\title{Internal Grid Functions}
\description{
Internal Grid functions
}
\details{
These are not to be called by the user (or in some cases are just
waiting for proper documentation to be written :).
}
\keyword{ internal }
こんな感じ
Rのパッケージを作る
先週はR版のpit(Config::Pit)をつくっていて、昨日今日でパッケージ化してみたので感じたことのメモ。
パッケージのスケルトンの作成はpackage.skeletonで行える。その際にRdフォーマットというTeXっぽい文書の雛形が作成される。パッケージのインストールはR CMD INSTALLで行うが、その際にドキュメントがちゃんと存在するかチェックされる。Rdフォーマットは R CMD Rd2dviでdviに出力できて、これはさらにpdfに容易に変換できる。Rのpdfがちゃんとしているのはそのため。
Rのドキュメントの充実度の秘密はこのあたりにあるのかと思った。Rはドキュメントのキチンと書かれていることに慣れきっていて、もしドキュメントがいい加減な(ましてや、Exampleのない)パッケージは使うときに相当イライラしますからな。
あと、R CMD checkすると、ドキュメントの存在しない関数は警告がでる。
* checking for missing documentation entries ... WARNING
Undocumented code objects:
pit_config pit_dir pit_profile_file pit.config pit.load
All user-level objects in a package should have documentation entries.
See the chapter 'Writing R documentation files' in manual 'Writing R
Rのrequireとlibrary
Loading and Listing of Packagesから。
library(package) and require(package) both load the package with name package. require is designed for use inside other functions; it returns FALSE and gives a warning (rather than an error as library() does by default) if the package does not exist.
面白そうなので、試してみる。
まずはRでパッケージを作る準備。package.skeletonは便利
f <- function(x,y) x+y
package.skeleton(list=c("f"), name="testpkg")
f.Rの先頭で存在しないパッケージを読み込む
library(dummypkg) # or require(dummypkg)
f <-
function(x,y) x+y
libraryのほうはエラーになってインストール出来ないが、
$R CMD INSTALL testpkg
* installing to library ‘/Library/Frameworks/R.framework/Resources/library’
* installing *source* package ‘testpkg’ ...
** R
** preparing package for lazy loading
Error in library(dummypkg) :
'dummypkg' という名前のパッケージが見当たりません
ERROR: lazy loading failed for package ‘testpkg’
* removing ‘/Library/Frameworks/R.framework/Resources/library/testpkg’
requireではwarningが出るが先に進む。
$ R CMD INSTALL testpkg
* installing to library ‘/Library/Frameworks/R.framework/Resources/library’
* installing *source* package ‘testpkg’ ...
** R
** preparing package for lazy loading
Loading required package: dummypkg
library(package, lib.loc = lib.loc, character.only = TRUE,\
logical.return = TRUE, 中で警告がありました:
'dummypkg' という名前のパッケージが見当たりません
ネットワーク分析 (Rで学ぶデータサイエンス 8)


Statistical Bioinformatics: with R
Rの基礎とプログラミング技法
確率微分方程式―入門から応用まで