はじめに
ケモインフォマティクスやマテリアルズインフォマティクスでは良く使用される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自体はRDKitやCDKに実装されおり、正確かつ高速に使うことができます。私の実装したものは、あくまで理解のためですので、アルゴリズムは多少間違っているかもしれません。
フィンガープリント
化学の分野で用いられるフィンガープリントは、分子構造を表現する方法の1つで2次元構造に基づいています。主に、分子の類似性検索や定量的構造活性・物性相関（QSAR, Quantitative Structure-Activity Relationship or QSPR, Quantitative Structure-Property Relationship）モデル構築に使用されます。
分子構造はSMILESと呼ばれる文字列情報で表現されます。また、分子はグラフ構造でも表せ、原子情報をノード、結合情報をエッジとして表現することもできます。原子＝ノード、結合＝エッジとして考えて構いません。
フィンガープリントで良く使用される方法はECFPです。また、2015年のニューラルネットワークによるフィンガープリント（Neural Graph Fingerprint, NFP）もあります。他にもグラフニューラルネットワークによる特徴量抽出など多数の方法がありますが、ECFPとNFPについて記載します。
図1はECFPとNFPの概念図です。左図はECFPの概略図で分子グラフのメッセージパッシングによりノード（原子の）情報が伝達され、最終的に固定長の0, 1のビットに変換されます。右図はNFPで情報伝達はニューラルネットワークで重み付けされ情報が伝達されます。
アルゴリズム
実際にアルゴリズムを見るのがわかりやすいと思います。
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)
自作版 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.Chem.rdmolops.GetAdjacencyMatrixだと上記の操作が行えます。
m = Chem.MolFromSmiles('c1ccccc1')
Chem.rdmolops.GetAdjacencyMatrix(m)
- 特徴量抽出
RDKitの便利な関数により、分子中の各原子ごとの情報を行列にしていきます。RDKitのECFPがどのような情報を用いているかわからなかったので、適当に選択しています。Module Hierarchyが参考になります。
def createNodeFeatures(mol) -> np.ndarray:
"""ノード特徴量を生成する
:param mol: rdkit mol object
:return: np.array
"""
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.GetAtomicNum(),
a.GetDegree(),
a.GetExplicitValence(),
a.GetImplicitValence(),
a.GetFormalCharge(),
a.GetTotalNumHs()] for a in mol.GetAtoms()], dtype=np.int32)
return features
特徴抽出の結果は、次のようになります。rowがノード数、columnが特徴数です。
>> print(createNodeFeatures(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 ECFP:
"""Extended Connectivity Fingerpritnsを計算する
:param mol: rdkit mol object
:param radius: 隣接する原子情報をどこまで見るか. メッセージパッシングの回数に相当する
:param nbits: 得られた部分構造の格納数するためのビット数
:param n_feat: 特徴量行列 (ノード数 × 特徴量数 )
"""
def __init__(self, mol, radius: int = 2, nbits: int = 2048, n_feat: np.array = None):
self.mol = mol
self.radius = radius
self.nbits = nbits
self.fps = np.zeros(shape=(self.nbits,), dtype=np.int32)
if n_feat is None:
n_feat = createNodeFeatures(mol)
n_feat = np.array(n_feat, dtype=np.int32)
n_atoms = n_feat.shape[0]
self.adj = Chem.GetAdjacencyMatrix(mol)
self.n_feat = np.array([''.join([str(f) for f in n_feat[atom]]) for atom in range(n_atoms)], dtype=str)
def _concat_neighbor(self, atom: int, n_feat: list) -> list:
"""隣接情報を付与する。メッセージパッシングとaggregationに相当する。concatenationで次元を減らす。
:param atom: 注目している原子のindex
:param n_feat: 更新前の特徴量情報
:return: ノード特徴量を返す
"""
nei_id = np.nonzero(self.adj[atom])[0]
vec = ''.join([str(ind) for ind in n_feat[nei_id]])
return vec
def calculate(self) -> np.ndarray:
"""ECFPを計算する
:return: result of fingerprint.
"""
n_atoms = self.n_feat.shape[0]
identifier = copy.deepcopy(self.n_feat)
for _ in range(0, self.radius):
for atom in range(n_atoms):
if _ == 0:
v = self.n_feat[atom]
else:
v = self._concat_neighbor(atom, identifier)
identifier[atom] = hash(v)
index = int(identifier[atom]) % self.nbits
self.fps[index] = 1
self.identifier = identifier
return self.fps
concat_neighbor()では、隣接する原子の情報を集めます。隣接のidと情報、結合次数をconcatenateしています。calculate()でECFPのアルゴリズム通り実装していきます。各原子ごとに隣接する特徴量を取得し、hash()を通し、これをその原子の識別子とします。これを固定長で割った余りが部分構造のindexになります。これを全ての原子に対して行います。メッセージパッシング操作を繰り返すほど、遠くの原子の情報を伝えることに相当します。
結果
ベンチマークで良く使用される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]
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
fp1 = np.vstack(list(map(lambda mol: AllChem.GetMorganFingerprintAsBitVect(mol, radius=2, nBits=2048, useFeatures=True), mols)))
run(fp1)
logPに関してはECFPだけでも予測できていそうです。すなわち、分子の部分構造から決まるということです。
自作 ECFP
def get_ecfp(mol, radius=2, nbits=2048):
n_feat = createNodeFeatures(mol)
obj = ECFP(mol, radius, nbits, n_feat)
return obj.calculate(conn_info)
fp2 = np.vstack(list(map(lambda dataset: get_ecfp(dataset, rad, nbits), mols)))
run(fp2)
精度はやはり自作の方が悪いです。また、radiusを増加させていくと精度が落ちていきます。また、RDKit版と自作版ではfingerprintの値は異なります。これは特徴量とハッシュ関数が異なるためと考えられます（もしくは、実装方法が異なるか・・・）。
最後に
やはりRDKitが便利なのですが、新しい手法は実装されていません。そのため、自分で実装していく必要があります。グラフネットワークを用いて実装したりしているのですが、実際の所、大抵の場合は求めたい精度に達しません。やはり楽な方法はなく、色々試行錯誤が必要となってきます。
NFPの方は試していませんが、DeepChemを用いると既に実装されています。Chainer Chemistryも気になっていますが、まだ試していません。こちらは、
- NFP
- GGNN
- Weave Net
- SchNet
が実装されています。試すだけなら良さそうですが、どうなんでしょうか。