Extended Connectivite Fingerprintsを実装する


はじめに

ケモインフォマティクスやマテリアルズインフォマティクスでは良く使用されるExtended Connectivity Circular Fingerprints (ECFPs, Circular Fingerprintsとも呼ばれます) を実装してみました。本方法は、Extended-Connectivity Fingerprints, David Rogers and Mathew Hahn, JCIM 2010の論文によるものです。ただし、実際にフィンガープリント実装の参考にしたのはConvolutional Networks on Graphs for Learning Molecular Fingerprints, David K. Duvenaus et al., NIPS2015です。

ECFP自体はRDKitCDKに実装されおり、正確かつ高速に使うことができます。私の実装したものは、あくまで理解のためですので、アルゴリズムは多少間違っているかもしれません。


フィンガープリント

化学の分野で用いられるフィンガープリントは、分子構造を表現する方法の1つで2次元構造に基づいています。主に、分子の類似性検索や定量的構造活性・物性相関(QSAR, Quantitative Structure-Activity Relationship or QSPR, Quantitative Structure-Property Relationship)モデル構築に使用されます。

分子構造はSMILESと呼ばれる文字列情報で表現されます。また、分子はグラフ構造でも表せ、原子情報をノード、結合情報をエッジとして表現することもできます。原子=ノード、結合=エッジとして考えて構いません。

フィンガープリントで良く使用される方法はECFPです。また、2015年のニューラルネットワークによるフィンガープリント(Neural Graph Fingerprint, NFP)もあります。他にもグラフニューラルネットワークによる特徴量抽出など多数の方法がありますが、ECFPNFPについて記載します。

ecfp_ngf_image.PNG

図1はECFPNFPの概念図です。左図はECFPの概略図で分子グラフのメッセージパッシングによりノード(原子の)情報が伝達され、最終的に固定長の0, 1のビットに変換されます。右図はNFPで情報伝達はニューラルネットワークで重み付けされ情報が伝達されます。


アルゴリズム

実際にアルゴリズムを見るのがわかりやすいと思います。

ecfp_ngf.PNG


Circular fingerprints (ECFP)

分子の半径Rとfingerprintの長さSを決めます。この半径がノードに情報を渡す回数です。半径2なら2つ隣まで、R=3なら3分子離れたノード情報を伝達していく操作になります。fを0で初期化します。長さSは1024や2048bitsとなることが多いです。分子中の原子ごとにAtom Propertiesである特徴量をリストに追加していきます。例えば、


  • 原子番号

  • 隣接重原子(水素を除いた原子)の数

  • 原子価

  • 環構造であるかどうか

などです。他にも追加していくことができます。注目している原子の隣接する情報を集め、(concatenate)結合していきます。これをハッシュに通し、fingerprintの長さSで余りを求めます。この余りのインデックスを1を書き込みます。これを全原子に対して行い(for each atom in molecule)、何度も繰り返す(for 1 to R)ことで情報を伝達させることができます。この操作はMorganアルゴリズムとも呼ばれています。


Neural graph fingerprints (NFP)

ECFPとの違いは、concatenate以降のところがニューラルネットワークになります。単純に結合し、ハッシュを通して、その余りのインデックスを1としていました。NFPでは学習によってノード情報にランク付けし、ソフトマックスで確率的にインデックスを決定し、フィンガープリントの値としています。そのため、0, 1のベクトルではなく連続値となります。論文名のConvolutional Networks on Graphともあるように、 これはグラフニューラルネットワークの考え方です。現在はより一般的なフレームワーク(DeepChem, GraphNetsやPyTorch geometricなど)があるので、そちらで実装した方が良いかもしれません。


pythonコード

分子を扱う場合はやはり、RDKitがとても便利です。


RDKit版 ECFP

from __future__ import print_function

from rdkit import Chem
from rdkit.Chem import AllChem, Draw
import numpy as np

smiles = 'c1ccccn1'

