『戦略分野１「予測する生命科学：医療および創薬基盤」』

2014-03-14

SiGN講習会資料：ベイジアンネットワークを用いた遺伝子ネットワークの推定と解析
土井 淳
atsushi_doi@cell-innovator.com
株式会社セルイノベーター
研究開発部
福岡市東区箱崎6-10-1
九州大学 産学連携棟I アントレプレナーシップ・センター 2階
http://www.cell-innovator.com

cell innovator

1. マネーボール：統計学の応用
2. 遺伝子発現とベイジアンネットワーク
3. 遺伝子ネットワーク

cell innovator

1. マネーボール：統計学の応用

cell innovator近年の統計学にまつわるトピック• マネーボール理論：経営論の参考にも。日経BP -http://special.nikkeibp.co.jp/ts/article/aaaa/114314/• ビッグデータ：Google、Facebook、Amazon などの企業によるイメージ。
• データアナリスト、データサイエンティストが25万人不足。
http://www.nikkei.com/article/DGXNZO57421630X10C13A7EA1000/「大量のデータを統計学を使って、なんとかしよう」というのがトレンド
cell innovatorマネーボール理論とは？
• 野球をアウトを取られないようにするゲームと定義。過去のデータをもとに導きだされた理論。
• バントをするな。
• フォアボールでいい。
• 初球に手を出すな。
• 盗塁もダメ。

• バントされても、2塁に投げるな。

安い選手で効率よく勝つための理論
cell innovatorここまでの話で、、、• 近年、統計学的なアプローチが、よく用いられるようになった。
• 統計学的なアプローチから得られたものが、必ずしも人間の直感に合わない。（裏、裏、裏と来たら、次は表と思いたいのが心情。）• 直感に合わなくても、役に立つかもしれない。（マネーボール理論のアスレチックス
は、シーズン中に20連勝。レッドソックスは、ワールドシリーズ優勝。）

cell innovator

2. 遺伝子発現とベイジアンネットワーク

cell innovator遺伝子発現と散布図• 遺伝子Aの発現量が、10のとき、遺伝子B• 遺伝子Bの発現量が、10なら、• 散布図に表すと、(x, y) = (10, 10)20• 同様に遺伝子Aの発現量が、20のとき、遺伝子Bの発現量が、20なら、 (x,y) = (20, 20)

1001020

遺伝子A
cell innovator

遺伝子の相関関係 (1)• つまり、遺伝子Aの発現量が低いとき、遺伝子Bの発現量も低い。

遺伝子B• また、遺伝子Aの発現量が高いとき、遺伝子Bの発現量も高い。
20• 遺伝子AとBの発現量には、正の相関が見られる。
1001020遺伝子AA

B
cell innovator

遺伝子の相関関係 (2)• その逆なら、遺伝子Aの発現量が低いとき、遺伝子Bの発現量は高い。

遺伝子B• また、遺伝子Aの発現量が高いとき、遺伝子Bの発現量は低い。
20• 遺伝子AとBの発現量には、負の相関が見られる。
1001020遺伝子AA

Bcell innovator

20遺伝子A

遺伝子Bどちらが上流？

101002010遺伝子A

2001020遺伝子B• X軸とY軸を入れ替えても同じなので、どちらが上流か分からない？？
cell innovatorデータを増やしていくと見えてくるもの???
• 上記は、400サンプル＝400個の点における関係を見たもの。
• サンプル数を増やしていくと、「関係の度合い」（＝確率）も見えそう。

遺伝子発現にも統計学的なアプローチを。
cell innovatorベイジアンネットワーク（モデル）• 遺伝子Aが、ある確率で、遺伝子Bを制御していて、• 遺伝子Bが、ある確率で、遺伝子Cを制御している。

ABC

条件付き確率で表されたネットワークが書ける。
cell innovatorベイジアンネットワーク（モデル）• ベイジアンネットワーク＝条件付き確率で表されたネットワークのうち、ループ構造がないもの。

ABC

ループは、なし

ABCどちらか？

ABC
cell innovatorベイジアンネットワーク（モデル）• Aが起こってから、Bが起こり、Cになるのか？
• Bが起こってから、AとCが起こるのか？
• 言い換えると、Aが原因なのか、Bが原因なのか？
• どちらのモデルか分かれば、どちらが原因か分かる。（因果推定）

ABC

原因はどちら？

AB

C
cell innovator例えば、雨とスプリンクラーと芝生の関係は？
• A: 雨が降る（降雨量）。
• B: スプリンクラーが作動する。

