グラフ理論
chemoinformatics

化学構造情報とグラフ理論

創薬 Advent Calendar 2018の14日目の記事です。

化学構造情報を扱うためのグラフ理論の概念やアルゴリズムについての備忘録です。

ちょっと類似化合物検索したことあるよ、くらいのケモインフォマティクスの知識がある読者を想定しています。


分子グラフ(molecular graph)モデル

原子をノード(node)、結合をエッジ(edge)として化合物をグラフで表現したモデルがよく使われます。水素は通常省略されます。

molgraph.png

ノードには原子の種類、電荷、多重度、質量など、エッジには結合の次数などの情報(ラベル)が格納されます。それぞれ芳香族性や立体異性に関する情報を持つこともあります。結合の次数についてはエッジではなくノードにπ電子数の形式で持たせた方が、実際の原子軌道と立体構造をより反映した表現になると思います。

分子グラフは一般的にエッジに向きが無い無向グラフ(undirected graph)で表現されます。ちなみにエッジに向きがある(一方通行の経路が存在する)グラフは有向グラフ(directed graph)と言います。

また、分子グラフは一般的に単純グラフ(simple graph)です。単純グラフは自己ループ(同一ノードをつなぐエッジ, self-loop)を持たず、また多重エッジ(2つのノードの間に2本以上あるエッジ, multi-edge)を持ちません。

あくまでこれは一つのモデルの例であり、分子グラフの表現方法は様々です。特にここ数十年で、三中心二電子結合など、単純グラフでは表せない非古典的な結合の存在が明らかになってきたので、分子をハイパーグラフ(エッジの概念を3つ以上のノードの集合に拡張することでグラフを一般化したもの)で表すアプローチが生まれています[1]


有機化合物の分子グラフの特徴


  • ノードの次数(degree)はおおよそ1-4

  • ほぼ全てが平面グラフ(planar graph)

  • その中でも外平面グラフ(outerplanar graph)が多い

次数(degree)とはノードに隣接するエッジの本数のことです。有機化合物の場合、一つの原子に結合する原子の数が4を超えることは稀です。遷移金属の配位や超原子価分子などで4を超えることもあります。

平面グラフ(planar graph)とは、平面上にグラフを配置した際にエッジが交差しないノードの配置が少なくとも一つ存在するグラフです。外平面グラフ(outerplanar graph)とは、平面グラフの中でも特に全てのノードがグラフの外縁に位置するようなグラフです。ちなみに、テトラへドランやフラーレンは立体的に見えますが、平面グラフです(Wikipediaのフラーレンの項目参照)

planar.png

分子グラフが比較的次数が低いグラフ(疎なグラフ、sparse graph)であることは重要です。疎グラフは行列(隣接行列)よりもマッピング(Pythonで言うところのDict)で実装した方が効率的です。また、汎用グラフアルゴリズムの中には特に疎グラフにおいて効率的に働くものがあります。同様に、非平面グラフではノード数の増加に応じて指数関数的に計算時間が増加するような問題でも、平面グラフ、外平面グラフであればより速く計算できるアルゴリズムが存在する場合があります。

次数が低く、立体障害など空間的な制約もあるので、非平面の分子グラフは極めて稀です。おそらく天然低分子としては存在せず、狙って合成する必要があります(なんと実際に非平面グラフ分子の合成報告があります) [2]


環構造の検出

化学構造における環は、グラフ理論で言うところの閉路(cycle)に相当します(より厳密には、同じノードを2回以上通らない閉路なので、単純閉路(simple cycle)です)。あるグラフが閉路を持つかどうかは、ある一点のノードから隣接ノードを順番に辿っていくとわかります。既に到達したノードに別の経路から到達することができればそのグラフには閉路があります。

spanning.png

このような探索を行うと、上図のような経路が生成されます。これをスパニングツリー(spanning tree)と呼びます。これは元のグラフの部分グラフのうち閉路を持たない最大のものともいえます。

