【RDKit】化合物中の特定の部分構造の数を数える①（カルボキシル基編）

## はじめに
**化合物中の特定の部分構造（官能基など）の数を数えていきたい**と思います。今回は例として、自ら用意した化合物中のカルボキシル基の数を数えてみたいと思います。例えばこれができれば、スクリーニングしてきた化合物の中からカルボキシル基を1つだけ持つ化合物を抽出したいといったシーンで役立つはずです。

※マテリアルズインフォマティクス関係の内容を他にも投稿していますので、よろしければ[こちら](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)のリンクにも詳しく記載があります。※なお、このリンクは[インフォマジシャン研究所の藤様](https://hideyoshifuji.com/)に教えていただきました

また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,0個持つ化合物を定義します
m0 = "CC(C1=CC=C(C(OC)=O)C=C1)=O"
m1 = "OC(C1=CC=C(C(OC)=O)C=C1)=O"
m2 = "OC(C1=CC=C(C(O)=O)C=C1)=O"
m3 = "CC(O)=O"
m4 = "OOCC"
```
それぞれの化合物を可視化します。
```python
Chem.MolFromSmiles(m0)
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1057114/a88866fe-2b52-f728-4008-caff9a092854.png)

```python
Chem.MolFromSmiles(m1)
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1057114/8bd9433c-1f8e-1d2b-6a6c-097373d00a17.png)

```python
Chem.MolFromSmiles(m2)
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1057114/b0b6be64-c964-b93c-7b90-6b1186483490.png)

```python
Chem.MolFromSmiles(m3)
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1057114/316409de-eb8a-3fc1-0142-25020750a88e.png)

```python
Chem.MolFromSmiles(m4)
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1057114/7ff0e3fb-6a75-4c65-898e-3b45d31ed128.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/b3888a58-dcbe-54eb-ce32-3fad182a15f8.png)
**実際にデータがある場合はこのようなデータフレームを読み込んでここから実行していくことが想定されます**。次にsmilesをMolに変換してリストで保持します。

```python
mol_list = [Chem.MolFromSmiles(x) for x in list(df["smiles"]) if x is not None]
```
### 数える部分構造の定義
今回はカルボキシル基をSMARTS表記で定義します。
```python
q = Chem.MolFromSmarts("*C(=O)-[OH1]")
q
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1057114/b44445ed-0b8a-c1da-80bd-43f2e74e9026.png)
### 部分構造の探索
HasSubstructMatchは指定した部分構造を持つか否かで返ってきます。[参考](https://future-chem.com/substruct-search-mcs/) 
```python
mol_list[0].HasSubstructMatch(q)
```
```:出力
False
```
```python
mol_list[1].HasSubstructMatch(q)
```
```:出力
True
```
0番目の分子はカルボキシル基を含み、1番目の分子はカルボキシル基を含まないのでこのような出力になります。

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

```python
import numpy as np

for i in np.arange(len(df)):
    print(len(mol_list[i].GetSubstructMatches(q)))
```
```:出力
0
1
2
1
0
```

### 実用を意識した例
実際の使用を想定して、部分構造を指定した個数分だけ持っている化合物のみを抽出していきます。ここでは、「カルボキシル基を1つだけ持つ化合物」を選抜していくことにします。
```python
# 着目する構造な数を指定する
num = 1

index_list = []

for i in np.arange(len(df)):
    count = len(mol_list[i].GetSubstructMatches(q))  # 着目する構造の数をカウント
    if count == num:
        index_list.append(i)  # 条件を満たす番号をリストに格納
        
new_df = df.iloc[index_list,:]
new_df
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1057114/faad17bd-2778-b078-913a-df6f7f950c78.png)
このようにカルボキシル基を1つだけ持つリストの1番目と3番目の化合物のみが抽出されました。

## まとめ
今回はカルボキシル基を例として化合物中の部分構造を数える方法を紹介いたしました。次は別の部分構造でトライしてみたいと思います。

### 参考
- [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/)