どの構造が寄与しているか解釈したい！ 〜Morgan Fingerprintの部分構造を特徴量に〜

こんにちはnissyです

定量的構造物性相関を取る上でどの構造が効いているか、直感的に解釈したい場面ありますよね。

そんな時に役に立つ手法を紹介しますので解説します！！！

分子内の特徴の重複回数を考慮してFingerprintを生成

分子内の特徴の重複回数を考慮してFingerprintを生成

早速、ライブラリーをインポートします。

import os
import pandas as pd
import numpy as np
import matplotlib.figure as figure
import matplotlib.pyplot as plt
from rdkit import Chem
from rdkit.Chem import AllChem, Draw
from rdkit.Chem.Draw import rdMolDraw2D
from IPython.display import SVG, display
import glob
from rdkit.Chem import AllChem
from collections import Counter

import os

import pandas as pd

import numpy as np

import matplotlib.figure as figure

import matplotlib.pyplot as plt

from rdkit import Chem

from rdkit.Chem import AllChem, Draw

from rdkit.Chem.Draw import rdMolDraw2D

from IPython.display import SVG, display

import glob

from rdkit.Chem import AllChem

from collections import Counter

データセットの準備

SMILESのフィンガープリントに変換するまでは、過去のブログで紹介してるのでそちらも参照してください。

#==========================
# 学習データのサンプル数が50になるように調整
#==========================
df = pd.read_csv('SMILES_logS_データセット.csv')
df_train, df_test= train_test_split(df, train_size=0.039, random_state=42)

#===========
# 目的変数の中心化　標準化をすると学習精度が悪くなる
#===========
y_name = "logS"
train_y = df_train[y_name]
train_y_ave = train_y.mean()
centered_train_y = train_y - train_y_ave

#================================
# SMILES⇒ECFP4への変換
#================================
radius= 3
bit =8192

mols_list = [Chem.MolFromSmiles(smiles) for smiles in df_train["SMILES"]]
ECFP = [AllChem.GetMorganFingerprintAsBitVect(mol, radius,bit) for mol in mols_list]
df_X = pd.DataFrame(np.array(ECFP, int))

#==========================

# 学習データのサンプル数が50になるように調整

#==========================

df = pd.read_csv('SMILES_logS_データセット.csv')

df_train, df_test= train_test_split(df, train_size=0.039, random_state=42)

#===========

# 目的変数の中心化　標準化をすると学習精度が悪くなる

#===========

y_name = "logS"

train_y = df_train[y_name]

train_y_ave = train_y.mean()

centered_train_y = train_y - train_y_ave

#================================

# SMILES⇒ECFP4への変換

#================================

radius= 3

bit =8192

mols_list = [Chem.MolFromSmiles(smiles) for smiles in df_train["SMILES"]]

ECFP = [AllChem.GetMorganFingerprintAsBitVect(mol, radius,bit) for mol in mols_list]

df_X = pd.DataFrame(np.array(ECFP, int))

radiusが大きい方が広く部分構造をとります。今回は3にしています。

カウントを有効にしたモルガンフィンガープリント

#==================================================
# カウントを有効にしたモルガンフィンガープリントの変換
#================================
# カウント付きのモルガンフィンガープリントを計算
fp_list = []
for mol in mols_list:
    fp = AllChem.GetMorganFingerprint(mol, radius, useCounts=True)
    fp_counts = Counter(fp.GetNonzeroElements())
    fp_list.append(fp_counts)
df_fp_count = pd.DataFrame(fp_list).fillna(0)  # fillna(0) で出現しないフィンガープリントを 0 にする
df_fp_count.head()

#==================================================

# カウントを有効にしたモルガンフィンガープリントの変換

#================================

# カウント付きのモルガンフィンガープリントを計算

fp_list = []

for mol in mols_list:

fp = AllChem.GetMorganFingerprint(mol, radius, useCounts=True)

fp_counts = Counter(fp.GetNonzeroElements())

fp_list.append(fp_counts)

df_fp_count = pd.DataFrame(fp_list).fillna(0) # fillna(0) で出現しないフィンガープリントを 0 にする

df_fp_count.head()

df_fp_countを開くと以下のようなデータフレームが得られてることが確認できます。

カラムにはハッシュ値が特徴量名として入力されています。

ある化学構造に対して、どのような部分構造をとっているかを確認します。

num = 0 #元データのインデックスナンバー