グラフにある閉路の数は、スパニングツリーに含まれないエッジの数と等しくなります(この数をcircuit rankと呼びます)。上図のグラフには4つの閉路があることになります。


Smallest set of smallest rings(SSSR)

これで環の個数は判明しましたが、環を構成する経路および環の大きさは何通りもあり得ます(下図参照)。

cyclebasis.png

このような場合、上図左のように環の個数と環の大きさが最小になるような組み合わせを選ぶことが多いです(smallest set of smallest rings; SSSR)。SSSRの決定にはグラフのminimum weight cycle basisを求めるアルゴリズムが利用できます。Hortonの方法や、Coelho de Pinaの方法が知られています [3][4]

どのような閉路の組み合わせを選んだとしても、先ほどのスパニングツリーに含まれない4つのエッジは、必ず全ての閉路に少なくとも一つずつは含まれます。つまり、これらのエッジは4つの閉路にそれぞれ一対一対応させることができます。このエッジ->閉路の組み合わせをfundamental cycle basisと呼び、ベクトルの集合で表現します(それぞれの閉路は総エッジ数と同じ次元のベクトルで、閉路を構成するエッジは1、そうでないエッジは0とするビット列)。これらのベクトルは、互いに排他的論理和をとると、それぞれの閉路を合成した閉路になるという特徴を持っています。SSSR決定アルゴリズムは、閉路ベクトル空間のこのような性質を利用しています。

fundamental.png


縮合環の検出

Bemis-Murcko scaffoldとして知られるような、縮合した環の構造(fused ring)は生理活性や易合成性の観点から創薬においてしばしば注目されます。このような部分構造は、グラフとしてみると全てのノードが2つ以上のエッジで繋がっているので、2-edge connected componentと呼ばれます。どのエッジを一か所だけ切断しても2つ以上のコンポーネントに分解されない部分グラフ、とも言えます。

scaffold.png

逆に、一つのエッジを切断すると2つ以上のコンポーネントに分解される場合、そのエッジを橋(bridge)と呼びます。橋はHopcroft-Tarjanアルゴリズムの応用で検出できます(下記リンク参照)。元のグラフから橋を取り除くことで2-edge connected componentが残ります。

Hopcroft-Tarjanのアルゴリズム(Biconnected component - Wikipedia)

https://en.wikipedia.org/wiki/Biconnected_component


部分グラフ同型性と構造検索

部分グラフ同型性(subgraph isomorphism)について、言葉で説明するのは難しいですが、化合物データを取り扱ったことのある人なら、部分構造一致と言えば直感的に理解しやすいと思います。

substr.png

部分グラフ同型性を判定する代表的なアルゴリズムとして、VF2アルゴリズムが知られています。これは比較的シンプルな深さ優先探索(DFS)ベースのアルゴリズムで、下図のように段階ごとに探索範囲を広げていきます。明らかに部分グラフ同型ではない場合、一旦前の段階に引き返して、別の可能性を探索していきます [5]

isomorph.png

部分グラフにはその生成方法によってノード誘導部分グラフ(node induced subgraph)やエッジ誘導部分グラフ(edge induced subgraph)などと呼ばれるものがあります。ノード誘導部分グラフは元グラフのノード集合の部分集合から一意に決まる部分グラフ、エッジ誘導部分グラフは元グラフのエッジ集合の部分集合から一意に決まる部分グラフです(下図)。

subgraph.png

