Transcript 第3部 - 京都大学

集中講義（九州大学数理学研究院）
バイオインフォマティクスにおける
カーネル法およびグラフ理論 (3)
構造のパターンマッチング
阿久津 達也
京都大学 化学研究所
バイオインフォマティクスセンター
内容





最大共通部分グラフ問題
最大共通部分木計算アルゴリズム
木の編集距離、木のアラインメント
ネットワークモチーフ
酵素反応式からの構造変換規則の推定
グラフの同型性、部分グラフ

G1(V1,E1) と G2(V2,E2)が同型



G2(V2,E2) が G1(V1,E1) の部分グラフ iff.
V2  V1 , E2  E1
部分同型性判定問題


V1 からV2 への１対１、上への写像（同型写像）φが存在し、
{u,v} E1 iff. {φ(u),φ(v)} E2
G2(V2,E2) は G1(V1,E1) の部分グラフに同型か？
例（図）： G1 と G2 は同型、 G3 は G1 の部分グラフと同型
G1
G2
G3
最大共通部分グラフ問題


入力： 無向グラフ G1(V1,E1), G2(V2,E2)
出力： G1 と G2 の両者の部分グラフと同
型で、かつ、頂点数最大の連結グラフ
G1
G2
最大共通部分グラフ問題のNP困難性


命題： 最大共通部分グラフ問題はNP困難
証明

ハミルトンパス問題（入力を G1 とする）を多項式時間還元


ハミルトンパス問題： 各頂点をちょうど一度づつたどるパスはあるか？
G2 として| V1 |個の頂点からなる鎖を考える

最大共通部分グラフのサイズが |V1| iff. ハミルトンパスが存在
G1
G2
最大共通部分木アルゴリズム(1)




最大共通部分問題はNP困難だが、木に対しては多項式時間
木： 閉路の無い連結なグラフ （RNA二次構造、進化系統樹、糖鎖など）
アルゴリズム： 動的計画法 ＋ 最大重み二部グラフマッチング
根つき木として取り扱う

O(nm)個の根のペアを考えれば良いから
（n,m: 木の頂点数）

根どうしはマッチするものと考える

子に順番は無いものとする
T1
T2
RNA二次構造
最大共通部分木アルゴリズム(2)


T(v): 頂点 v とその子孫から
構成される T の部分グラフ
D(u,v)： T1(u) と T2(v) の最大
共通部分木のサイズ（頂点
数）


ただし、 u と v は対応すること
が必要
動的計画法の適用


例
D(2,b) = 1, D(2,c)=1, D(2,d)=1
D(3,b) = 3, D(3,c)=3, D(3,d)=2
D(4,b) = 3, D(4,c)=2, D(4,d)=2
D(u,v) を葉から根へ
と計算
D(u,v) を計算する時
には、その子のペア
(u’,v’) の D(u’,v’) は
すべて計算済み
2
b
4
3
5
a
T1
1
6
7 e
8
9
T2
c
f g
j
d
h
i
最大共通部分木アルゴリズム(3)

D(u,v) の計算



2
5


1
U と W の間にしか辺が無いよ
うに頂点集合を (U,W) に分割
可能なグラフ
マッチング
端点を共有しない辺の集合
最大重みの場合も含め、多項
式時間で計算可能
5
4
二部グラフ


1
子のペアに対する D(…) から二
部グラフを構成
二部グラフの最大重みマッチン
グを計算
7
2
U
W
マッチング
最大共通部分木アルゴリズム(4)

D(1,a) の計算
 D(1,a) = 最大重みマッチン
グの重み ＋ 1 = 8
 (2,d), (3,c), (4,b) が対応
2
b
4
3
5
a
T1
1
6
8
D(2,b) = 1, D(2,c)=1, D(2,d)=1
D(3,b) = 3, D(3,c)=3, D(3,d)=2
D(4,b) = 3, D(4,c)=2, D(4,d)=2
c
f g
7 e
9
T2
j
h
d
2
i
3
1
1
b
3
1
3
c
2
3
4
2
2
d
木の編集距離 (1)


文字列間の編集距離の
根つき木への拡張
３種類の編集操作


a
c
頂点ラベルの置換、頂点の削除、
頂点の挿入
T1, T2 間の
編集距離

a
T1 を T2 に変え
るのに必要な
最小の編集
操作回数
c
挿入
置換
削除
a
a
c
b
c
b
c
c
a
c
b
a
a
b
c
b
木の編集距離 (2)

T1, T2 間の mapping のコスト




１対１、親子関係を保存
コスト： ラベルが異なる頂点数 ＋ mapping に含まれない頂点数
補題： 編集距離 ＝ 最小コストmappingのコスト
編集距離の計算


順序木（子の順序を保存）： 動的計画法により、多項式時間
無順序木： NP困難
a
c
a
d
b
e
c
c
コスト＝２
b
d
c
e
木のアラインメント

