Skip to content

Latest commit

 

History

History
336 lines (250 loc) · 19.8 KB

ch07_graph.asciidoc

File metadata and controls

336 lines (250 loc) · 19.8 KB
Raw

7章: グラフ構造を利用した類似性の評価

jupyter

グラフとはノード（頂点）群とノード間の連結関係を示すエッジ（枝）群で構成されるデータのことを指します。化学構造はこのグラフで表現できます。つまり原子をノード、結合をエッジとしたグラフ構造で表せます。

通常、6章で紹介したようなフィンガープリントを使い分子同士の類似性を評価することが多いですが、グラフ構造を利用して類似性を評価する手法もあります。次に紹介するMCS（Maximum Common Substructure）は対象となる分子集合の共通部分構造のことを指します。共通部分構造が多いほとそれらの分子はより似ていると考えます。

主要な骨格による分類(MCS)

最大共通部分構造Maximum Common Substructure(MCS)とは与えられた化学構造群において共通する最大の部分構造のことです。RDKitではMCS探索のためにrdFMCSというモジュールが用意されています。

今回はMCS探索のサンプルデータとしてrdkitに用意されているcdk2.sdfというファイルを利用します。RDConfig.RDDocsDirが、サンプルデータのディレクトリを表す変数で、そのディレクトリ以下のBooks/data/にcdk2.sdfというファイルが存在するので、os.path.joinメソッドでファイルパスを設定します。尚、os.path.joinはosのパスの違いを吸収するためのpythonの組み込みモジュールです。

import os
from rdkit import Chem
from rdkit.Chem import RDConfig
from rdkit.Chem import rdFMCS
from rdkit.Chem.Draw import IPythonConsole
from rdkit.Chem import Draw
filepath = os.path.join(RDConfig.RDDocsDir, 'Book', 'data', 'cdk2.sdf')
mols = [mol for mol in Chem.SDMolSupplier(filepath)]
# 構造を確認します
Draw.MolsToGridImage(mols[:7], molsPerRow=5)
compounds

読み込んだ分子を使ってMCSを取得します。RDKitではMCSの取得方法に複数のオプションが指定できます。 以下にそれぞれのオプションでの例を示します。

  1. デフォルト

  2. 原子がなんであっても良い（構造とボンドの次数があっていれば良い）

  3. 結合次数がなんでも良い（例えば、ベンゼンとシクロヘキサンは同じMCSとなる）

result1 = rdFMCS.FindMCS(mols[:7])
mcs1 = Chem.MolFromSmarts(result1.smartsString)
mcs1
print(result1.smartsString)
#[#6]1:[#7]:[#6](:[#7]:[#6]2:[#6]:1:[#7]:[#6]:[#7]:2)-[#7]
MCS01 
result2 = rdFMCS.FindMCS(mols[:7], atomCompare=rdFMCS.AtomCompare.CompareAny)
mcs2 = Chem.MolFromSmarts(result2.smartsString)
mcs2
print(result2.smartsString)
#[#6]-,:[#6]-,:[#6]-[#6]-[#8,#7]-[#6]1:[#7]:[#6](:[#7]:[#6]2:[#6]:1:[#7]:[#6]:[#7]:2)-[#7]
MCS02 
result3 = rdFMCS.FindMCS(mols[:7], bondCompare=rdFMCS.BondCompare.CompareAny)
mcs3 = Chem.MolFromSmarts(result3.smartsString)
mcs3
print(result3.smartsString)
#[#6]1:[#7]:[#6](:[#7]:[#6]2:[#6]:1:[#7]:[#6]:[#7]:2)-[#7]
MCS03

RDKitではMCSに基づく類似性を数値化するアルゴリズムのひとつにFraggle Similarityが実装されています。これを利用することでクラスタリングや、類似性に基づいた解析が行なえます。

from rdkit.Chem.Fraggle import FraggleSim
sim, match = FraggleSim.GetFraggleSimilarity(mols[0], mols[1])
print(sim, match)
#0.925764192139738 *C(C)C.*COc1nc(N)nc2[nH]cnc12
match_st = Chem.MolFromSmiles(match)
match_st
FraggleSimilarity

このようにFraggleSimilarityは類似性及びマッチした部分構造を返します。ECFPを利用した類似性よりもケミストの感覚に近いことが多いです。詳しくは参考リンクを参照してください。