我々のよく知っている部分構造検索は、おおよそedge induced subgraphの同型性を判別しています(例外として、SMARTSで[#16;X2;!R]のような一原子のクエリはエッジを含まないので別途対応)。

VF2アルゴリズムはnode induced subgraphの同型性を判定する手法なので、edge induced subgraphに適用する場合は、分子グラフのline graph(Wikipedia)を生成し、line graphのnode induced subgraphの同型性を判定します。その場合、line graphへの変換によってdelta-Y exchangeが生じていないことを確認する必要があります(詳しくは文献 [6])。

部分グラフ同型性の判定問題はNP完全なので、分子グラフのサイズが大きくなると指数関数的に計算時間が増加する可能性があります。が、drug-likeくらいのサイズであれば計算時間が問題になることはまずありません。実際のライブラリ検索ではVF2を適用する前にノード数、エッジ数、原子の種類、環の数・大きさなどで明らかに部分グラフ同型にならないものをあらかじめフィルタリングすることで高速化が可能です。


最大共通部分構造 (Maximum common substructure)

mcs.png

クエリ分子がデータベース分子と完全に一致しない場合でも、どれくらい共通構造があるかその度合いを知りたい場合があると思います。その度合いが最大共通部分構造 (Maximum common substructure; MCS)で、部分グラフ同型性と同様の手法で計算することが可能です。共通する結合(エッジ)数をそのまま閾値としたり、あるいはJaccard/Tanimoto係数のような類似度指標に変換することも可能です。

MCSの決定はNP困難であることが知られています。部分構造一致の場合、構造が一致した瞬間(あるいは一致しないことが確定した瞬間)に探索を打ち切ることができますが、MCSの場合は全ての可能性を探索するまで最適解を出力することができません。よって、部分構造一致と異なりdrug-likeくらいのサイズでもコンピューターをハングアップさせることがしばしばあります。そのため、計算時間の上限を決める、共通エッジ数が閾値を超えた時点で探索を打ち切る、高速な近似解アルゴリズムを使用するなどの工夫が必要です。

VF2アルゴリズムでもMCSを計算することが一応可能ですが、より効率化・用途特化されたアルゴリズムが数多く開発されています。オープンソースのMCSの実装としては、VF2同様DFSベースのアルゴリズムであるFMCSがRDKitにあります。他に、最大クリーク問題を解くアプローチもあり、こちらの方が近似アルゴリズムなどを実装しやすいと思います [7][8]

RDKit FMCS

http://rdkit.org/docs/source/rdkit.Chem.fmcs.fmcs.html


参考文献


  1. Konstantinova, E. V., & Skorobogatov, V. A. (1995). Molecular Hypergraphs: The New Representation of Nonclassical Molecular Structureswith Polycentric Delocalized Bonds. Journal of Chemical Information and Computer Sciences, 35(3), 472–478. https://doi.org/10.1021/ci00025a015

  2. Simmons, H. E., & Maggio, J. E. (1981). Synthesis of the first topologically non-planar molecule. Tetrahedron Letters, 22(4), 287–290. https://doi.org/10.1016/0040-4039(81)80077-9

  3. https://dare.uva.nl/search?identifier=93573ea1-c3ea-4321-a479-294c74b7f0bb

  4. Kavitha, T., Mehlhorn, K., Michail, D., & Paluch, K. E. (2008). An $\tilde{O}(m^{2}n)$ Algorithm for Minimum Cycle Basis of Graphs. Algorithmica, 52(3), 333–349. https://doi.org/10.1007/s00453-007-9064-z

  5. Cordella, L. P., Foggia, P., Sansone, C., & Vento, M. (2001). An improved algorithm for matching large graphs. Proceedings of the 3rd IAPR Workshop on Graph-Based Representations in Pattern Recognition. https://doi.org/10.1.1.101.5342

  6. Raymond, J. W., & Willett, P. (2002). Maximum common subgraph isomorphism algorithms for the matching of chemical structures. Journal of Computer-Aided Molecular Design, 16(7), 521–533. https://doi.org/10.1023/A:1021271615909

  7. Englert, P., & Kovács, P. (2015). Efficient Heuristics for Maximum Common Substructure Search. Journal of Chemical Information and Modeling, 55(5), 941–955. https://doi.org/10.1021/acs.jcim.5b00036

  8. https://www.jstage.jst.go.jp/article/ciqs/2017/0/2017_P4/_article/-char/ja/