あらかじめクラスタ数を決めないでクラスタリングする方法(Affinity Propagation)

K-meansのように予めクラスタ数を指定すると、「そのクラスタ数は正しいの?」っていう疑問が浮かぶと思う。

「なんらかの統計値に基づいて適切なクラスタに分割して欲しい」そんな願いを叶えるのがAffinity Propagationというクラスタリングアルゴリズムである

exemplara(セントロイドとかクラスタ中心)になるべきパラメータ(responsibility)とクラスタメンバに属しやすさ(availability)を交互に更新していって収束させる手法なので、K-meansのような初期値依存性がないらしい。

クラスタ数は類似度行列の対角要素(自分との類似度)に依存する(デフォルトはmedian)のでここを変更するとクラスタ数も変わるんだけどね。

Scikit-learnではAffinity Propagationが実装されているのでsykのケミストリースペースを作ってクラスタリングしてみた。ちなみにスライドのPCAの説明は間違っていた(pca.fit(fps).transform(fps)としなければいけなかった)。

from rdkit import Chem
from rdkit.Chem import AllChem, DataStructs
from sklearn.decomposition import PCA
from sklearn.cluster import AffinityPropagation
from ggplot import *
import numpy as np
import pandas as pd

suppl = Chem.SDMolSupplier('syk.sdf')

fps = []
for mol in suppl:
    fp = AllChem.GetMorganFingerprintAsBitVect(mol, 2)
    arr = np.zeros((1,))
    DataStructs.ConvertToNumpyArray(fp, arr)
    fps.append(arr)

print len(fps)
pca = PCA(n_components=2)
x = pca.fit(fps).transform(fps)
af = AffinityPropagation().fit(x)

d = pd.DataFrame(x)
d.columns = ["PCA1", "PCA2"]
d["labels"] = pd.Series(af.labels_)
g = ggplot(aes(x="PCA1", y="PCA2", color="labels"), data=d) + geom_point() + xlab("PCA1") + ylab("PCA2")
ggsave("ap.png", g)

1404818907

RDKitでFatal Python error: PyThreadState_Get: no current thread

homebrewのRDKit入れたらこんなエラーに遭遇

$ python mkdesc.py 
Fatal Python error: PyThreadState_Get: no current thread
Abort trap: 6

調べたらboost関連らしい

brew rm boost
brew install --with-icu --build-from-source boost
brew uninstall rdkit
brew install rdkit

で動くようになったけど brew installで再度boostのインストールを始めたので2番目のboostのインストールは必要ないかもしれない

土曜日はpython-ggplotの話もします

ChEMBLからデータを取ってきて、活性が強いほど色が濃くなってドットのサイズも大きくなるようにしてみた。python-ggplotを使うと数行で描くことができるので調子いいですみたいな話をします。

あとはscikit-learnとrdkitで予測モデルを作ったりクラスタリングをしたりといった話もできればなーと考えているので、 興味があれば参加して下さい。

import pandas as pd
from ggplot import *
import numpy as np

d = pd.read_csv("syk.csv")
d["pIC50"] = 9 - np.log10(d["IC50"])
print d

p = ggplot(aes(x='MWT', y='ALOGP', color="pIC50", size="pIC50"), data=d) + geom_point()
#print p
ggsave("2dplot.png", p)

1404120520

scikit-learnで交差検定が地味に便利

train_test_splitっていうメソッドが便利。これでテストデータと訓練データを分けてくれる。

SVC with linear kernelでちょっと検定

>>> import numpy as np
>>> from sklearn import cross_validation
>>> from sklearn import datasets
>>> from sklearn import svm
>>> iris = datasets.load_iris()
>>> X_train, X_test, y_train, y_test = cross_validation.train_test_split(
...     iris.data, iris.target, test_size=0.4, random_state=0)
>>> X_train.shape, y_train.shape
((90, 4), (90,))
>>> X_test.shape, y_test.shape
((60, 4), (60,))
>>> clf = svm.SVC(kernel='linear').fit(X_train, y_train)
>>> clf.score(X_test, y_test)
0.96666666666666667

