MCSを求めるには比較したい二つの化合物のModular Productを求めて、それの最大クリーク探索をすればいい。

Modular Productとは

• 二つのグラフの属性の対応が同じノードが両方共隣接していてエッジの属性が同じ場合
• 二つのグラフの属性の対応が同じノードが両方共隣接していない場合

にエッジを張ったグラフ。

これを作ってあとはigraphの最大クリーク探索メソッドで

import openbabel as ob
from igraph import *

sm1 = 'OCC(N)=O'
sm2 = 'O=CC(N)=O'

obc = ob.OBConversion()
obc.SetInAndOutFormats('smi','smi')

source = ob.OBMol()      
obc.ReadString(source,sm1)
source_atoms = [a for a in ob.OBMolAtomIter(source)]

target = ob.OBMol()      
obc.ReadString(target,sm2)
target_atoms = [a for a in ob.OBMolAtomIter(target)]

pairs = []
for i in range(len(source_atoms)):
    for j in range(len(target_atoms)):
        if source_atoms[i].GetAtomicNum() == target_atoms[j].GetAtomicNum():
            pairs.append((source_atoms[i].GetIdx(),target_atoms[j].GetIdx()))

p = []
for x in range(len(pairs)-1):
    for y in range(x+1,len(pairs)):
        u = pairs[x]
        v = pairs[y]
        if u[0] != v[0] and u[1] != v[1]:
            sb = source.GetBond(u[0],v[0])
            tb = target.GetBond(u[1],v[1])
            if sb != None: sbo = sb.GetBondOrder()
            if tb != None: tbo = tb.GetBondOrder()

            if (sb == None and tb == None) or (sb != None and tb != None and sbo == tbo):
                p.append((pairs.index(u),pairs.index(v)))

g = Graph(p)
mc = g.largest_cliques()
for c in mc:
    print [pairs[i] for i in c]

実際にはline graph(エッジとノードを入れ替えたグラフ)のModular Productをもとめてクリーク探索したほうが計算量がかなり減るらしい。

が、今回はやっていない。

Handbook of Chemoinformatics Algorithms (Chapman & Hall/Crc Mathematical and Computational Biology)

Chapman & Hall / ￥ 8,895 ()
通常1～3週間以内に発送

Chemoinformatics and Computational Chemical Biology (Methods in Molecular Biology)

Humana Press / ￥ 14,386 ()
近日発売　予約可

MCSSTanimotoで組んだネットワークのメリットはなんといってもエッジの属性にMCSの情報(smilesとかInChi)をのっけられることであろう。

import sys
from Gaston import MCSSTanimoto,Gaston
import openbabel as ob

obc = ob.OBConversion()
obc.SetInAndOutFormats("sdf", "smi")

input = "pc_sample.sdf"
mol = ob.OBMol()
next = obc.ReadFile(mol,input)
mols = [mol]

while next:
    mol.StripSalts(14)
    mols.append(mol)
    mol = ob.OBMol()
    next = obc.Read(mol)

mols = mols[1:11]

for i in range(len(mols)-1):
    for j in range(i+1,len(mols)):
        sys.stdout.flush()
        mcs = MCSSTanimoto(mols[i],mols[j])
        if mcs.score() > 0.2:
            print "%s\t%s\t%s\t%s\t%2.3f" % (mols[i].GetTitle(), "mcs", \
            mols[j].GetTitle(), obc.WriteString(mcs.mcs).split()[0], mcs.score())

でも、実際に組んでみたら、解釈が意外に難しそうなことに気がついた。MCSSTanimotoのアルゴリズム自体が問題なのかもしれない。

MCSよりはmolblaster的に細切れにしてみて最頻出のフラグメントからネットワークを組んでいくか、ありがち置換基の定義済みテーブルを使って類似性の高い置換基動どうしをつなぐという経験ベースのアプローチのほうが直感的かもしれないなぁ。

GastonでMCSを探索するのは色々無理があったようで、どうするかな困った困ったと探していたら、SMSDを見つけた。cdkには既に組み込まれているらしくてさっきビルドしてみたんだけど、よく考えたらCDKの使い方あんまよくわかんないや、、、ってことでGUIで遊んだ。

ちなみに、java6でビルドしないといけないので、osx(10.5)の場合は環境変数を設定した。

とりあえず、vildagliptinとsitagliptinのMCSをもとめた。

なんで5員環と6員環がマッチしてるんだろうか？あとでペーパーちゃんと読もうっと。

追記openbabelでのdiscussもあった

