RDKit:PandasのDataFrameのSMILES列にapplyを適用してcanonicalizeする
背景
化学者のDr_Oです。有機合成化学と呼ばれる分野で仕事しています。
pythonやRを使ったデータ解析について、自分も勉強しながらですが、化学に登場しがちな場面で気づいたこと、経験したことを共有できたらと思い、記事を書き始めました。
はじめに
前回に引き続き、分子の構造を表すSMILESに関する記事です。化合物群をPandasのDataFrameで管理する際、SMILESの列を一気にcanonicalizeする場合について、記載します。
１つのSMILESをrdkitでcanonicalizeすることや、Pandasのデータフレームの特定の列に関数を適用して一気に変換することは、それぞれ色々なところで解説されていますが、これを組み合わせた場合の記述はあまりないように思ったので、記事にしました。
概要
lambda式と三項演算子を組み合わせることにより、より少ない行数で、欠損値を持つpandasのDataFrameのSMILES列に対し、rdkitの関数をapplyを用いて適用し、canonicalizeすることができる。
SMILES列に欠損値がない場合に、SMILES列をmolオブジェクトに変換する場合には、あえてlambda式とする必要はないが、欠損値がある場合には、やはりlambda式と三項演算子の組み合わせが便利。
内容
1. 事例用のDataFrame作成
化合物６種類に対して、化合物番号（comID）、構造（SMILES）、特性値（今回は分子量；MW）の３列を持つpandasのDataFrameを作ります。構造不明な化合物がある場合も想定して、comID 6の構造は空白にしました。構造が未定でも、特性値がわかる場合もあるので、MW列には値を入れています。
import numpy as np
import pandas as pd
r1 = ['c1', 'c2', 'c3', 'c4', 'c5', 'c6']
r2 = ['C1CC(C)CCC1', 'C=CC(=O)OC1CCCC1', 'C=Cc1ccccc1', 'CC(=O)CC', 'C=Cc1ccncc1', '']
r3 = [98.19, 140.18, 104.15, 72.11, 105.14, 18.01]
df = pd.DataFrame(data={'comID':np.array(r1),
'SMILES':np.array(r2),
'MW':np.array(r3)})
df
| comID | SMILES | MW | |
|---|---|---|---|
| 0 | c1 | C1CC(C)CCC1 | 98.19 |
| 1 | c2 | C=CC(=O)OC1CCCC1 | 140.18 |
| 2 | c3 | C=Cc1ccccc1 | 104.15 |
| 3 | c4 | CC(=O)CC | 72.11 |
| 4 | c5 | C=Cc1ccncc1 | 105.14 |
| 5 | c6 | 18.01 |
2. やりたいこと
rdkit.Chem.rdmolfiles moduleを使い、SMILES列のSMILESを、全てcanonicalizeします。公式チュートリアルにもあるように、SMILESの文字列を一旦molオブジェクトとし、再度SMILESにすることにより、canocicalizeできます。
from rdkit import Chem
m=Chem.MolToSmiles(Chem.MolFromSmiles('C1CC(C)CCC1'))
m
[出力結果]: 'CC1CCCCC1'
これを、SMILES列の各要素に適用してやるのですが、データセットが大きくなった時などに備え、for文ではなく、applyで処理していきたいと思います。
3. 自作関数の利用
しかし、単純に、上記コードを使い、例えば下のようにするとerrorになります。PandasのSeriesにChem.MolFromSmilesを適用しようとしたためです。
#失敗事例
df['SMILES'] = df['SMILES'].apply(Chem.MolToSmiles(Chem.MolFromSmiles))
そこで、項目2. でのコードの内容を、自作関数にします。
def cansmi(m):
m=Chem.MolFromSmiles(m)
m=Chem.MolToSmiles(m)
return m
この関数cansmiを、dfのSMILES列に適用します。
df['SMILES'] = df['SMILES'].apply(cansmi)
df
これでSMILES列の内容を、canonical SMILESに書き換えることができます。
| comID | SMILES | MW | |
|---|---|---|---|
| 0 | c1 | CC1CCCCC1 | 98.19 |
| 1 | c2 | C=CC(=O)OC1CCCC1 | 140.18 |
| 2 | c3 | C=Cc1ccccc1 | 104.15 |
| 3 | c4 | CCC(C)=O | 72.11 |
| 4 | c5 | C=Cc1ccncc1 | 105.14 |
| 5 | c6 | 18.01 |
今回の例では、comID = c1, c4 のSMILESが書き変わっています。
4. lambda式の利用
自作関数を利用する方法は便利ですが、どうしてもコードの行数が増えてしまいます。同じコードの中で、自作した関数を何回か利用する機会があれば良いのですが、今回のような場合では、最初にSMILESをcanonicalizeすれば、その後同様の操作を行うケースは少ないと思います。
そこで、lambda式を使い、1行ですっきりまとめることとしました。結果として、6行分のコードを1行に収めることができました。
df['SMILES'] = df['SMILES'].apply(lambda m : Chem.MolToSmiles(Chem.MolFromSmiles(m)))
pythonの解説記事などで、lambda式が取り上げられていても、今ひとつ有り難み、使い途が判らなかったのですが、今回の事例で、便利さを実感できました。
5. csvファイルから読み込んでDataFrameとする場合
さて、データ解析の場面では、項目1のようにpythonのリストなどからDataFrameを作成するのではなく、別途csvファイル（や表計算ソフト）にデータをまとめておき、これをread_csv()で読み込んでDataFrameとする場合が多いと思います。この状況を模倣するため、項目1のDataFrame dfを一旦csvで書き出し、改めてDataFrame df2として読み込んでみます。
df.to_csv('test.csv', columns=['comID', 'SMILES', 'MW'], index=False)
df2 = pd.read_csv('test.csv')
df2
| comID | SMILES | MW | |
|---|---|---|---|
| 0 | c1 | C1CC(C)CCC1 | 98.19 |
| 1 | c2 | C=CC(=O)OC1CCCC1 | 140.18 |
| 2 | c3 | C=Cc1ccccc1 | 104.15 |
| 3 | c4 | CC(=O)CC | 72.11 |
| 4 | c5 | C=Cc1ccncc1 | 105.14 |
| 5 | c6 | NaN | 18.01 |
comID = 6のSMILES列が、空白からNaNに変わっています。read_csv()では、空白の要素が欠損値として取り扱われるためです。このDataFrame df2に、項目4のコードを適用すると、errorが返されます。欠損値に対して、Chem.MolFromSmilesを適用したためです。そこで、lambda式とif文を組み合わせることにします。
具体的には、
- DataFrameの要素が欠損値ならば（pd.isna()）、欠損値（np.nan）を返す。
- そうでなければ、canonical SMILESに変換する。
という動きにします。
df2['SMILES'] = df2['SMILES'].apply(lambda m : np.nan if pd.isna(m) else Chem.MolToSmiles(Chem.MolFromSmiles(m)))
もちろん、逆の動作、つまり
- DataFrameの要素が欠損値でなければ（pd.notna()）、canonical SMILESに変換する。
- そうでなければ、欠損値（np.nan）を返す。
df2['SMILES'] = df2['SMILES'].apply(lambda m : Chem.MolToSmiles(Chem.MolFromSmiles(m)) if pd.notna(m) else np.nan)
としても同じ動きをします。
| comID | SMILES | MW | |
|---|---|---|---|
| 0 | c1 | CC1CCCCC1 | 98.19 |
| 1 | c2 | C=CC(=O)OC1CCCC1 | 140.18 |
| 2 | c3 | C=Cc1ccccc1 | 104.15 |
| 3 | c4 | CCC(C)=O | 72.11 |
| 4 | c5 | C=Cc1ccncc1 | 105.14 |
| 5 | c6 | NaN | 18.01 |
このやり方は、欠損値を持つSMILES列から変換したMolオブジェクトの列を、DataFrameに追加する際にも使えます。
df2['ROMol'] = df2['SMILES'].apply(lambda m : Chem.MolFromSmiles(m) if pd.notna(m) else np.nan)
こうやって発生させたmolオブジェクトに対して、種々操作を加え、化学にまつわる解析を行なって行くことができます。
6. 化学の視点から検証
項目1のcomID = c6のSMILESに、空白（''）の代わりに、あえて手で'NaN'と入れてみました。
r2 = ['C1CC(C)CCC1', 'C=CC(=O)OC1CCCC1', 'C=Cc1ccccc1', 'CC(=O)CC', 'C=Cc1ccncc1', 'NaN']
Naはナトリウム、Nは窒素で、SMILESでは水素を省略するため、ナトリウムアミド
として誤認識することを懸念しました。
しかし、NaNはNaNのようで、やはり項目4の書き方ではerror、項目5の書き方では、欠損値として取り扱われました。
ナトリウムアミドをSMILES記法で書くと、[NH2-].[Na+]で、欠損値とナトリウムアミドを取り違えてしまう、ということは起こらないことも確認できました。
まとめ
lambda式の便利さを、実感として理解できました。特に条件分岐を持つlambda式は、python初心者にとってなかなか取っ付きにくい感じがありましたが、こうやって便利さを実感できると、利用していこう、という気持ちになります。
Discussion