SSSRを使って環状化合物と鎖状化合物を区別しようとする話

以前の記事で、創薬レイドバトルの共闘用シェアデータを２つにクラスタリングした後、環状ペプチドとそれ以外の分子の集まりとに手作業で分類しました。これでは数増えた時に大変ですし、見落としもありそうです。そこで今回はRDKit オンラインドキュメンテーションに記載のあったSSSRという考え方を頼りに同様の分類できるかどうか試してみたいと思います。

RDKitオンラインドキュメンテーションの該当箇所は下記です。
1. Ring Information
2. Ring Finding and SSSR

SSSR（smallest set of smallest rings; SSSR）については @mojaie さんのこちらの記事 化学構造情報とグラフ理論の説明が非常にわかりやいのでご参照ください。

とにかく、分子に含まれる環構造のうち「環の個数と環の大きさが最小になるような組み合わせ」だそうです。（ふわっとした理解）

環の情報がわかるなら、大員環とそれ以外の分子と区別することもできるのでは？？？　というのが今回の記事の目的です。

まずはデータの準備

from rdkit import rdBase, Chem
from rdkit.Chem import AllChem, Draw, Descriptors, PandasTools
import pandas as pd
 

df = pd.read_csv('../pd1_inhibitor_dataset-master/PD1_inhibitor_dataset.csv')
PandasTools.AddMoleculeColumnToFrame(df, 'smiles')

GetSSSR関数

RDKitのドキュメンテーションによると、SSSRに関してはGetSSSRとGetSymmSSSRの２つの関数があるそうです。

まずはGetSSSR関数を用いて見ます。こちらはSSSRの数を取得することができる関数です。

データセットの最初の分子をテスト化合物としてどのような結果が出るのか眺めて見ます。

最初の分子の構造はこのようなものです。

test_mol = df.loc[0, 'ROMol']
Draw.MolToImage(test_mol)

この分子のSSSRの数は？？？

test_SSSR = Chem.GetSSSR(test_mol)
print(test_SSSR)
# 8

8個だそうです！

・・・わからん。

SSSRの数だけではよくわかりません。

GetSymmSSSR関数

もう一方の関数、GetSymmSSSRを使うと、認識された各環構造の情報を取得することもできるそうです。

test_SymmSSSR = Chem.GetSymmSSSR(test_mol)
print(type(test_SymmSSSR))
# <class 'rdkit.rdBase._vectNSt3__16vectorIiNS_9allocatorIiEEEE'>

print(type(test_SSSR))
# <class 'int'>

GetSSSRで得た結果は整数値（int）ですが、GetSymmSSSRで得た結果はrdkitのオブジェクトでSSSRについてより詳細な情報をもっていそうです。

GetSymmSSSRの結果から、SSSRの数を取得するには長さを求めれば良いそうです。

print(len(test_SymmSSSR))
#8

先ほどの結果と一致しました。

最大の環構造を探す

ドキュメンテーションの説明をみる限り、GetSymmSSSRで得られるものは、認識された環構造のリストとなっており、リストの要素は各環構造に含まれる原子（atom index）のリストとなっているようなオブジェクトのようです（リストのリスト）。

8個の各環構造のサイズを求めてみます。

for i in range(len(test_SymmSSSR)):
    print('{}番目の環のサイズ:{}'.format(i, len(test_SymmSSSR[i])))

"""
出力結果
0番目の環のサイズ:45
1番目の環のサイズ:5
2番目の環のサイズ:5
3番目の環のサイズ:6
4番目の環のサイズ:5
5番目の環のサイズ:6
6番目の環のサイズ:5
7番目の環のサイズ:6
"""

0番目が環員数 45となっています。ぼんやりと眺めていましたが環状ペプチドの長さはこんなにも長いものなんですね。具体的な数値を見るとあらためてびっくりしました。

この最大の環構造を構成する原子の番号（atom index）を取得して見ます。

largest_ring_member = list(test_SymmSSSR[0])
print(largest_ring_member)

# [4, 5, 7, 8, 9, 11, 12, 13, 14, 15, 16, 18, 19, 20, 22, 23, 24, 26, 27, 28, 30, 34, 35, 37, 38, 39, 41, 42, 43, 45, 46, 47, 49, 50, 51, 53, 54, 55, 57, 58, 59, 61, 62, 63, 65]

構成要素の原子がわかったので、この情報を用いて分子のどの原子が該当の環構造に含まれるのか、ハイライトして描画して見ました。

期待通り、環状ペプチドの骨格となる構造を認識できていそうです。

環構造を順番に眺める

ipywidgetsを利用して各環構造を順番に眺めて見ます。

# 描画用
from rdkit.Chem.Draw import rdMolDraw2D
from IPython.display import SVG
from matplotlib.colors import ColorConverter

# インタラクティブなビューワー
from ipywidgets import interact,fixed,IntSlider
import ipywidgets