あるディレクトリにあるファイルのなかで一番新しく作られたものを調べる(Python)

globとosモジュールを使う

>>> import os, glob
>>> max(glob.iglob('*.txt'), key=os.path.getctime)
'list_sample.txt'

参考

ProductName みんなのPython 第3版
柴田 淳
SBクリエイティブ株式会社 / ?円 ( 2013-12-04 )


CheckiOはPython3の学習サイト

遊んでみたけど結構楽しい

checkio

pandasとggplotを学ぶための本

最近グラフはPandas+ggplotで描いていて、深夜のバッチ処理で毎日沢山のグラフを生成させている。

rpy2?

使ってないなーみたいな。

Python版ggplotはR版のAPIと同じものを提供することをゴールにしているのでggplotのクックブックを読んでおくといいことがあります。六角形の密度プロットとか早く実装されないかなー。

ProductName Rグラフィックスクックブック ―ggplot2によるグラフ作成のレシピ集
Winston Chang
オライリージャパン / 3570円 ( 2013-11-30 )


尚、下の本も参考になります。

ProductName グラフィックスのためのRプログラミング
H. ウィッカム
丸善出版 / 4200円 ( 2012-02-29 )


そして、データの操作はPandasで。

ProductName Pythonによるデータ分析入門 ―NumPy、pandasを使ったデータ処理
Wes McKinney
オライリージャパン / 3780円 ( 2013-12-26 )


機械学習はscikit-learnを使えばよい。

pychembldbを使うための手順

今週末の勉強会の内容です。ChEMBLダンプのimportに時間がかかるので、pythonラッパー弄りたい人はChEMBLを予め入れてから参加するといいでしょう。

当日行う手順を記しておきます。

Pythonのインストール

2.7.xを入れます。OSXは標準でインストール済みです。WIndowsの場合はここからインストーラをダウンロードしてインストールしてください。

pipのインストール

pipはPythonのパッケージ管理ツールです。

ez_setup.pyをダウンロードして実行します

sudo python ez_setup.py

その後get-pip.pyをダウンロードして実行

sudo python get-pip.py

pychembldbのインストール

最初にPostgreSQLのアダプタであるpsycopg2をインストールし、続いてORマッパーをインストールし、最後にChEMBLdbインターフェースをインストールします。

sudo pip install psycopg2
sudo pip install sqlalchemy
sudo pip install pychembldb

pychembldbの設定

pychembldbはデフォルトでMySQLを使うようになっているので.bashrcに以下のように記述しておきます。

export CHEMBL_URI="postgresql+psycopg2://localhost/chembl_17"

編集したら、新しいタブを開くかsource .bashrcをしてください(source効くのかな?)

(注)

例えばuserがpostgresでパスワードがsakeと設定している場合は

export CHEMBL_URI="postgresql+psycopg2://postgres:sake@localhost/chembl_17"

と記述します

使ってみる

python 対話環境を起動します

python

pythonコードを実行します(select pref_name from target_dictionary をpythonで実行)。

>>> from pychembldb import *
>>> for target in chembldb.query(Target).all():
...   print target.pref_name

ターゲット名がずらずらと表示されればOKです。

OSXにChEMBLを入れる

今週末の勉強会の内容です。ChEMBLダンプのimportに時間がかかるので、pythonラッパー弄りたい人はChEMBLを予め入れてから参加するといいでしょう。

当日行う手順を記しておきます。

Homebrewのインストール

パッケージ管理システムです(WIndowsユーザーは飛ばしてください)。

ruby -e "$(curl -fsSL https://raw.github.com/Homebrew/homebrew/go/install)"

こんだけ。

PostgreSQLのインストール

brewを入れてあればインストールは下のコマンドを打つだけ。(WIndowsユーザーはここから最新版のインストーラをダウンロードしてください)

brew install postgresql

インストール後にデータベースを初期化します

initdb /usr/local/var/postgres

サーバー起動

postgres -D /usr/local/var/postgres

別のターミナルから動いているか確認をします

