バイオインフォマティクス理論 第9回: RNA配列情報解析

本日の内容
バイオインフォマティクス理論
RNA (RiboNucleic Acid) とは?
Gene
 RNA2次構造予測
Transcription
 RNAシュードノット構造予測
第9回: RNA配列情報解析
Gene
ncRNA
 RNA構造アラインメント
Gene
Genomic DNA sequence
ncRNA
Interaction
mRNA
Translation
Protein
11月24日(火)4時限目
加藤 有己
 Non‐coding RNA (ncRNA)
 Functional
 Non‐functional
京都大学 iPS細胞研究所
講義資料
http://tcs.cira.kyoto‐u.ac.jp/~ykato/lecture/bioinfo15‐2/
Correlation
Function
Structure
3
2
non‐coding RNA (ncRNA) 解析の重要性
RNAの構造
RNA配列情報解析
 ゲノム全体でncRNAに転
 1次構造(塩基配列)
 立体構造解析は高コストかつ低スループット
→コンピューターを用いて塩基配列から2次構造を予測
CAAAAGUCUGGGCUAAGCCCACUGAUGAGUCUCUGAAAUGAGACGAAACUUUUG
写される領域の割合は
高等生物になるほど高
い
 2次構造(塩基対の集合)
 モデル化が容易
 3次構造(立体構造)
 実験的解析が困難
CAAAAGUCUGGGCUA
AGCCCACUGAUGAGU
CUCUGAAAUGAGACG
AAACUUUUG
 ncRNA機能解析
→高等生物の複雑な生体
内制御機構の解明
 どんな複雑な構造も扱いたい
 塩基対が複雑に交差する
Watson‐Crick相補塩基対
 2つのRNAが結合する、etc.
A-U -: 水素結合
G-C
[Shabalina & Spiridonov 2004]
4
 ゲノム網羅的な解析のためになるべく速く、正確に
6
5
最も基本的なRNA2次構造
RNA2次構造予測
Nussinovのアルゴリズム
 2次構造
 入力: RNA塩基配列 x = x1x2…xn
 動的計画法(dynamic programming, DP)
 部分問題の解からより大きな問題の解を導き出す
 出力: x の2次構造 y
ステム
 部分問題の最適解の定義
ループ
 塩基配列(1次構造)
 S(i, j): 部分配列 x[i..j] 上の塩基対の個数の最大値
 ある目的関数 f の下での最適化問題と考える
(1  i  j  n; n は配列全体の長さ)
(f(y) が最大(最小)となる構造 y を計算)
 目的関数の例
 塩基対の個数
 構造の自由エネルギー
塩基配列上では塩基対が
入れ子になって出現
[Eddy04]
より転載
 入れ子型塩基対の個数が最大となる2次構造を計算
CAAAAGUCUGGGCUAAGCCCACUGAUGAGUCUCUGAAAUGAGACGAAACUUUUG
7
→Nussinovのアルゴリズム [Nussinov et al. 1978]
i..j 上で入れ子型塩基対を含む構造を作れるパターンはこの4種類のみ
8
9
DPによる S(i, j) の計算
DPの表(初期化)
DPの表(反復)
 初期化
 1個の塩基は塩基対を作らない
 すでに埋まっているマス目の値から隣接するマス目
 S(i, i) = 0
1個の塩基は塩基対を作らない
 S(i, i+1) = 0
隣接する塩基は塩基対を作らない
1
G G G A A A U C C
 反復(漸化式)
左塩基生成
右塩基生成
(i と j が塩基対)
の値を漸化式に従って計算
 隣接する塩基は塩基対を作らない
j
塩基対生成
分岐
G
G
G
A
A
i A
U
C
C
0
0
0
1
0
0
2
0
0
3
0
0
0
0
0
4
S(i, i) = 0
S(i, i+1) = 0
0
5
6
0
0
2
3
4
5
6
7
8
9
j
G G G A A A U C C
7
0
8
i
0
9
G
G
G
A
A
A
U
C
C
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
11
10
12
DPの表(反復)
DPの表(終了)
DPの表(トレースバック)
 すでに埋まっているマス目の値から隣接するマス目
 DPの最終結果は S(1, n):
 実際の2次構造は S(1, n) からトレースバックして計算
