はじめに
RDKitの部分構造について解説していきます。具体的には以下の2点に焦点を当てます。
- ある共通の部分骨格を持つ分子を抜き出す方法
- 最大共通部分骨格を求める方法
内容に関しては、RDKitの公式ページを参考にさせていただきました。
※マテリアルズインフォマティクス関係の内容を他にも投稿していますので、よろしければこちらの一覧から他の投稿も見て頂けますと幸いです。
環境
- windows10
- conda 4.9.2
- python 3.7.1
- rdkit 2020.03.2.0
データの読み込み
RDKitのGitHubページから取得した"cdk2.sdf"のデータを読み込んでいきます。
from rdkit import Chem
# データの読み込み
# suppleのデータ型はrdkit.Chem.rdmolfiles.SDMolSupplier
suppl = Chem.SDMolSupplier('data/cdk2.sdf')
# if文で条件を付けて、最初に空のリスト生成せずに一気に処理
# rdkit.Chem.rdchem.Molのデータ型のリストの作成
ms = [x for x in suppl if x is not None]
読み込んだデータはmsとしてrdkit.Chem.rdchem.Molのデータ型のリストで取得します。読み込んだデータを少し分析してみます。
# データセットに含まれる分子の数と分子量を調べてみる
from rdkit import Chem
from rdkit.Chem import rdMolDescriptors
weights = []
for i in ms:
weight = rdMolDescriptors._CalcMolWt(i)
weights.append(weight)
print("number_of_dataset:" + str(len(ms)))
print("Max_Molecularweight: " + str('{:.1f}'.format(max(weights))))
print("Min_Molecularweight: " + str('{:.1f}'.format(min(weights))))
number_of_dataset:47
Max_Molecularweight: 449.5
Min_Molecularweight: 235.2
読み込んだデータには47個の分子が含まれており、最大分子量は449.5であり最小分子量は235.2であることが分かります。分子量の算出はこちらを参考にさせて頂きました。分子量のヒストグラムを書いてみます。
import matplotlib.pyplot as plt
plt.hist(weights, bins=10)
plt.xlim(225, 475)
plt.xlabel("Molecula_rweight")
plt.ylabel("Frequency")
plt.grid()
plt.show()
共通の部分骨格を持つ分子を抜き出す
まずは取得したデータの確認として最初の8分子を可視化してみます。
from rdkit.Chem import Draw
# msの要素x(class：rdkit.Chem.rdchem.Mol)は.GetProp"_Name"関数で名前を習得できる
img=Draw.MolsToGridImage(ms[:8],molsPerRow=4,subImgSize=(200,200),legends=[x.GetProp("_Name") for x in ms[:8]])
img
次に共通の部分骨格をもつデータのみを抜き出していきます。今回は部分骨格として次の構造pを指定します。
p = Chem.MolFromSmiles('[nH]1cnc2cncnc21')
p
ちなみに、このような構造の分子あったっけな？と気になって調べてみましたが該当なしでした。こういったときは東京化成工業の構造式検索などが使えます。
取得したデータ（47個の分子が含まれる）の中で部分骨格として構造pを持つものを抜き出し、その数を確認します。
subms = [x for x in ms if x.HasSubstructMatch(p)]
len(subms)
14
47個の分子のうち、14個の分子が構造pを持つことが分かりました。共通骨格pでアライメントをとって、この14個を可視化していきます。
#共通骨格pでアライメントを取って可視化
from rdkit.Chem import AllChem
AllChem.Compute2DCoords(p)
for m in subms: AllChem.GenerateDepictionMatching2DStructure(m,p)
img=Draw.MolsToGridImage(subms,molsPerRow=4,subImgSize=(200,200),legends=[x.GetProp("_Name") for x in subms])
img
このようにデータセットの中からある共通の部分骨格を持つ分子を抜き出し、可視化することができます。
最大共通部分骨格を求める
続いて、分子間の最大共通部分骨格を求める方法を紹介いたします。上記で可視化した"subms"の5,6,7番目（上から2列目の左の3分子）の3分子の最大共通部分骨格(MCS：Maximum Common Substructure)を求めていきます。複数分子から最大共通部分骨格を求める際はFindMCS関数を使います。
from rdkit.Chem import rdFMCS
# MCSを求める3つの分子を定義
mol1 = subms[5]
mol2 = subms[6]
mol3 = subms[7]
mols = [mol1, mol2, mol3]
# MCSを求める（MCSはSMARTS文字列で返ってくる）
res = rdFMCS.FindMCS(mols)
#smart→molに変換
mcs = Chem.MolFromSmarts(res.smartsString)
mcs
この構造が3つの分子の最大共通部分骨格です。最後にこの最大共通部分骨格が各分子の中でどこに位置するかを可視化します。コードに関してですが、こちらを参照すると、「分子中の特定構造はGetSubstructMatchもしくはGetSubstructMatchesで取り出すことができ、GetSubstructMatchでは最初にマッチした部位の原子インデックスが、GetSubstructMatchesでは複数マッチした場合にも全ての原子インデックスがタプルとして取得できること」が分かりました。
mol1_match = mol1.GetSubstructMatch(mcs)
mol2_match = mol2.GetSubstructMatch(mcs)
mol3_match = mol3.GetSubstructMatch(mcs)
match = [mol1_match, mol2_match, mol3_match]
Draw.MolsToGridImage(mols, molsPerRow=len(mols), subImgSize=(400,400), highlightAtomLists=match)
最大共通部分骨格は"subms"の5番目の分子そのものであったことが分かります。今回検討した3つの分子は構造的にかなり近しい印象です。この3分子が数値的にどれぐらい類似しているのかについてこちらの投稿で議論しましたので是非ご参考にして下さい。
まとめ
RDKitの部分構造について、以下2点に焦点を当てて解説いたしました。
- ある共通の部分骨格を持つ分子を抜き出す方法
- 最大共通部分骨格を求める方法