HGC スパコン Web サービス利用法講習会
2011年10月25日 @ ヒトゲノム解析センター

大規模遺伝子ネットワーク
推定ソフトウェア SiGN
玉 田 嘉 紀tamada@is.s.u-tokyo.ac.jp
東京大学 大学院 情報理工学系研究科 コンピュータ科学専攻
理化学研究所 次世代計算科学研究開発プログラム
RIKEN Computational Science Research Program
Human Genome Center�
Institute of Medical Science, University of Tokyo	

本日の内容• 遺伝子ネットワーク推定ソフトウェア
SiGN（サイン）の紹介– SiGN-SSM, SiGN-BN, SiGN-L1• ベイジアンネットワークを用いた SiGNBN の機能紹介

Gene Network Estimation• Gene network: model for transcriptome level
gene–gene regulation using directed graphs.
• Gene network estimation is to estimate genenetworks from high-throughput biological datae.g. gene expression data.
Gene	Directed Graph:

Gene

Mathematical model consistingof nodes and directed edges
connecting them.

Regulation

KD1	

KD2	

KD3	 …	Gene	  1	 1.45	 -­‐1.54	 1.23	 …	
Gene	  2	 3.21	
…	…	-­‐2.1	

1.44	 …	

…	

…	

…	

Gene

Estimation
Gene Network

Gene Expression Data

3	

SiGN （サイン）• SiGN: A collection of large-scale gene
network estimation software designed forutilizing super computers.
– SiGN-BN: Bayesian networks

本日ベータ版公開（ベイジアンネットワーク）– SiGN-SSM: State space models（状態空間モデル）

オープンソース
公開中– SiGN-L1: L1-regularization based models 準備中（L1正則化）SiGN Web Site: http://sign.hgc.jp/

SiGN （サイン）• HGCスパコンで動くスパコン用専用ソフトウェア• 実行にはHGCのアカウントが必要• 実行はターミナルアプリから手でコマンド入力– 多少の Unix 操作の知識が必要

SiGNでできること（まとめ）• 遺伝子発現データから遺伝子間の発現の依存関係を予測・推定する– 遺伝子発現データ（マイクロアレイデータ）• 患者サンプルから得られる細胞• ノックダウン実験• 薬剤投与などの時系列に観測したデータ• ある程度のサンプル数が必要– 必要な量はモデルや実験デザインにより様々

SiGN （サイン）• SiGN: A collection of large-scale gene
network estimation software designed for
utilizing super computers.
– SiGN-BN: Bayesian networks

本日ベータ版公開（ベイジアンネットワーク）– SiGN-SSM: State space models（状態空間モデル）

オープンソース
公開中– SiGN-L1: L1-regularization based models 準備中（L1正則化）

SiGN Web Site: http://sign.hgc.jp/

SiGN-SSM
Dynamic gene network estimation software using a State Space Model (SSM).
Open source software distributed at http://sign.hgc.jp/signssm/（状態空間モデル）

Suitable for modeling time series gene expression data.

xn = Fxn−1 + vn , vn ~ N(0,Q)

Definition:[System model]

yn = Hxn + wn , wn ~ N(0, R) [Observation model]
System Model (= modules) xn = F xn–1 + vnx0(k)-vector

Fx1
H

Fx2

F...

HFxN
HParameter estimation =EM algorithm( x0 , F,Q, H, R and k )# dimensionsk << p# genesy1(p)-vector

y2

...

yN

Observation Model (= input data) yn = H xn + wn

Needs to repeat
estimation with various
initial values.

SiGN-SSM: 特徴• 時系列マイクロアレイデータからの遺伝子
ネットワーク推定• 不等間隔・繰り返し計測データに対応• SSMのパラメータ推定．
• 推定結果に基づきネットワーク構造の決定．
• オープンソース– 一部機能はHGCスパコン限定

SiGN-SSM: 詳細はウェブサイトへ
http://sign.hgc.jp/signssm/

SiGN （サイン）• SiGN: A collection of large-scale genenetwork estimation software designed forutilizing super computers.
– SiGN-BN: Bayesian networks本日ベータ版公開（ベイジアンネットワーク）– SiGN-SSM: State space models（状態空間モデル）

オープンソース
公開中– SiGN-L1: L1-regularization based models 準備中（L1正則化）SiGN Web Site: http://sign.hgc.jp/

Nonparametric Bayesian Network Model
We use the nonparametric Bayesian network as models for gene networks
Joint Probability by a DAG (Directed Acyclic Graph)
Node = Gene
X1	X2	f (X1, X 2 ,…, X 6 )= f1 (X1 ) f 2 (X 2 ) f 3 (X 3 | X1 ) f 6 (X 6 | X 3 , X 4 )Network score = Posterior ProbabilitynX5	i =1X4	

€G :  gene networkX :  expression data

Nonparametric regression by B-spline
X6	