の値を漸化式に従って計算
配列全体に対する塩基対の個数の最大値
j
G
G
G
A
A
i A
U
C
C
0
j
j
G G G A A A U C C
0
0
0
0
0
0
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
G G G A A A U C C
G G G A A A U C C
G
G
G
A
A
i A
U
C
C
0
0
0
0
0
0
0
1
2
3
0
0
0
0
0
1
2
3
0
0
0
0
1
2
2
0
0
0
1
1
1
0
0
1
1
1
0
0
0
0
0
0
0
0
0
G
G
G
A
A
i A
U
C
C
0
13
0
0
0
0
0
0
1
2
3
0
0
0
0
0
1
2
3
予測2次構造
0
0
0
0
1
2
2
0
0
0
1
1
1
GGGAAAUCC
.((.(.)))
0
1
1
1
0
0
0
0
0
0
0
0
0
0
ドット・ブラケット表現
(Vienna形式):
対応する括弧で
塩基対を表現
0
20
19
Nussinovのアルゴリズムの計算量の評価
RNA2次構造予測(再掲)
自由エネルギー最小化法の限界
 S(i, j) を計算するために必要な時間と領域
 入力: RNA塩基配列 x = x1x2…xn
 現在の利用可能なエネルギーモデルは不完全なため予
測精度は限定的
 出力: x の2次構造
[Mathews06]
が多い
 ある目的関数の下での
最適化問題と考える
 目的関数を構成する要素の例
 スタック塩基対・各種ループ構造のエネルギー
 自由エネルギーが最小となる2次構造を求める
←動的計画法(DP)(O(n3) 時間)
 時間計算量
j の動く範囲: 1  i  j  n
の動く範囲: i < k < j
⇒ O(n3)
 i,
 漸化式の分岐における k
 領域計算量
j の動く範囲: 1  i  j  n
⇒ O(n2)
 i,
The energy landscape of
structures for Agrobacterium
tumefaciens 5S rRNA
[Ding et al. 2005]
 mfold [Zuker & Stiegler 1981]
 RNAfold [Hofacker et al. 1994]
21
 最小自由エネルギー構造が正解の構造に似ていないこと
22
23
最近のRNA2次構造予測法の傾向
分配関数
塩基対確率
 RNA構造形成に対して統計力学の考え方を適用する
 可能な2次構造の集合(アンサンブル) S(x) 上で、ある2次
 アンサンブル S(x) 上で、塩基 xi と塩基 xj が塩
構造 y の事後形成確率 P(y | x) はBoltzmann因子に比例
 P(y | x)  e-βF(y)
 β: (温度×Boltzmann定数)の逆数
 F(y): y の自由エネルギー
 可能な2次構造の集合(アンサンブル)を考慮した
分配関数、塩基対確率を利用し、アンサンブルの
重心などを求める
 シュードノットを予測できるようにする
 複数のRNAのファミリーにおいてはシュードノットを
構成する塩基対が少なからず(0.2~14.4%)存在
[Mathews et al. 1999]
基対を形成する確率(塩基対確率)
 分配関数のDP表を利用して、O(n3) 時間の動的計
 絶対確率を得るためにBoltzmann因子を分配関数 Z で
画法で計算可能 [McCaskill 1990]
正規化
 McCaskillのアルゴリズムはRNAfold [Hofacker et al. 
 相同配列の共通2次構造を予測する
1994] で実装されている
 配列が長くなれば可能な2次構造の数は指数的に増加す
 マルチプルアラインメントを利用
るため、Z は直接計算できない
→O(n3) 時間の動的計画法(DP) [McCaskill 1990]
24
25
期待精度最大化法
本日の内容
 正しい予測に対するゲイン関数
 RNA2次構造予測
: 正解2次構造
Crossing
A
5’—UC
 RNAシュードノット構造予測
Connect base
AGACG
pairs with arcs
A
GC—3’
A
5’—UCAGCAGAAAGC—3’
-
 RNA構造アラインメント
: 予測2次構造
正解塩基対をなるべく多く予測したい
RNA2次構造: シュードノット
-
真の塩基対の個数 真の非塩基対の個数
26
 アンサンブル S(x) 上の確率分布 P(y | x) が与えられ
H型シュードノット
たとき、期待精度を最大化する2次構造 を計算
 シュードノットはRNAの機能において様々な役割を果たす
 翻訳・スプライシングの調節、など
 立体構造形成を補助
 CentroidFold [Hamada et al. 2009]
 塩基対確率を用いたNussinov型DPを用いる