二つの木から一つの木を構成
それぞれの木に空白ラベル頂点（ギャップに相当）を挿入することにより同
一の木を構成


木アラインメントの計算
順序木、および、次数限定の無順序木： 動的計画法により、多項式時間
次数に制約の無い無順序木： NP困難


a
e
b
d
c
(a,a)
a
b
(e,-)
f
c
d
(b,b)
(c,-)
(-,f)
(-,c)
(d,d)
ネットワークモチーフ

モチーフ



ネットワークモチーフ



（ランダムなネットワークと比べて）実際のネットワークにおい
て頻出する（統計的に有意に）ネットワークのパターン
ネットワークのパターン： 部分グラフ
ネットワークモチーフの検出


配列解析において現れる機能と関連した配列パターン
例： L-x(6)-L-x(6)-L-x(6)-L (ロイシンジッパーモチーフ)
（多数のグラフからの）最大共通部分グラフ問題を含むため、
一般にNP困難 ⇒ 多くのヒューリスティックなアルゴリズム
ネットワークモチーフの例



フィードフォワード制御
Single Input Module
Dense Overlapping regulons
ネットワークモチーフの例
feedforward
loop
single input
module
X
dense overlapping
regulons
X
Y
Z
Z1 Z2
Zn
crp
argR
arg I
argF
argE
argD
araBAD
argCBH
araC
X1 X2 X3
Xn
Z1 Z2 Z3
Zn
内容





最大共通部分グラフ問題
最大共通部分木計算アルゴリズム
木の編集距離、木のアラインメント
ネットワークモチーフ
酵素反応式からの構造変換規則の推定
構造変換規則の推定問題



化学反応式から、どの部分構造が置き換わったのかを同定
（可能な組み合わせの数え上げ）
X-A + Y-B → X-B + Y-A から、A,B,X,Y を同定
応用分野：トレーサー実験のシミュレーション、薬剤設計、代謝
データベースのチェックなど
O
H
O
O
O
H
A
NH2
O
H
O
B
O
O
H
O
O
O
H
O
O
H
B
O
O
H
O
O
NH2
A
O
H
問題の例


解は一種類とは限らない⇒数え上げも必要
以下の例では、(A,B)の組として３種類が可能
O
H
O
O
O
H
NH2
O
H
O
O
O
H
O
O
O
H
O
O
H
O
O
H
O
O
NH2
O
H
酵素反応データベース

KEGG (Kanehisa et al., Kyoto Univ.)



http://www.genome.ad.jp/kegg/
Kyoto Encyclopedia of Genes and Genomes
LIGAND (KEGGの一部)


数千以上の酵素反応式
ARM system (Arita, Univ. of Tokyo)




http://www.metabolome.jp/arm.html
反応式 + 解析ソフト
トレーサー実験のシミュレーションなど
構造変換規則抽出の機能も備える
(最大共通部分グラフに基づく方法を利用)
従来法との比較

従来法（最大共通部分グラフに基づく方法）




数え上げが困難
NP困難
多数の反応式からの規則獲得は難しい
提案手法





グラフ分割とグラフ同型性判定の組み合わせ
数え上げも比較的容易
実用上のほとんどの場合に多項式時間
多数の反応式からの規則獲得も可能
実用バージョンは簡単にインプリ可能
諸定義

化学構造

次数D以下の連結無向グラフ




頂点はラベルを必ず持ち、辺にラベルがあっても可
化学反応式 I 1 + … + I p ⇔ O 1 + … + O q


頂点数をnとすると辺の個数はO(n)
同型性判定、正規形計算が多項式時間で可能
左辺のラベルの多重集合
＝右辺のラベルの多重集合
化学構造 G に対する
大きさ C の chemical cut


高々C 個の辺を削除して、
G をいくつかの連結成分に分解
連結成分と削除した辺はstarを構成
C=5のchemical cut の例
問題の定義

入力



化学反応式
 I 1 + I 2 + … + I p ⇔ O 1 + O2 + … + O q
整数 C >0
出力

I 1,…I p,O 1,…,O q のそれぞれに対する大きさ C の
chemical cut で以下を満たすもの
 右辺からcutにより得られる連結成分の多重集合
＝左辺からcutにより得られる連結成分の多重集合
（ただし、同型な化学構造は同一の要素とみなす）
結果
C,p,q 全てに制約⇒多項式時間
 C=2 もしくは p=q=2 でも、NP困難
 C=1,p=q=2,入力が木 ⇒ O(n1.5)時間



（いくつかの拡張が可能）
実用的アルゴリズム
実データ（KEGG/LIGAND）を用いた計算
機実験
多項式時間アルゴリズム