参考リンク

Matched Molecular PairとMatched Molecular Series

jupyter

創薬研究の構造最適化ステージにおいて、起点となる化合物（リード化合物）をどのように構造変換していくかは非常に重要ですが、ステージが進んだ場合どの構造変換が活性や物性に影響を及ぼしたかというレトロスペクティブな解析することもまた非常に大切です。

Tip
 興味があればhttps://sar.pharm.or.jp/wp-content/uploads/2018/09/SARNews_19.pdfを読むとよいです。

Matched Molecular Pair(MMP)は、二つの分子のうち一部の部分構造だけが異なりそれ以外は同一な分子のペアのことです。例として、クロロベンゼンとフルオロベンゼンはCl基とF基のみが異なるのでMMPです。このようなペアの特性の変化を大量に解析することで、置換基変換のトレンドを掴むことができます。これをMatched Molecular Pair Analyisis（MMPA)と呼びます。大規模なデータでMMPAを行うことにより、置換基の変化がもたらす特性変化の普遍的なルールを抽出できます。このようなルールを理解していれば構造最適化を効率的に進められることになります。

ここではRDKitのContribに提供されているRDKit/Contrib/MMPA[mmpa]を使ってMMP解析を行います。

RDKitインストール先の下にあるContrib/mmpaに移動し、pythonスクリプトを順次実行します。

python rfrag.py <MMPAを実施したいFileの名前 >フラグメント化したデータの保存ファイル名
# 例えば
# python rfrag.py <data/sample.smi >data/sample_fragmented.txt

python indexing.py <先のコマンドでできたフラグメントのファイル >MMP_アウトプットファイル.CSV
# 例えば
# python indexing.py <data/sample_fragmented.txt >data/mmp.csv

以上のコマンドを実行すると分子A, 分子B, 分子AのID, 分子BのID, 変換された構造のSMIRKS, 共通部分構造（context）のcsvファイルが生成されます。このデータに基づき活性や物性などを紐つけることでMMPAが行えます。

Note
 SMIRKSは分子の変換をSMILESのように文字列表記によって表現する手法です。

MMPの拡張としてMatched Molecular Series(MMS)という手法も提案されています。MMPは分子のペアですが、このペアを共通構造をもった3つ以上の集団としてリスト化したものがMMSです。

実際にMMSを作ってみましょう。以下の例ではChEMBLよりFactor Xaのデータをダウンロードし、一例として使いました。MMSの実装に関してはNoel O’Boyle氏のRDKit UGMでのプレゼンテーションのコードを利用しています。

まず利用するライブラリの読み込みと、データの読み込みを行い、SaltRemoverを使い脱塩しています。

import sys
import os
import pandas as pd
from rdkit import Chem
from rdkit.Chem import rdMMPA
from rdkit.Chem import RDConfig
from rdkit.Chem import rdBase
from rdkit.Chem.Draw import IPythonConsole
from rdkit.Chem import Draw
from rdkit.Chem import SaltRemover

# RDKit contrib/mmpaをPythonのSystem Pathに追加して次のパートで分子をフラグメント化するためのrfragをインポートできるようにします
mmpapath = os.path.join(RDConfig.RDContribDir, 'mmpa')
sys.path.append(mmpapath)

# データをPandasのDataFrameとして読み込み、SmilesをMolオブジェクトに変換、CHEMBLのIDをプロパティとして登録します。
df = pd.read_csv('Chembl_FXa.txt', sep='\t')
remover = SaltRemover.SaltRemover()
mols = []
for smi, chembl_id in zip(df.CANONICAL_SMILES, df.CMPD_CHEMBLID):
    try:
        mol = Chem.MolFromSmiles(smi)
        mol.SetProp('CMPD_CHEMBLID', chembl_id)
        mol = remover.StripMol(mol)
        mols.append(mol)
    except:
        print(smi)