→精緻な構造解析ではシュードノットを考慮すべき
28
27
シュードノットを考慮したRNA2次構造予測
シュードノット付き2次構造予測の既存手法
 入力: RNA配列 x = x1x2…xn
Approach
 出力: シュードノットを含む x の2次構造
 (注) y は2値の上三角行列
DP
yij = 1
5’
x
Input
3’
5’
i
Restricted
Heuristic
j
3’
Class of pseudoknots†
Wide
Time complexity
O(n4)
~
O(n6)
~O(n3)
Optimality (MFE‡)
Available software
Yes
PKNOTS, NUPACK, pknotsRG
N/A
ILM,
HotKnots, FlexStem, ProbKnot
x
3’
Input
Approximate probability distribution
P(y | x)
base‐pairing probabilities
Maximum expected
accuracy structure
Maximize expected accuracy
[Akutsu 2000; Lyngsø & Pedersen 2000]
 ‡MFE: Minimum Free Energy
ある目的関数の下での最適化問題(エネルギー最小化問題など)
IPknot: Integer Programming‐based prediction of RNA pseudoKNOTs
5’
 †任意のシュードノットを予測する問題はNP困難
Output
29
Integer programming
5’
3’
Output
IPknot
高速かつ高精度の予測法が求められている
[Sato et al., 2011]
30
31
32
確率分布の近似
整数計画法(Integer Programming, IP)
比較手法
 シュードノット構造上で定義された確率分布を因数分
 整数計画問題
 シュードノットを考慮したRNA2次構造予測
 ProbKnot [Bellaousov&Mathews, 2010]
 FlexStem [Chen et al., 2008]
 HotKnots [Ren et al., 2005; Andronescu et al., 2010]
 pknotsRG [Reeder&Giegerich, 2004]
 ILM [Ruan et al., 2004]
解により近似
 整数値を取る変数に関する線形の不等式制約の下
で、線形の目的関数を最大化(最小化)する問題
Time‐consuming to compute
Level 1
Pseudoknot
free
Level 2
Pseudoknot
free
Level 3
Pseudoknot
free
Pseudoknotted
 シュードノットを考慮しないRNA2次構造予測
 CentroidFold [Hamada et al., 2009]
 RNAfold [Hofacker et al., 1994]
 主なアルゴリズムは分枝限定法と分枝カット法
 実際にはソルバー(CPLEX, GLPK)を利用して解く
33
予測精度評価
35
34
予測精度(388本のRNA配列上)
実行時間
10000
SEN =
TP
TP + FN
 Positive predictive value
PPV =
TP
TP + FP
TP: # of correctly predicted base pairs
1000
IPknot
100
FN: # of true base pairs that were not predicted
ProbKnot
FlexStem
Time (s)
 Sensitivity
 refinement:
10
FP: # of incorrectly predicted base pairs
HotKnots
pknotsRG
refine base‐
pairing probabilities defined over pseudoknot‐free structures
ILM
1
0
200
1000
1200
NUPACK
Sequence length (nt)
38
DAFS: Dual decomposition for Aligning & Folding RNA sequences Simultaneously
比較に基づくRNA構造予測
 RNA2次構造予測
 複数配列上での比較解析は1配列解析より予測精度
の点で優れている
5’
 RNA構造アラインメント
800
37
本日の内容
 RNAシュードノット構造予測
600
0.1
0.01
36
400
良いアラインメント
5’
良い構造
a
b
3’
3’
Approximate probability distribution P(θ | a, b)
Input
Maximum expected
accuracy structural
alignment
 RNA配列のアラインメントと構造予測を同時に行う
 Sankoffのアルゴリズム(1985)
 O(L3N) 時間 & O(L2N) 領域
 L: アラインメントの長さ
Maximize expected accuracy
時間がか
かりすぎる
 N: 配列の本数
Posterior base‐pairing/
alignment‐matching probabilities
Integer programming
5’
5’
3’
3’
Output
DAFS
[Sato et al., 2012]
39
40
41
確率分布の近似
期待精度最大化
期待精度最大化(続き)
 構造アラインメント上の確率分布を因数分解により近似
 期待精度を最大にする構造アラインメント を求める
 目的関数
maximize
Time‐consuming to compute
: 構造アラインメント空間
: 参照構造アラインメント
: 予測構造アラインメント
 α, τ, σ: バランスパラメータ
真の予測に対するゲイン関数
# of true positives
 塩基対・マッチング確率
 pij(a): 塩基 ai が塩基 aj と対になる