p=q=2, C=1 の時のみを示す（他も同様）
forallpartition ( X , A ) of I 1 do
forallpartition (Y , B ) of I 2 do
forallpartition ( X , B ) of O 1 do
forallpartition (Y , A ) of O 2 do
X
ifX
and Y
Y
and A  A and B
B
then output ( A , B , B , A )
O
O
O
H
H
O
NH2
Ai
O
H
O
O
Bj
O
O
O
H
O
H
H
O
O
Bk
O
H
O
O
Ah
O
NH2
O
H
NP困難性

C が2でも、p,q に
制約が無いならNP
困難


３DMを還元
p=q=2でも、C に
制約がないならNP
困難
{{A,1,a},{A,1,c},{A,2,b},
{B,2,c},{C,2,a},{C,3,b}} からの還元
A
A
A
B
C
C
ββββββ
αααααααααααα
ββββββ
1 a 1 c 2 b 2 c 2 a 3 b
A
A
A
B
C
C
α α ββ ββ α α α α ββ
ββββββαααααα
1 2 3 a b c 1 2 2 a b c
NP困難性からの示唆

構造を木に限定してもNP困難



木の最大共通部分構造検出は多項式時間で実行
可能
最大共通部分構造検出法では対応できない場合
が存在することを示唆
カットのサイズおよび入出力の個数が多きくな
るにつれ、計算量が増大するのはほとんど不
可避
グラフの正規形
正規形：同型なグラフは、そのデータ構造が全く同一
となるように、頂点を順番づけ
（正規型のサイズ＝グラフのサイズ）
 本手法で用いるID付け（IDは定数サイズ）
 アルゴリズム中で生成されるグラフにのみ計算




同型なグラフ ⇔ 同型なID
木の場合(ただし、C=1): O(n)時間
平面的グラフ(ただし、C=1): O(n2)時間
（グラフ１個あたりは O(n) 時間）
木の場合のIDづけ


頂点にボトムアップに
IDをつける
サイズ１のchemical
cut で生成される木に
対するIDづけ



すべての辺（両方向）
に対してIDをボトム
アップに計算
木の個数: O(n)
計算時間: O(n)
4
1
5
2
1
3
特殊な場合の高速化



p=q=2、C=1、入力構造は木：X-A + Y-B ⇔ X-B + Y-A
長さ４の有向閉路検出問題に変換
O(n1.5)時間(← cycle finding algorithm by Alon et al.)
forallpartition ( X i , Ai ) , ( Y j , B j ) , ( X k , Bk ) , (Yh , Ah ) of I1 , I 2 , O1 , O2 do
add edges {( X i , Ai ) , ( Y j , B j ) , ( X k , Bk ) , (Yh , Ah )};
find a directed cycle of length 4
X
X A
Y B
Na O H
H Cl
X B
Y A
Na Cl
H O H
A
B A’
Y
X’=NaO A’=H
Y’=OH
X’
Y’
長さ４の閉路の見つけ方

対象とするグラフ（有向）


高次数頂点：次数がΔ以上



A,B,X,Y はdisjoint、|V|=n, |E|=O(n)
そのような頂点の個数は O(n/Δ)
各頂点からDFSで閉路を探索
→ O(n2/Δ)時間
X’
X
A
以下のパスを数え上げ




A → X → B 型のパス
B → Y → A 型のパス
パスの総数は O(nΔ)
A,Bが一致するパスを radix ソート
で検出 → O(nΔ)時間
Δ=0.5n とおいて、O(n1.5)時間
X’’
B
Y
高次数頂点を除去したグラフ


B’’
B’
Y’
A’
Y’’
A→X→B
A→X→B’
A→X’→B’’
A’→X’→B’’
B→Y→A
B→Y→A’
B→Y’→A’
B’→Y’’→A’
複数の反応式からの規則の抽出：
特殊な場合の高速化

p=q=2、C=1、入力構造は木、反

応式は２個：
X-A + Y-B ⇔ X-B + Y-A
X’-A + Y’-B ⇔ X’-B + Y’-A
図aの部分グラフを検出（→単純な
方法ではO(n 2)時間）
図bへの変換を一部の頂点（次数
Δ=n 1/5以下の頂点）にのみ適用す
ると、O(n 2-1/5)時間




高次数頂点数 O(n/Δ)⇒高次数頂
点の処理にO(n2/Δ)
高次数頂点を除いた後、変換した
グラフのサイズはO(nΔ)⇒変換後
のグラフの処理にO((nΔ)1.5)時間
(nΔ)1.5=n2/Δを解く
X
(a)
X’
A
B
Y’
Y
<X,Y> (b)
X <X,Y>
A
A
Y
B
A
<X’,Y’>
反応式が K=2L 個の場合

p=q=2、C=1、入力構造は木、反応式が K=2L 個：