# 分子例（mols_listの一部を使って例示）
mol = mols_list[num]  # 任意の分子を使います
print(df_fp_count.iloc[num,:])
# フィンガープリントと情報を同時に取得
info = {}
fp = AllChem.GetMorganFingerprint(mol, radius, bitInfo=info)  # bitInfoに情報を保存

# コンテナの準備
view = rdMolDraw2D.MolDraw2DSVG(500,300)
tm = rdMolDraw2D.PrepareMolForDrawing(mol)
view.DrawMolecule(tm)
view.FinishDrawing()
svg = view.GetDrawingText()
display(SVG(svg)) か

num = 0 #元データのインデックスナンバー

# 分子例（mols_listの一部を使って例示）

mol = mols_list[num] # 任意の分子を使います

print(df_fp_count.iloc[num,:])

# フィンガープリントと情報を同時に取得

info = {}

fp = AllChem.GetMorganFingerprint(mol, radius, bitInfo=info) # bitInfoに情報を保存

# コンテナの準備

view = rdMolDraw2D.MolDraw2DSVG(500,300)

tm = rdMolDraw2D.PrepareMolForDrawing(mol)

view.DrawMolecule(tm)

view.FinishDrawing()

svg = view.GetDrawingText()

display(SVG(svg)) か

確認する構造は以下を選定しています。

どのような部分構造が取られたかを確認して行きましょう！

# ハッシュ値と対応する部分構造を取得して表示
substructures = []
legends = []  # ハッシュ値をラベルとして保持するリスト

for bit_id, atom_info in info.items():
    for atom_idx, rad in atom_info:
        if rad == radius:
            env = Chem.FindAtomEnvironmentOfRadiusN(mol, radius, atom_idx)
            atoms = set()
            for bond_id in env:
                bond = mol.GetBondWithIdx(bond_id)
                atoms.add(bond.GetBeginAtomIdx())
                atoms.add(bond.GetEndAtomIdx())
            sub_mol = Chem.PathToSubmol(mol, env)
            substructures.append(sub_mol)
            legends.append(f"BitID: {bit_id}")  # ハッシュ値をラベルに追加
print(len(legends))

# 部分構造を描画して、ハッシュ値とともに表示
img = Draw.MolsToGridImage(substructures, molsPerRow=4, subImgSize=(200, 200), legends=legends)

# Jupyter Notebookで画像を表示
display(img)

# ハッシュ値と対応する部分構造を取得して表示

substructures = []

legends = [] # ハッシュ値をラベルとして保持するリスト

for bit_id, atom_info in info.items():

for atom_idx, rad in atom_info:

if rad == radius:

env = Chem.FindAtomEnvironmentOfRadiusN(mol, radius, atom_idx)

atoms = set()

for bond_id in env:

bond = mol.GetBondWithIdx(bond_id)

atoms.add(bond.GetBeginAtomIdx())

atoms.add(bond.GetEndAtomIdx())

sub_mol = Chem.PathToSubmol(mol, env)

substructures.append(sub_mol)

legends.append(f"BitID: {bit_id}") # ハッシュ値をラベルに追加

print(len(legends))

# 部分構造を描画して、ハッシュ値とともに表示

img = Draw.MolsToGridImage(substructures, molsPerRow=4, subImgSize=(200, 200), legends=legends)

# Jupyter Notebookで画像を表示

display(img)

BitIDのファイル名でデータセット内の部分構造の画像を保存

下記コードを走らせることで、データセットないの部分構造の画像が保存されます。

#=====================================
# 画像の保存先のフォルダーを指定
#=====================================
output_dir = r" " #保存先のフォルダパスを入れる
os.makedirs(output_dir, exist_ok=True)

for bit_id, atom_info in info.items():
    for atom_idx, rad in atom_info:
        if rad == radius:
            # 部分構造の環境を取得
            env = Chem.FindAtomEnvironmentOfRadiusN(mol, radius, atom_idx)
            atoms = set()
            for bond_id in env:
                bond = mol.GetBondWithIdx(bond_id)
                atoms.add(bond.GetBeginAtomIdx())
                atoms.add(bond.GetEndAtomIdx())
            
            # 部分構造のサブ分子を取得
            sub_mol = Chem.PathToSubmol(mol, env)
            
            # 部分構造を画像として保存
            img = Draw.MolToImage(sub_mol, size=(200, 200), legend=f"BitID: {bit_id}")
            img.save(os.path.join(output_dir, f"{bit_id}.png"))

            # 画像を保存するごとに確認表示（Jupyter Notebook用）
            display(img)