• Somebody interested in a VF2 based structure matcher for OB?
• gist

ノード間で類似性が最大となるような関係にひとつだけエッジを張るようにしてみた。

#!/usr/bin/env python
# -*- encoding:utf-8 -*-

import pybel

mols = list(pybel.readfile("sdf", "pc_sample.sdf"))
fps = [x.calcfp() for x in mols] 

for i in range(len(fps)-1):
    max_sim = 0
    max_name = ""
    for j in range(i+1,len(fps)):
        sim = fps[i] | fps[j]
        if sim > max_sim: 
            max_sim = sim
            max_name = mols[j].title
    print "%s\t%s\t%s\t%2.3f" % (mols[i].title, 'sim', max_name, max_sim)

ついでに類似度属性で色が変わるようにしてみた。

わからん。さっきのエッジに加えて類似性が高い化合物同士のエッジは加えてみた。

for i in range(len(fps)-1):
    max_sim = 0
    max_name = ""
    targets = []
    for j in range(i+1,len(fps)):
        sim = fps[i] | fps[j]
        if sim > max_sim: 
            max_sim = sim
            max_name = mols[j].title
        if sim > 0.7:
            print "%s\t%s\t%s\t%2.3f" % (mols[i].title, 'sim', mols[j].title, sim)
    if len(targets) == 0:
        print "%s\t%s\t%s\t%2.3f" % (mols[i].title, 'sim', max_name, max_sim)

ちょっとまとまりが出てきたけど、あれだなぁ。もとのsdfがそもそもあかんのかなぁ。

で、最小全域木(MST)で描くことも考えたけど化合物間の類似性の最小経路を求めてもなぁ、、、とか思ってヤル気がおきない。類似性でつないだネットワークのMSTの意義ってなんなんだろうなぁと。

なんやろかー、なんやろかー。論文もう一回ちゃんと読んでみようかなぁ。

• Finding key members in compound libraries by analyzing networks of molecules assembled by structural similarity.

openbabelにはMaximum Common Substructure(MCS)を求めるメソッドがないのでGastonを使った。

似たようなのは昔書いたけど扱いにくいのでクラスにしといた。ついでにMCSSTanimotoも求めるようにした。

MCSSTanimoto(A,B) = 　MCS /(A + B - MCS) #全てヘテロ原子数

という、MCS版のタニモト距離。例えばOc1ccccc1 とCc1ccccc1のMCSTanimotoは

6 /(7 + 7 - 6) = 0.75

スキャフォールドが決まっている場合には割と手軽に使えんじゃないかなぁ。まぁヘテロの性質無視して同一性だけで判断しているのでハロゲンの違いは割と似ているとかそういう補正は入れたほうがいいのかもしれないが。

import openbabel as ob
import os,re
from tempfile import mkstemp

class Gaston:
    def __init__(self,mols=[]):
        self.mols = mols

    def writefile(self,filename):
        """ convert file to gaston format"""
        f = open(filename, 'w')
        for (molnum, mol) in enumerate(self.mols):
            f.write("t # %d\n" % molnum)
            for i,atom in enumerate(ob.OBMolAtomIter(mol)):
                f.write("v %d %d\n" % (i,atom.GetAtomicNum()))
            for i,bond in enumerate(ob.OBMolBondIter(mol)):
                f.write("e %d %d %d\n" % (bond.GetBeginAtomIdx()-1,bond.GetEndAtomIdx()-1,bond.GetBondOrder()))

        return filename

    def readfile(self,file):
        txt = open(file).read()
        p = re.compile('#.+?(?=(#|$))',re.S)
        mols = p.finditer(txt)

        result = [self._gas2ob(gas.group()) for gas in mols]
        return result

    def _gas2ob(self,gf):
        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 mol

    def search(self,freq=None):
        if freq == None: freq = len(self.mols)

        (m,gasfile) = mkstemp()
        (n,output) = mkstemp()
        self.writefile(gasfile)
        os.system("/opt/local/bin/gaston %d %s %s > /dev/null 2>&1" % (freq,gasfile,output))
        mols = self.readfile(output)
        os.unlink(gasfile)
        os.unlink(output)
        return mols

    def mcs(self):
        substructures = self.search()
        mcs = substructures[0]
        for substructure in substructures[1:]:
            if substructure.NumAtoms() > mcs.NumAtoms():
                mcs = substructure
            elif substructure.NumAtoms() == mcs.NumAtoms():
                if substructure.NumBonds() > mcs.NumBonds():
                    mcs = substructure
        return mcs

