🌊

RDKitの3次元座標を使いこなす（移動＆回転）

2023/04/23に公開

Python

RDKit

tech

RDKitの３次元座標操作

分子系のケモインフォマティクスの代表ツールともいえるRDKitですが、その機能は実に多彩です。日本語で紹介されている記事では単分子での扱いがほとんどですが、実は複数分子の取り扱いやScipyなどと組み合わせることで座標操作も可能です。この記事では複数分子の実際の取り扱いや、分子の移動・回転操作の一例を書いてみました。

これらの分子の操作を使いこなすと、二量体のGaussianのインプットファイルを自動生成したり、
錯体のような複雑なモデルも自動で生成したりすることができます。

複数分子を埋め込む

系の中に2つの分子を埋め込むこともできます。やり方はシンプルでChem.CombineMols()を用いて、2つの分子を指定するだけです。

smiles = "CC(=O)OC1=CC=CC=C1C(=O)O"
mol = Chem.MolFromSmiles(smiles)
mol = Chem.AddHs(mol)
mols = Chem.CombineMols(mol, mol)
AllChem.EmbedMolecule(mols)
v = py3Dmol.view(width=400, height=400)
mb = Chem.MolToMolBlock(mols)
v.addModel(mb, 'sdf')
v.setBackgroundColor('0xeeeeee')
v.setStyle({'stick':{}})
v.zoomTo()
v.show()

思いっきり分子が重なっていますが、確かに2つの分子がいますね。

分子を平行移動させる

重なっているのも気持ち悪いので平行移動してみましょう。RDKitの3次元座標はConformer→Positionという階層で記録されているので、mols.GetConformers()→conformer.GetPositions()の順番に座標を参照可能です。

下記では系内の一つ目の分子の座標を10Å平行移動しています。

from rdkit.Geometry import Point3D
mols = Chem.CombineMols(mol, mol)
AllChem.EmbedMolecule(mols)
conf = mols.GetConformer()

for c in mols.GetConformers():
    for i, p in enumerate(c.GetPositions()):
        if i < len(mol.GetAtoms()):
            conf.SetAtomPosition(i, Point3D(p[0]+10, p[1], p[2]))

v = py3Dmol.view(width=400, height=400)
mb = Chem.MolToMolBlock(mols)
v.addModel(mb, 'sdf')
v.setBackgroundColor('0xeeeeee')
v.setStyle({'stick':{}})
v.zoomTo()
v.show()

こんな感じでできました。

分子を回転させる

分子の回転も平行移動と似たようなやり方でできます。3次元座標の回転はいろいろなやり方がありますが、例えば下記の記事に記載されているような方法を用いることで分子を回転させることができます。ここでは上で作成したアスピリンの二量体に対してカルボキシル基を対向させてみようと思います

RDKitで分子の回転を取り扱う場合、一旦各分子の座標をnumpyやpandasで取り出してScipyで回転させ、その後RDKitの３次元座標を書き換えます。おそらくもうちょっとクールなやり方があると思いますが、著者のRDKit力ではこれが限界です・・。

RDKitのMolオブジェクトからの３次元座標の取り出しは下記のスクリプトで実施できます。

from rdkit.Geometry import Point3D
import pandas as pd
import numpy as np

mols = Chem.CombineMols(mol, mol)
AllChem.EmbedMolecule(mols)
conf = mols.GetConformer()

for c in mols.GetConformers():
    position = pd.DataFrame(c.GetPositions(), columns=['x', 'y', 'z'])

print(position.head())

          x         y         z 
0  3.456982 -0.672642 -0.487726 
1  2.027099 -0.640773 -0.108061 
2  1.757157 -0.847229  1.108713 
3  1.005543 -0.399226 -1.015894 
4 -0.303193 -0.401190 -0.500701

各分子の回転は下記のように実施できます。下記のスクリプトでは2つ目の分子の0番目の原子と分子の重心のベクトルが(-1, 0, 0)方向になるように分子を回転しています。

def rotate_molecule(position_df, target_atom_id, vector):

    def rotation_matrix(frm, to):
        if np.array_equal(frm, to):
            return np.identity(3)
        if np.array_equal(frm, -to):
            return -np.identity(3)
        
        s = frm + to
        return 2.0 * np.outer(s, s) / np.dot(s, s) - np.identity(3)
    
    position_df[['x', 'y', 'z']] -= position_df[['x', 'y', 'z']].mean()

    atom_vector = position_df.iloc[target_atom_id].values
    print(atom_vector)
    atom_vector = atom_vector/np.linalg.norm(atom_vector)
    rot_matrix = rotation_matrix(atom_vector[0], vector)

    for i, value in enumerate(position_df[['x', 'y', 'z']].values):
        position_df.at[i, 'x'] = np.dot(rot_matrix, value)[0]
        position_df.at[i, 'y'] = np.dot(rot_matrix, value)[1]
        position_df.at[i, 'z'] = np.dot(rot_matrix, value)[2]

    return position_df

mol2_pos = position.iloc[len(mol.GetAtoms()):].copy().reset_index(drop=True)
mol2_pos = rotate_molecule(mol2_pos, 0, np.array([1, 0, 0]))

rot_position = pd.concat([mol1_pos, mol2_pos]).reset_index(drop=True)

回転後の原子座標を書き換えて、可視化してみると・・・

for c in mols.GetConformers():
    for i, p in enumerate(c.GetPositions()):
        rot_p = rot_position.iloc[i].values
        conf.SetAtomPosition(i, Point3D(rot_p[0], rot_p[1], rot_p[2]))

v = py3Dmol.view(width=400, height=400)
mb = Chem.MolToMolBlock(mols)
v.addModel(mb, 'sdf')
v.setBackgroundColor('0xeeeeee')
v.setStyle({'stick':{}})
v.zoomTo()
v.show()

これが・・・

こうなります。