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
canonicalSMILES | |
---|---|
PubChem | C1=CC=C2C(=C1)C(=C3C=CC(=CC3=N2)Br)NCCCCCCCCNC4=C5C=CC(=CC5=NC6=CC=CC=C64)Br |
pybel | Brc1ccc2c(c1)nc1c(c2NCCCCCCCCNc2c3ccccc3nc3c2ccc(c3)Br)cccc1 |
rdkit | Brc1ccc2c(NCCCCCCCCNc3c4ccccc4nc4cc(Br)ccc34)c3ccccc3nc2c1 |
・・・いや、3つとも、バラバラやん。
芳香環がケクレ構造か否か、とか以前に、文字の並び方が、三者三様です。準拠しているルールが、それぞれ異なるようです。もちろん、それぞれのやり方の中では、繰り返し試しても、再現性よく同じ文字列が発生します。ルールの出典は、それぞれの公式マニュアル等に記載があります。
まとめ
化合物群を含む二つのデータセットを比較して、重複する化合物だけ抜き出す、とか、重複を除く、という操作は、化学のデータ解析の場面で、よく登場します。その際、SMILESをキーとして比較、抽出することも、多いです。
ところが、二つの化合物群が、それぞれ別々の方法で canonicalize されていたら、正確な比較ができません。また、例えば同じrdkitでcanonical SMILESを発生させていたとしても、片方はデフォルトオプションで、もう一方は別のオプションを指定して canonicalize していたとしたら、やはり同一化合物が、異なるSMILESで表現されてしまう、ということが起こります。
従って、化合物のSMILESを含むデータセットを取り扱う時は、自分の手元で一旦 canonicalize してから、解析を始める習慣にしています。
色々工夫して、解析が散々進んでから、「あーこれとこれは一緒だった😵」とか、「抜けてた〜」、となるのは、ダメージ大きいですもんね。
Discussion