class MCSSTanimoto:
    def __init__(self,mol1,mol2):
        self.mol1 = mol1
        self.mol2 = mol2

    def score(self):
        mcs = Gaston(mols=[self.mol1,self.mol2]).mcs()
        return float(mcs.NumAtoms()) / (mol1.NumAtoms() + mol2.NumAtoms() - mcs.NumAtoms())


if __name__ == '__main__':
    import sys
    obc = ob.OBConversion()
    obc.SetInAndOutFormats("smi", "smi")

    mol1 = ob.OBMol()
    obc.ReadString(mol1, "CCC1=CC=CS1")
    mol2 = ob.OBMol()
    obc.ReadString(mol2, "OCC1=CC=CS1")
    mol3 = ob.OBMol()
    obc.ReadString(mol3, "CCCC1=CC=CS1")

    gaston = Gaston(mols=[mol1,mol2,mol3])
    print ":::Frequent Structures:::"
    for m in gaston.search():
        sys.stdout.write(obc.WriteString(m))
    print ":::MCS:::"
    sys.stdout.write(obc.WriteString(gaston.mcs()))
    print ":::MCSSTanimoto:::"
    score = MCSSTanimoto(mol1,mol2).score()
    print score

実行結果

:::Frequent Structures:::
CC  3
CC=C    3
C(=C)C=C    3
C(=C)C=CS   3
CC(=C)S 3
CC=CC   3
CC(=CC)S    3
CC=CC=C 3
CC(=CC=C)S  3
Cc1cccs1    3
CC=CC=CS    3
CC=CS   3
CC=CSC  3
c1ccsc1 3
CC=C(SC)C   3
CC=CSCC 3
CCS 3
CCSC    3
CC(=C)SC    3
CCSC=C  3
CC(=C)SC=C  3
C=C 3
C=CS    3
C=CSC   3
C=CSC=C 3
CS  3
CSC 3
:::MCS:::
Cc1cccs1    3
:::MCSSTanimoto:::
0.75

あと実際にモジュール化する場合にはコマンドがどこにあるかわからないし、もしかしたたらインストールされてないかもしれないけど、そういう場合にどういう感じで処理すればいいんだろうか？whichで調べるとかすんのかな？

昨日の化合物の類似性ネットワークのやつは何がメリットなのかなぁと考えてみたけど、一つにはedgeの数を自由に決められるあたりなんじゃないかと。あんまり低い閾値だとなんでもかんでもつながって何見てるかわからんし、逆に厳しい閾値にしちゃうとシングルトンばかりになっちゃうし。

というわけで、プロット欠かせてそれ見て閾値決める必要はあるんだろうなと。最終的にはeddgesの数とか、ノードあたりに生えているedgeの数の平均とかを統計的に処理して自動で閾値を決めるのが楽なんだろうけど、、、（一般的な方法あるのかな？）。X-meansに倣って、BICとかAIC使ってみようかなぁ。

で、プロット描いた。最初に書いたコードがナイーブすぎて遅くていらっときたので書き直した。まだかなり遅いけど、耐えられるのでこれでよしとした。

import pybel

mols = list(pybel.readfile("sdf", "pc_sample.sdf"))
fps = [x.calcfp() for x in mols] 

similarity_matrix = []
for i in range(len(fps)):
    for j in range(i+1,len(fps)):
        similarity = fps[i] | fps[j]
        similarity_matrix.append(similarity)

threshold = 0.0
print "threshold\tnodes"
while threshold < 100:
    t = threshold / 100.0
    new_matrix = [e for e in similarity_matrix if e > t]
    nodes = len(new_matrix)
    print "%2.1f\t%d" % (threshold, nodes)

    similarity_matrix = new_matrix
    threshold += 0.1

あとはRで

library(ggplot2)
nodes <- read.table("xxxx",header=T)
qplot(threshold,edgess, data=nodes, log="y")
dev.copy(device=png,filename="xxxxxx")
dev.off()

node数が1158なのでedgeの数が10000から100あたりの間でちょうどいい閾値を探したいんだけど、これだけじゃわからん。やっぱシングルトンの数とかノードあたりに生えているエッジ数の平均とか数えないとだめかな。

chemoinformaticsでもネットワークアナリシスは重要だと思う。今はマイナーだけど今後精力的に研究されんじゃないかなぁと。

pybel使うとシミラリティベースのネットワークは簡単に作れる。

import pybel

mols = list(pybel.readfile("sdf", "pc_sample.sdf"))
fps = [x.calcfp() for x in mols] 

