異種ゲノムデータの統合による 遺伝子ネットワーク推定手法

第

54

2

333–356 2006 c

［研究詳解］

異種ゲノムデータの統合による 遺伝子ネットワーク推定手法

玉田 嘉紀

¹

²

²

（受付

2006

4

10

2006

6

28

要 旨

細胞内で生成されるタンパクは生物の主要な構成要素であり，遺伝子はその設計図に相当す る．遺伝子がタンパクに変換される時期・量の制御も遺伝子の働きによるものであり，生物は 遺伝子同士が協調して作用することによって生命を維持している．このような遺伝子間の依存 関係を，頂点と枝から構成されるグラフを用いて表現したものを遺伝子ネットワークという．

近年のマイクロアレイ技術の発展により，細胞内の遺伝子の活動状態を網羅的に観測できるよ うになり，遺伝子発現データとして蓄積されている．遺伝子発現データに基づく遺伝子ネット ワークの推定問題は，バイオインフォマティクスにおいて最も重要な課題の

1

遺伝子ネットワークの推定問題は，遺伝子の発現量を確率変数として見なすことにより，グラ フィカルモデルの推定問題として定式化される．しかし，ネットワークに含まれる遺伝子は一 般に数百以上と多く，そのためモデルに含まれるパラメータの数は膨大となる．したがって，

発現データへの過適合を避けるためのモデリングの方法論を構築することが必要不可欠といえ る．本稿ではこのような問題を解決するための方法として，著者らが開発した

2

1

DNA

2

両手法は，ベイジアンネットワークを遺伝子ネットワークのモデルとして用い，ネットワーク をグラフ構造の事後確率最大化に基づいて推定する．その際，配列情報および進化情報をネッ トワークの事前確率を構成するために用いることが特徴となっている．開発した手法はシミュ レーションおよび実データへの適用を通してその有効性を確認した．

キーワード： 遺伝子ネットワーク，遺伝子発現データ解析，制御配列，進化情報，ベ イジアンネットワーク．

1.

遺伝子は細胞内で生成されるタンパクの設計図に相当するもので，細胞の核内に存在する

DNA

DNA

1

GNI

3–8–33–608）

2

4–6–1

て

mRNA

この様に，遺伝子が読み取られタンパクに変換されることを「遺伝子が発現する」と表現する．

遺伝子の発現を制御している物質は転写因子と呼ばれるが，転写因子もまた遺伝子から合成さ れたタンパクの一種である．そのため，ある転写因子が発現すると，その被制御遺伝子の発現 制御を促すことにつながる．このように遺伝子の発現は互いに依存しており，協調し複雑に影 響しあうことによって生命活動の維持に必要なタンパクを適切に制御している．遺伝子ネット

ワーク（

gene network

たものである．遺伝子ネットワークでは，各遺伝子はグラフ中の頂点として表され，遺伝子間 の発現の依存関係は頂点を結ぶ有向枝で表される．

近年，DNAマイクロアレイなどのツールの登場によって，様々な実験的状況下において細 胞内の遺伝子の発現状態を網羅的に観測することができるようになった（DeRisi et al., 1997）．

それらの遺伝子発現データを用いて，遺伝子ネットワークを推定する研究が盛んにおこなわれ ている（van Someren, 2002）．遺伝子ネットワーク推定の問題は，遺伝子の発現量を確率変数 としてとらえたグラフィカルモデルの推定問題として定式化される．従来のグラフィカルモデ ル推定問題と比較して，きわめて多くの頂点（遺伝子）からなる構造未知のネットワークを数限 られたサンプル（マイクロアレイデータ）から推定することがその特徴であると言える．また，

ネットワークの制御や挙動予測よりもネットワークの構造推定に重点が置かれるのもその特徴 である．実際の生物における具体的な遺伝子数としては，出芽酵母だと約

6000，ヒトに至って

は

20000

1000

同一細胞組織（肝細胞，血管内皮細胞など）のマイクロアレイデータは，コストなどの面からそ の数は数百のオーダーでしか集まらないのが現状である．つまり，モデルに含まれるパラメー タ数がサンプル数よりも多くなる状況が頻繁に生じ，そのためデータへ過適合し，高精度な遺 伝子ネットワークを推定することは非常に困難な問題となっている．

著者らは，このような問題に対し

2

1

DNA

Tamada et al., 2003），他方は異なる生物種間に進化的に保存されている情報（Tamada et al., 2005）をそれぞ

同一の転写因子から直接の制御を受ける遺伝子の

DNA

2

本論文の構成は次の通りである．2章では，ベイジアンネットワークによる遺伝子ネットワー ク推定法を説明し，関連研究としてこれまでに行われてきた遺伝子ネットワーク推定研究の概 要を説明する．3章では制御配列の情報を利用した手法，4章では進化情報を利用した手法を 解説する．

