化合物から水素原子を引き抜いたラジカルの組み合わせ
ある化合物がCYPにより水酸化されて代謝されると分かっている場合に、全ての水素原子の乖離エネルギーを求めてラジカル生成エネルギーの低い順に代謝部位のあたりをつけたりしますが、毎回手で水素を消去した組み合わせを作るのは面倒なのでopenbabelでやってみます。
molオブジェクトのコピーをとるメソッドが見つからなかったのでcloneという関数を定義してます。
import openbabel as ob import sys def clone(mol): new_mol = ob.OBMol() for atom in ob.OBMolAtomIter(mol): new_mol.AddAtom(atom) for bond in ob.OBMolBondIter(mol): new_mol.AddBond(bond) return new_mol def abs_h(mol): radicals = [] mol.AddHydrogens() h_indexs = [atom.GetIdx() for atom in ob.OBMolAtomIter(mol) if atom.GetType() == 'H'] for h_index in h_indexs: new_mol = clone(mol) new_mol.DeleteAtom(new_mol.GetAtom(h_index)) new_mol.SetTotalSpinMultiplicity(2) new_mol.SetTotalCharge(0) radicals.append(new_mol) return radicals if __name__ == '__main__': input = sys.argv[1] obc = ob.OBConversion() obc.SetInAndOutFormats("smi", "mol") mol = ob.OBMol() next = obc.ReadFile(mol,input) while next: radicals = abs_h(mol) for r in radicals: print obc.WriteString(r) mol = ob.OBMol() next = obc.Read(mol)
あとはGamessでエネルギー計算すればいいですね。
GAMESSでさくっと構造最適化をしたい
forcefieldを用いた構造最適化の限界は、特に芳香環での部分電荷なんかに顕著に現れます。
そのため、必要最低限の量子化学計算のスキルはケミストには必須だと思います。しかし、量子化学計算ソフトウェアは歴史的経緯からか、ライトユーザーにとってはどうでもいい中間ファイルを出力したり、インプット自体が複雑怪奇だったりしますね？
そこで、そこら辺をまとめて面倒見つつ、必要最低限の作業をサクっとするようなものを書いてみました。
これを使うと、段々基底関数の精度をあげつつ最適化をするというルーチンがこんなかんじで書けます
import openbabel as ob from gamess import Gamess g = Gamess() obc = ob.OBConversion() obc.SetInFormat("mol") methane = ob.OBMol() next = obc.ReadFile(methane,"methane.mol") # STO-3Gで構造最適化 methane_sto3g = g.optimize(methane, basis="sto3g") # STO-3Gで最適化された構造を初期値にして、さらに6-31Gで最適化 methane_631g = g.optimize(methane_sto3g, basis="631g") for a in ob.OBMolAtomIter(methane_631g): print "%s x: %1.3f, y: %1.3f, z: %1.3f charge: %1.3f" % \ (a.GetType(), a.x(), a.y(), a.z(), a.GetPartialCharge())
実行結果
C3 x: 0.937, y: 0.067, z: 0.028 charge: -0.660 HC x: 2.019, y: 0.101, z: 0.086 charge: 0.165 HC x: 0.522, y: 0.935, z: 0.527 charge: 0.165 HC x: 0.578, y: -0.834, z: 0.512 charge: 0.165 HC x: 0.630, y: 0.066, z: -1.011 charge: 0.165
実際にユーザーが使いやすいようにするのはCGIのようなWEBから実行出来るようにしたほうがいいかもしれませんね。
アロマティックなヘテロ環を漏れなくさがす
MORPHというアロマティックなヘテロ環を探索するツールの論文があります。
結局、総当りで探していってバリデータで確認すればいいっぽいので、単環で探すようなものを書いてみました。
戦略としては単純で、まず総当りで評価するSMILES Stringを生成します。原子のタイプ(atomchoices)とボンドのタイプ(bondchoices)を別々に生成してからコンビネーションをとれば、全ての組み合わせがつくれます(もれなくだぶりありなのでミッシーではないですが)。
アロマティックなリングかどうかを判定するのに、アロマティックなリングは小文字になるという性質を使います。要するにopenbabelでSMILES文字列を読んで、SMILES文字列を出力させてみて小文字だったら、芳香環が生成できたと判断します。find_aromatic_ringsという部分です。
from pybel import * import re def atomchoices(n): if n == 1: return ["C","N","O"] else: return [atom+atoms for atoms in atomchoices(n-1) for atom in ["C","N","O"]] def bondchoices(n): if n == 1: return ["-","="] else: return [bond+bonds for bonds in bondchoices(n-1) for bond in ["-","="]] def make_smiles(n): atoms = atomchoices(n) bonds = bondchoices(n) atoms = [a for a in atoms if "C" in a] bonds = [b for b in bonds if b.count("=") / float(n) <= 0.5] for a in atoms: for b in bonds: smiles = "" for i in range(n): smiles = smiles + a[i] if i == 0: smiles += "1" smiles = smiles + b[i] if i == n-1: smiles += "1" if "=O" not in smiles and "O=" not in smiles and "O-O" not in smiles: mymol = str(readstring("smi",smiles))[:-2] yield mymol def find_aromatic_rings(n): rearoma = re.compile('^[a-z1-9]+$') aromatic_compounds = list(set([smi for smi in make_smiles(n) if rearoma.search(smi)])) print aromatic_compounds
縮合環でやるならそれ用のテンプレート文字列を渡せるようにするか内部で辞書的に持っておいてそれを置き換えるようにすればいいんじゃないかと思います。
ライブラリ中の特定の構造を削除する
サブストラクチャーサーチで集めたライブラリや、テンプレート構造を使って作ったバーチャルライブラリ群のプロパティを分析ししたいときに、テンプレート構造やサーチクエリに使った構造は要らない（むしろ邪魔）ことが多いと思います。
そういう場合プロパティの全平均からの差分と評価したりするのですが、構造自体を消してしまいたいこともあります。
import openbabel as ob import sys input = sys.argv[1] obc = ob.OBConversion() obc.SetInAndOutFormats("smi", "sdf") mol = ob.OBMol() next = obc.ReadFile(mol,input) sp = ob.OBSmartsPattern() sp.Init("c1ccccc1") while next: sp.Match(mol); maplist = sp.GetUMapList() if len(maplist) > 0: del_atoms = [] for i,atom in enumerate(ob.OBMolAtomIter(mol),1): if i in maplist[0]: del_atoms.append(atom) for a in del_atoms: mol.DeleteAtom(a) print obc.WriteString(mol) mol = ob.OBMol() next = obc.Read(mol)
これで例えばCOCc1ccccc1からベンゼン環を削除してCOCという構造だけを残すことができます。
GASTONで頻度よく現れる部分構造を探る
環構造を抽出する
PerlMolにはRing::Findというモジュールがあって、環のサーチをしてくれるのですが、返ってくるものが、アトムとボンドのIDの詰め合わせなのでちょっとあれです。本来ならば、見つかったリングの部分構造の分子オブジェクトのリストが返ってくると嬉しい場合も結構あります。
use warnings; use strict; use Chemistry::File::SMILES; use Chemistry::Ring 'aromatize_mol'; use Chemistry::Ring::Find ':all'; use Perl6::Say; sub get_rings { my $origin = shift; my @mols = (); my @rings = find_rings($origin, all => 1); for my $fg (@rings){ my $mol = Chemistry::Mol->new(id => $fg->{id}, name => $fg->{id}); for my $a (@{$fg->{atoms}}) { $mol->new_atom(name => $a->name, id => $a->id, symbol => $a->symbol) } for my $b (@{$fg->{bonds}}) { $mol->new_bond(name => $b->name, type => $b->type, atoms => [$mol->by_id($b->{atoms}->[0]->{id}), $mol->by_id($b->{atoms}->[1]->{id})], order => $b->order) } $mol->add_implicit_hydrogens; aromatize_mol($mol); push @mols, $mol; } return @mols; } my $smi = "c1c(CC=O(O))cccc1CC2CC2CCC3CCCCCCC3O"; my $origin = Chemistry::Mol->parse($smi, format => 'smiles'); aromatize_mol($origin); for my $mol (get_rings($origin)) {say $mol->print(format => 'smiles', unique => 1, name => 1);}
