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

化学

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

化学は分子を中心に扱う学問ですので,Molオブジェクトの扱いに習熟することは大変重要です.

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

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

from rdkit import rdBase, Chem
from rdkit.Chem import AllChem, Draw
print('rdkit version: {}'.format(rdBase.rdkitVersion))
# rdkit version: 2019.09.3
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オブジェクトとBondオブジェクトは互いにリンクしており,

  • Atomオブジェクトはどの結合を形成しているかと
  • Bondオブジェクトはどの原子から結合が形成されているか

といった情報が相互に記録されています.

Atomオブジェクト

Mol.GetAtoms()
Mol.GetAtomWithIdx(idx)

Atomオブジェクトには

  • 全原子をまとめて取得
  • 原子のインデックスを使って選択的にアクセス

することが可能です.

このようにして得られたAtomオブジェクトには,

原子の情報 メソッド
インデックス atom.GetIdx()
原子番号 atom.GetAtomicNum()
元素記号 atom.GetSymbol()
質量 atom.GetMass()
混成 atom.GetHybridization()
環内原子か否か atom.IsInRing()
芳香族か否か atom.GetIsAromatic()
結合リスト atom.GetBonds()

などのように,様々な要素を取得するメソッドが用意されています.

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

for mol in mols[:3]:
     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: 11
###
Dimethyl 4-Acetoxyisophthalate
index for sp3 carbon: 12
index for sp3 carbon: 13
index for sp3 carbon: 16
###

Bondオブジェクト

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

Bondオブジェクトについても,

  • 全結合をまとめて取得
  • 結合のインデックスを使ってアクセス

によって取得が可能です.

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

Bondオブジェクトについても下表のように

結合の情報 メソッド
結合のインデックス bond.GetIdx()
結合の種類 bond.GetBondType()
芳香族結合か否か bond.GetIsAromatic()
環内結合か否か bond.IsInRing()
始点原子を取得 bond.GetBeginAtom()
始点原子をIdで取得 bond.GetBeginAtomIdx()
終点原子を取得 bond.GetEndAtom()
終点原子をIdで取得 bond.GetEndAtomIdx()

様々な要素に対して対応するメソッドが設定されており,これらの情報が容易に取得可能になっています.

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

  1. bond.GetIsAromaticで芳香族結合か否かを判断
  2. 結合始点と終点の原子をGetBeginAtom, GetEndAtomで取得
  3. 元素記号とインデックスを取得

という手順で処理をしています.

for mol in mols[:2]:
    print(mol.GetProp('PRODUCT_NAME'))
    for b in mol.GetBonds():
        if b.GetIsAromatic():
            atom_a = b.GetBeginAtom()
            atom_b = b.GetEndAtom()
            print('bond between {}{}-{}{} is aromatic.'.
            format(atom_a.GetSymbol(), atom_a.GetIdx(),
                   atom_b.GetSymbol(), atom_b.GetIdx()))
    print('---')
Methyl Acetylsalicylate
bond between C0-C1 is aromatic.
bond between C1-C2 is aromatic.
bond between C2-C3 is aromatic.
bond between C3-C4 is aromatic.
bond between C4-C5 is aromatic.
bond between C0-C5 is aromatic.
---
Methyl o-Anisate
bond between C0-C1 is aromatic.
bond between C1-C2 is aromatic.
bond between C2-C3 is aromatic.
bond between C3-C4 is aromatic.
bond between C4-C5 is aromatic.
bond between C0-C5 is aromatic.
---

RingInfoオブジェクト

Mol.GetRingInfo()

環状構造は分子において特別な部分構造ですので,RDkitでは環構造情報に特化したRingInfoオブジェクトが存在しています.

ここでは以下に示す多環性スピロ化合物を例にRingInfoオブジェクトについて見ていきたいと思います.

spiro = Chem.MolFromSmiles('C1CCC12CCCC23CC3')
spiro_ring = spiro.GetRingInfo()

この分子のAtomとBondのインデックスは以下のようになっています.黒数字がAtom,青数字がBondオブジェクトのインデックスに対応しています.

環を構成するAtomやBondの情報

RingInfo.NumRings()
RingInfo.AtomRings()
RingInfo.BondRings()

分子内に存在する環構造について,それぞれ構成するAtomやBondの集合をタプルのタプルで返します.