Directed Edge =Regulatory
Relationships

child geneX3	π (G | X ) ∝ π (G ) ∫ ∏ f ( xi1 ,..., xip | θG )π (θG | λ)dθG

xij = m j1 ( pi(1j ) ) + ... + m jq j ( piq( jj) ) + ε ij
ε ij ~ N (0, σ 2j )m jk ( pik( j ) ) = ∑l =1jk γ lk blk( j ) ( pik( j ) )Mparent gene

Dynamic Bayesian Networks
DBN (Dynamic Bayesian Network): model for time course data.
DBN assumes the dependency between consecutive time points.
Time course gene expression data: X1, ..., XT
X1

X2

f (X2 | X1 )

...

X3

f (X3 | X2 )XTf (XT | XT −1 )Gene Network: Bipartite Graph	x1(T)	

x2(2)	...	

x2(T)	...	

...	

xp(1)	

xp(2)	

...	

...	

€

...	

x1(1)	

x1(2)	

x2(1)	€€...	

xp(T)	
13	Gene Network Estimation by
Dynamic Bayesian Networks
Estimated Bipartite Graph	x1(t-1)	

x1(t)	x2(t-1)	

x2(t)	x3(t-1)	

x3(t)	

x4(t-1)	

x4(t)	x5(t-1)	

x5(t)	

Corresponding Gene Network	Self-loop	x1	x2	x3	

Feedback
regulation	

x4	x5	
14	

Difficulty in Bayesian Network Estimation
A huge number of possible DAGs : impossible to search the optimal one#	  of	  nodes	#	  of	  DAGs	#	  of	  nodes	#	  of	  DAGs	

1	

1	16	

≈ 8.37 × 1046	

2	3	17	

≈ 6.26 × 1052	

3	25	

18	

≈ 9.93 × 1058	

4	543	

19	

≈ 3.32 × 1065	

5	29,281	

20	

≈ 2.34 × 1072	

6	3,781,503	

21	

≈ 3.46 × 1079	

7	1,138,779,265	

22	

≈ 1.07 × 1087	

8	783,702,329,343	

23	

≈ 6.97 × 1094	

9	1,213,442,454,842,881	

24	

≈ 9.43 × 10102	

10	

≈ 4.17 × 1018	

25	

≈ 1.86 × 10111	

11	

≈ 3.15 × 1022	

・・・	・・・	

12	

≈ 5.21 × 1026	

30	

≈ 2.71 × 10158	

13	

≈ 1.86 × 1031	

・・・	・・・	

14	

≈ 1.43 × 1036	

・・・	・・・	

15	

≈ 2.37 × 1041	

40	

≈ 1.12 × 10276	

Exceeds thenumber of atoms
in the universe

Network Size and Algorithms
Different search algorithms are developed depending on the size of networks#	  of	  Genes

2,000〜20,000
Neighbor Node Sampling & Repeat (NNSR) algorithm(Tamada et al., 2011)	
〜2,000
Greedy Hill-climbing Algorithm (HC) +Bootstrap (Imoto et al. 2002)	
〜500
Extended COS (Kojima et al. 2010)	
Constrained Optimal Search algorithm (COS) (Perrier et al. 2008)	

〜30

Parallel OS (Para-OS) (Tamada et al. 2011) 	
Optimal Search algorithm (OS) by Dynamic
Programming (Ott et al. 2004)	

Greedy Hill-Climbing Algorithm (HC)
Algorithm for searching the local optimal DAG structure	
Heuristics algorithm applicable to estimate gene networks for 〜 100 genes.
StartX3	

X5	

X1	

X2	

Loop I

X3	

X4	

X6	

X1	

X5	

X2	

Loop II

X3	

X4	

X6	

1. Begins with an empty graph.
2. Visits nodes in a random order.

X5	

X1	

X2	

Loop III

X4	

X3	

X6	

X5	

X1	

X2	

X4	

X6	

4. Employs the best operation that
improves the score.

3. Calculates local scores for all
possible candidate parents.

Add/Delete/Reverse

5. Repeats until any operation
can improve the score.

※ Need to check every time whether a cyclic path is made or not.

17	

※ Repeats this many times, then employs the best structure because they are local optimal.

HC + Bootstrap〜 1000 genes

Bootstrapping is required for calculating the reliability of edges.	
Networkestimation

X1	

Randomsampling withreplacement

X3	

X5	

Original expression
dataset-Re-sampleddataset

X2	

X4	

X6	

Estimated network

Estimate networks many times for re-sampled dataset. (1,000 times 〜)
The final structure is determined by the frequencies of edges during therepeated estimation.
We can perform each network estimation independently for the re-sampleddatasets in parallel.
-Parallelization is easy for Bootstrap HC.

Estimated Networks
Bootstrap probabilities

0.8	

X	

0.6	

0.5	

Y	• Edges in estimatednetworks have “Bootstrap
Probability.”