確率
 pik(a,b):塩基 ai が塩基 bk とマッチする
確率
# of true negatives
なるべく多くの正しい塩基対・マッチングを予測したい
42
制約条件
高々1塩基とのみ塩基対を構成可能
43
44
制約条件(続き)
双対分解
 マッチング塩基対
 マッチング塩基対に関する制約式を目的関数に移動
シュードノットはない
Nussinov型DPで解ける
Needleman−Wunsch型
DPで解ける
高々1塩基とのみマッチする
塩基対は保存されなければ
ならない
交差マッチはない
問題の複雑化
45
 λ, μ, ν: Lagrange乗数
46
47
双対分解(続き)
Lagrange乗数の決定
マルチプルアラインメントへの拡張と計算機実験
 マッチング塩基対に関する制約式を目的関数に移動
 λ, μ, ν に関して L(λ, μ, ν) を最小化
 (注) L は微分不可能
 平均確率を用いた累進法によるマルチプルアラインメ
ントへの拡張:
 平均塩基対確率
 平均アラインメントマッチング確率
 劣勾配法
シュードノット予測に対
してIPknotで解ける
 初期化: λ = 0, μ = 0, ν = 0
アラインメントをアラインする
 反復 (T回):
 Rfam 11.0 [Burge+2013] からのデータを用いた実験
 現在の λ, μ, ν を用いて x, y, z, w を計算
Needleman−Wunsch型
DPで解ける
 x, y, z が移動した制約式を満たす
 691個のシュードノットを含まないアラインメント
⇒ 計算終了
(PKfree dataset)
 82個のシュードノットを含むアラインメント
 そうでなければ
 λ, μ, ν: Lagrange乗数
⇒ 劣勾配を用いて λ, μ, ν を更新
48
(PK dataset)
49
50
評価尺度
PKfreeデータセット上での実行時間
PKfreeデータセット上での予測精度
 SPS: sum‐of‐pairs score
 SCI: structure conservation index
Best MCC
0.9
アラインメント評価用
70000
0.8
60000
0.7
50000
0.6
 SEN: sensitivity
 PPV: positive predictive value
2次構造評価用
40000
0.5
SPS
0.4
 MCC: Matthews correlation coefficient
Time (s)
SCI
MCC
30000
0.3
20000
Very fast
0.2
10000
0.1
0
0
DAFS
CentroidAlign
RAF
LARA
LocARNA
DAFS
ProbConsRNA
51
PKデータセット上での予測精度
0.9
CentroidAlign
RAF
LARA
LocARNA
ProbConsRNA
52
53
まとめ
PKデータセット上での実行時間
180,000
 RNA2次構造予測は機能推定のための鍵
0.8
160,000
 RNA2次構造予測
0.7
140,000
 Nussinov型DPはRNA2次構造予測法の基盤
0.6
120,000
 主流は自由エネルギー最小化法からアンサンブル
Best MCC
0.5
SPS
SCI
0.4
MCC
100,000
80,000
0.3
60,000
0.2
40,000
0.1
20,000
0
0
DAFS (IPknot decoding)
DAFS (Nussinov decoding)
RNASampler
Time (s)
Very fast
 分割した問題を独立に解き結果を重ね合わせる手
法は速度のみならず精度の観点からも有効
DAFS (IPknot decoding)
DAFS (Nussinov decoding)
RNASampler
54
55
参考文献
レポート課題
 阿久津達也、バイオインフォマティクスの数理とアル
 RNA配列解析について、情報科学技術の有用
ゴリズム、共立出版、2007年
性をまとめ、構造解析の今後の展開を自由に
発想し、A4用紙1枚以内で論じなさい
 Eddy, S., Nat. Biotechnol., 22, 1457−1458, 2004
 Mathews, D.H. and Turner, D.H., Curr. Opin. Struct.  提出期限
 2015年12月7日(月)
Biol., 16, 270–278, 2006
 加藤有己、二次構造に基づくRNA配列解析ソフトウェ
アの進展(有田正規編、使えるデータベース・ウェブツ
ール第4章第8節)、羊土社、2011年
 提出先(電子メール)
 kyomu (at) bic.kyoto‐u.ac.jp
 講義資料
http://tcs.cira.kyoto‐u.ac.jp/~ykato/lecture/bioinfo15‐2/
57
を考慮する手法へ
 複雑なRNA2次構造解析
 RNAシュードノット構造予測
 RNA構造アラインメント
58
56