A 70%• C: 芝生が濡れる。
• 芝生が濡れるのは、雨が降ったか、または、スプリンクラーが作動したから。

CB 30%

cell innovator

濡れた芝生の面積

芝生が濡れたら、雨が降る？

A?
?

CB
降雨量
• 雨が降ったから、芝生が濡れたのか？ A --> C• 芝生が濡れたから、雨が降ったのか？ C --> Acell innovatorスプリンクラーの影響を考慮• もし、芝生が濡れたから、雨が降ったのなら、B --> C --> A• つまり、スプリンクラーが作動すると、雨に何らかの影響があることになる。

• これは調べれば分かる。スプリンクラーが作動しても、天気に影響はない。

A?

CB

すべてのパターンを調べれば、どちらのモデルが適切か分かる！
cell innovator実際は、、、• Bスプラインによるノンパラ回帰• DAG 探索問題
• Greedy Hill Climbing アルゴリズム• BNRC スコア、オーバーフィッティング• 、、、、（詳細は玉田さんの資料をご覧ください）

イメージ的には、とにかく総当たりで、すべてのネットワークのパターンをチェックして、もっともらしいネットワークの状態を推定
cell innovator補論
簡単な数式で紹介する
ベイジアンネットワーク

cell innovatorこの補論で取り扱うベイジアンネットワーク

ベイジアンネットワークとは
１．大量の遺伝子発現解析データから求められる同時確率分布 P(G1,G2,…,GN) を用いて、２．遺伝子 G1, G2, …, GN 間の統計的依存関係を見出し，
３．非循環型の有向グラフでネットワーク構造を表現する方法
注１） ここでは、個々の実験で得られるN個の遺伝子の発現量 G1, G2, …, GN を確率変数とみなし，細胞の状態をそれら確率変数の同時確率分布 P(G1, G2, …, GN) で表現している。
注２） 非循環型とは、「フィードバックループなどの循環構造がネットワークの中には存在しない」という制限がこの方法にはあらかじめ課せられていることを示している。
ただし、この制限はダイナミック・ベイジアンネットワークでは取り除かれる。

cell innovator条件付き確率と有向グラフ
ベイジアンネットワークでは１．遺伝子間（例：G1 と G2 ）の条件付確率 P(G1|G2) がグラフに付与され，２．遺伝子 G2 から遺伝子 G1 への統計的依存関係をG2 → G1と方向性を持って，すなわち因果関係として表現する。

P(G1|G2) ⇔ G2 → G1
cell innovator条件付き確率での子遺伝子 G と親遺伝子 π(G)P(G1|G2) ⇔G 2 → G1

１．ここで G1 は子と呼ばれ，G2 は G1 の親 π(G1)と呼ばれ、π(G1)=(G2) とリスト表記される。
２．遺伝子 G1 の親遺伝子が複数，例えば G2 と G3 である場合、π(G1) は π(G1) = (G2, G3) とリスト表記される。
３．なお，遺伝子 G1 に親がいない場合には便宜的に，
π(G1) = (ø) と空集合 ø でリスト表記される．

cell innovator

遺伝子ネットワークの表現

ベイジアンネットワークによる遺伝子ネットワークの表現は、同時確率分布 P(G1, G2, …, GN) を条件付確率 P(Gi|π(Gi)) の積として展開することと同等となる。

cell innovator

仮想的な遺伝子発現の実験結果 D
ケース G1（実験）123456789102212112121G2 G312121221211222122121注）1: 発現していない 2: 発現している
cell innovator

仮想的な遺伝子発現の実験結果 D
マイクロアレイのデータで考えると, 下の表のようなイメージ.
（10サンプルぶんの実験結果のうち, 3つの遺伝子だけを見ている.）

exp.1

exp.2

exp.3

exp.4

exp.5

exp.6

exp.7

exp.8

exp.9exp.10

G12212112121G21212122121G31222122121注）1: 発現していない 2: 発現している
cell innovator遺伝子が３つの場合のグラフ表現
遺伝子が３つの場合には、全ての有向グラフ、つまり考えられる全ての統計的因果関係を数えきれる。

グラフ番号

構造Ni

同時確率分布P(Ni)=P(G1,G2,G3)の表現

G1

P(G1,G2,G3) = P(G1)P(G2)P(G3)

1

G2

G3

G1

2-7

P(G1,G2,G3) = P(G3|G2)P(G2)P(G1)
G2

G3

8-13

G1

G2

G3