本論文で使用したシミュレーションや実データへの適用の際のデータは，公開可能なものを 著者らのウェブページ（http://bonsai.ims.u-tokyo.ac.jp/˜tamada/supplement.html）に掲載して

いる．

2.

2.1

本節では，Imoto et al.（2002）によって提案された

B -スプラインを用いたノンパラメトリッ

ベイジアンネットワークは，確率変数間の条件付き独立性を非循環有向グラフを用いて表し たものである．グラフの各頂点に確率変数が

1

1

1 , . . ., X p

p

各変数間の依存関係を表すグラフ構造を

G

1 , . . ., X p

Pr( X 1 , . . . , X p ) =

p

j=1

Pr( X j |P a ( X j )) . (2.1)

ただし，P a

( X j )

X j

P a ( X j ) = ∅

1

X j

X j

今，p個の遺伝子

X 1 , . . . , X p

N

N × p

X

X

( i, j )

x ij

i

X j

G = ( V, E )

= {X 1 , . . . , X p }

⊆ V × V

2.1

f (X|G,Θ) =

^N

i p

j

f j ( x ij |pa( x ij ) , θ j ) .

ただし，pa

( x ij )

i

X j

Θ = (θ ₁ , . . . , θ _p )

j (1 ≤ j ≤ p )

X j

図

1.

発現データ

x ij

pa( x ij ) = ( p ^(j) _i1 , . . . , p ^(j) _{iq j} )

x ij = m j1 ( p ^(j) _i1 ) + ··· + m jq j ( p ^(j) _{iq j} ) + ε ij .

jk ( · )

B

m jk ( p ^(j) _ik ) =

M _jk

m=1

γ _mk ^(j) b ^(j) _mk ( p ^(j) _ik )

と表せる．ここで

{b ^(j) _1k ( · ) , . . . , b ^(j) _{M jk ,k} ( · ) }

M jk

B

{γ _1k ^(j) , . . ., γ _M ^(j) _{jk ,k} } ∈ R ^{M jk}

B

ε ij

0

σ ² _j

x ij

f j ( x ij |pa( x ij ) , θ j ) = 1

2 πσ _j ² exp

− {x ij −

^q _k=1 ^j m jk ( p ^(j) _ik )} ² 2 σ _j ²

.

発現データからの遺伝子ネットワークの推定問題は，遺伝子ネットワークのグラフ構造

G

G

G

発現データ

X

G

π ( G|X ) ∝ π ( G ) · π ( X|G )

= π ( G )

f (X|G,Θ) π (Θ|λ) dΘ.

(2.2)

ただし

π ( Θ|λ )

λ

Θ

従って，遺伝子ネットワークの構造は，

MAP

G ˆ = argmax _G π ( G|X )

20

DAG

2 . 34 × 10 ⁷²

そのため発見的な探索法である

greedy hill-climbing

Friedman et al., 1998; Imoto et al., 2002

Ott et al.

2004

greedy hill-climbing

greedy hill-climbing

Greedy hill-climbing

S ( E )

Input:

X，頂点（変数）集合 V = {X 1 , . . . , X p }，親候補数 m， 反復回数 T Output:

G = ( V, E )

Step 1. E := ∅

:= 1.

Step 2.

i , j

1 ≤ i, j ≤ p , i = j）に対し s ij := S ( { ( X i , X j ) } )

. Step 3. C j :={s _α (j)

1 , j , . . ., s _α (j)

m , j },

{α ^(j) ₁ , . . . , α ^(j) m }

s 1j , . . . , s pj

m

Step 4. 1 , . . . , p

π 1 , . . ., π p

Step 5.

X π j (1 ≤ j ≤ q )

X π j

Y ∈ C π j ∪ {X : ( X, X π j ) ∈ E}

を考え，以下の操作で得られるスコア

S ( E )

（

a

:= E ∪ { ( Y, X π j ) } if ( Y, X π j ) ∈ E

（b）E

:= E \ {( Y, X π j )} ∪ {( X π j , Y )} if ( Y, X π j ) ∈ E

（c）E

:= E \ { ( Y, X π _j ) } if ( Y, X π _j ) ∈ E

（a）

,（b） ,（c）のうち最もスコアが向上する操作を採用し，E := E

操作は

E

Step 6.

Step 4 ∼ 5

た場合，その時点のグラフを

G t

:= t + 1 .

Step 7. t ≤ T

E = ∅

Step 4 ∼ 6

Step 8. G 1 , . . . , G T

実際にここで示した

greedy hill-climbing

S ( E )

S ( E ) =

X _j ∈V

S j ( X j , P a ( X j )) .

事後確率に基づくスコア（式（2.2））はこのような分割が可能である．

2.2

遺伝子ネットワークを推定するために，これまで様々な数理モデルとその推定のためのアル ゴリズムが提案されてきた．ブーリアンネットワークは遺伝子の発現量を

