【RDKit】化合物中の特定の部分構造の数を数える②（エステル結合編）

## はじめに
**化合物中の特定の部分構造（官能基など）の数を数えていきたい**と思います。[その①（カルボキシル基編）](https://qiita.com/oki_kosuke/items/45952ae5299f404d2e3a)の続編となりますので、その①をお読みになっていな方はその①も参考にしてください。今回は自ら用意した化合物中のエステル結合の数を数えてみたいと思います。
この記事に関しては[インフォマジシャン研究所の藤様](https://hideyoshifuji.com/)にいろいろとアドバイスを頂いた内容を投稿しています。（掲載の許可はいただいております）。インフォマジシャン研究所ではインフォマティクスを専門としている方と議論することができます！

※マテリアルズインフォマティクス関係の内容を他にも投稿していますので、よろしければ[こちら](https://qiita.com/oki_kosuke/items/34b1071224953776c211#%E6%8A%95%E7%A8%BF%E4%B8%80%E8%A6%A7)の一覧から他の投稿も見て頂けますと幸いです。

## SMILESとSMARTS
本題に入る前に「SMILES」と「SMARTS」について整理していきましょう。一般的な内容は次のサイトなどを参考にしてください。
- [RDKitの公式](https://www.rdkit.org/docs_jp/Supplementary_Text_jp.html)
- [MOLSIS社のHP](https://www.molsis.co.jp/datamanagement/daylight/smiles/)
- [化学の新しいカタチ](https://future-chem.com/smiles-smarts/#toc2)

**SMILES**は分子の表記方法でDr. David Weiningerが考案した分子記述言語です。
![sample.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1057114/66986092-2165-f5c5-083c-9282fa0a335c.png)

上記の分子は次のSMILESで表現されます。
```
ClC(C)(C)OC(/C=C/C1=CN=CC=C1)=O
```
**SMARTS**はSMILES言語を検索クエリ用に拡張した構造検索のための表記方法です。例えば[C,N]は脂肪族炭素または脂肪族窒素を表し、SMILESよりもより曖昧な表現が可能です。SMARTSに関する詳細は[DAYLIGHT社のHP](https://www.daylight.com/dayhtml_tutorials/languages/smarts/)が参考になります。[こちら](https://www.daylight.com/dayhtml/doc/theory/theory.smarts.html)のリンクにも詳しく記載があります。

またSMARTSは[SMARTS PLUS](https://smarts.plus/)で確認することが可能です。SMARTSはSMILESよりも曖昧に表記できて検索性が高い分、使いこなすのはかなり難しい印象です。

### 環境
- python 3.7.13
- rdkit 2022.03.4

## 化合物中のエステル結合を数えていく
### 化合物の準備（実際にデータがある場合は不要）
まずは自分でサンプルとなる化合物を作っていきます。
```python
from rdkit import Chem

# それぞれのエステル結合は0,1,2,1,2個
m0 ="CC(OC)OC"
m1 ="CC(OCC(C)=O)=O"
m2 ="COC(C1=CC=CC=C1C(OC)=O)=O"
m3 ="O=C(OC)OC"
m4 ="O=C(OO1)C2=CC=CC=C2C1=O"
```
それぞれの化合物を可視化します。
```python
Chem.MolFromSmiles(m0)
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1057114/9f46375e-c714-dd76-f919-9493324198fe.png)

```python
Chem.MolFromSmiles(m1)
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1057114/20820f6c-3b67-274c-372a-e1e51d48edc5.png)


```python
Chem.MolFromSmiles(m2)
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1057114/3b991eab-e569-536a-c4c1-f52d8086e5dc.png)


```python
Chem.MolFromSmiles(m3)
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1057114/9da45554-446a-9457-47f5-05abb1d877b1.png)


```python
Chem.MolFromSmiles(m4)
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1057114/e165a4a7-e8f1-085e-36c4-7bf4587b640d.png)


それぞれの化合物を一括で処理できるようにデータフレームに格納します。
```python
import pandas as pd

df = pd.DataFrame(
    data={"compound": ["m0", "m1", "m2", "m3", "m4"], 
          "smiles": [m0, m1, m2, m3, m4]}
)
df
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1057114/f2b016a9-1e2e-ebc4-c0d5-e9b9c90bd729.png)

**実際にデータがある場合はこのようなデータフレームを読み込んでここから実行していくことが想定されます**。次にsmilesをMolに変換してリストで保持します。

```python
mol_list = [Chem.MolFromSmiles(x) for x in list(df["smiles"]) if x is not None]
```
### 部分構造の探索（その１）
検索する部分構造をSMARTS表記で定義します。
```python
q1 = Chem.MolFromSmarts("*[CX3](=O)[OX2H0]*")
q1
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1057114/bb4e49b0-4b62-f45b-7f34-a3261ed68413.png)

GetSubstructMatchesは指定した部分構造がマッチする原子インデックスが返ってくるので、これをlenで取り出して部分構造の数を数えます。今回は一括で処理します。

```python
[len(x.GetSubstructMatches(q1)) for x in mol_list]
```
```:出力
[0, 1, 2, 2, 0]
```
m3のエステル構造が2、m4のエステル構造が0と認識されている点が問題です。m3についてはC=Oの左右にOがあるので、ダブルカウントされていることが考えられます。まずはm4でエステル構造が0と認識れている問題について考察していきます。

### 芳香族性の認識
m4のSMILESを確認してみます。
```python
Chem.MolToSmiles(mol_list[4])
```
```:出力
'O=c1ooc(=O)c2ccccc12'
```
このことからm4の左右の環状構造はともに芳香族として認識されていることが分かります。左側の環状構造はベンゼン環ではないなの、なぜ芳香族として認識されるのでしょうか？これはRDKitの芳香族の認識のルールに起因します。[こちら](https://www.rdkit.org/docs_jp/The_RDKit_Book_jp.html)に詳細は載っていますが、「芳香族性は単純であると同時にどうしようもなく複雑な不愉快な話題の一つです。」との記載があり、化学者の直感的な理解とは異なるようです。

### 部分構造の探索（その２）
上記の内容を踏まえて、例えば炭素原子はC（脂肪族中の炭素）ではなく#6で定義することで芳香族の問題をクリアしてみようと思います。
```python
q2 = Chem.MolFromSmarts("*[#6;X3](=O)[#8]*")
q2
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1057114/41388355-dba1-b7a0-9d3a-01c6b332cf34.png)

```python
[len(x.GetSubstructMatches(q2)) for x in mol_list]
```
```:出力
[0, 1, 2, 2, 2]
```
これでm4のエステル構造が認識されないという問題はクリアしました。次はm3のエステル結合のダブルカウントに関してです。

### 部分構造の探索（その２）
m3のエステル結合を1つとしてカウントするためには、C=Oだけを数えるとケトンもカウントされてしまうので、-O-C(=O)-O-もカウント対象として対処します。
```python
Chem.MolFromSmarts("[#6;X3](=O)([#6])([#8])")
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1057114/52da2acd-e901-f852-9db0-9a58f673e201.png)
```python
Chem.MolFromSmarts("[#6;X3](=O)([#8])([#8])")
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1057114/968c6f4e-815c-134f-aef9-d26e25bf7d3d.png)
この2つをorでカウントできるように、以下のように部分構造を定義します。
```python
q3 = Chem.MolFromSmarts("[$([#6;X3](=O)([#6])([#8])),$([#6;X3](=O)([#8])([#8]))]")
[len(x.GetSubstructMatches(q3)) for x in mol_list]
```
```:出力
[0, 1, 2, 1, 2]
```
これで上手く化合物中のエステル結合の数をカウントすることができました。

## まとめ
今回はエステル結合を例として化合物中の部分構造を数える方法を紹介いたしました。

### 参考
- [RDKitの公式](https://www.rdkit.org/docs_jp/Supplementary_Text_jp.html)
- [MOLSIS社のHP](https://www.molsis.co.jp/datamanagement/daylight/smiles/)
- [化学の新しいカタチ](https://future-chem.com/smiles-smarts/#toc2)
- [DAYLIGHT社のHP](https://www.daylight.com/dayhtml_tutorials/languages/smarts/)