Direction reliabilities (Static model only)
X	
0.8	

0.5

0.3

Y	
Direction reliability = 0.625• Edges also have reliability ofthe direction.
– Dynamic model can haveedges in both directions.

Hubness & Hub Index
Hubness: simply counts the number of children regardless of the confidence of theedge direction.
Hub genes may have important roles as master regulators in the networks.
Useful for analyzing the estimated networks.	
Direc:on	  reliability	
0.8	

X	

Hub	  Index	  of	  X	  =	  1.9	

0.6	

0.5	

Hubness	  of	  X	  =	  3	Y	

Hubness	  of	  Y	  =	  0	Hub	  Index	  of	  Y	  =	  0.4	

Hub Index: takes the confidence of edge direction into account. Simply the sum of“BS.Direction” for all edges connected to the target node.	
Note: Hub Index does not take bootstrap probabilities into account.

実行・解析の流れ
1. 入力となるアレイデータの準備– テキストファイルの入力ファイル• EDF, タブ区切り独自形式

2. HGCスパコンへのジョブの投入–

SGE のアレイジョブによりBootstrapの繰り返し計算を並列実行

3. 出力結果のまとめ処理–

結果がファイルに出力されるので，それを１つのネットワークファイルにまとめる•CSML, タブ区切りテキストファイル

4. Cell Illustrator Online などによる解析–

ハブ遺伝子の分析やターゲット遺伝子関連サブネットワークの抽出・分析

SiGN-BN の公開状況• ベータ版 (rel. 0.9.0) が ~tamada/sign 以
下にインストールされている• 近日中に正式版およびドキュメントを公開予定– HC+Bootstrap 法以外の SiGN-BN のアルゴリズ
ムや SiGN-L1 なども準備が整い次第順次HGC
スパコンユーザ向けに公開予定

入力ファイル• EDF 形式によるアレイデータ• 独自タブ区切りテキストファイル– 詳細は後日公開するドキュメントを参照
タブ区切りテキストファイルによる独自フォーマット
Probe_1	 Gene_1	 0.5	

0.3	

0.52	

...	

Probe_2	 Gene_2	 1.5	

1.44	

1.55	

....	

Probe_3	 Gene_3	 ...	

...	

...	

....	

...	

Probe_p	 Gene_p	...	

実行
• HGCスパコンにログインし，SGEのqsubコマンドによりジョブを投入・実行
ブートストラップの実行
qsub -t 1-X [SGE options] /usr/local/bin/signbn-hc.sh--bs -o file_prefix [SiGN options] input_file※ X はブートストラップ回数．1000以上を推奨．
結果のまとめ処理
qsub [SGE options] /usr/local/bin/signbn-hc.sh --bin
signproc --bs prefix=file_prefix [,other options] --output
file=output_file,type=file_type※ SGE の -t オプションは付けない
file_type は CSML か TXT

Cell Illustrator Online による解析推定したネットワーク (CSML) は CIO で閲覧・編集・解析が可能

Edge Properties
The following information is generated and assigned to each estimated edge. You can
see them on CI Online.
BS.Prob

Bootstrap probability

edgeScore

Same as BS.Prob.

BS.Gain
BS.up/downAverage of the edge gain
One of “up”, “down”, or “unknown” that represents theestimated type of regulation.

BS.edgeTypeRatio

The ratio of up, down, and unknown regulation.

BS.Direction

Frequency (confidence) of the edge direction.

Edge properties shown in CI Online.

サンプルデータによる実習（ベータ版）
サンプルファイル sample003.txt が ~tamada/sign/samples にあります．
(1) 出力用ディレクトリの準備と移動
mkdir ~/test
cd test(2) SGEへのジョブの投入

赤字はベータ版でのみの指定

qsub -t 1-10 ~tamada/sign/signbn-hc.sh --dir ~tamada/sign --bs -oresult ~tamada/sign/samples/sample003.txt(3) SGEのジョブの確認
qstat(4) 出力ファイルのまとめ
qsub ~tamada/sign/signbn-hc.sh --dir ~tamada/sign --bin signproc --bs
prefix=result --output file=result.csml,type=csml~/test に result.csml が作られます．ローカルPCに転送しCIOで開いてください．

Acknowledgments
Computational time was provided by the Super
Computer System, Human Genome Center,
Institute of Medical Science, The University ofTokyo, and RIKEN Integrated Cluster of Clusters(RICC) system.

Members
Teppei Shimamura, Rui Yamaguchi, Masao Nagasaki, Seiya Imoto, Satoru Miyano

Collaborators
Cristin Print, Hiromitsu ArakiD. Steve Charnock-Jones
Osamu Hirose
Ryo Yoshida, Tomoyuki Higuchi

The University of Auckland, NZ
The University of Cambridge, UKThe University of Tokyo
The Institute of Statistical Mathematics