matplotlibでPCA

web.pyを使ってサービスを提供できるようになったので
何か使えないかな〜と考えて見た。
任意のSDFがあったらPCAやって結果を返すとかあると
PCAってなに？みたいな人もspotfireで図が見えるしいいかもと思った。
PCRの方がより良いと教えていただいたのでこちらの実装も考えよう。
さて、自分がやるならRを使うのだが、web.pyに入れるとなると全部pythonがいい。
rpy2入れたしこれで。と思っていたら動かないで肩すかし、、どうも社内の環境のwinで構築したR2.10はだめなようだ
RDKitのモジュールはもうメンテしてないし、matplotlibがいいじゃないと言われたのでそれで取りあえず書いてみました。

rom rdkit import Chem
from rdkit.Chem import AllChem
from rdkit.ML.Data import Stats
from matplotlib import mlab
from pylab import *
import numpy
import sys

mols = [mol for mol in Chem.SDMolSupplier(sys.argv[1])]
fps = [AllChem.GetMorganFingerprintAsBitVect(mol,2) for mol in mols]

mat = []
for fp in fps:
    bits = fp.ToBitString()
    bitsvec = [int(bit) for bit in bits]
    mat.append(bitsvec)
mat=numpy.array(mat)
nrows = len(mat)
ncols = len(mat[0])

if nrows <= ncols:
    res = mlab.PCA(mat.T)
    proj = res.Y.T
else:
    res = mlab.PCA(mat)
    proj = res.Y

print "-"*20
print "ID","PC1","PC2","PC3"
for i in range(len(proj)):
    print i+1,proj[i][0], proj[i][1], proj[i][2]

print "Finish"
scatter(proj[:,0],proj[:,1])
show()

化合物数＞ビット（2048）の場合はそのままでいいみたいですが、
化合物数＜ビットの場合はそのままではエラーになるので転置します。
で最後にもう一回スコアを転置して第三成分までを出力。
手ものとのデータでやったらあまりぱっとしない感じ。
もちっとデータ集めてやってみる。

Advertisement