mol = Chem.MolFromSmiles(smiles)
bit_morgan1 = {}
fp1 = AllChem.GetMorganFingerprintAsBitVect(mol, radius=2, nBits=12, bitInfo=bit_morgan1)
bit1 = list(fp1)
print(bit1)

RDKitを用いれば、このように簡単にSMILESからECFPを行うことができます。このフィンガープリントを説明変数として目的変数を回帰したものが、QSAR / QSPRと呼ばれるものになります。

>> [0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0]

参考までに分子を描画します。

def mols2grid_image(mols, molsPerRow=1):

mols = [e if e is not None else Chem.RWMol() for e in mols]
for mol in mols:
AllChem.Compute2DCoords(mol)

return Draw.MolsToGridImage(mols, molsPerRow=molsPerRow, subImgSize=(150,150))

mols2grid_image([mol], molsPerRow=2)

mol.png


自作版 ECFP

自作するためには、SMILESから原子の情報とエッジの情報を表す特徴を抽出しなければなりません。複数分子に対応できるようにmolオブジェクトをリストで渡すようにしていますが、その後の処理は行っていません。


  • 隣接行列

def GetAdjacencyMatrix(smols: list, connected=True):

max_length = max(mol.GetNumAtoms() for mol in mols)
bond_labels = [Chem.rdchem.BondType.ZERO] + list(sorted(set(bond.GetBondType() for mol in mols for bond in mol.GetBonds())))
bond_encoder_m = {l: i for i, l in enumerate(bond_labels)}

A = np.zeros(shape=(max_length, max_length), dtype=np.int32)
begin, end = [b.GetBeginAtomIdx() for b in mol.GetBonds()], [b.GetEndAtomIdx() for b in mol.GetBonds()]
bond_type = [bond_encoder_m[b.GetBondType()] for b in mol.GetBonds()]
A[begin, end] = bond_type
A[end, begin] = bond_type
return A

GetAdjacencyMatrixで隣接行列を取得します。smilesからChem.MolFromSmilesによりmolオブジェクトにした後、結合種類によってラベル付けをし、分子のインデックス、結合情報が得られるのでを書き込んでいきます。


  • 特徴量抽出

RDKitの便利な関数により、分子中の各原子ごとの情報を行列にしていきます。RDKitのECFPがどのような情報を用いているかわからなかったので、適当に選択しています。Module Hierarchyが参考になります。

def GetFeatures(mols):

max_length = max(mol.GetNumAtoms() for mol in mols)
features = np.array([[
*[a.GetDegree() == i for i in range(5)],
*[a.GetExplicitValence() == i for i in range(9)],
*[int(a.GetHybridization()) == i for i in range(1, 7)],
*[a.GetImplicitValence() == i for i in range(9)],
a.GetIsAromatic(),
a.GetNoImplicit(),
*[a.GetNumExplicitHs() == i for i in range(5)],
*[a.GetNumImplicitHs() == i for i in range(5)],
*[a.GetNumRadicalElectrons() == i for i in range(5)],
a.IsInRing(),
*[a.IsInRingSize(i) for i in range(2, 9)],
# Add Features
a.GetAtomicNum(),
a.GetDegree(),
a.GetExplicitValence(),
a.GetImplicitValence(),
a.GetFormalCharge(),
a.GetTotalNumHs()] for a in mol.GetAtoms()], dtype=np.int32)

return np.vstack((features, np.zeros((max_length - features.shape[0], features.shape[1]))))

特徴抽出の結果は、次のようになります。rowがノード数、columnが特徴数です。

>> print(GetFeatures([mol]))

[[0. 1. 0. ... 0. 3. 1.]
[0. 0. 1. ... 0. 2. 2.]
[0. 0. 0. ... 0. 0. 3.]
...
[0. 0. 1. ... 0. 1. 2.]
[0. 0. 0. ... 0. 0. 3.]
[0. 1. 0. ... 0. 0. 1.]]


  • MyECFPクラス

上記の隣接行列と原子の特徴行列からECFPを実装します。

class MyECFP:

def __init__(self):
self.A = None
self.F = None

