12082010 chemoinformatics Python
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で調べるとかすんのかな？
ハイ・コンセプト「新しいこと」を考え出す人の時代
Emacsテクニックバイブル ～作業効率をカイゼンする200の技～
Emacs 辞典 (DESKTOP REFERENCE)
ツイッターで会社をPRする本 Twitter企業の活用例100
電子書籍と出版─デジタル／ネットワーク化するメディア
ポケットモンスター ホワイト
ポケットモンスター ブラック
おつまみ横丁―すぐにおいしい酒の肴185