X-A + Y-B ⇔ X-B + Y-A
X’-A + Y’-B ⇔ X’-B + Y’-A
X’’-A + Y’’-B ⇔ X’’-B + Y’’-A
…
2個の場合を繰り返し実行(2L個 → 2L-1 個 → …)
K  2 で O(n
L
n
1
3 f ( L )1
2  f (1L )
) 2
L 1
で O(n / )  O(( n)
とおくと上式は O(n
2
2  3 f (1L )1
2  f (1L )
)。
)となるので、
f ( L  1)  3 f ( L)  1 。
f (0)  2 より、f ( L)  32  3L  12 。よって、O(n
2 3
1
log
3 1
K
2
2
)。
実用的アルゴリズム



正規形計算にMorganアルゴリズム(1965)を
利用。正当性は必ずしも保証されないが、ほ
とんどの場合に正しく動作。O(n2)時間
長さ４の閉路発見のために配列を利用
C=1,p=q=2の場合、O(n 3)時間


O(n)個の部分構造に対して正規形を計算するた
め
平面的グラフに対しては線形時間の正規形計算
アルゴリズムを用いることにより、O(n2)時間で常
に正しく動作
Morganアルゴリズム
①
②
③
最初は頂点の種類に応じて頂点に番号を割り当て
自分の番号と隣接する頂点の番号を加算。これを
全ての頂点に対し同時に実行。
②を異なる番号数が増加しなくなるまで繰り返す。
H
H
C
H
O
C*
H
*
*
*
H H C C O
1 1 2 2 3
3 3 7 8 5
10 11 24 23 13
計算機実験(1)

KEGG/LIGANDデータベースを利用


以下の二つの場合をテスト



各種生物の生体内で利用される化学反応式などを格納した
データベース（京大化研金久研で開発）。国際的に有名。
X-A + Y-B ⇔ X-B + Y-A
X-Y ⇔ X + Y
（反応式数＝2346)
（反応式数＝367)
最初の場合については以下の３種類をテスト



3a: 水素原子、辺の種類の両者を考慮
3b: 水素原子、辺の種類を無視。
ただし、切断辺に関する情報は考慮
3c: 水素原子、辺の種類、切断辺に
関する情報を無視
H O H
H
O H
H
O H
計算機実験(2)

実用アルゴリズムの詳細




O
3a: 水素原子、辺の種類の両者を考慮
3b: 水素原子、辺の種類を無視。ただ
し、切断辺に関する情報は考慮
3c: 水素原子、辺の種類、切断辺に関
する情報を無視
HO
P

O
O
H
H
O
NH2
Ai
O
H
O
Bj
H
O
P
OH
HO
HO
O
H
H
O
O
Bk
O
P
OH
HO
O
O
H
H
OH
HO
O
O
O
OH
P
O
3a は下の変換の検出に失敗
3b は右の変換の検出に失敗
O
O
HO
例

O
H
O
O
Ah
O
NH2
O
H
計算機実験(3)

総計算時間、解が計算できた式の
個数、解の平均個数を計測







3a: 水素原子、辺の種類の両者を考慮
3b: 水素原子、辺の種類を無視。ただ
し、切断辺に関する情報は考慮
3c: 水素原子、辺の種類、切断辺に関
する情報を無視
C=1,p=q=2については3cがBEST
X-A + Y-B ⇔ X-B + Y-A
計算時間

ワイルドカードが反応式に存在
環の切断が反応式に存在
解の平均個数
3a
150.4 sec.
1104
2.50
3b
68.0 sec.
1071
2.18
3c
66.7 sec.
1912
2.15
1反応式あたり平均で2個程度
1反応式あたりの処理時間は十分
高速（0.1sec以下）
解無しとなった場合の主な理由

解あり反応式数
(全反応式数=2346)
X-Y ⇔ X + Y
計算時間
4c
13.4 sec.
解あり反応式数
330
(全反応式数=367)
結論と課題

結論




化学反応式から構造変換規則を抽出する新しい手法を提案
一般の場合には、NP困難であるが、多くの実用的な場合には多項式時間
実用的アルゴリズムも開発
課題

C=1,p=q=2以外の場合の効率化




最小カット数の性質が有用？
部分構造に対する正規形計算法の改良
解に順位をつけて出力するための
ヒューリスティクスの開発、組み込み
環構造の変換、不斉炭素などへの対処
O
NH2
HO
NH2
O
NH3
HO
HN
参考文献




G. Valiente: Algorithms on Trees and Graphs, Spriner, 2002
D. Shasha & K. Zhang: Fast algorithms for the unit cost editing distance
between trees, J. Algorithms, 11:581-621, 1990.
T. Jiang, L. Wang & K. Zhang: Alignment of trees – an alternative to tree
edit, Theoretical Computer Science, 143:137-148, 1995.
T. Akutsu: Efficient extraction of mapping rules of atoms from enzymatic
reaction data, J. Computational Biology, 11:449-462, 2004.