P(G1,G2,G3) = P(G3|G2)P(G2|G1)P(G1)

14-16

G1

G2

G3

P(G1,G2,G3) = P(G3|G1)P(G2|G1)P(G1)

G1

G2

G3

P(G1,G2,G3) = P(G2|G1,G3)P(G3)P(G1)

17-19
20-25

G1
G2

G3

P(G1,G2,G3) = P(G3|G1,G2)P(G2|G1)P(G1)cell innovatorモデル選択の評価関数「データD が与えられたとき、遺伝子が３つの場合には25個のうちどのネットワーク構造が最も確からしい構造として推定されるか」＝「データD からのモデル選択の問題」＝「一般的にはある評価関数を設定し、その値が最大値（あるいは最小値）をとるネットワーク構造 Ni が選択される。」⇒
例えば、「データDが与えられたときネットワーク構造 Ni が実現する事後確率，すなわち条件付同時確率 P(Ni|D) を評価関数に用いる。」

cell innovator

Bayes の定理から事後確率 P(Ni|D) を求める

注） ネットワーク構造 Ni の出現確率 P(Ni) は、データ D 以外の情報から先験的に推測される事前確率で，もしなんらの付加情報がない場合には P(Ni) = P(Nj) となる。

cell innovator

実験データ D から P(D|Ni) を求める
Cooperらは以下の４つの条件
条件１：遺伝子の発現量は離散値（例：1=発現していない、2=発現している）で表現される。
条件２：ベイジアンネットワークのモデルが与えられると，実験データDの各行は相互に独立して現れる。
条件３：実験データDに、欠損値データはない。
条件４：実験データを得る前は、ネットワーク構造 Ni にあらかじめ付与する情報について完全に無知である。
が仮定できる場合には，条件付同時確率 P(D|Ni) が次式で与えられることを示した。

ただし、
Cooper,G.F. et al.(1992) A Bayesian method for the induction of probabilistic networks from data, Mach.Learn., 9, 309

cell innovator

Cooperらの式の説明

１） n は遺伝子の数で、今回の場合、n = 3。
２） qi は遺伝子 Gi の親 π(Gi) が取りうる状態の数。
今回の場合、親が１つの遺伝子 Gm の場合，すなわちπ(Gi) = (Gm) の場合には qi = 2。
親が２つ遺伝子 Gm と Gn の場合，すなわち π(Gi) = (Gm, Gn) の場合には qi =4。
３） ri は遺伝子 Gi が取りうる状態の数で，今回の場合、ri = 2．
４） Nijk は，遺伝子 Gi が「k =1:発現していない」あるいは「k=2:発現している」のどちらかの値をとり，Giの親 π(Gi) が1 ≦ j ≦ qi の j 番目の状態を取っている数（ケースの数）。

cell innovator

実験データDから事後確率 P(Ni|D) を計算する
P(N8)=P(G3|G2)P(G2|G1)P(G1)
G1

G2

(なし)

遺伝子G1の状態
1:非発現 2:発現

N111=5

G3

G1

G2

G3

ネットワーク構造N14：

ネットワーク構造N8：π(G1)P(N8)=P(G3|G1)P(G2|G1)P(G1)N112=5計π(G1)

遺伝子G1の状態
1:非発現 2:発現N11=10(なし)N111=5計π(G2)

遺伝子G2の状態
1:非発現 2:発現計

N21=5N22=5(G1=1)(G1=2)

N211=4
N221=1

N21=5N22=5計

N31=5N32=5

N112=5π(G2)

遺伝子G2の状態
1:非発現 2:発現(G1=1)(G1=2)

N211=4
N221=1π(G3)

遺伝子G3の状態
1:非発現 2:発現計

π(G3)

遺伝子G3の状態
1:非発現 2:発現(G2=1)(G2=2)

N311=4
N321=0

N31=5
N32=5(G1=1)(G1=2)

N311=3
N321=1

N212=1
N222=4

N312=1
N322=5N212=1
N222=4

N312=2
N322=4計

N11=10

Cooperらの式からP(D|N8) = 2.23×10-9

P(D|N14) = 2.23×10-10

もし P(Ni) = P(Nj) とすると、Bayesの定理から
P(N8|D) = 0.109
P(N14|D) = 0.011
ただし、最大の事後確率を与えるネットワーク構造は
N13 : G3 → G2 → G1 で、P(N13|D) = 0.112

cell innovator

3. 遺伝子ネットワーク

cell innovator