psql -l

(エラーになる場合)

Lionだと標準でPostgreSQLがインストールされているようで、/usr/bin/psqlが存在します。homebrewでインストールする場合には/usr/local/bin/psqlなので、

which psql
/usr/local/bin/psql

とならない場合はPATHを通すか

/usr/local/bin/psql -l

と打ってください

環境変数に設定する場合は~/.bashrcに以下の行を追加します

export PATH=/usr/local/bin:$PATH

その後新しいターミナルを立ち上げた後、

psql -l

ChEMBLを入れる

ChEMBLのftpサイトから最新版(chembl_17_postgresql.tar.gz)をダウンロードします。

展開します

tar xvfz chembl_17_postgresql.tar.gz
cd chembl_17_postgresql

PostgreSQLにChEMBL用のデータベースを用意します

psql postgres
postgres=# create database chembl_17;
#Ctrl-D で抜けます

ダンプをインポートします

psql chembl_17 < chembl_17.pgdump.sql

数時間待ちます。(MBAの2013バージョンだと1.5hくらいかかったそうです)

僕のMBAではこのような結果でした。(ERRORは特に問題ありません)

ALTER TABLE
ALTER TABLE
REVOKE
ERROR:  role "postgres" does not exist
ERROR:  role "postgres" does not exist
GRANT

real    44m46.001s
user    0m10.153s
sys 0m13.130s

動作確認

$ psql chembl_17
psql (9.2.4)
Type "help" for help.

chembl_17=# select pref_name from target_dictionary;

と打ってターゲット名が表示されればきちんと動いています。

逆引きPandas (DataFrame編)

R逆引きハンドブックのMatrixの章をPandasで。

ProductName R言語逆引きハンドブック
石田基広
シーアンドアール研究所 / 4830円 ( 2012-01-26 )


作者の書いた本はオススメ

ProductName Pythonによるデータ分析入門 ―NumPy、pandasを使ったデータ処理
Wes McKinney
オライリージャパン / 3780円 ( 2013-12-26 )


データフレームを作成する

>>> pd.DataFrame([[1,2,3],[4,5,6]])
   0  1  2
0  1  2  3
1  4  5  6

データフレームに列名と行名を追加

>>> x = pd.DataFrame([[1,2,3],[4,5,6]])
>>> pd.DataFrame([[1,2,3],[4,5,6]])
   0  1  2
0  1  2  3
1  4  5  6
>>> x = pd.DataFrame([[1,2,3],[4,5,6]])
>>> x
   0  1  2
0  1  2  3
1  4  5  6
>>> x.index = list("ab")
>>> x.columns = list("cde")
>>> x
   c  d  e
a  1  2  3
b  4  5  6

データフレームの列にアクセスする

>>> x = pd.DataFrame([[1,2,3],[4,5,6]], index=list("ab"), columns=list("cde"))
>>> x
   c  d  e
a  1  2  3
b  4  5  6
>>> x["c"]
a    1
b    4
Name: c, dtype: int64
>>> x.c
a    1
b    4
Name: c, dtype: int64

Seriesが返ってくるのでメソッドもよべる

>>> x.c.mean()
2.5
>>> x.c.sum()
5

データフレームに要素を追加する

>>> x = pd.DataFrame([[1,2,3],[4,5,6]], index=list("ab"), columns=list("cde"))
>>> y = pd.Series([7,8,9], index=list("cde"))
>>> x.append(y, ignore_index=True)
   c  d  e
0  1  2  3
1  4  5  6
2  7  8  9

データフレームを結合する

rbind

::: python
>>> x
   0  1
0  1  0
1 -2  3
>>> x.append(x)
   0  1
0  1  0
1 -2  3
0  1  0
1 -2  3
>>> pd.concat([x, x], axis=0)
   0  1
0  1  0
1 -2  3
0  1  0
1 -2  3

cbind

::: python
>>> pd.concat([x, x], axis=1)
   0  1  0  1
0  1  0  1  0
1 -2  3 -2  3

データフレームから一部を抽出する

