RDKitの分子Molオブジェクトを扱う

RDKitでケモインフォマティクスに入門」という記事では,RDKitの分子であるMolオブジェクトの作成方法を中心に学習しました.今回はそのMolオブジェクトに対してどのような操作が行えるかを学んでいきます.

必要なライブラリのインポート

まずは必要なライブラリをimportし,分子を読み込んでおきましょう.SDFとしては前回と同じく,東京化成から入手したサリチル酸誘導体を使います.

from rdkit import rdBase, Chem
from rdkit.Chem import AllChem, Draw
print('rdkit version: {}'.format(rdBase.rdkitVersion))
# rdkit version: 2017.09.1
suppl = Chem.SDMolSupplier('./sdf_20180716131641.sdf')
mols = [x for x in suppl if x is not None]
len(mols)
# 73

Molオブジェクトの構成要素

Molオブジェクトは分子ですので,当然原子から構成されています.また分子は原子と原子の間に結合が形成されることによりできあがります.こういった化学者の感覚通りに,RDKitのMolオブジェクトにもAtomオブジェクトやBondオブジェクトが用意されています

結合は原子間に存在することから各原子がどの結合を形成しているか,また逆に各結合を形成する原子は何かといった情報も記録されています.

Atomオブジェクト

Mol.GetAtoms()
Mol.GetAtomWithIndex(idx)

Atomオブジェクトには全原子をまとめて取得するか,インデックスを使ってアクセスが可能です.Atomオブジェクトには,インデックス,原子番号,原子記号,質量,混成,価数,をはじめとして様々な要素が設定されており,いずれの情報も簡単に取得可能です.

ここでは73分子のうち最初の5個について,sp3混成の炭素原子のインデックスを表示してみましょう.注意点としてはGetHybridization()で得られる結果はChem.rdchem.HybridizationType.SP3というオブジェクトなので,このままでは’SP3’という文字列とは比較できません.そこでstringオブジェクトに変換する作業が必要になります.

for mol in mols[:5]:
     print(mol.GetProp('PRODUCT_NAME'))
     for a in mol.GetAtoms():
         if a.GetSymbol() == 'C' and str(a.GetHybridization()) == 'SP3':
             print('index for sp3 carbon: {}'.format(a.GetIdx()))
     print('###')
Methyl Acetylsalicylate
index for sp3 carbon: 9
index for sp3 carbon: 12
###
Methyl o-Anisate
index for sp3 carbon: 9
・・省略・・
index for sp3 carbon: 16
index for sp3 carbon: 25
###

Bondオブジェクト

Mol.GetBonds()
Mol.GetBondWithIdx(idx)
Mol.GetBondBetweenAtoms()

結合についても全結合をまとめて取得するか,インデックスを使ってアクセスが可能です.また結合を形成する原子のインデックスがわかっている場合にはその結合を直接指定することも可能です.

Bondオブジェクトについても開始・終了原子のインデックス,結合の種類,結合次数,環内の結合かどうか,芳香族かどうか,共役しているかどうか,など様々な要素が設定されており容易に取得可能です.

ここでは73分子の最後の3分子について,芳香族結合を構成する原子のインデックスを出力してみます.

for mol in mols[-3:]:
     print(mol.GetProp('PRODUCT_NAME'))
     for b in mol.GetBonds():
         if b.GetIsAromatic():
             print('bond between {}-{} is aromatic.'.
             format(b.GetBeginAtomIdx(),b.GetEndAtomIdx()))
     print('###')
2-Ethylhexyl Salicylate
bond between 3-4 is aromatic.
bond between 3-8 is aromatic.
bond between 4-5 is aromatic.
bond between 5-6 is aromatic.
bond between 6-7 is aromatic.
bond between 7-8 is aromatic.
###
Propyl Salicylate
bond between 3-4 is aromatic.
bond between 3-8 is aromatic.
・・省略・・

プロパティ

GetProp(name)
SetProp(name, value)
GetPropNames()

Molオブジェクトにはプロパティと呼ばれる任意の値を付加可能です.設定にはSetProp,参照にはGetPropを使います.もし元のSDFに最初からプロパティが設定してある場合にはオブジェクト作成と同時にそれらの値も設定されています.

先ほどのコードでもPRODUCT_NAMEというプロパティを参照している箇所がありますね.既にMolオブジェクトにどのようなプロパティが存在するかを知りたい場合にはGetPropNamesで取得可能です.

他に特別なプロパティとしてSDFのタイトル行に相当する_Nameというものもあります.こちらはマジックプロパティと呼ばれる隠れプロパティ扱いになっており,GetPropNamesでは取得できません.

下のコードでは設定されているプロパティとその値の一覧を表示してます.

for p in mols[0].GetPropNames():
     print('{}: {}'.format(p, mols[0].GetProp(p)))
PRODUCT_NUMBER: A0114
PRODUCT_NAME: Methyl Acetylsalicylate
MOLECULAR_FORMULA: C10H10O4
MOLECULAR_WEIGHT: 194.19
CAS_NUMBER: 580-02-9
MDL_NUMBER: MFCD00014978
マジックプロパティの一覧については「The RDKit Book:“Magic” Property Values」を参照してください.