for i in range(len(fps)):
    ids = []
    for j in range(i,len(fps)): 
        if (fps[i] | fps[j]) > 0.5 and i != j: 
            ids.append(mols[j].title)
    print "%s %s %s" % (mols[i].title, "sim", " ".join(ids))

SIF(simple interaction format)形式なので、拡張子sifでcytoscapeに読み込まさせればOK。

あとでケモインフォクックブックにも載せておこう。

昨日のワークショップでは、面白い話が聞けてよかった。いくつかアイデアが浮かんだのでメモっておく。特に二次会の他愛のない話はブレスト的でよかった。

文献中の画像から分類

どういう解析をしたらいいかわからないとか、どういうグラフとかプロットが好まれるのかわからないという話があって、解析結果を含む（当該ソフトウェア名を含むとかで絞るのか？）論文を集めてきて、pdfから画像を抜き出し、プロット、チャートを画像類似検索とか組み合わせて、検索できるサイトがあれば便利なんじゃなかろうか？

R Graphical Manualの論文版みたいな感じ。

Patentと文献を自動的に解析

製薬企業が文献を出すのは、色々制約があって(故に制約企業)。たいていはプロジェクトがポシャったとき（またはゴールしたとき）と相場が決まっているので、patent中の実施例化合物名を会社名と紐付けて収集しておいてデータベース化しておいて、論文が出たらそっちのデータベースと照会するっていうのはどうかなと。パテントは結構そのプロジェクトのトーンがわかるでしょう(例えば実施例がやたらと多いパテントとかやっつけ仕事か？とか)。あと論文とかも学位対策か？とかもauthorと特許の出願者照らし合わせながらみればわかったりしませんかね。BMCLとか、「実験記録じゃねーのこれは？」というのが結構多くて、なんでこんなもん投稿すんのかなぁとか思ったりして背景が知りたかったりするのですよね。学位対策度87%とか出てくればdiscussionのとこだけ読んでたいしたことなければ読み飛ばして終了とかできるのでありがたかったりする。

というように、文献検索系のマイニングをビジネスというか組織ハックとか業界ハック的な側面から考えていくと色々面白そうなこともやれるよなぁと。

Rによるテキストマイニング入門
石田 基広
森北出版 / ￥ 2,940 ()
在庫あり。

昨日bayonでどんだけいけるかって話になって、そういえばあんまり大きい数のクラスタリングしたことないなと思ったので早速やってみた(白macbook)。

データはPubchemから最新の25万件分(Compound_45925001_45950000からCompound_46150001_46175000まで)をダウンロードしてきた。bayon用のデータセットを作るためのpythonスクリプトは以前書いたものを使った。

for i in *.sdf.gz ; do babel -imol $i -ofpt ${i%.sdf.gz}.fpt  -xh -xfFP2; done
for i in *.fpt ; do ../python/f2bayon.py $i > ${i%.fpt}.tsv; done

とりあえず25000件くらい

$ wc Compound_45925001_45950000.tsv
   24770 6696478 20005168 Compound_45925001_45950000.tsv

$ time bayon -n 100 Compound_45925001_45950000.tsv > Compound_45925001_45950000.cls

real    0m37.312s
user    0m36.378s
sys     0m0.402s

1分かからず終了。続いて25万件くらい

$ wc all.tsv 
  248232 74728370 222905963 all.tsv

$ time bayon -n 10000 all.tsv > all.cls

real    9m49.447s
user    9m4.833s
sys     0m8.037s

これは10分かからずに終了。

「スライドは全て英語にしといたほうが無難だな」と思った。日本のユーザー会といっても、日本で行われるユーザー会であって、日本語のユーザー会とは異なるのですな。でも、そのほうがいいのかも。

• chemoinformaticsにおいてaromacityはruleにすぎない! ってのは名言だと思った。
• javaはサポートしているからjruby使えとか、groovy使えとかScala使えというのは商用ソフトウェアのサポート体制として利にかなっているなと。
• 同様に.NETサポートするからIronRuby使えとかIronPython使えとかF#使えとか。でもF#はあなたの趣味でしょ？と思った。

あとはAmazon EC2でガバッと計算したときにbillionオーダーの化合物の計算が安価に速く出来ればとか。そういう事例はないかなと思ったり。中小製薬企業で中途半端な数のクラスタマシンをもつより、そういう方法で計算させる方向にシフトしていけば費用対効果もいいんじゃないかなと思うんだけど、どうなんでしょうね？