続いてRDKit contribに登録されているmmpaのrfragをインポートして、分子をフラグメントに分割します。 以下のコードで用いているrfrag.fragment_molで呼ばれている関数の実態はrdMMPA.FragmentMolです。この関数は、与えたMolオブジェクトを回転可能ボンドにて切断したフラグメントを生成します。 出力の形式はfragment(SMILES, COMPOUND_ID)とした場合、SMILES, COMPOUND_ID, CORE, CHAINS というカンマ区切りのテキストのリストになります。COREというのは例えばA-B-Cという分子（ーが回転可能ボンドを意味します）を考えると、切断の仕方がA/B-C, A/B/C, A-B/Cという三つのケースがあります。二番目のケースの時にBがCoreにA, CがCHAINSになります。一番目と、三番目はどちらも (A,B-C), (A-B,C) というフラグメントはCHAINSにアサインされ、COREはNONEになります。

import rfrag
remover = SaltRemover.SaltRemover()
rfragdata = []
for smi, chembl_id in zip(df.CANONICAL_SMILES, df.CMPD_CHEMBLID):
    try:
        mol = Chem.MolFromSmiles(smi)
        mol = remover.StripMol(mol)
        smi = Chem.MolToSmiles(mol)
        out = rfrag.fragment_mol(smi, chembl_id)
        rfragdata.append(out)
    except:
        print(smi, chembl_id)

MMSを作成する関数を定義します。コードはUGMの資料に記載されているものをほぼそのまま利用しますが、Jupyter上で全ての処理をおこなうため読み込み先をファイルからリストに変更しました。

以下MMSの作成プロセスの概要です。

  1. 各分子を一定のルール（回転可能結合で切断など）でカット

  2. カットしたフラグメントがkeyの辞書を作成、同じキーを持つ分子のフラグメントを辞書のvalueに格納

上記の作業を繰り返すとことで共通のスキャフォールドを持つ分子をまとめられます。共通のスキャフォールドでまとまった分子は、スキャフォールド以外の置換基が異なる分子となります。

スキャフォールドとは？

創薬において、前臨床試験の前のステージに構造最適化というステージがあり、そこでは化合物の主要骨格以外の部分をちょこまかと変換して薬にふさわしいバランスの取れたプロパティにします。

この主要骨格のことをスキャフォールドと呼びます。例えばこの特許ではRを除いた部分は固定されておりこの主要骨格をスキャフォールドと呼びます。

scaffold

以下のコードの簡単な説明をします。

Fragというnamedtupleはid, scaffold, rgroupを属性に持ちます。Namedtupleを利用することで簡単なクラスと同じ働きをするオブジェクトを作ることができます。

getFrags関数は、フラグメント化してあるデータを元にFragクラスにデータを登録します。登録の際に、フラグメントの原子数を比較しどちらをScaffoldにするかを決めています。この関数を呼ぶことで各分子のid, scaffold, rgroupのペアを作ります。

次のgetSeries関数はFragオブジェクト内に格納されたscaffoldと与えられたデータを照合し、同じscaffoldを持つrgroupのデータをまとめます。これによってペアからシリーズにしています。Seriesクラスはscaffoldが文字列（SMILES)、rgroupsはリストで複数のフラグメントを持てるようになっています。

from collections import namedtuple

Frag = namedtuple( 'Frag', ['id', 'scaffold', 'rgroup'] )

class Series():
    def __init__( self ):
        self.rgroups = []
        self.scaffold = ""

def getFrags(rfrags):
    frags = []
    for lines in rfrags:
        for line in lines:
            broken = line.rstrip().split(",")
            if broken[2]: # single cut
                continue
            smiles = broken[-1].split(".")
            mols = [Chem.MolFromSmiles(smi) for smi in smiles]
            numAtoms = [mol.GetNumAtoms() for mol in mols]
            if len(numAtoms) < 2:
                continue
            if numAtoms[0] > 5 and numAtoms[1] < 12:
                frags.append(Frag(broken[1], smiles[0], smiles[1]))
            if numAtoms[1] > 5 and numAtoms[0] < 12:
                frags.append(Frag(broken[1], smiles[1], smiles[0]))
    frags.sort(key=lambda x:(x.scaffold, x.rgroup))
    return frags

def getSeries(frags):
    oldfrag = Frag(None, None, None)
    series = Series()
    for frag in frags:
        if frag.scaffold != oldfrag.scaffold:
            if len(series.rgroups) >= 2:
                series.scaffold = oldfrag.scaffold
                yield series
            series = Series()
        series.rgroups.append((frag.rgroup, frag.id))
        oldfrag = frag
    if len(series.rgroups) >= 2:
        series.scaffold = oldfrag.scaffold
        yield series

