RDKitで原子の特徴を計算しPandasに取り込んで使い倒そう

#はじめに
以前、原子の特徴を計算するプログラムを作成してみた という記事を書いたがそこではCDKを使った。今回はRDKitでそれをやってみたい。

#環境

• RDKit 2021.09.2

#対象記述子

• MorganFingerprint
• AtomPairFingerprint

対象記述子とはいいながらもフィンガープリントである。
MrogainFingerprintはある原子を中心とした一定の半径内の環境を表現するものであり、AtomPairFingerprintは原子同士のペアを表現するものである。これらは通常は分子を指定して計算するが、引数を指定することで特定の原子に関連するものだけを計算することもできる。

今回は原子単位でそれらのフィンガープリントを計算し、結果をPandasのデータフレームに取り込む方法を解説する。

#ソース

まず、いつものようにライブラリをインポートしよう。

import numpy as np
import pandas as pd
from rdkit import rdBase, Chem
from rdkit import DataStructs
from rdkit.DataStructs import ExplicitBitVect
from rdkit.Chem import AllChem, Descriptors
from rdkit.ML.Descriptors import MoleculeDescriptors                  
from collections import defaultdict

続いてsmilesの格納されたcsvを読み込もう。

# ファイルの読み込み
import codecs
with codecs.open('./datas/bace.csv', 'r', 'utf-8', 'ignore') as f:
    df = pd.read_csv(f)

続いてSMILESと名前を取り出そう

# SMILESのデータを取り出す
smiles_list = df["mol"].values
# 名前を取り出す
names_list = df["CID"].values

RDKitのMOLオブジェクトに変換しよう。

# SMILESのリストをMOLに変換
mols = []

for i, smiles in enumerate(smiles_list):
    mol = Chem.MolFromSmiles(smiles)
    if mol:
        mols.append(mol)
    else:
        print("{0} th mol is not valid".format(i))

さぁ、いよいよFingerPrintの計算だ。やり方は、分子内の原子をGetAtomsメソッドで順次取り出し、fromAtomsに原子のインデックスを指定して、フィンガープリント計算関数を呼び出すだけだ。

mg_fps_tmp = []
ap_fps_tmp = []
mg_keys = set()
ap_keys = set()
atom_names = []

# 分子でループ
for i, (mol, name) in enumerate(zip(mols, names_list)):   
    # 分子内の原子でループ
    for atom in mol.GetAtoms():
        # 一意な原子名を付ける
        atom_name = str(name) + "_" + str(atom.GetIdx())
        atom_names.append(atom_name)
        # MorganFingerprintの計算
        mg_fp = AllChem.GetMorganFingerprint(mol, 2, fromAtoms=[atom.GetIdx()])
        mg_fps_tmp.append(mg_fp.GetNonzeroElements())
        for key in mg_fp.GetNonzeroElements():
            mg_keys.add(key)
        # AtomPairFingerprintの計算
        ap_fp = AllChem.GetAtomPairFingerprint(mol, fromAtoms=[atom.GetIdx()])
        ap_fps_tmp.append(ap_fp.GetNonzeroElements())
        for key in ap_fp.GetNonzeroElements():
            ap_keys.add(key)

#結果をまとめる
fp_dict = defaultdict(list)
for i, (mg_fp, ap_fp) in enumerate(zip(mg_fps_tmp, ap_fps_tmp)):
    print(i, len(mg_fps_tmp))
    for mg_key in mg_keys:
        if mg_key in mg_fp:
            fp_dict["mg." + str(mg_key)].append(mg_fp[mg_key])
        else:
            fp_dict["mg." + str(mg_key)].append(0)

    for ap_key in ap_keys:
        if ap_key in ap_fp:
            fp_dict["ap." + str(ap_key)].append(ap_fp[ap_key])
        else:
            fp_dict["ap." + str(ap_key)].append(0)
            
#データフレームの構築
df_fp = pd.DataFrame(data=fp_dict, index=atom_names)

このデータフレームを表示すると、こんな感じで原子毎にフィンガープリントが得られる。

さぁ、これで原子の性質を予測するモデルを作ってみよう！

#おわりに
今回は2種類のフィンガープリントで試したが、他のフィンガープリントもできるかもしれない。