Bioinformatics版のProject Euler

RosalindというProject Eulerっぽいサイトがあったので少し遊んでみた。

楽しくbioinformaticsの問題を解きながらアルゴリズムを覚えていく感じなので、BioPythonを覚えるのにちょうどいいかなと思った。

GC Contentを計算するメソッドくらいあるだろうとは思うんだけど知らないので素朴に計算。リスト内包表記でmaxで取り出すのがいいかなと思ったけど、面倒くさいのでこれまたループで素朴に。

import sys
from Bio import SeqIO

def count_gc(seq):
    return float((seq_record.seq.count("G") + seq_record.seq.count("C"))) * 100 / len(seq_record.seq)

if __name__ == '__main__':
    fasta = sys.argv[1]
    highest_id = ""
    highest_content = 0
    for seq_record in SeqIO.parse(fasta, "fasta"):
        c = count_gc(seq_record)
        if c > highest_content:
            highest_id = seq_record.id
            highest_content = c

    print highest_id
    print "{0:.6f}".format(highest_content)

僕はBioPerlでバイオインフォマティクス関連のプログラミングをしてたので、BioPythonで配列解析をしたことは殆ど無いんですよね。

あらかじめクラスタ数を決めないでクラスタリングする方法(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;

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