遺伝子ネットワーク• 遺伝子発現レベルのデータから推定されたベイジアンネットワークが、遺伝子ネットワーク。
• ただ、相関係数を調べて、線で結んだわけではない。
• 矢印（エッジ）には方向がある。

cell innovator

遺伝子ネットワークの意味するもの• 遺伝子ネットワークは、いわゆる「パスウェイ」ではない。
• いわゆる「パスウェイ」は、下記の情報のいずれか。
• タンパク間相互作用 = Protein-Protein Interaction (PPI) network。
• 遺伝子発現制御 = 転写因子と、その転写制御領域を持つ遺伝子の関係。
• 共発現 = ともに発現している遺伝子の関係。
• 文献情報 = 文献に、「制御関係あり」と報告された関係。

• 遺伝子ネットワークは、パスウェイとは異なる、新たな相互作用の情報。

cell innovatorパスウェイ解析と遺伝子ネットワーク解析の違い？

• パスウェイ解析は、「どの遺伝子が増加、減少した遺伝子した」のか、既知の情報をもとに結果を表示するもの。
• 遺伝子ネットワーク解析は、「どの遺伝子の影響が強い」のか、原因を予想するもの。また、未知の情報を含む。

cell innovator

遺伝子ネットワークの利点と欠点• 利点• 純粋にマイクロアレイデータのみから推定できるため、文献情報や、配列情報などのアノテーション情報を必要としない。（データドリブン）• lincRNAなど、機能が不明な遺伝子であっても、制御関係を推定できる。
• これまでに未知の制御関係を発見できる可能性がある。

• 欠点• 数十から数百個のマイクロアレイデータが必要。＝高いコスト• 高レベルの計算機環境が必要。（スーパーコンピューターなど）

cell innovator現在では、推定時の問題を回避可能• NCBI の Gene Expression Omnibus (GEO) に公開されているマイクロアレイデータを用いて推定を行う。 --> 高コストの問題を回避。
• 例えば、 Cancer Cell Line Encyclopedia (CCLE) には、およそ 1000 サンプル分のマイクロアレイデータが公開されている。[GSE36133]• 計算には、「京（SCLS）」などのスーパーコンピューターを利用。 --> 計算機環境の問題をクリア。

cell innovator

遺伝子ネットワークのグラフ論的な解釈• 数学的には、丸を「ノード」、矢印を「エッジ」と呼ぶ。

ノード＝遺伝子• エッジの始点になるノードが「親」

エッジ＝制御関係• エッジの終点になるノードが「子」• ネットワークの構造としては、一部の親に多数の子が集中するという構造になることが多い。（スケールフリー）• 特に「子が多いノード」は、「ハブ」と呼ばれる。

親子

cell innovator

遺伝子ネットワークの利用方法• 「ハブ」を探す＝ネットワーク中で影響力の強い遺伝子を見つける。（ハブの発現レベルが変化す
ると、子の発現レベルが変化するはず。）• 遺伝子ネットワークのノードを、logFCなどで色づけ。（パスウェイと同様、マイクロアレイデータの解析に利用。）• 上流解析：発現変動遺伝子を制御するのは、どの遺伝子か？（原因はどれか？）

cell innovator

解析事例（ハブをノックダウン）

Affara, M., Dunmore, B., Savoie, C., Imoto, S., Tamada, Y., Araki, H., Charnock-Jones, D. S., et al. (2007). Understanding
endothelial cell apoptosis: what can the transcriptome, glycome and proteome reveal? Philosophical transactions of the Royal
Society of London Series B, Biological sciences, 362(1484), 1469–1487. doi:10.1098/rstb.2007.2129

cell innovatorよくある質問、疑問• エッジの何パーセントが当たっているのか？エッジの何パーセントが既知で、何パーセントが未知の情報か？
• シグナル伝達系の活性化される順序は、分からないのか？
• レセプターが、リガンドを活性化しているように見えるが？
• 「ハブ」といっても、ただのキナーゼでは？転写因子でないから、転写は制御できな
いハズ。

データからはそう見える（バントしないほうがいい）といっているにすぎない。
cell innovator

Gene Network of apoptosis related genes

cell innovator

ANXA5 (top gene)

TP53(bottom gene)cell innovator

ANXA5

TP53 regulatinggenes

TP53

cell innovator• 統計学的に得られた結論は、感覚的には合わないかもしれませんが、参考にするのはどうでしょうか？

• 絶対、バントしてはダメだとか、言うつもりではありません。

• 昨年も、レッドソックスが優勝しましたね。。。。

cell innovator