def _neighbors_concat(self, a):
conc = []
idn = np.nonzero(self.A[a])[0]
for i, x in enumerate(self.r[idn]):
conc.append(str(i) + str(x) + str(self.A[a][idn][i]))

vec = ""
for ind in conc:
vec += ind

return vec

def _features_concat(self, F):
v=[]
F = np.array(F, dtype=np.int32)
for i in range(self.A.shape[0]):
tmp = ""
for j in F[i]:
tmp += str(j)
v.append(tmp)
v = np.array(v, dtype="str")
return v

def GetECFP(self, mols, radius=2, nBits=1024):
self.A = GetAdjacencyMatrix(mols)
self.F = GetFeatures(mols)

# add connectivity info
conn = np.array([len(np.nonzero(self.A[i])[0]) for i in range(self.A.shape[0])])
self.F = np.c_[self.F, conn]

self.r = self._features_concat(self.F)
#self.r = np.sum(self.F, axis=1)
fp = np.zeros(shape=(nBits,), dtype=np.int32)
for ite in range(0, radius):
for a in range(self.A.shape[0]):
if ite==1:
v = self.r[a]
else:
v = self._neighbors_concat(a)
self.r[a] = hash(v)
index = int(self.r[a]) % nBits
fp[index] = 1
return fp

_neighbors_concat()では、隣接する原子の情報を集めます。隣接のidと情報、結合次数をconcatenateしています。_features_concat()では、特徴量行列の複数の特徴を1つにまとめています。このconcatの部分をニューラルネットワークにした方が良いのは明らかです。GetECFP()でアルゴリズム通り実装していきます。


結果

ベンチマークで良く使用されるzincデータを用います。3500万の化合物でデータセットです。データ読込みし、logP(溶解性の指標)に対してRandomForest()で回帰してみます。

import pandas as pd

from sklearn.ensemble.forest import RandomForestRegressor
from sklearn.model_selection import train_test_split
import matplotlib.pyplot as plt

df = pd.read_csv("./data/zinc.csv")
logP = df.logP[0:5000]
smiles = df.smiles[0:5000]

mols = [Chem.MolFromSmiles(smi) for smi in smiles]
fp1 = np.array([list(AllChem.GetMorganFingerprintAsBitVect(mol, radius=1, nBits=1024)) for mol in mols])

def run(fp):
X_train, X_test, y_train, y_test = train_test_split(fp, logP, test_size=0.1, random_state=42)
rf = RandomForestRegressor(n_estimators=10)
rf.fit(X_train, y_train)
y_train_pred = rf.predict(X_train)
y_test_pred = rf.predict(X_test)

plt.plot(y_train, y_train_pred, "ro", label="train")
plt.plot(y_test, y_test_pred, "bo", label="test")
plt.plot(logP, logP, "k")
plt.legend()

plt.show()


RDKit版 ECFP

run(fp)

rdkit_ecfp1.png

logPに関してはECFPだけでも予測できていそうです。すなわち、分子の部分構造から決まるということです。


自作 ECFP

test = ECFP()

fp2 = np.array([test.GetECFP(mol, radius=2, nBits=1024) for mol im mols])
run(fp2)

my_ecfp.png

精度はやはり自作の方が悪いです。また、radiusを増加させていくと精度が落ちていきます。また、RDKit版と自作版ではfingerprintの値は異なります。これは特徴量とハッシュ関数が異なるためと考えられます(もしくは、実装方法が異なるか・・・)。


最後に

やはりRDKitが便利なのですが、新しい手法は実装されていません。そのため、自分で実装していく必要があります。グラフネットワークを用いて実装したりしているのですが、実際の所、大抵の場合は求めたい精度に達しません。やはり楽な方法はなく、色々試行錯誤が必要となってきます。

NFPの方は試していませんが、DeepChemを用いると既に実装されています。Chainer Chemistryも気になっていますが、まだ試していません。こちらは、


  • NFP

  • GGNN

  • Weave Net

  • SchNet

が実装されています。試すだけなら良さそうですが、どうなんでしょうか。


参考文献