🐰

構成元素を取り出し化合物を分類する

2023/01/06に公開

Python

chemistry

tech

はじめに

以前、こちらの記事で、「分子式から元素記号を取り出すロジックを思いつけない」と記載しましたが、正規表現を使って達成できたので、記事にしました。
「特定の元素を含む化合物を抽出する」際にも使用できますが、今回は、「候補化合物のリストから、指定した元素以外の元素を含む化合物を除く、というケースを取り上げます。
「C, H, N, Oの他、Sは含んで良いけど、SiやPを含む化合物は除きたい」といった場合を想定しています。本記事では、SMILES記法での構造式から選別しますが、分子式の場合は、前処理が不要のため、内容の2. 以降を参照ください。

概要

含んで良い元素をポジティブリストとして与え、SMILESや分子式に含まれる元素と比較した。
元素記号が「英字大文字一文字」または「英字大文字一文字＋小文字一文字（合計二文字）」であることを利用し、reモジュールでのパターンマッチングにより、SMILES表記の構造式（分子式でも可）から構成元素を取り出した。

内容

1. SMILESの場合（前処理）

芳香族化合物をSMILES記法で表す場合、元素記号を小文字で表記する場合があり、この際は前処理が必要です。RDkitで、デフォルトの設定でcanonicalizeした場合などがこれに該当し、芳香性を持つ該当元素が小文字になってしまいます（例：ベンゼン'C1=CC=CC=C1'は'c1ccccc1'と表記される）。
PubChemのcanonical SMILESなどは、ケクレ構造で表記されており（ベンゼンは'C1=CC=CC=C1'）、元素記号どおりの表記となっているため、本項の処理は不要です。
RDkitでも、Chem.MolToSmilesのオプションに'kekuleSmiles=True'を指定することで、ケクレ構造で出力できます。

from rdkit import Chem
m = 'C1=CC=CC=C1'
Chem.MolToSmiles(Chem.MolFromSmiles(m), kekuleSmiles=True)

[出力]: 'C1=CC=CC=C1'

化合物として、ギ酸、酢酸メチル、酢酸エチル、酢酸ナトリウム、酢酸トリメチルシリル、トリス(3-メチルフェニル)リン酸エステルを選び、データフレームとして下記を与えました。

import numpy as np
import pandas as pd
from rdkit import Chem
df = pd.DataFrame(['O=CO', 'COC(C)=O', 'CCOC(C)=O', 'CC(=O)[O-].[Na+]', 'CC(=O)O[Si](C)(C)C', 
                   'Cc1cccc(OP(=O)(Oc2cccc(C)c2)Oc2cccc(C)c2)c1'], columns=['compound',])
df

	compound
0	O=CO
1	COC(C)=O
2	CCOC(C)=O
3	CC(=O)[O-].[Na+]
4	CC(=O)O[Si](C)(C)C
5	Cc1cccc(OP(=O)(Oc2cccc(C)c2)Oc2cccc(C)c2)c1

このデータフレーム中のSMILESを、RDkitを使い、ケクレ構造のcanonical SMILESにしました。

df['compound'] = df['compound'].apply(lambda m : Chem.MolToSmiles(Chem.MolFromSmiles(m),kekuleSmiles=True))
df

	compound
0	O=CO
1	COC(C)=O
2	CCOC(C)=O
3	CC(=O)[O-].[Na+]
4	CC(=O)O[Si](C)(C)C
5	CC1=CC=CC(OP(=O)(OC2=CC=CC(C)=C2)OC2=CC=CC(C)=C2)=C1

No5の行の表記が変わっていることが確認できます。ここら辺の手順は、必要に応じて、前出の記事（RDKit:canonical SMILESを比較する、RDKit:PandasのDataFrameのSMILES列にapplyを適用してcanonicalizeする）も参照ください。
以降、このデータフレームに含まれている化合物を選別していきます。今回、この前処理以外では、RDkitは使用しません。

2. 手順の設定

以下の手順で実施しました。