分子の2次元構造描画

今回取り扱うの分子は全て原子の座標を持っていますが,もしも初期設定されていないようでしたら,以前見たようにAllChem.Compute2DCoordinateで分子を描画する前に計算しておく必要があります.

1) 1分子の描画
a. Chem.Draw.MolToImage(mol)
b. Chem.Draw.MolToFile(mol, file_name)
2) 複数分子の描画
a. Chem.Draw.MolsToImage(mols)
b. Chem.Draw.MolsToGridImage(mols, molsPerRow, subImgSize, legends)

ファイルに保存したい場合にはMolToFileを,Jupyter Notebookにそのまま打ち出したい場合にはMolToImageを使います.複数分子を描画したい場合にはMolsと複数形にします.凡例を分子毎につけたい場合などはMolsToGridImageを使います.他にも一列あたりに並べるグリッドの数や,グリッド1つの大きさなど,細かく設定可能です.MolsToImageの使いどころは正直よくわかってません.

# 1分子の描画
Draw.MolToImage(mols[0])
# 複数分子の描画(凡例付き)
Draw.MolsToGridImage(mols[:9], molsPerRow=3, subImgSize=(300,200),
                     legends=[x.GetProp('PRODUCT_NUMBER') for x in mols[:9]])


構造式の描画については「RDKitでの構造式描画を詳しく解説」という記事で詳細に解説しています.参照してみてください.

テンプレート構造を用いた構造の並列化

先ほどのグリッドイメージでは化合物の並び方が一致していないため,類似構造がぱっと見ではわかりにくいです.そこでサリチル酸をテンプレートとして分子を回転させてみます.Chem.AllChemにあるGenerateDepictionMatching2DStructureを使います.テンプレートのMolオブジェクトを作成して2次元座標の計算を行うだけです.

# テンプレート作成
import pubchempy as pcp
tmp = pcp.get_compounds('salicylic acid', 'name')
tmp = tmp[0]
tmp_smiles = tmp.canonical_smiles
template = Chem.MolFromSmiles(tmp_smiles)
AllChem.Compute2DCoords(template)
# テンプレートをもとに2次元座標の再計算
for mol in mols:
    if mol.HasSubstructMatch(template):
        AllChem.GenerateDepictionMatching2DStructure(mol, template)
# グリッドイメージで再び描画
Draw.MolsToGridImage(mols[:9], molsPerRow=3, subImgSize=(300,200),
                    legends=[x.GetProp('PRODUCT_NUMBER') for x in mols[:9]])

今回のケースではもともとの描画の仕方がきれいなので,5個目のメチレンジオキシ部位や9個目の糖の部分がオリジナルの方が美しい気もしますが,作業のポイントはわかっていただけると思います.

部分構造の取り扱い

Mol.HasSubstructMatch(substruct)

先ほどのテンプレートを用いた描画でさりげなく使いましたが,Molオブジェクトには部分構造を引数に取るメソッドが用意されています.例えば部分構造を有するか否かを調べることで,分子のフィルタリングに使えます.

ここでは例としてアスピリン(アセチルサリチル酸)を部分構造に持つ分子の数を調べ,描画してみます.73個の分子うち,6個の分子のフェノールがアシル化されているとわかりました.

# テンプレート作成
aspirin = pcp.get_compounds('aspirin', 'name')
aspirin = aspirin[0]
aspirin_sm = aspirin.canonical_smiles
aspirin_mol = Chem.MolFromSmiles(aspirin_sm)
AllChem.Compute2DCoords(aspirin_mol)
# 部分構造を有する分子の描画
match = []
for mol in mols:
    if mol.HasSubstructMatch(aspirin_mol):
        match.append(mol)
print(len(match)) # 6
Draw.MolsToGridImage(match, molsPerRow=3, subImgSize=(300,200),
                     legends=[x.GetProp('PRODUCT_NAME') for x in match])

部分構造検索や複数分子間での共通構造の探索については「RDKitを用いた部分構造検索とMCSアルゴリズム」という記事で詳しく解説しています.参照してみてください.

終わりに

RDKitでケモインフォマティクスに入門」ではRDKitにおける分子の表現方法であるMolオブジェクトの作成方法をいろいろと扱いました.今回は分子を構成する原子結合に関する情報の取得からはじまり,分子の2次元表記について構造の描画方法や部分構造の扱い方などを学びました.

次回からは分子を2次元から3次元へと展開していくことで,コンフォマーの生成や分子力学法を用いた構造最適化エネルギー計算といった内容を扱います..ここまでを学びきることができれば記述子を使ったケモインフォマティクス分析から,Gaussianを始めとした外部プログラムとの連携までいろいろな可能性が広がります.是非頑張っていきましょう.

>>>次の記事:「RDKitによる3次元構造の生成

シェアする

  • このエントリーをはてなブックマークに追加

フォローする