ggplotの組込みデータを使う。 wt>3のname,disp列を抽出

>>> from ggplot import mtcars
>>> mtcars[mtcars.wt > 3][["name","disp"]]
                   name   disp
3        Hornet 4 Drive  258.0
4     Hornet Sportabout  360.0
5               Valiant  225.0
6            Duster 360  360.0
7             Merc 240D  146.7
8              Merc 230  140.8
9              Merc 280  167.6
10            Merc 280C  167.6
11           Merc 450SE  275.8
12           Merc 450SL  275.8
13          Merc 450SLC  275.8
14   Cadillac Fleetwood  472.0
15  Lincoln Continental  460.0
16    Chrysler Imperial  440.0
21     Dodge Challenger  318.0
22          AMC Javelin  304.0
23           Camaro Z28  350.0
24     Pontiac Firebird  400.0
28       Ford Pantera L  351.0
30        Maserati Bora  301.0

データフレームの列を変換する

transform

>>> x = pd.DataFrame([(1,2),(3,4),(5,6),(7,8)])
>>> x
   0  1
0  1  2
1  3  4
2  5  6
3  7  8
>>> x[2] = x[0]+2
>>> x
   0  1  2
0  1  2  3
1  3  4  5
2  5  6  7
3  7  8  9

因子ごとに組み合わせを作成する

R expand.grid() function in Python

データフレームの列や行ごとの合計を求める

::: python
>>> x
   0  1  2
0  1  2  3
1  4  5  6
2  7  8  9
>>> x.sum(axis=0)
0    12
1    15
2    18
dtype: int64
>>> x.sum(axis=1)
0     6
1    15
2    24
dtype: int64

データフレームの列ごとに関数を適用する

>>> x = pd.DataFrame([(1,2),(3,4),(5,6),(7,8)])
>>> x.apply(sum)
0    16
1    20
dtype: int64

データフレームから欠損値を取り除く

>>> x = pd.DataFrame([(1,2),(3,np.nan),(np.nan,6),(7,8)])
>>> x
    0   1
0   1   2
1   3 NaN
2 NaN   6
3   7   8
>>> x.dropna()
   0  1
0  1  2
3  7  8

データフレームの列を分解、統合する

multiindex便利

>>> x = pd.DataFrame([(1,"Y"),(3,"N"),(6,"N"),(7,"Y")])
>>> x
   0  1
0  1  Y
1  3  N
2  6  N
3  7  Y
>>> y = x.groupby(1)
>>> y.head()
     0  1
1        
N 1  3  N
  2  6  N
Y 0  1  Y
  3  7  Y
>>> y.mean()
     0
1     
N  4.5
Y  4.0

縦長形式のデータフレームを横長形式に変換する

>>> x = pd.DataFrame([(1,0,1),(1,1,2),(1,2,0),(2,0,1),(2,1,1),(2,2,3),(3,0,0),(3,1,1),(3,2,2)], columns=["subject","time","conc"])
>>> x
   subject  time  conc
0        1     0     1
1        1     1     2
2        1     2     0
3        2     0     1
4        2     1     1
5        2     2     3
6        3     0     0
7        3     1     1
8        3     2     2
>>> x.pivot(index="time", columns="subject")
         conc      
subject     1  2  3
time               
0           1  1  0
1           2  1  1
2           0  3  2

データフレームを分割する

groupby

データフレームを並べ替える

>>> x = pd.DataFrame([(1,0,1),(1,1,2),(1,2,0),(2,0,1),(2,1,1),(2,2,3),(3,0,0),(3,1,1),(3,2,2)], columns=["subject","time","conc"])
>>> x
   subject  time  conc
0        1     0     1
1        1     1     2
2        1     2     0
3        2     0     1
4        2     1     1
5        2     2     3
6        3     0     0
7        3     1     1
8        3     2     2
>>> x.sort("time")
   subject  time  conc
0        1     0     1
3        2     0     1
6        3     0     0
1        1     1     2
4        2     1     1
7        3     1     1
2        1     2     0
5        2     2     3
8        3     2     2