fingerprintの由来を知りたい
openbabelのFP4フィンガープリントは部分構造由来のビット列なので、どのパターンに由来するビットが立っているのか知りたいときがよくあります。(MACCSなんかも同様に)
pybelを使えば簡単にサーチできるようです。
import pybel def readsmartsfile(filename="/usr/local/share/openbabel/2.2.0/SMARTS_InteLigand.txt"): patterns = [] inputfile = open(filename, "r") for line in inputfile: line = line.strip() if line and line[0]!="#": colon = line.find(":") name = line[:colon] smarts = line[colon+1:].strip() patterns.append([pybel.Smarts(smarts), name]) return patterns if __name__=="__main__": import sys patterns = readsmartsfile() file = sys.argv[1] for molecule in pybel.readfile("sdf", file): print "Looking at molecule %s" % molecule.title print "It has the following functional groups:", for smarts, name in patterns: if smarts.findall(molecule): print name, print "\n"
openbabelのfingerprintをビット列にする
openbabelではxhオプションをつけることでフィンガープリントを16進表記で書き出すことができます。
babel -ismi test.smi -ofpt test.fpt -xh -xfFP2 less test.fpt >test1 18 bits set 00002000 00000000 00000100 00010200 00000000 00000000 00000000 00000000 00000000 00000840 00000000 40008000 00000000 00000000 00000000 00000000 00020000 02000000 04000000 00000000 00000000 08000000 10000000 00400000 00000000 80000000 00000000 00000040 00004000 00020000 00000000 00000000 >test2 Tanimoto from test1 = 0.027027 00010000 01000000 00000000 00000000 00080000 00000000
ただしこれだと、解析するときに不便なのでビット列に直したい時があります。なのでpythonでビット列に変換するコードを書いてみました
def hex2bin(fingerprint): bf = "" h2b = {"0":"0000","1":"0001","2":"0010","3":"0011", "4":"0100","5":"0101","6":"0110","7":"0111", "8":"1000","9":"1001","a":"1010","b":"1011", "c":"1100","d":"1101","e":"1110","f":"1111", } for l in fingerprint: for c in l: b = h2b.get(c) if b: bf += b return bf def convert(file): result = "" for data in open(file,"r").read().split("\n>"): fp = "" for list in data.split("\n")[1:]: fp += hex2bin(list) result += data.split("\n")[0].split(" ")[0] + " " + fp + "\n" return result if __name__ == "__main__": import sys file = sys.argv[1] sys.stdout.write(convert(file))
フィンガープリントをさらにsplitしてcsvで書き出せばRであれとかこれとかQSAR的な解析が出来るようになりますね。
RECAPルールで構造を切断したい
molblasterのように構造をランダムに切断するのではなくて、もう少し合成を考慮した方法としてRECAPルールが使えます。ルールは全部で11あります。
RECAPを簡単に使える方法がなかったので、Chemistry::RECAPというPerlのモジュールを作ってみました。PerlMolとUNIVERSAL::requireが必要です。
use Chemistry::RECAP; use Perl6::Say; my $recap = Chemistry::RECAP->new(); my @frags = $recap->cleave("CC(=O)NCCn1cccc1-c2ccccc2"); say join ",", @frags;
実行結果
$ perl recap.pl CC=O,CCN,n1cccc1,c1ccccc1
RECAPは切断部位のみで切断後の構造に関しては特に指定してないのでこのモジュールは切断した位置には全て水素を付加しています。実際にフラグメントを反応させたい場合には単に水素を埋めるだけではなく、酸クロライドみたいな試薬としても適切な形でキャッピングするようにしたほうがよいでしょう。
構造をランダムに切断したい
現在の創薬シーンにおいてはフラグメントベースのドラッグデザイン(FBDD)やそれをcomputational的に行うのは主流に近いと思います。
ランダムに結合を切る方法としてhttp://www.ncbi.nlm.nih.gov/pubmed/16995724というのがあるのでopenbabelを使って簡単に実装してみます。
molblasterのアルゴリズムは単純で、ボンドをサンプリングしてカットしていくというのを何回か繰り返してフラグメントの頻度を数えるというものです。というわけで、一回にどのくらいの数のボンドを切断するかと、何回繰り返すかで精度が変化するようです。
import openbabel as ob from random import sample def randomsplit(mol,cutnum=5): cutlist = sample(range(mol.NumBonds()),cutnum) delbonds = [] for i,bond in enumerate(ob.OBMolBondIter(mol)): if i in cutlist: delbonds.append(bond) for b in delbonds: mol.DeleteBond(b) return mol def molblaster(smi,iter=100): obc = ob.OBConversion() obc.SetInAndOutFormats('smi','smi') for i in range(iter): mol = ob.OBMol() obc.ReadString(mol,smi) yield obc.WriteString(randomsplit(mol))[:-1] if __name__ == "__main__": smiles = 'C(C)CCc1ccccc1Cl' for smi_string in molblaster(smiles): print smi_string
自分専用お気に入りフラグメントデータベースを持つことがモデラーとしてのキャリアの第一歩だという噂です。
脱塩したい(pybelで)
化合物の物性などを予測する時には、イオン化の状態を揃えて中性にしたり脱塩して、単一の化合物にするというルーチンワークが必須です。
openbabelにはそのためのメソッドとしてStripSaltというものが用意されています。例として以下の化合物の脱塩をしてみます。
http://pubchem.ncbi.nlm.nih.gov/summary/summary.cgi?cid=24847969&loc=ec_rcs:aspirin
import pybel mol = pybel.readstring('smi','CCC(C)C(=O)OC1CC(C=C2C1C(C(C=C2)C)CCC(CC(CC(=O) \ [O-])O)O)O.CC(=O)OC1=CC=CC=C1C(=O)O.[Na]') mol.OBMol.StripSalts(14) mol.write('smi') # 'CCC(C)C(=O)OC1CC(C=C2C1C(C(C=C2)C)CCC(CC(CC(=O)[O-])O)O)O\t\n'
引数として与えた数未満のヘテロアトム数のcompoundが除かれるようです。
実際には、化合物本体よりも大きい塩があったり、R体とS体の両方が一つの化合物として入力されていたりと解析者泣かせのデータを扱うことが多いので、もう少し複雑怪奇な脱塩ルーチンが使われることが多いようですが、、、（謎）
pybelで特性を計算する
pybelはcalcdescメソッドでlogP,MR,TPSAを計算することができます
import pybel mol = pybel.readstring('smi','COC=C') mol.calcdesc() #{'LogP': 0.77629999999999999, # 'MR': 17.146000000000001, # 'TPSA': 9.2300000000000004} mol.OBMol.GetExactMass() #58.041864812
さらに、OBMol属性からopenbabelのメソッドを呼ぶことができます。今回はExaxtMassを計算してみました。
perlでもOASAライブラリを利用して構造描画したい
Inline::Pythonというモジュールを使うとPythonのコードをperlに埋め込むことができます。
use Inline Python => <<'END'; import pybel def draw_png(smiles,file): mymol = pybel.readstring('smi',smiles) mymol.draw(filename=file, show=False) END draw_png("CCCc1ccccc1OC","sample.png");
Inline系のモジュールは色々あります。参考までにInline::JavaでCDKなどを使った例もあげておきます。
Structure-CDKの絵は奇麗なので気に入っていますが、もっとこうアートな感じの描画ライブラリがあってもいいのになとは思います。
