🐼

RDKit:PandasのDataFrameのSMILES列にapplyを適用してcanonicalizeする

2022/05/11に公開

背景

化学者のDr_Oです。有機合成化学と呼ばれる分野で仕事しています。
pythonやRを使ったデータ解析について、自分も勉強しながらですが、化学に登場しがちな場面で気づいたこと、経験したことを共有できたらと思い、記事を書き始めました。

はじめに

前回に引き続き、分子の構造を表すSMILESに関する記事です。化合物群をPandasのDataFrameで管理する際、SMILESの列を一気にcanonicalizeする場合について、記載します。
1つのSMILESをrdkitでcanonicalizeすることや、Pandasのデータフレームの特定の列に関数を適用して一気に変換することは、それぞれ色々なところで解説されていますが、これを組み合わせた場合の記述はあまりないように思ったので、記事にしました。

概要

lambda式と三項演算子を組み合わせることにより、より少ない行数で、欠損値を持つpandasのDataFrameのSMILES列に対し、rdkitの関数をapplyを用いて適用し、canonicalizeすることができる。
SMILES列に欠損値がない場合に、SMILES列をmolオブジェクトに変換する場合には、あえてlambda式とする必要はないが、欠損値がある場合には、やはりlambda式と三項演算子の組み合わせが便利。

内容

1. 事例用のDataFrame作成

化合物6種類に対して、化合物番号(comID)、構造(SMILES)、特性値(今回は分子量;MW)の3列を持つ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文を組み合わせることにします。
具体的には、

  1. DataFrameの要素が欠損値ならば(pd.isna())、欠損値(np.nan)を返す。
  2. そうでなければ、canonical SMILESに変換する。
    という動きにします。
df2['SMILES'] = df2['SMILES'].apply(lambda m : np.nan if pd.isna(m) else Chem.MolToSmiles(Chem.MolFromSmiles(m)))

もちろん、逆の動作、つまり

  1. DataFrameの要素が欠損値でなければ(pd.notna())、canonical SMILESに変換する。
  2. そうでなければ、欠損値(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では水素を省略するため、ナトリウムアミド

Na^+NH_2^-

として誤認識することを懸念しました。
しかし、NaNはNaNのようで、やはり項目4の書き方ではerror、項目5の書き方では、欠損値として取り扱われました。
ナトリウムアミドをSMILES記法で書くと、[NH2-].[Na+]で、欠損値とナトリウムアミドを取り違えてしまう、ということは起こらないことも確認できました。

まとめ

lambda式の便利さを、実感として理解できました。特に条件分岐を持つlambda式は、python初心者にとってなかなか取っ付きにくい感じがありましたが、こうやって便利さを実感できると、利用していこう、という気持ちになります。

Discussion

ログインするとコメントできます