MMSを作る準備ができたので実行します。同じスキャフォールドに対して４つ以上置換基の変換があったデータのみを可視化します。

frags = getFrags(rfragdata)
series = getSeries(frags)
series =[i for i in series]
from IPython.display import display
for s in series[:50]:
    mols = [Chem.MolFromSmiles(s.scaffold)]
    ids = ['scaffold']
    for r in s.rgroups:
        rg = Chem.MolFromSmiles(r[0])
        mols.append(rg)
        ids.append(r[1])
    if len(mols) > 5:
        display(Draw.MolsToGridImage(mols, molsPerRow=5, legends=ids))
        print("########")
MMS

スキャフォールドに対して5つの置換基のMMSが表示されました。

Note
 このMMSを利用して活性予測を行うこともできます。

Cytoscapeを使ってMMPネットワークを可視化する

Warning
 この内容は入門の内容を超えるので興味がなければ飛ばしてください

MMPは変換前、変換後の情報をノード、変換ルールをエッジとするグラフ構造と考えることができます。Cytoscapeなどのネットワーク可視化ツールを利用するとこのグラフ構造を直感的に把握できます。

RDKitには先に紹介したMMPAの他にmmpdbという別プロジェクトがあります。 こちらはコマンドラインのツール群とデータベースシステムとして提供されているため、長期的な管理がしやすいという特徴があります。本セクションではこのmmpdbとCytoscapeを利用したMMPの可視化を紹介します。

Cytoscapeのインストール

Cytoscapeはオープンソースのネットワーク可視化ソフトで色々なシーンで広く使われています。化合物の構造表示用プラグインを使うことで構造のネットワークを表示することができます。

インストールは簡単でダウンロードサイトから対応するOSのインストーラをダウンロードして指示のとおりにインストールするだけです。

インストールが完了したらCytoscapeを起動して化合物構造描画用のChemviz2プラグインをインストールします。手順は簡単でApps→App Managerからchemviz2を選択してインストールします。

AppManager

mmpdbからgmlファイルを作成する

今回利用するデータは<Inhibition of recombinant GSK3-beta> J. Med. Chem. (2008) 51:2062-2077の151化合物です。MMPAを行うにはHTSのような探索データではなくて構造最適化のようにスキャフォールドが決まっているものを使うのが原則です。

コマンドの流れを載せておきます。SMILESのtextと活性や物性値のデータは別々にデータベースに登録する必要があります。

$ mmpdb fragment smiles.txt -o CHEMBL930273.fragments     # fragmentation
$ mmpdb index CHEMBL930273.fragments -o CHEMBL930273.db   # make db
$ mmpdb loadprops -p act.txt CHEMBL930273.db              # load properties

そのあとCytoscapeで読み込むためのgmlファイルを作成しますが、これは本書の範囲を超えるので割愛します。もし興味があるのであればコードを直接読んでもらうといいのですが 流れは以下のとおりです。

  1. mmpdbからpython-igraphを使ってgmlファイルを作る

  2. gmlファイルをCytoscapeで読み込む

  3. Cytoscapeで属性を各パラメータにアサインして視覚的に理解しやすくする

    1. ノードの大きさを物性値に対応

    2. エッジの色を活性差に対応

    3. chemviz2 pluginで構造を描画してノードに貼り付ける

解釈する

MMPネットワークを見てみましょう。あまり活性差のないMMPが左上の方に固まっています。右下の方にはエッジが赤い（活性差が大きい）ものが観測されます。このような小さな置換基変化が大きな活性差を生むもMMPをActivity Cliffと呼びます。一般的にActivity Cliffは創薬プロジェクトにおいてブレークスルーとなることが多いため、このような活性変化を見逃さないことは大切です。

MMPN MMPN

実際にどういう置換が行われたのかを確認すると、OH基がMeO基に置換されることで活性の消失が起こっています。

MMPだけではこのように単純に事実しかわからないので、もう少し深く考察するために類似体の複合体結晶構造を探してみました。するとPDBID:5OY4というGSK3βと類似化合物の複合体が見つかりました。

MMPN MMPN

OH基をMeO基に置換するとポケットの壁にぶつかりそうですね。つまりこのActivity Cliffはリガンドと蛋白質の立体障害 により引き起こされたと考えられます。