RPy2を使ってPythonからRを呼び出してデータ解析
Pythonで書いたコードの中でRを呼び出す必要に迫られたので、備忘録を兼ねてやり方を記録しておきます。
この前から取り組んでいるPPiプロジェクトの一環で、タンパク質相互作用ネットワークを超幾何分析してスコアを記録しよう、みたいな話になった。
超幾何分析自体の話はbeyond the scopeすぎるので専門家に譲るとして、分析そのものはphyper
というR言語の関数でちゃちゃっとできる。
phyper(q, m, n, k, lower.tail = FALSE, log.p = FALSE)
q
:調べる蛋白質のモジュール内のリンク数
m
:モジュール内のリンクの総数
n
:モジュール外のリンクの総数
k
:調べるタンパク質の全リンク数
それで、この前Pythonで生成したファイルにRを走らせるのが一番単純な手法かもしれないけど、面倒だしスマートじゃないので、やはりPythonコードのなかでRを動かしたいところ。
少し調べると、PythonとRの連携はRPy2
というパッケージが人気の模様。
ちゃんとPython 3.xでも動作するということで、試してみる。
まずは$ sudo pip install rpy2
みたいな感じでローカルにRPy2をインストール。
それからPythonコードで読み込みます。
from rpy2.robjects import r
Rの関数の引数としてPython側の変数を使いたいものの、どうもうまくいかない。 調べてみると直接使わずにR側の新しい変数にポートして活用するらしい。
r.assign('RVariable', PyVariable)
という感じ。 これを引数に使う変数全部に対して行って、いよいよ関数にぶちこみます。
r('関数名(引数1,引数2,...)')
実際のコードでは、まずnamelist1
とwholelist
からそれぞれモジュール内のリンクの総数moduleedges
とタンパク質相互採用ネットワーク全体のリンク数wholeedges
を計算し、wholeedges
からmoduleedges
を引くことでモジュール外のリンクの総数nonmoduleedges
を得ました。
それから、モジュール内のタンパク質の名簿であるpnamelist
をループで回して、各タンパク質のモジュール内リンク数と総リンク数をそれぞれRにポートし、最後にそれをもとに超幾何分析を行ってdict型に収納、という手順を踏みました。
from rpy2.robjects import r moduleedges = len(namelist1) wholeedges = len(wholelist) nonmoduleedges = wholeedges - moduleedges from rpy2.robjects import r #超幾何分析結果を収納するdictを用意 hypergeo = {} #pythonの変数をRの変数にポート r.assign('module', moduleedges) r.assign('nonmodule', nonmoduleedges) #モジュール内の各タンパクについて超幾何分析 for x in pnamelist: r.assign('intralink', intralink[x]) r.assign('wholelink', wholelink[x]) p = r('phyper(intralink, module, nonmodule, wholelink, lower.tail = FALSE, log.p = FALSE)') #dictに分析結果を収納 hypergeo[x] = p[0]
あとは前回のようにcsvに出力して終わり。
僕以外にも統計分析にはRを使うけどそれ以外は全部Python、みたいな人は多いらしく、今回使ったRPy2のようにPythonの中でネイティブなRを走らせる以外にも、Rの大人気グラフ作成パッケージggplot2をPythonで再現する試みとかいろいろあるみたいです。
やっぱりデータ分析分野はR言語で書かれた遺産が多いので、うまく使っていきたいところです。
今回参考にしたウェブサイト
R: The Hypergeometric Distribution
http://rpy.sourceforge.net/rpy2/doc-2.1/html/