RDKit:canonical SMILESを比較する

背景

化学者のDr_Oです。有機合成化学と呼ばれる分野で仕事しています。
この分野でも、色々な方が、今まであまり使わなかったデータ解析の方法を取り入れたりしておられます。自分も勉強しながらですが、化学に登場しがちな場面で気づいたこと、経験したことを共有できたらと思い、記事を書き始めました。

はじめに

ということで記念すべき（？）第一回です。化学構造式を表現するSMILES記法や、それを用いた化合物の選別について気付いた情報を共有します。

概要

確認したこと

「canonical SMILES」と書いてあれば、常に同一化合物は、同一の文字列で表わされているか。

結果

canonical SMILESと書いてあっても、「何（ナニ）で」正規化処理（canonicalize）したかによって、文字（原子と結合）の順番が異なることがある。SMILESを基準にして、化合物の重複をチェックする場合などでは、注意が必要。

内容

1. SMILESってなんや

化学の人は、「構造式」で化合物を認識してます。こんなやつ。

しかし ! データ解析の場面では、この絵（構造式）ではコンピュータが理解してくれないので、構造式を文字列にして入力します。ざっくり、化合物中の（水素以外の）原子の並びを、元素記号で並べて表現してる感じです。ここら辺は、良い記事がたくさんあるので、「SMILES記法」などで検索して下さい 😅
例えば上の化合物（スチレンね）は、「C=CC1=CC=CC=C1」となったりします。

2. ではcanonical SMILESとは

単に原子の並びを文字（元素記号） の並びで表現する、だと、同じ化合物を、複数の異なる文字列で表現でできてしまいます。例えば、上のスチレンだと、

• C=CC1=CC=CC=C1
• C1=CC=CC=C1C=C
• C1=C(C=C)C=CC=C1

など、いずれでも書けます。SMILES的にはOKなんだけど、これでは、文字列を見ても同じ化合物かどうか判りにくいですよね。
ということで、ルールに従って並び順を決め（正規化 = canonicalize）、それに従って記述しましょう、というのがcanonical SMILESです。
つまり、化合物の分子を、一義的に決まる文字列に対応させて表現しましょう、ということです。

3. やっと本題

元となるSMILESをどうやって準備するか、というのは、「化学の話」なのでちょっと置いておいて、SMILESをcanonicalizeする方法（ツールといっても良いかも）は、いくつかあります。
現時点、化学系でのdata drivenな取り組み（ざっくりケモインフォとかマテリアルズインフォマティクス (MI) のイメージ）では、pythonベースで物事を取り扱う方が多いと思います。とすると、rdkit、openbabel（またはpybel）あたりが選択肢になってくるでしょう。また、PubChemで化合物検索をして、Canonical SMILESの項目から持ってくることもできます。で、canonicalなんだから、どのやり方でも答えは同じなの？というのが今回の記事の趣旨です。

4. やってみた

PubChemから無作為に選んだ、CID 10024の化合物のCanonical SMILESを、文字列 s1 としました。
こんな構造です。

でかいけど対称性はいいな。

#CID 100224  
s1='C1=CC=C2C(=C1)C(=C3C=CC(=CC3=N2)Br)NCCCCCCCCNC4=C5C=CC(=CC5=NC6=CC=CC=C64)Br'

pybel（検証環境ではopen babel 2.4.1訂正 3.1.1）で、s1をSMILESとして読み込み、モルオブジェクト m としました。その後、m をcanonical SMILESに変換して、s2として出力しました。なお、なぜかタブ区切りと、改行の不可視コードが含まれていたので、これは除きました。

from openbabel import pybel
mp=pybel.readstring("smi", s1)
s2 = mp.write('can').replace('\t\n', '')

次に、同じくs1を、rdkitでモルオブジェクトにした後、そのままSMILESに変換し、s3 として出力しました。rdkitでは、これでcanonical SMILESとして出力されます。

from rdkit import Chem
s3 = Chem.MolToSmiles(Chem.MolFromSmiles(s1))

s1、s2、s3をデータフレームに並べて、比較しやすくしました。

import pandas as pd
df = pd.DataFrame(pd.Series([s1, s2, s3],index=['PubChem', 'pybel', 'rdkit']), 
                  columns=['canonicalSMILES'])
df

・・・いや、３つとも、バラバラやん。
芳香環がケクレ構造か否か、とか以前に、文字の並び方が、三者三様です。準拠しているルールが、それぞれ異なるようです。もちろん、それぞれのやり方の中では、繰り返し試しても、再現性よく同じ文字列が発生します。ルールの出典は、それぞれの公式マニュアル等に記載があります。

まとめ

化合物群を含む二つのデータセットを比較して、重複する化合物だけ抜き出す、とか、重複を除く、という操作は、化学のデータ解析の場面で、よく登場します。その際、SMILESをキーとして比較、抽出することも、多いです。
ところが、二つの化合物群が、それぞれ別々の方法で canonicalize されていたら、正確な比較ができません。また、例えば同じrdkitでcanonical SMILESを発生させていたとしても、片方はデフォルトオプションで、もう一方は別のオプションを指定して canonicalize していたとしたら、やはり同一化合物が、異なるSMILESで表現されてしまう、ということが起こります。
従って、化合物のSMILESを含むデータセットを取り扱う時は、自分の手元で一旦 canonicalize してから、解析を始める習慣にしています。
色々工夫して、解析が散々進んでから、「あーこれとこれは一緒だった😵」とか、「抜けてた〜」、となるのは、ダメージ大きいですもんね。