jythonがmacbookにインストールできなかったのでとりあえずlinuxで。

といってもgcjだとエラーを吐くので、sunのjavaをインストール。 ここみて設定。alternativeコマンドを使うとjavaの共存ができるのね。いままで、シンボリックリンクを上書きしてたのでめんどいなーと思ってたけど、これだとらくちん。

jythonをインストールしたら、あとはopsinのjarを落としてきてクラスパスに通す。

>>> import uk.ac.cam.ch.wwmm.opsin as opsin
>>> opsin.NameToStructure().parseToCML("4-iodobenzoic acid").toXML() 

二行でIUPAC名がCMLに。 すばらしい。

jrubyの例もある。

ちなみにjrubyはmacbookにさくっと入って、この例の通りにやれば動いた。

perlで同じ事をやる場合にはInline::Javaを使ってやればいいけど、Javaのライブラリを有効に利用するのはJavaで実装された言語処理系がやっぱ楽だ。

元ネタはRadar / Spider Chars

五角形にしたかったのでrule of fiveにPolar Surface Areaを加えておいた。

#!/usr/bin/env python
from matplotlib.projections.polar import PolarAxes
from matplotlib.projections import register_projection
from pylab import *

class RadarAxes(PolarAxes):
    """Class for creating a radar chart (a.k.a. a spider or star chart)        
    http://en.wikipedia.org/wiki/Radar_chart
    """
    name = 'radar'
    # use 1 line segment to connect specified points
    RESOLUTION = 1

    def draw_frame(self, x0, y0, r):
        verts = [(r*cos(t) + x0, r*sin(t) + y0) for t in theta]
        return Polygon(verts, closed=True)

    def set_varlabels(self, labels):
        self.set_thetagrids(theta * 180/pi, labels)

    def get_axes_patch(self):
        x0, y0 = (0.5, 0.5)
        r = 0.5
        return self.draw_frame(x0, y0, r)

if __name__ == '__main__':

    register_projection(RadarAxes)
    N = 5

    theta = 2*pi * linspace(0, 1, N+1)[:-1]
    theta += pi/2
    labels = ['HBA', 'HBD', 'cLogP', 'MWT', 'PSA']
    rule_of_five = [10, 5, 5, 500, 140]
    desc = [12, 3, 3.6, 532, 160]
    desc_rate = [100*desc[i]/float(v) for (i,v) in enumerate(rule_of_five)]

    ax = subplot(111, projection='radar')

    ax.fill(theta, [100]*N)
    ax.fill(theta, desc_rate)

    for patch in ax.patches:
        patch.set_alpha(0.5)

    ax.set_varlabels(labels)
    rgrids((20, 40, 60, 80, 100))

    grid(True)
    show()

rule of fiveのようにある範囲内に収まっていること(超えるとリスク)というような指標を表すのにレーダーチャートは適しているんだろうか。つまり充足している事を示すような面積の表現はいいのかなぁ。あと、レンジが負になったりするのでそれもどうかと思う。

再考の余地はあるな。

TODO:多変量がある決まったレンジ内に収まっているかどうかを視覚的に捉えやすい表現手段を探す。

XYZでreadした分子(CC)もAddHydrogensできない

>>> from openbabel import *
>>> conv = OBConversion()
>>> conv.SetInFormat("xyz")
True
>>> mol = OBMol()
>>> conv.ReadFile(mol, "test.xyz")
True
>>> mol.NumAtoms()
2  
>>> mol.AddHydrogens()
True
>>> mol.NumAtoms()
2

でも分子中のそれぞれの原子を指定して

mol.AddHydrogens(atom)

としてやればきちんと水素が付加されることがわかったので、とりあえずスクリプト中でループをまわせばいいかな

mixiでpubchemからSMILESを抜き出すのは？みたいなエントリがたってたのだけど消えちゃったみたい。

なんか、ここは宿題まるなげちゃうわーみたいな厳しいコメントついてたからかな。まぁああいったコメント書くヒトはすっきりするのかもしれんけど、見るほうにとっては情報量ゼロのゴミなんだよなー。まだ、課題まるなげっていう情報のほうが情報量的に有意義。

ちょっと考えたのでここに書いておく。

XMLをBeautifulSoupで

import urllib2,sys
from BeautifulSoup import BeautifulSoup

cid = sys.argv[1]
url = 'http://pubchem.ncbi.nlm.nih.gov/summary/summary.cgi?cid=%s&disopt=DisplayXML' % cid

opener = urllib2.build_opener()
html = opener.open(url).read()
soup = BeautifulSoup(html)

print soup.findAll('pc-infodata_value_sval')[-2].string 

xmlがなんかfindしにくいので配列の要素指定に-2とかやって良くない香り。この程度のだったらxmlじゃなくてSDFを正規表現でいじるな。

import urllib2,sys,re

cid = sys.argv[1]
url = 'http://pubchem.ncbi.nlm.nih.gov/summary/summary.cgi?cid=%s&disopt=DisplaySDF' % cid

html = urllib2.build_opener().open(url).read()
p = re.compile('<PUBCHEM_OPENEYE_CAN_SMILES>\n(.+)\n')
m = p.search(html)
print  m.group(1)

biopythonとかでもいけそうな気がするし、urllib2+openbabelの組み合わせでも良いかもしれない。

perl,python,Rでオームな解析をするための本。ツールやデータベースの説明が主なので、ハウツーな感じではなくて、インフォマティクス側からみたバイオロジーとかケミストリーのサービスとかツールのレビューに近い感じかも。