化合物に含まれても良い元素のリスト（＝ポジティブリスト）を作成した
reモジュールで、正規表現での置換subを用い、SMILES表記から、大文字および小文字の英字以外の文字を除いた。
```
re.sub(r'[^a-zA-Z]', '', 'SMILES表記の構造式等')
```
pythonドキュメント re---正規表現操作によれば「集合の最初の文字が '^' なら、集合に含まれない全ての文字にマッチします。」とある。これにより、英大小文字以外を「入力なし」と置換することができた。結合を表す'='や'#'、数字を除くことにより、元素記号のみを取り出すことができる、という発想である
「大文字一文字＋小文字0文字または一文字（'?'で表現）」となるよう、re.findall関数と正規表現を使ってマッチング処理し、元素記号に相当する文字列をリストの要素として取得した。例えば
```
import re
re.findall('[A-Z][a-z]?', 'CC(=O)O[Si](C)(C)C')
```
[出力]: ['C', 'C', 'O', 'O', 'Si', 'C', 'C', 'C']
となる。ここで、元素記号が「英字大文字一文字」または「英字大文字一文字＋小文字一文字（合計二文字）」であることを利用している（将来三文字の新元素が発見？されると修正しなきゃ）。前回はこの部分が実現できなかった。
list型のままでは演算しにくいので、ポジティブリスト、候補化合物から取得した元素のリストともset型に変換し、両者の和をとったあと、ポジティブリストに含まれる要素を差し引いた集合を作った。ポジティブリストに含まれない元素が候補化合物に含まれる場合、その元素が新たなリストに残り、含まれない場合、新たなリストの要素の数は0になる。
新たなリストの長さが0かそれ以外かを判定し、bool値を返した。ポジティブリストに含まれる元素のみで候補化合物が構成される場合にはTrue、ポジティブリストに含まれない化合物を有する場合にはFalseを返す。

3. 自作関数化

項目2.の動作を自作関数にしたものが、下記elemchk。if文の部分以外は、一箇条が1行に対応しています。
ポジティブリストに含まれる元素はCHNOの他、Sとハロゲンにしました（自分の趣味^^）。このリストを書き換えれば、他の元素にも対応できます。

def elemchk(m):
    poslst = ['H','C','N','O','S','F','Cl','Br','I']
    m2 = re.sub(r'[^a-zA-Z]', '', m)
    melem = re.findall('[A-Z][a-z]?', m2)
    l = set(poslst+melem)-set(poslst)
    if len(l)==0:
        a=True
    else:
        a=False
    return a

4. 実行

項目1.で準備したデータフレームの'compound'列に、関数elemchkを適用し、ポジティブリストに含まれる元素だけで構成される候補化合物を抜き出しました。

df[df['compound'].apply(elemchk)]

	compound
0	O=CO
1	COC(C)=O
2	CCOC(C)=O

それぞれの化合物に対する判定結果を書き出すと以下のとおり。意図どおり判定できています。

df['class'] = df['compound'].apply(elemchk)
df

	compound	class
0	O=CO	True
1	COC(C)=O	True
2	CCOC(C)=O	True
3	CC(=O)[O-].[Na+]	False
4	CC(=O)O[Si](C)(C)C	False
5	CC1=CC=CC(OP(=O)(OC2=CC=CC(C)=C2)OC2=CC=CC(C)=...	False

まとめ

分子式から、pyrhonの正規表現を使い、構成元素を取り出しました。SMILES表記の構造式の場合、表記法によっては、RDkitを使わざるを得ませんが、単に分子式から構成元素を取り出すのであれば、標準モジュールのみでも実施可能な方法です。
ただし、分子式から、元素の原子数を合わせて取り出すには、もう一工夫必要で、前記事の方法の方が良いかもしれません。

余談

箇条書きの連番の途中にコードチャンクを入れる場合、行頭に半角スペースを４つ程度入れると、番号が続き、インデントするようになりました。この半角がないと、インデントしませんでした。内容の2.の部分です。