# ビューワーに与える関数を作成
def SSSR_viwer(idx):
    ring_member_atoms = list(test_SymmSSSR[idx])
    prepro_test_mol = rdMolDraw2D.PrepareMolForDrawing(test_mol)
    target_ring_atoms = list(test_SymmSSSR[idx])
    
    #デフォルトの設定は見づらいので色を変える
    color_dict = {}
    radius_dict = {}
    for i in target_ring_atoms:
        color_dict[i] = ColorConverter().to_rgb('khaki')
        radius_dict[i] = 0.7
    
    #コンテナとなるオブジェクトの作成
    view = rdMolDraw2D.MolDraw2DSVG(600,600)
    
    #コンテナの描画設定
    option = view.drawOptions()
    option.circleAtoms=False
    option.continuousHighlights=False
    
    #最大の環構造に含まれる原子をハイライトに設定
    view.DrawMolecule(prepro_test_mol, 
                      highlightAtoms=target_ring_atoms,
                      highlightAtomColors= color_dict,
                      highlightAtomRadii= radius_dict,
                      highlightBonds=[])

    #コンテナをファイナライズ
    view.FinishDrawing()

    #コンテナに書き込んだデータを取り出す
    svg = view.GetDrawingText()

    #データを描画
    return(SVG(svg.replace('svg:,','')))

#ビューワーを実行
interact(SSSR_viwer, idx=ipywidgets.IntSlider(min=0,max=len(test_SymmSSSR)-1, step=1));

順番に表示させることでSSSRがどのような環構造を認識するのか、非常にわかりやすくなりました。

SSSRを利用して大員環を区別

大体SSSRでどのような情報が得られるのかわかったので、これを用いてマクロサイクル型の化合物とそれ以外を区別できるか検証して見たいと思います。

Wikipediaによると大員環化合物（macrocyclic compound)は「概ね10個以上の原子からなる環状構造を持つ有機化合物」ということなので、これを基準として用います。SSSRで得た環構造のうち、最大の大きさの環構造を構成する原子数が10以上の場合、大員環化合物とすることとしたいと思います。

まずは、各化合物のSSSRをもとめ、最大の環構造の情報を取得します。

関数を作成し、Molオブジェクトを含むカラムに適用し、得た値を新しいカラムへと格納します。

def largest_ring_member_count(mol):
    SymmSSSR = Chem.GetSymmSSSR(mol)

    size_list = []
    
    # 各環の大きさを順番にリストに追加
    for i in range(len(SymmSSSR)):
        size_list.append(len(SymmSSSR[i]))
    
    # リストの最大値をとりだす
    max_size = max(size_list)
                   
    return(max_size)

# 最大の環の情報をもつ新しいカラムを作成
df['largest_ring_size'] = df['ROMol'].map(largest_ring_member_count)

環のサイズが10以上の時はマクロサイクル型（macro）、それより小さい場合は鎖型（chain）として分類してみます。

#classifierとか一回言ってみたかった
def macro_classifier(size):
    if size >= 10:
        return('macro')
    elif size < 10:
        return('chain')

df['macro_chain'] = df['largest_ring_size'].map(macro_classifier)

どのような分類になったか数を確認してみます。

df['macro_chain'].value_counts()

"""
chain    39
macro    14
Name: macro_chain, dtype: int64
"""

マクロサイクル型 14個、鎖状型39個、となりました。

前回の分類では、外れ値としたマクロサイクル型化合物１つの除いた、合計52個を、マクロサイクル型13個と鎖状型39個とに分けました。数を見る限り同様の結果が得られていそうです。

それぞれ別のDataFrameに分けたあと、構造を確認します。

マクロサイクル型・・・

df_macro = df[df['macro_chain']== 'macro']
macro_list = list(df_macro['ROMol'])

def macro_viewer(idx):
    mol = macro_list[idx]
    return(display(Draw.MolToImage(mol)))
    
interact(class_1_viewer, idx=ipywidgets.IntSlider(min=0,max=len(macro_list)-1, step=1));

鎖状型・・・

df_chain = df[df['macro_chain']== 'chain']
chain_list = list(df_chain['ROMol'])

def chain_viewer(idx):
    mol = chain_list[idx]
    return(display(Draw.MolToImage(mol)))
    
interact(chain_viewer, idx=ipywidgets.IntSlider(min=0,max=len(chain_list)-1, step=1));

しっかりと望み通りの分類ができていそうです。

まとめ

今回はSSSRを使って、大員環化合物とそれ以外の化合物を分類できるか？ということを行ってみました。
RDKitのドキュメンテーションを読んでいる時は、どのような目的で使えば良いものかさっぱりわからず、どうして The SSSR Problem という節を設けるほど解説してあるのか疑問に思っていたのですが、実際に使って見るととても使い勝手の良い重要な考え方ということがわかりました。

かなり適当に使ってしまったので、概念の理解の誤り等ございましたらご指摘いただければ幸いです。