オープンソースで始めるゲノム・プロテオーム・メタボローム解析
樋口 千洋
オーム社 / ￥ 6,090 ()
在庫あり。

個人的には5章のケモインフォのとこと、9章のRの章が面白かった。

SDFからGASTON実行できるようにしたので、結果を解釈するためにsmilesに変換するコードを用意した。

import re
import openbabel as ob

def convert(gf):
    obc = ob.OBConversion()
    obc.SetOutFormat('smi')
    mol = ob.OBMol()

    for l in gf.split('\n'):
        if len(l) > 0 and l[0] == 'v':
            a = mol.NewAtom()
            atomic_num = int(l.split(' ')[2])
            a.SetAtomicNum(atomic_num)

        elif len(l) > 0 and l[0] == 'e':
            begin_atom_idx = int(l.split(' ')[1]) + 1
            end_atom_idx   = int(l.split(' ')[2]) + 1
            bond_order     = int(l.split(' ')[3])
            b = mol.AddBond(begin_atom_idx, end_atom_idx, bond_order)

        elif len(l) > 0 and l[0] == '#':
            title = l.split(' ')[1]
            mol.SetTitle(title)

    return obc.WriteString(mol)

if __name__ == '__main__':
    txt = open('gaston.out').read()
    p = re.compile('#.+?(?=(#|$))',re.S)
    m = p.finditer(txt)

    for ss in m:
        print convert(ss.group())[:-1]

実行結果

NC  11
NCC 11
NCCC    11
NCC(=C)C    11
NCCC=C  11
NCC(=C)C=C  11
NCCC=CC 11
NCC(=C)C=CC 11
NCC=C   11

頻度の高いフラグメントのSMILESと、その分子のIDが出力される

それにしてもpythonの正規表現はややこしいなぁと、みんなのpythonを読み返して復習

みんなのPython
柴田 淳
ソフトバンククリエイティブ / ￥ 2,940 ()

正規表現文字列 -> コンパイル -> 正規表現オブジェクト -> マッチ処理 -> マッチオブジェクト

gSpanやGASTONの入力はラベル付きのvertexとedgeなので、openbabelでsdfを読み込んで、GASTONのinputに変換するものを作ってみた。

import openbabel as ob

def convert(sdf):
    obc = ob.OBConversion()
    obc.SetInAndOutFormats('sdf','smi')

    mol = ob.OBMol()
    next = obc.ReadFile(mol,sdf)
    molnum = 0
    while next:
        # mol
        print "t # %d" % molnum
        # atom
        for i,atom in enumerate(ob.OBMolAtomIter(mol)): 
            print "v %d %d " % (i,atom.GetAtomicNum())

        # bond
        for i,bond in enumerate(ob.OBMolBondIter(mol)): 
            print "e %d %d %d" % (bond.GetBeginAtomIdx()-1,bond.GetEndAtomIdx()-1,bond.GetBondOrder())

        mol = ob.OBMol()
        next = obc.Read(mol)

        molnum += 1
    return True

if __name__ == "__main__":
    sdffile = 'pubchem_sample.sdf'

    convert(sdffile)

SDFはpubchemから適当に選んだがCID: 16757835は外しといた。

./gaston 11 pubchem_data pubchem_out

GASTONを実行した出力の一部

# 14
t 163
v 0 6
v 1 6
v 2 1
e 0 1 2
e 1 2 1
# 12
t 164
v 0 6
v 1 6
v 2 1
v 3 1
e 0 1 2
e 0 3 1
e 1 2 1

アロマティックなボンドは3みたいに別のラベルを与える必要がある気はするが、、、

後は、GASTONの出力ファイルから構造を構築するものを用意すれば、頻出する部分構造を扱えるようになる。

今度の発表の掴みを考えてみた。

htsはアッセイという関数を引数にとりコンパウンドリストを返す高階関数と考えることができる。

def hts(assay,library):
    criteria = 50 # %inh
    return [c for c in library if assay(c) > criteria]

うーん、マニアすぎるか

def chemist(comp):
    return modify_compound(c)

これまたあれだな。

おー、このタイムライン分かりやすい。

インフラとしてはディスクリプターをどう発生させとくか（毎度毎度モデル構築のたびに発生させるのはだるいし、基本的に一回作っとけば良いのでコンパウンドをどっかに登録するタイミングでバックグラウンドでやればいい仕事）とデータへのアクセス方法（いまはRESTっぽいのがいいんじゃないかと思ってる）をきちんとやっておけば、僕らはモデル構築の試行錯誤にそのリソースのほとんどを投入できていいよね。

RESTful Webサービス
Leonard Richardson,Sam Ruby
オライリー・ジャパン / ￥ 3,990 ()
在庫あり。

あと、こういう状況だと、モデル構築能力的には常に最先端のアプローチを試せるような状況にしとかないとまずいから、自然とデフォルトの環境がRになるような気がするし、立ち位置もデータマイニング寄りになると思うんだよね。

そして、立ち位置がマイニング寄りになっちゃうと、ケモインフォマティクスとかバイオインフォマティクスってのは、新しいモデル化の方法を試すための一つのターゲットに過ぎなくて、プロジェクトのモデル化とか創薬ビジネスモデル構築とかのほうに興味が移っちゃうのだよなぁ。

というエントリをケモインフォクックブックにのせておきました。