28 08 2012 chemoinformatics Python Tweet
論文読んでたら、コンベンショナルなMMPの求め方が載ってたので、似たようなやり方を考えてみた。
要するに網羅的に結合を切断していって
{Scaffold; [fragment, fragment, fragment]}
っていうハッシュを作っていく。fragmentは共通のScaffoldを持つので配列中の任意のフラグメントのペアがMMPということになる。
僕の切断ルール
僕の切断ルール
- シングルボンド
- 環の一部ではない
- ボンドのどちらかの端がCであること
CHEMBL327743の結果
('CC[C@@H]([C@H](NC(=O)[C@@H]1CCCCN1CC(=O)*)/C=C/c1ccccc1)C', 'COc1cccc(c1*)OC', 'CHEMBL327743') ('CC[C@@H]([C@H](NC(=O)[C@@H]1CCCCN1C*)/C=C/c1ccccc1)C', 'COc1cccc(c1C(=O)*)OC', 'CHEMBL327743') ('CC[C@@H]([C@@H](/C=C/c1ccccc1)N*)C', 'COc1cccc(c1C(=O)CN1CCCC[C@H]1C(=O)*)OC', 'CHEMBL327743') ('CC[C@@H]([C@H](NC(=O)*)/C=C/c1ccccc1)C', 'COc1cccc(c1C(=O)CN1CCCCC1*)OC', 'CHEMBL327743') ('CC[C@@H](C(NC(=O)[C@@H]1CCCCN1CC(=O)c1c(OC)cccc1OC)*)C', '*/C=C/c1ccccc1', 'CHEMBL327743') ('COc1cccc(c1C(=O)CN1CCCC[C@H]1C(=O)N[C@H](/C=C/c1ccccc1)*)OC', 'CC[C@H](C)*', 'CHEMBL327743') ('*c1ccccc1', '*/C=C\\[C@H]([C@H](CC)C)NC(=O)[C@@H]1CCCCN1CC(=O)c1c(OC)cccc1OC', 'CHEMBL327743') ('*OC', 'CC[C@@H]([C@H](NC(=O)[C@@H]1CCCCN1CC(=O)c1c(*)cccc1OC)/C=C/c1ccccc1)C', 'CHEMBL327743') ('*OC', 'CC[C@@H]([C@H](NC(=O)[C@@H]1CCCCN1CC(=O)c1c(*)cccc1OC)/C=C/c1ccccc1)C', 'CHEMBL327743') ('*CC', 'COc1cccc(c1C(=O)CN1CCCC[C@H]1C(=O)N[C@@H](C(C)*)/C=C/c1ccccc1)OC', 'CHEMBL327743') ('CC[C@@H]([C@H](NC(=O)[C@@H]1CCCCN1CC(=O)c1c(OC)cccc1OC)/C=C/c1ccccc1)*', 'C*', 'CHEMBL327743') ('*C[C@@H]([C@H](NC(=O)[C@@H]1CCCCN1CC(=O)c1c(OC)cccc1OC)/C=C/c1ccccc1)C', 'C*', 'CHEMBL327743')
コード
どこで切断されたかどうかを明確にするために切断されたボンドにダミーアトム(atomicnumが0のアトム)をおいた。
canonicalSMILESで出力して.の文字でsplitした文字列をそのままキーに使っているが、openbabelのcanonicalizationアルゴリズムはMorgan Algorithmを改変したものを使っているのでおそらく大丈夫。
下のコードでは標準出力に書きだしたが、実際はデータベース化することになる。多分Mongo使う。
smi_string,id = clone.write(format="can")[:-1].split('\t') frag1, frag2 = smi_string.split('.') dbh.mmp.update({"core": frag1}, {"$push": {"frags": (frag2, id)}}, True) dbh.mmp.update({"core": frag2}, {"$push": {"frags": (frag1, id)}}, True)
こんな感じだろうか。
以下MMPを求めるのに使ったコード
import pybel import openbabel import sys input = sys.argv[1] def mkclone(mol): obc = openbabel.OBConversion() obc.SetInAndOutFormats("mol", "mol") molstring = obc.WriteString(mol) new_mol = openbabel.OBMol() obc.ReadString(new_mol, molstring) return new_mol def ring_member(rings, atom): for i, ring in enumerate(rings): if ring.IsMember(atom): return i return -1 def is_same_ring_member(rings, a1, a2): a1_idx = ring_member(rings, a1) a2_idx = ring_member(rings, a2) if a1_idx == -1 or a2_idx == -1: return False return a1_idx == a2_idx def detect_del_bonds(mol): rings = mol.OBMol.GetSSSR() del_bonds = [] for a in mol: if a.atomicnum == 6: for na in mol: n_a = na.OBAtom b = a.OBAtom.GetBond(n_a) if b != None and b.GetBondOrder() == 1 and not is_same_ring_member(rings, a.OBAtom, n_a): if a.OBAtom.GetIndex() < n_a.GetIndex(): del_bonds.append((a.OBAtom.GetIndex(), n_a.GetIndex())) return del_bonds for mol in pybel.readfile("sdf", input): chembl_id = mol.data["CHEMBL ID"] print chembl_id del_bonds = detect_del_bonds(mol) for a1, a2 in del_bonds: clone = pybel.Molecule(mkclone(mol.OBMol)) # !!! index starts 0 but GetAtom starts 1 !!! ob_a1 = clone.OBMol.GetAtom(a1 + 1) ob_a2 = clone.OBMol.GetAtom(a2 + 1) bond = ob_a1.GetBond(ob_a2) clone.OBMol.DeleteBond(bond) new_atom1 = clone.OBMol.NewAtom() new_atom1.SetAtomicNum(0) new_bond1 = clone.OBMol.NewBond() new_bond1.SetBegin(ob_a1) new_bond1.SetEnd(new_atom1) new_bond1.SetBondOrder(1) new_atom2 = clone.OBMol.NewAtom() new_atom2.SetAtomicNum(0) new_bond2 = clone.OBMol.NewBond() new_bond2.SetBegin(ob_a2) new_bond2.SetEnd(new_atom2) new_bond2.SetBondOrder(1) smi_string, id = clone.write(format="can")[:-1].split('\t') frag1, frag2 = smi_string.split('.') #print (frag1, chembl_id) #print (frag2, chembl_id) print (frag1, frag2, chembl_id)