応用生命科学・情報生命学 第2回: RNA配列情報解析

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