On

Off

2

AND, OR, NOT

Liang et al.

1998

Akutsu et al.

1998

Akutsu et al.

2

Akutsu et al.

2000

Shmulevich et al.

Murphy and Mian

1998

Friedman et al.

D’haeseleer et al.

1999

Chen et al.

（1999），de Hoon et al.（2003）は線形微分方程式を利用した方法をそれぞれ提案している．ベ イジアンネットワークを連続値データに適用したものとしては

Friedman et al.

2000

Murphy and Mian, 1999; Kim et al., 2004

1998; Higuchi and Kitagawa, 2000）に基づく遺伝子ネットワーク推定手法（Rangel et al., 2004;

Yamaguchi and Higuchi, 2006; Yamaguchi et al., 2006; Hirose et al., 2006

Toh and Horimoto

Basso et al.

Imoto et al., 2002

マイクロアレイにより細胞内で発現している

mRNA

マイクロアレイ実験には金銭的コストが必ずかかりデータをいくらでも多く観測するわけには いかない．また，（1）マイクロアレイによる

mRNA

2

mRNA

mRNA

RNA

さらに遺伝子ネットワーク推定を生物学における知識発見にまで結びつけるようなより高精度 なものにするために，マイクロアレイデータに加えて遺伝子発現制御に関わる生物学的データ の併用が研究されている．

Hartemink et al.

2002

た．

Imoto et al.

2004

な枠組みを提案している．本稿で解説する

Tamada et al.

ネットワーク推定と配列探索の

2

の方法は，ネットワーク中の各枝に対して生物学的な裏付けがあるか・ないかの離散的な情報 をネットワーク推定に使用するものであるが，一般に生物学的な情報は連続値をとる場合が多 い．この弱点を克服する方法として

Bernardo and Hartemink

2005

Imoto et al.

Bernardo and Hartemink

1

ないが，

Imoto et al.

2006

組みを示し，それを利用して薬剤反応経路の同定を行った．本稿で解説するもう一つの方法で ある

Tamada et al.

2005

2

3.

本章では，Tamada et al.（2003）が開発した，遺伝子の

DNA

DNA

DNA

DNA

このような共通に存在する制御配列の情報を活用することによって，直接の制御か否かを識別 できる可能性がある．図

2

FKH2

3

ACE2，CLB2，SWI5

FKH2

3

ACE2，CLB2，SWI5

GTAAACA

開発した手法の実行手順を図

3

図

2.

図

3.

情報と発現データを組み合わせて，遺伝子ネットワークを再推定する．制御配列の情報はネッ トワークの事前確率を構成するために利用される．以上のステップを推定されるネットワーク が変化しなくなるまで繰り返すことにより共通配列の探索と遺伝子ネットワークの推定を同時 に行う．

遺伝子の

DNA

3.1

本節では推定した遺伝子ネットワークの情報を利用して，DNA配列中に共通に存在する配 列を探索する共通配列探索手法について述べる．共通配列探索の対象となる遺伝子は，推定さ れた遺伝子ネットワークの構造に基づいて選択される．まず，転写因子の候補となる遺伝子を 選び出し，その転写因子から制御を受けている可能性のある遺伝子を選び出す．具体的には，

ある転写因子の候補に対してその直接の子供，および孫からなる遺伝子群を制御配列探索の対 象とする．ここで転写因子の候補と書いたのは，実際には転写因子として働くが，転写因子で あることが分かっていない遺伝子を考慮したためである．つまり，ネットワーク推定の結果か ら多くの遺伝子（ここでは

4

そのような遺伝子を転写因子の候補とした．

一般的にクラスター分析に基づいた共通配列探索では，配列探索の対象となる遺伝子すべて に共通の配列が存在することが仮定されている．一方，推定されたベイジアンネットワークに 基づき共通配列探索の対象となる遺伝子を選択する本手法では，各遺伝子毎に転写因子候補に 対してそれらがどの程度被制御遺伝子として尤もらしいかを尤度として計算することができる．

この尤度が高い遺伝子には転写因子候補に対して共通配列が存在する確率も高く，逆に尤度の 低い遺伝子には共通配列が存在する確率は低いと考えるのは自然と言えよう．このような考え に基づいた共通配列探索を行えば，推定された遺伝子ネットワークの情報をより有効に活用し た生物学的にもっともらしい共通配列を探索できることが期待できる．Bannai et al.（2002）が 開発した

string pattern regression

X TF

X j

S ({( X TF , X j )})

⊂ Σ ^∗ × R

Σ ^∗

は文字列の集合，Rは実数全体の集合とする．つまり

D

DNA

String pattern regression

v

MSE ( D, v ) =

(s,r)∈D v ( µ ( D v ) − r ) ² +

