Extended Connectivite Fingerprintsを実装する

## はじめに

ケモインフォマティクスやマテリアルズインフォマティクスでは良く使用されるExtended Connectivity Circular Fingerprints (ECFPs, Circular Fingerprintsとも呼ばれます) を実装してみました。本方法は、[Extended-Connectivity Fingerprints, David Rogers and Mathew Hahn, JCIM 2010](<https://pubs.acs.org/doi/10.1021/ci100050t>)の論文によるものです。ただし、実際にフィンガープリント実装の参考にしたのは[Convolutional Networks on Graphs for Learning Molecular Fingerprints, David K. Duvenaus et al., NIPS2015](<https://papers.nips.cc/paper/5954-convolutional-networks-on-graphs-for-learning-molecular-fingerprints>)です。



ECFP自体は[RDKit](<https://www.rdkit.org/>)や[CDK](<https://cdk.github.io/>)に実装されおり、正確かつ高速に使うことができます。私の実装したものは、あくまで理解のためですので、アルゴリズムは多少間違っているかもしれません。



## フィンガープリント

化学の分野で用いられるフィンガープリントは、分子構造を表現する方法の1つで2次元構造に基づいています。主に、分子の類似性検索や定量的構造活性・物性相関（**QSAR**, **Q**uantitative **S**tructure-**A**ctivity **R**elationship or **QSPR**, **Q**uantitative **S**tructure-**P**roperty **R**elationship）モデル構築に使用されます。



分子構造は[SMILES](<https://en.wikipedia.org/wiki/Simplified_molecular-input_line-entry_system>)と呼ばれる文字列情報で表現されます。また、分子はグラフ構造でも表せ、原子情報をノード、結合情報をエッジとして表現することもできます。原子＝ノード、結合＝エッジとして考えて構いません。



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

![ecfp_ngf_image.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/230659/28f5d6b9-79ae-7cc1-0fa0-99a31ab34226.png)


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



## アルゴリズム

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

![ecfp_ngf.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/230659/afbcf314-a13b-d97d-1813-452ff7d7d302.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

```python
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と呼ばれるものになります。

```python
>> [0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0]
```



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

```python
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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/230659/34d709bd-19b4-854f-840f-c32bfe5b8762.png)


### 自作版 ECFP

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

- 隣接行列

```python
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`だと上記の操作が行えます。

```python

m = Chem.MolFromSmiles('c1ccccc1')
Chem.rdmolops.GetAdjacencyMatrix(m)
```



- 特徴量抽出

RDKitの便利な関数により、分子中の各原子ごとの情報を行列にしていきます。RDKitのECFPがどのような情報を用いているかわからなかったので、適当に選択しています。[Module Hierarchy](<http://www.rdkit.org/docs-beta/api/rdkit.Chem.rdMolDescriptors-module.html#GetMorganFingerprintAsBitVect>)が参考になります。

```python
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が特徴数です。

```python
>> 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を実装します。

```python

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データ](<https://zinc.docking.org/>)を用います。3500万の化合物でデータセットです。データ読込みし、`logP(溶解性の指標)`に対して`RandomForest()`で回帰してみます。


```python
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

```python
fp1 = np.vstack(list(map(lambda mol: AllChem.GetMorganFingerprintAsBitVect(mol, radius=2, nBits=2048, useFeatures=True), mols)))
run(fp1)
```



![rdkit_ecfp1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/230659/c94c77db-69c4-a0b4-1d5a-3da263a93492.png)

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


### 自作 ECFP

```python
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)
```



![my_ecfp.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/230659/7f2cf40f-8894-8072-61c1-1e0e329115ef.png)


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



## 最後に

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



NFPの方は試していませんが、DeepChemを用いると既に実装されています。[Chainer Chemistry](<https://research.preferred.jp/2017/12/chainer-chemistry-beta-release/>)も気になっていますが、まだ試していません。こちらは、

- NFP
- GGNN
- Weave Net
- SchNet

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


## 参考文献

- [CheｍAxon ECFP](<https://docs.chemaxon.com/display/docs/Extended+Connectivity+Fingerprint+ECFP)
- [RDKit Fingerprints Document](<https://www.rdkit.org/UGM/2012/Landrum_RDKit_UGM.Fingerprints.Final.pptx.pdf>)

- [化合物でもDeep Learningがしたい](<https://kivantium.net/deep-for-chem>)