print(spiro_ring.AtomRings())
print(spiro_ring.BondRings())

上で示した番号と一致しています.

((0, 3, 2, 1), (4, 3, 7, 6, 5), (8, 7, 9))
((9, 2, 1, 0), (3, 10, 6, 5, 4), (7, 11, 8))

環サイズや環数に関する情報

NumAtomRings(atomIdx)
IsAtomInRingOfSize(atomIdx, ringSize)
MinAtomRingSize(atomIdx)

それぞれ

  • 原子が何個の環に含まれているか
  • 原子が指定サイズを環を形成しているか
  • 原子を含む最小の環サイズ

について返すメソッドになります.

同様にBondオブジェクトに関しても対応するメソッドが存在しています.

print(spiro_ring.NumAtomRings(3)) # 2
print(spiro_ring.IsAtomInRingOfSize(0,4)) # True
print(spiro_ring.MinAtomRingSize(3)) # 4

Molオブジェクトのプロパティ

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

Molオブジェクトにはプロパティと呼ばれる任意の値を付加可能です.

  • プロパティの設定にはSetProp
  • プロパティの参照にはGetProp

を使います.

もし元のSDFに最初からプロパティが設定してある場合にはオブジェクト作成と同時にそれらの値も設定されています.先ほどのコードでもPRODUCT_NAMEというプロパティを参照している箇所がありました.これはMolオブジェクト作成に利用したSDFに最初から設定されていた値になります.

既にMolオブジェクトにどのようなプロパティが存在するかを知りたい場合にはGetPropNamesで取得可能です.

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

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

for prop in mols[0].GetPropNames():
     print('{}: {}'.format(prop, mols[0].GetProp(prop)))
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」を参照してください.
MOLファイルやSDFについては,「MOLファイル・SDFとはどんな化学情報ファイルなのか?」という記事で解説しています.参照してみてください.
MOLファイル・SDFとはどんな化学情報ファイルなのか?
我々の目はモノの形やパターンを読み取る能力に優れています.化学者が用いる構造式は人間の目で理解しやすいように設計された分子の表現です.例えば下のアスピリンの構造式を見るだけで,化学者は多くの情報を瞬時に理解できます.ところがコンピュータにとってはこのような構造式の画像フ...

分子の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だと思います.このメソッドはオプション指定で,

オプション名 内容
molsPerRow 一列あたりに描画する分子の数
subImgSize 分子1つの大きの設定
legends 分子毎に凡例をつける

など,細かく設定可能になっています.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]])


この例で明らかですが,MolToImageMolsToGridImageでは画像の質が大きく違います.

構造式の描画については「RDKitでの構造式描画を詳しく解説」という記事で詳細に解説しています.参照してみてください.
RDKitでの構造式描画を詳しく解説
構造式を2次元に描画することは人間が分子の形・性質を理解する第一歩です.これまで「RDKitの分子Molオブジェクトを扱う」という記事ではRDKitにおける分子の扱い方や描画方法を学びました.また「RDKitを用いた部分構造検索とMCSアルゴリズム」という記事では複数分子の間...

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

先ほどのグリッドイメージでは化合物の並び方が一致していないため,類似構造がぱっと見ではわかりにくいです.そこでサリチル酸をテンプレートとして分子を回転させてみます.

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オブジェクトの作成方法をいろいろと扱いました.

今回は「RDKitの分子Molオブジェクトを扱う」という話題について,

  • 分子を構成する原子結合に関する情報を取得する方法
  • 分子内の環構造に関する情報を有するRingInfoオブジェクトの扱い方
  • 分子の2次元表記について構造の描画方法や部分構造の扱い方

などを学びました.

次回からは分子を2次元から3次元へと展開していくことで,コンフォマーの生成や分子力学法を用いた構造最適化エネルギー計算といった内容を扱います.

ここまでを学びきることができれば記述子を使ったケモインフォマティクスにおける予測モデルの作成から,Gaussianを始めとした外部の量子化学計算プログラムとの連携までいろいろな可能性が広がります.是非頑張っていきましょう.

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

コメント

  1. S より:

    「Atomオブジェクト」の節 (#Atom) で (i) の黄色枠に Mol.GetAtomWithIndex(idx) としておりますが Mol.GetAtomWithIdx(idx) の誤りではないでしょうか、確認お願いいたします。

タイトルとURLをコピーしました