(s,r)∈D v ¯ ( µ ( D ¯ v ) − r ) ²

|D| .

ただし，|D|は集合

D

v

v

D

¯ v = D \ D v , µ ( D ) =

(s,r)∈D r

|D |

D ⊂ D

v

substring pattern

Bannai et al., 2002）．

3.2

次に共通配列探索の結果に基づきグラフ構造

G

π ( G )

X TF

Y 1 , . . . , Y n

DNA

( X TF , Y 1 ) , . . . , ( X TF , Y n )

G

π ( G )

G = ( V, E )

X i → X j

ζ 1

ζ 2

< ζ 1 < ζ 2

2

π ( G )

π ( G ) = Z ⁻¹ exp

−

(X i ,X j )∈E

ζ α(i,j)

.

ただし

α ( i, j )

X i → X j

1

2

3.3

開発した手法のアルゴリズムを以下にまとめる．

Step 1.

遺伝子ネットワークを推定する．推定したネットワークを

G = ( V, E )

Step 2.

X ∈ V

X

G

X

る集合とする．|D

X | ≥ 4

X

X

X

Step 3.

X TF

てネットワーク

G

G

G

（a）共通配列探索によって発見された共通配列が存在する

X TF

Y

X TF → Y

G

（b）X

TF

W

DNA

Step 2

X TF → W

G

Step 4.

π ( G )

( G )

greedy hill-climbing

G

Step 5. Step 2

Step 4

図

4.

このアルゴリズムでは得られるネットワークに変化がなくなるまで遺伝子ネットワークの推 定と共通配列探索を繰り返し行う．この方法ではアルゴリズムの収束性は保証されないが，現 実的には比較的少ない回数で終了する．

図

4

3.4

開発した手法の有効性を検証するために，モンテカルロシミュレーションと実データへの適 用を行った．

3.4.1

データ. 図

5

100

DNA

signal-to-noise ratio

DNA

A，T，C，G

4

gene7

TATAT

1000

結果. モンテカルロシミュレーションの結果を表

1

Sp

“specificity”

Sn

“sensitivity”

で正解の枝のうち推定することのできた枝の割合を表す．事前確率

π ( G )

specificity

sensitivity

1 = 1 . 0 , ζ 2 = 7 . 0）．発現データと共通配列の情報を組み合わせることによって，特に Sp

0.384 → 0.540

12,727 → 4,934

TATAT

1,000

433

6

図

5.

表

1.

トワークの典型的な例を示す．以上のことからモンテカルロシミュレーションにおいては，共 通配列の情報を用いることにより，より正確な遺伝子ネットワークが推定できることが確認で きた．

3.4.2

データ. 次に開発した手法を実データへ適用した結果を示す．対象とする生物種は出芽酵母 であり，遺伝子発現データは遺伝子破壊実験によって得られたマイクロアレイ

120

図

6.

図

7.

（

Aburatani et al., 2003

Imoto et al.

2002

CHA4

GAL11

SWI6

3

124

DNA

GenBank

（http://www.ncbi.nlm.nih.gov/Genbank/index.html）データベースから上流

800 bp

結果. 本手法をデータに適用した結果，

4

CHA4

7

8

CHA4

MIPS

functional category

GAL11

GAL2

GAL2

GAL2

ARG2

GAL11

GAL2

図

8.

表

2.

MCM1-FKH2

REB1

ARG2） .

4

AAAGA，AAACG（2

このうち

TAAAC

FKH2

FKH2

MCM1

ACE2

SWI5

MCM1

FKH2

TAAAC

MCM1

MCM1-FKH2

実際に，それらの遺伝子の

DNA

MCM1

REB1

ARG2

2

FKH2

MCM1

FKH2

FKH2

2

CHA4

MCM1-FKH2

REB1

と

ARG2

MCM1-FKH2

MCM1-FKH2

4.

本章では

Tamada et al.

2005

法について解説する．細胞周期など，細胞が生きていくために必要な基本的な機能は，異なる 生物種においてもよく保存されており，それらに関わる遺伝子から生成されるタンパク配列も よく類似している．このような遺伝子は，進化の過程において同一の遺伝子を起源としてもつ と考えられており，そのような遺伝子同士を

ortholog

A

B

X a , X b

A

a

b

B

a

Y a , X b

Y b

ortholog

X a

Y a , X b

Y b

A

X a → X b

生物種

B

Y a → Y b

X a

X b

a

Y b

9

KEGG

http://www.genome.ad.jp/kegg/

ortholog

BLAST E-value

BLAST

1990）によって計算されるタンパク配列の類似度を表す指標で，値が小さいほど配列が類似し

4.1

今，2つの生物種

A

B

G A = ( V A , E A )，G B = ( V B , E B )

V A , V B

A, B

A , E B

A

B

H AB

図

9.

KEGG

value）．ただし MCM

ORC

間の

E-value