Traitement en cours

Veuillez attendre...

Paramétrages

Paramétrages

Aller à Demande

1. WO2011062166 - PROCÉDÉ D'ESTIMATION D'ÉTAT DYNAMIQUE DE STRUCTURE SECONDAIRE D'ARN

Document

明 細 書

発明の名称 RNAの二次構造の動態を予測する方法

技術分野

0001   0002   0003   0004  

先行技術文献

特許文献

0005  

非特許文献

0006  

発明の概要

発明が解決しようとする課題

0007   0008  

課題を解決するための手段

0009   0010  

発明の効果

0011  

図面の簡単な説明

0012   0013   0014   0015   0016   0017   0018   0019   0020   0021   0022   0023   0024   0025   0026   0027   0028   0029   0030   0031   0032   0033   0034   0035   0036   0037   0038   0039   0040   0041   0042   0043   0044   0045   0046   0047   0048   0049   0050   0051   0052   0053   0054   0055   0056   0057   0058   0059   0060   0061   0062   0063   0064   0065   0066   0067   0068   0069   0070   0071   0072   0073   0074   0075   0076   0077   0078   0079   0080   0081   0082   0083   0084   0085   0086   0087   0088   0089   0090   0091   0092   0093   0094   0095   0096  

実施例

0097   0098   0099   0100   0101   0102   0103   0104   0105   0106   0107   0108   0109  

産業上の利用可能性

0110  

請求の範囲

1   2   3   4   5   6   7   8   9   10   11   12  

図面

1   2   3   4   5   6   7   8   9   10   11  

明 細 書

発明の名称 : RNAの二次構造の動態を予測する方法

技術分野

[0001]
 本発明は、RNAの二次構造が平衡状態となる前の時点、すなわち、該二次構造がマクスウェル・ボルツマンの分布則に従った存在確率の分布を示すよりも前の時点における、該RNAの二次構造の動態を予測する方法、ならびにそのためのプログラムおよび装置に関する。
[0002]
(発明の背景)
 近年、アンチセンスオリゴヌクレオチドやsiRNAを用いて、mRNAが創薬のターゲットになり得るという認識が広がり、ターゲット分子としてのmRNAについて構造生物学的な研究が徐々に注目を集めつつある。
 構造生物学は、生体内の分子の特定の構造とその構造に起因する機能を研究することで発展してきた。この分野では、タンパク質のみならず、機能を持つRNAであるtRNAやrRNAの構造について、従来からよく研究が進んでいる。しかしながら、遺伝子であるDNAの遺伝情報をタンパク質に翻訳するための「伝令」としての役割としての認識が強かったmRNAについては、これまでその高次構造についてはほとんど注目されてこなかった。これは、mRNAが遺伝情報を伝達する働き以外に、注目すべき機能を持たないと考えられていたからである(機能を持たないのであれば、それを引き起こす構造もないであろうと漠然と推察されていた)。
 mRNAは、細胞内のイベントに応じてその全体として形態(一次・二次・三次構造)を大きく変える。核内でmRNAが合成されたとき、スプライシングを受けているとき、核外に運搬されたとき、タンパク質の翻訳が行われているとき、さらにその後など、mRNAの全体的な構造は一意には論じることができない(本明細書では、核外のmRNAの挙動について述べるので、以後mRNAの一次構造の変化については注目しない)。
 一方RNAは自己相互作用しやすい高分子であるため、溶液中で何らかの高次構造を取り得ると推定される。しかし、長いRNAは最安定構造と準最安定構造群とがエネルギー的に差異が少ないので、結果として溶液中の長いRNAは、複数の高次構造の混合物として存在しているはずである(すなわち、細胞質内において一意の構造をとることが保証されていない)。
 こうした二つの要因から、mRNAの全体の高次構造を研究するための方法が存在せず、ひいてはこのためmRNA構造が生命現象に影響を与えていると意識されることが少なかったのである。
[0003]
 ところが、昨今の研究から、mRNA上の局所的な二次構造が、RNAの編集、riboswitch、IRES、フレームシフトなどの、分子生物学的イベントに関与していることが次第にわかってきた。RNAiやアンチセンス医薬品の結合部位による効率の差、miRNAの結合部位にも、そうした局所的な二次構造の関与が示唆され得る。つまり、全体の構造を把握不可能なmRNAの高次構造であっても、対象をある細胞内イベントに対応したmRNAに限定すると、その局所的高次構造については生物学的に機能があることがわかってきた。このように機能を持つ構造体に対しては、医薬品の開発ができることが示唆され、事実、例えば、ピリチアミンはTPP riboswitchに対して作用することで抗生物質として働いていることが近年分かってきた。但し、その発見と解析の難しさから、このように「構造を持つ・機能を持つ」ことが報告されているmRNAはごく少数であり、したがって創薬ターゲットとなり得るmRNAの二次構造は非常に限定されていた。
[0004]
 本発明者は以前、〔自己相互作用しやすいというRNA分子の特性に基づけば、あらゆるmRNAは確率的に何らかの高次構造をとることがあり、それが生体内で特別な機能を持たないかまたは機能を発揮し得るほど安定でないがために(或いは単純にその機能が未知であるがために)、これまで構造として認識されてこなかっただけである〕との発想の下で、mRNA分子内の局所的な二次構造を、コンピュータを用いてcDNAデータベースから網羅的に検索し得るソフトウェアを開発した(特許文献1)。さらに、該検索により抽出された局所的な二次構造が、生体内で比較的安定に存在し得る(すなわち、ある物質(例えば、低分子化合物)との特異的結合により容易に安定化される)ことを、インビトロの実験系を用いて実証した。
 本発明者はさらに、対象とするRNAが、抽出された局所的な二次構造を、実際に生体内でとるかどうかを精度よく予測し得る方法を開発した(特許文献2)。該方法は、mRNAのような長い配列の対象区間上のある位置に、ある短い長さのフレームを想定し、そのフレームで規定される小区間内の各特定位置において或る局所的な二次構造が何程の蓋然性で存在するか(第1の蓋然性)を求める、という操作を、フレームを微小長さ(フレーム長よりも十分短い長さ)だけ移動させながら各小区間で行った後、各小区間で判明した特定位置が、全ての小区間の間でどれだけ互いに関連しているかの指標、すなわち各特定位置自体がどの程度の信憑性を持つかを示す指標(第2の蓋然性)を求める、というものである。第1の蓋然性が高く、かつ、第2の蓋然性が高い特定位置は、その特定位置自体が該局所的な二次構造の位置である信憑性が高く、しかも、そこに該局所的な二次構造がより確かに存在するということになる。従って、該方法によって、目的とする局所的な二次構造が有意義な確実性で存在するかどうかを定量的に予測することが可能となった。
 該方法は、mfold(GCG Software;非特許文献1参照)などの従来の構造予測ソフトウェアが好ましい予測結果を与えない比較的長い区間(例、200塩基以上)に対して、対象区間がどれだけ長くても信頼性の高い予測結果を与えることができる。該方法によって、例えば、アンチセンス核酸、RNAi、低分子化合物などの、二次構造部分もしくはその近傍に結合して該構造を安定化することにより、mRNAからタンパク質への翻訳を抑制する(あるいはmRNAの分解安定性を改善して翻訳効率を増強する)化合物の開発・設計をより効率よく行うことが可能となった。

先行技術文献

特許文献

[0005]
特許文献1 : 国際公開2006/054788号パンフレット
特許文献2 : 国際公開2008/072713号パンフレット

非特許文献

[0006]
非特許文献1 : Proc. Natl. Acad. Sci. USA, 86: 7706-10 (1989)

発明の概要

発明が解決しようとする課題

[0007]
 上記特許文献2で開示された方法では、各小区間における局所的な二次構造の抽出および存在確率の算出は、各小区間について得られる予測構造群のうち、該局所的な二次構造が発見される予測構造にわたって、ギッブス自由エネルギーから算出されるそれらの存在確率を総和することに基づいて行われる。これは、エネルギーのより低い二次構造の方が、その存在確率はより高い、ということを前提としている。ところでRNAは、生体内では、二次構造を形成するために十分な時間を与えられたときに、エネルギーの観点からより安定な二次構造として存在する可能性が高くなり、十分な時間経過後の各二次構造の存在確率の分布は、マクスウェル・ボルツマンの分布則に従うものとなると考えられる。それゆえ、上記前提はもっともであり、実際、上記特許文献2で開示された方法は、生物学的に妥当な結果を与え得ることが確認されてきた。さらに、局所的な二次構造の存在確率を論じるにあたり、それが発見される予測構造自体のエネルギーを用いることで、局所的な二次構造間の存在確率の大小を一般的に論じることができるようになった(それまで、長さと配列の違う局所的な二次構造間ではその存在確率の大小を論じることができなかった)。
[0008]
 しかしながら、本発明者によるその後の研究により、RNAの二次構造群の存在確率の分布を論じるにあたり、マクスウェル・ボルツマンの分布則に従うものとなるより前(すなわち、系が平衡状態に達する前)の時点の二次構造を予測することが重要であり得ることが明らかとなってきた。
 すなわち、mRNAなどのRNAが構造を形成するために与えられた時間が非常に短かった場合に、いわゆる動力学的トラップ(kinetic trap)の状態(すなわち、熱力学的に安定な二次構造を取るまでの、速度論的に安定な構造を取る状態)となり、動力学的トラップに陥っていた構造は、さらに時間を与えられたときに、カノニカルな塩基対相互作用だけを考慮した熱力学的に安定な状態(ここでは、「熱力学的トラップ(thermodynaic trap)の状態」と呼ぶ)となると考えられる。
 このように、生体内でのRNAの高次構造は段階的に変遷し、それぞれの時間スケールとしては、シングルストランドのRNAの状態から開始して、動力学的トラップの状態については1秒以内、熱力学的トラップの状態については6~10秒程度の時間で達するものと非常におおまかに想定される。さらにその後、熱力学的トラップの構造間で平衡状態に達するであろう。上記特許文献2の方法は、熱力学的トラップの平衡状態での二次構造を暗に想定するものである。該方法が、上述したように、生物学的に妥当な結果を与え得ることは、生体内でのRNAの状態が、熱力学的トラップの平衡状態に少なくとも近似できることを示唆しているものと考えられる。
 しかしながら、生体内でのRNA(特にmRNA)が構造を形成するために与えられた時間は、系が平衡状態に達するのに十分なものであると結論付ける理由がない。従って、生体内でのRNAの機能発現を検討する際に、平衡状態に達する前を想定し、とり得る二次構造として動力学的トラップの状態と熱力学的トラップの状態との間の二次構造を取り扱う方がより現実に近いと考えられる。ところが、従来、そのような時間的区間における二次構造を予測する方法はおろか、RNAの二次構造をそのような時間的区間において解析することの重要性の認識すら存在しなかった。

課題を解決するための手段

[0009]
 本発明者は、上記課題を解決するために鋭意検討した結果、構造が動力学的トラップから熱力学的トラップへ遷移する時間帯であっても、構造を組み始めるスタート時の影響を強く受ける時点を脱した時間帯を扱うのであれば、たとえ構造が遷移する時間帯においてであっても熱力学的トラップの状態のみを扱えばよいという考えに至った。その場合、動力学的トラップの状態から熱力学的トラップの状態への遷移はほぼ不可逆と見なすことができるから、熱力学的トラップの状態で存在する構造群が、動力学的トラップの状態で存在する構造群を無視できるほど優勢に存在していると考えることができる。そのとき、RNAの二次構造を解析すべき時間的区間においては、熱力学的トラップの状態での二次構造のみが相互に変換していると考えてよい。それゆえ、熱力学的トラップのRNA二次構造を、従来公知の二次構造予測方法を用いて決定しかつそれらの存在確率の分布を適当な初期条件として与え、それらの二次構造同士の相互変換を連立微分方程式で記述し、かつ該連立微分方程式を上記初期条件の下で解くことによって、上記時間的区間内の任意の時点における各二次構造の存在確率を予測することが可能となることを見出し、本発明を完成させた。
[0010]
 すなわち、本発明は、次の特徴を有するものである。
[1]RNAの二次構造の動態を予測する方法であって、
 (A)RNAの塩基配列から長さL1を有する区間を選択し、該区間において、平衡状態における2種以上の特定二次構造を決定し、かつ、各特定二次構造のギッブス自由エネルギーと平衡状態での存在確率とを算出するステップと、
 (B)前記特定二次構造同士の間の相互変換における遷移状態の二次構造を決定するとともに、各遷移状態の二次構造のギッブス自由エネルギーを求めるステップと、
 (C)各遷移状態の二次構造およびそれらのギッブス自由エネルギーから、各特定二次構造同士の間の相互変換の経路の活性化エネルギーを算出するステップと、
 (D)各活性化エネルギーを基にして、平衡状態に達するまでの時間的区間内の各時点における各特定二次構造の存在確率を決定するための連立微分方程式を形成するステップと、
 (E)下記(i)~(iv)のうちのいずれか1つの解法によって、各特定二次構造が、前記時間的区間内の各時点において何程の確率で存在し得るかを決定するステップと、
を有することを特徴とする、RNAの二次構造の動態を予測する方法。
  (i)前記平衡状態における各特定二次構造の存在確率の逆数を、前記時間的区間の或る時点t 0における各特定二次構造の存在比率とし、そこから導かれる該時点t 0での各特定二次構造の存在確率を初期値として前記連立微分方程式を解き、求める存在確率とする解法。
  (ii)前記時間的区間の或る時点t 0において各特定二次構造の存在比率が全て等しいものであるとして、各特定二次構造の該時点t 0での存在確率を決定し、これを初期値として、前記連立微分方程式を解き、求める存在確率とする解法。
  (iii)前記時間的区間の或る時点t 0における各特定二次構造の存在確率を、それぞれ正規分布として与えるものとし、それぞれの正規分布は、前記平衡状態における各特定二次構造の存在確率の逆数を、前記時間的区間の或る時点t 0における各特定二次構造の存在比率とすることで導かれる存在確率を平均値とし、かつ、該平均値と、各特定二次構造が等しい比率で存在した場合の存在確率との差の絶対値を標準偏差とするものとし、
 各特定二次構造毎に前記正規分布に従った1個の乱数を与えて対応する存在比率を求め、そこから導かれる各存在確率を初期値として前記連立微分方程式を解き、該初期値を用いた場合の存在確率を求める、という操作をN回繰り返した上で、各時点毎に、該N回の操作毎の存在確率の総和をNで除して平均し、求める存在確率とする解法。
  (iv)前記時間的区間の或る時点t 0における各特定二次構造の存在確率を、それぞれ正規分布として与えるものとし、それぞれの正規分布は、前記平衡状態における各特定二次構造が等しい比率で存在した場合の存在確率を平均値とし、かつ、該平均値と、前記平衡状態における各特定二次構造の存在確率の逆数を、前記時間的区間の或る時点t 0における各特定二次構造の存在比率とすることで導かれる存在確率との差の絶対値を標準偏差とするものとし、
 各特定二次構造毎に前記正規分布に従った1個の乱数を与えて対応する存在比率を求め、そこから導かれる各存在確率を初期値とし、前記連立微分方程式を解き、該初期値を用いた場合の存在確率を求める、という操作をN回繰り返した上で、各時点毎に、該N回の操作毎の存在確率の総和をNで除して平均し、求める存在確率とする解法。
[2]以下の(I)~(VI)のステップをさらに有するものである、上記[1]に記載の方法。
 (I)前記長さL1よりも大きい長さLを有する対象区間Sを与えるステップ。
 (II)前記対象区間S内の位置xの一端を位置x=x 1として、位置x 1と位置x 1+L1とを両端とする小区間をS(x 1)とし、
 上記のステップ(A)で選択される区間が前記小区間S(x 1)であるとして、上記のステップ(A)~(E)を行って、該小区間S(x 1)での前記各特定二次構造が、前記時間的区間内の各時点で何程の確率で存在し得るかを決定するステップ。
 (III)前記位置x 1を、所定の微小長さrだけ他端側へ移動させた位置をx 2 =x 1+rとし、位置x 2と位置x 2+L1とを両端とする小区間をS(x 2)とし、前記小区間S(x 1)の代わりに該小区間S(x 2)を用いる以外は前記(II)のステップと同様にして、該小区間S(x 2)での前記各特定二次構造が、前記時間的区間内の各時点で何程の確率で存在し得るかを決定する、という操作を、前記位置x=x 1 、x 2に続いて、小区間の一端が対象区間Sの他端に達するまで、位置x=x 3 、...、x Mについて繰り返し行い、それによって、最初の小区間S(x 1)から最後の小区間S(x M)までの各小区間毎の前記各特定二次構造が、前記時間的区間内の各時点で何程の確率で存在し得るかを決定するステップ。
 (IV)前記(II)~(III)のステップの結果に基づいて、前記の小区間S(x 1)からS(x M)のうちの任意の小区間S(x i)において、任意の時間tにおける前記各特定二次構造の存在確率を決定し、
 前記各特定二次構造において、位置および構造を表す所定の特徴パラメータ値を有する任意のモチーフM k が存在するか否かを調べ、その結果と上記で得られた前記時間tでの各特定二次構造の存在確率とに基づいて、該小区間S(x i)内において、前記時間tにおいて該モチーフM k が何程の蓋然性で存在するかを決定し、これを、小区間S(X i)における該モチーフM k に関する第1の蓋然性とする、
 という操作を、前記の小区間S(x 1)からS(x M)の全てに対して行なうステップ。
 (V)前記(IV)のステップにおいて各小区間内の各特定二次構造における前記モチーフM k の存在および非存在を調べた結果に基づいて、各小区間において、モチーフM k を有する特定二次構造が少なくとも1つ存在するか否かを決定する、
 という操作を、前記の小区間S(x 1)からS(x M)の全てに対して行い、そのような特定二次構造が少なくとも1つ存在する小区間の頻度を、モチーフM k に関する第2の蓋然性とするステップ。
 (VI)前記(IV)および(V)のステップで求められた、第1の蓋然性と、第2の蓋然性とに基づいて、前記時間tにおける前記モチーフM k の存在を決定するステップ。
[3]前記モチーフM k が、ステム・ループ構造を含むものである、上記[2]に記載の方法。
[4]前記モチーフM k が、2つのステム・ループ構造の間にあって他の塩基配列と相互作用していない間隙部分と、その両側にある前記ステム・ループ構造のそれぞれの内側のステム部分とを有してなる凹形構造を含むものである、上記[2]または[3]に記載の方法。
[5]RNAの二次構造の動態を予測するためのコンピュータプログラムであって、
 コンピュータを、
 (P1)RNAの塩基配列から長さL1を有する区間を選択し、該区間において、平衡状態における2種以上の特定二次構造を決定し、かつ各特定二次構造のギッブス自由エネルギーおよび平衡状態での存在確率を算出する手段、
 (P2)前記特定二次構造同士の間の相互変換における遷移状態の二次構造を決定するとともに、各遷移状態の二次構造のギッブス自由エネルギーを求める手段、
 (P3)各遷移状態の二次構造およびそれらのギッブス自由エネルギーから、各特定二次構造同士の間の相互変換の経路の活性化エネルギーを算出する手段、
 (P4)各活性化エネルギーを基にして、平衡状態に達するまでの時間的区間内の各時点における各特定二次構造の存在確率を決定するための連立微分方程式を形成する手段、
 (P5)下記(i)~(iv)のうちのいずれか1つの解法によって、各特定二次構造が、前記時間的区間内の各時点において何程の確率で存在し得るかを決定する手段、
として機能させるための前記コンピュータプログラム。
  (i)前記平衡状態における各特定二次構造の存在確率の逆数を、前記時間的区間の或る時点t 0における各特定二次構造の存在比率とし、そこから導かれる該時点t 0での各特定二次構造の存在確率を初期値として前記連立微分方程式を解き、求める存在確率とする解法。
  (ii)前記時間的区間の或る時点t 0において各特定二次構造の存在比率が全て等しいものであるとして、各特定二次構造の該時点t 0での存在確率を決定し、これを初期値として、前記連立微分方程式を解き、求める存在確率とする解法。
  (iii)前記時間的区間の或る時点t 0における各特定二次構造の存在確率を、それぞれ正規分布として与えるものとし、それぞれの正規分布は、前記平衡状態における各特定二次構造の存在確率の逆数を、前記時間的区間の或る時点t 0における各特定二次構造の存在比率とすることで導かれる存在確率を平均値とし、かつ、該平均値と、各特定二次構造が等しい比率で存在した場合の存在確率との差の絶対値を標準偏差とするものとし、
 各特定二次構造毎に前記正規分布に従った1個の乱数を与えて対応する存在比率を求め、そこから導かれる各存在確率を初期値として前記連立微分方程式を解き、該初期値を用いた場合の存在確率を求める、という操作をN回繰り返した上で、各時点毎に、該N回の操作毎の存在確率の総和をNで除して平均し、求める存在確率とする解法。
  (iv)前記時間的区間の或る時点t 0における各特定二次構造の存在確率を、それぞれ正規分布として与えるものとし、それぞれの正規分布は、前記平衡状態における各特定二次構造が等しい比率で存在した場合の存在確率を平均値とし、かつ、該平均値と、前記平衡状態における各特定二次構造の存在確率の逆数を、前記時間的区間の或る時点t 0における各特定二次構造の存在比率とすることで導かれる存在確率との差の絶対値を標準偏差とするものとし、
 各特定二次構造毎に前記正規分布に従った1個の乱数を与えて対応する存在比率を求め、そこから導かれる各存在確率を初期値とし、前記連立微分方程式を解き、該初期値を用いた場合の存在確率を求める、という操作をN回繰り返した上で、各時点毎に、該N回の操作毎の存在確率の総和をNで除して平均し、求める存在確率とする解法。
[6]コンピュータを、さらに以下の(p1)~(p6)の手段として機能させるためのものである、上記[5]に記載のコンピュータプログラム。
 (p1)前記長さL1よりも大きい長さLを有する対象区間Sを入力する手段。
 (p2)前記対象区間S内の位置xの一端を位置x=x 1として、位置x 1と位置x 1+L1とを両端とする小区間をS(x 1)とし、
 上記の手段(P1)で選択される区間が前記小区間S(x 1)であるとして、上記の手段(P1)~(P5)を用いて、該小区間S(x 1)での前記各特定二次構造が、前記時間的区間内の各時点で何程の確率で存在し得るかを決定する手段。
 (p3)前記位置x 1を、所定の微小長さrだけ他端側へ移動させた位置をx 2 =x 1+rとし、位置x 2と位置x 2+L1とを両端とする小区間をS(x 2)とし、前記小区間S(x 1)の代わりに該小区間S(x 2)を用いる以外は前記(p2)の手段と同様にして、該小区間S(x 2)での前記各特定二次構造が、前記時間的区間内の各時点で何程の確率で存在し得るかを決定する、という操作を、前記位置x=x 1 、x 2 に続いて、小区間の一端が対象区間Sの他端に達するまで、位置x=x 3 、...、x M について繰り返し行い、それによって、最初の小区間S(x 1)から最後の小区間S(x M)までの各小区間毎の前記各特定二次構造が、前記時間的区間内の各時点で何程の確率で存在し得るかを決定する手段。
 (p4)前記(p2)~(p3)の手段を用いて得られた結果に基づいて、前記の小区間S(x 1)からS(x M)のうちの任意の小区間S(x i)において、任意の時間tにおける前記各特定二次構造の存在確率を決定し、
 前記各特定二次構造において、位置および構造を表す所定の特徴パラメータ値を有する任意のモチーフM k が存在するか否かを調べ、その結果と上記で得られた前記時間tでの各特定二次構造の存在確率とに基づいて、該小区間S(x i)内において、前記時間tにおいて該モチーフM k が何程の蓋然性で存在するかを決定し、これを、小区間S(X i)における該モチーフM k に関する第1の蓋然性とする、
 という操作を、前記の小区間S(x 1)からS(x M)の全てに対して行なう手段。
 (p5)前記(p4)の手段を用いて各小区間内の各特定二次構造における前記モチーフM k の存在および非存在を調べた結果に基づいて、各小区間において、モチーフM k を有する特定二次構造が少なくとも1つ存在するか否かを決定する、
 という操作を、前記の小区間S(x 1)からS(x M)の全てに対して行い、そのような特定二次構造が少なくとも1つ存在する小区間の頻度を、モチーフM k に関する第2の蓋然性とする手段。
 (p6)前記(p4)および(p5)の手段で求められた、第1の蓋然性と、第2の蓋然性とに基づいて、前記時間tにおける前記モチーフM k の存在を決定するための出力を行う手段。
[7]前記モチーフM k が、ステム・ループ構造を含むものである、上記[6]に記載のコンピュータプログラム。
[8]前記モチーフM k が、2つのステム・ループ構造の間にあって他の塩基配列と相互作用していない間隙部分と、その両側にある前記ステム・ループ構造のそれぞれの内側のステム部分とを有してなる凹形構造を含むものである、上記[6]または[7]に記載のコンピュータプログラム。
[9]RNAの二次構造の動態を予測するための装置であって、
 (Y1)RNAの塩基配列から長さL1を有する区間を選択し、該区間において、平衡状態における2種以上の特定二次構造およびそれらの存在確率を算出する第一演算部と、
 (Y2)前記特定二次構造同士の間の相互変換における遷移状態の二次構造を決定するとともに、各遷移状態の二次構造のギッブス自由エネルギーを求める第二演算部と、
 (Y3)各遷移状態の二次構造およびそれらのギッブス自由エネルギーから、各特定二次構造同士の間の相互変換の経路の活性化エネルギーを算出する第三演算部と、
 (Y4)各活性化エネルギーを基にして、平衡状態に達するまでの時間的区間内の各時点における各特定二次構造の存在確率を決定するための連立微分方程式を形成する第四演算部と、
 (Y5)下記(i)~(iv)のうちのいずれか1つの解法によって、各特定二次構造が、前記時間的区間内の各時点において何程の確率で存在し得るかを決定する第五演算部と、
を有する、前記装置。
  (i)前記平衡状態における各特定二次構造の存在確率の逆数を、前記時間的区間の或る時点t 0における各特定二次構造の存在比率とし、そこから導かれる該時点t 0での各特定二次構造の存在確率を初期値として前記連立微分方程式を解き、求める存在確率とする解法。
  (ii)前記時間的区間の或る時点t 0において各特定二次構造の存在比率が全て等しいものであるとして、各特定二次構造の該時点t 0での存在確率を決定し、これを初期値として、前記連立微分方程式を解き、求める存在確率とする解法。
  (iii)前記時間的区間の或る時点t 0における各特定二次構造の存在確率を、それぞれ正規分布として与えるものとし、それぞれの正規分布は、前記平衡状態における各特定二次構造の存在確率の逆数を、前記時間的区間の或る時点t 0における各特定二次構造の存在比率とすることで導かれる存在確率を平均値とし、かつ、該平均値と、各特定二次構造が等しい比率で存在した場合の存在確率との差の絶対値を標準偏差とするものとし、
 各特定二次構造毎に前記正規分布に従った1個の乱数を与えて対応する存在比率を求め、そこから導かれる各存在確率を初期値として前記連立微分方程式を解き、該初期値を用いた場合の存在確率を求める、という操作をN回繰り返した上で、各時点毎に、該N回の操作毎の存在確率の総和をNで除して平均し、求める存在確率とする解法。
  (iv)前記時間的区間の或る時点t 0における各特定二次構造の存在確率を、それぞれ正規分布として与えるものとし、それぞれの正規分布は、前記平衡状態における各特定二次構造が等しい比率で存在した場合の存在確率を平均値とし、かつ、該平均値と、前記平衡状態における各特定二次構造の存在確率の逆数を、前記時間的区間の或る時点t 0における各特定二次構造の存在比率とすることで導かれる存在確率との差の絶対値を標準偏差とするものとし、
 各特定二次構造毎に前記正規分布に従った1個の乱数を与えて対応する存在比率を求め、そこから導かれる各存在確率を初期値とし、前記連立微分方程式を解き、該初期値を用いた場合の存在確率を求める、という操作をN回繰り返した上で、各時点毎に、該N回の操作毎の存在確率の総和をNで除して平均し、求める存在確率とする解法。
[10]さらに以下の(y1)~(y6)の手段を有するものである、上記[9]に記載の装置。
 (y1)前記長さL1よりも大きい長さLを有する対象区間Sを入力する手段。
 (y2)前記対象区間S内の位置xの一端を位置x=x 1として、位置x 1と位置x 1+L1とを両端とする小区間をS(x 1)とし、
 上記の手段(Y1)で選択される区間が前記小区間S(x 1)であるとして、上記の手段(Y1)~(Y5)を用いて、該小区間S(x 1)での前記各特定二次構造が、前記時間的区間内の各時点で何程の確率で存在し得るかを決定する手段。
 (y3)前記位置x 1を、所定の微小長さrだけ他端側へ移動させた位置をx 2 =x 1+rとし、位置x 2と位置x 2+L1とを両端とする小区間をS(x 2)とし、前記小区間S(x 1)の代わりに該小区間S(x 2)を用いる以外は前記(y2)の手段と同様にして、該小区間S(x 2)での前記各特定二次構造が、前記時間的区間内の各時点で何程の確率で存在し得るかを決定する、という操作を、前記位置x=x 1 、x 2 に続いて、小区間の一端が対象区間Sの他端に達するまで、位置x=x 3 、...、x M について繰り返し行い、それによって、最初の小区間S(x 1)から最後の小区間S(x M)までの各小区間毎の前記各特定二次構造が、前記時間的区間内の各時点で何程の確率で存在し得るかを決定する手段。
 (y4)前記(y2)~(y3)の手段によって決定された結果に基づいて、前記の小区間S(x 1)からS(x M)のうちの任意の小区間S(x i)において、任意の時間tにおける前記各特定二次構造の存在確率を決定し、
 前記各特定二次構造において、位置および構造を表す所定の特徴パラメータ値を有する任意のモチーフM k が存在するか否かを調べ、その結果と上記で得られた前記時間tでの各特定二次構造の存在確率とに基づいて、該小区間S(x i)内において、前記時間tにおいて該モチーフM k が何程の蓋然性で存在するかを決定し、これを、小区間S(X i)における該モチーフM k に関する第1の蓋然性とする、
 という操作を、前記の小区間S(x 1)からS(x M)の全てに対して行なう手段。
 (y5)前記(y4)の手段を用いて各小区間内の各特定二次構造における前記モチーフM の存在および非存在を調べた結果に基づいて、各小区間において、モチーフM を有する特定二次構造が少なくとも1つ存在するか否かを決定する、
 という操作を、前記の小区間S(x 1)からS(x M)の全てに対して行い、そのような特定二次構造が少なくとも1つ存在する小区間の頻度を、モチーフM k に関する第2の蓋然性とする手段。
 (y6)前記(y4)および(y5)の手段で求められた、第1の蓋然性と、第2の蓋然性とに基づいて、前記時間tにおける前記モチーフM k の存在を決定するための出力を行う手段。
[11]前記モチーフM k が、ステム・ループ構造を含むものである、上記[10]に記載の装置。
[12]前記モチーフM k が、2つのステム・ループ構造の間にあって他の塩基配列と相互作用していない間隙部分と、その両側にある前記ステム・ループ構造のそれぞれの内側のステム部分とを有してなる凹形構造を含むものである、上記[10]または[11]に記載の装置。

発明の効果

[0011]
 本発明の方法により、従来の方法では扱うことができなかった、RNAの二次構造がマクスウェル・ボルツマンの分布則に従う存在確率の分布を示すよりも前の時点における該RNAの二次構造の動態を予測することが可能となる。それにより、生体内でRNAが実際に取っている二次構造をより反映する形で、RNAの二次構造を予測することができる。本発明はまた、該方法をコンピュータに実行させるためのプログラム、および該方法を実施するための装置を提供する。

図面の簡単な説明

[0012]
[図1] 本発明の一つの実施形態による、装置またはプログラムの構成の例を模式的に示すブロック図である。
[図2] 本発明の一つの実施形態による、装置またはプログラムの構成の例を模式的に示すブロック図である。
[図3] 本発明のプログラムにおける、パラメータの入力のためのパネルの一例を示す。
[図4] 本発明のプログラムにおける、ユーザーに計算ステップの進捗状況を示すためのアイコンの一例を示す。
[図5] 本発明のプログラムにおける、ユーザーに微分方程式の計算の進捗状況を示すためのアイコンの一例を示す。
[図6] 実施例1からの3つのスペクトルを示す図である。(a)はスペクトル1、(b)はスペクトル2、(c)はスペクトル3に対応している。各スペクトルについて、横軸は塩基アドレス、縦軸は第1の蓋然性、ピークの色の濃さは第2の蓋然性を表している。閾値を設定して、第3の蓋然性が90%以上のピークのみを表示させている。
[図7] 実施例2からの3つのスペクトルを示す図である。(a)はスペクトル4、(b)はスペクトル5、(c)はスペクトル6に対応している。各スペクトルについて、横軸は塩基アドレス、縦軸は第1の蓋然性、ピークの色の濃さは第2の蓋然性を表している。閾値を設定して、第3の蓋然性が90%以上のピークのみを表示させている。
[図8] 実施例3における、T fold =3.0秒以内で、ほぼ平衡状態に達する系の2つの例を示す図である。
[図9] 実施例3における、平衡状態に向かって進みつつも、T fold =3.0秒では平衡に達していない系の2つの例を示す図である。図中では二次構造の動態のおおよその様子を見るために、0.001秒の刻み幅を明記するとともに6秒までプロットしている。
[図10] 実施例3における、0.001秒以内(最小刻み幅Δt内)に平衡状態に達してしまう系の例を示す図である。(a)は各特定二次構造の存在確率に基づくプロットであり、(b)は構造間の活性化エネルギーと速度定数のマトリックスである。このマトリックスにおいては各特定二次構造を1~4とし、最上段にそれぞれの二次構造のその時点でのポピュレーションを表示するとともに、矢印“=>”の記載の後に、マクスウェル・ボルツマン分布した場合のポピュレーション(平衡状態でのポピュレーション)を記載している。さらに、下部のマトリックスにおいては、1=>2のように、どの構造からどの構造への遷移であるかの記載の下に、その遷移の活性化エネルギー(kcal/mol)と遷移の速度定数Kobsと刻み幅Δtをかけた実際の遷移の量が記載されている。また本計算では、ある数種の構造間で遷移の量が一定以上である場合には、その関係が刻み幅内の時間ですでに平衡に達しているとしているが、それはマトリックス中では対応する構造の遷移のマスを同じ色で塗りつぶすことによって明記している。
[図11] 実施例3における、全く平衡状態に向かって進んで行かない系の例を示す図である。(a)は各特定二次構造の存在確率に基づくプロットであり、(b)は構造間の活性化エネルギーと速度定数のマトリックスである。このマトリックスにおいては各特定二次構造を1~6とし、最上段にそれぞれの二次構造のその時点でのポピュレーションを表示するとともに、矢印“=>”の記載の後に、マクスウェル・ボルツマン分布した場合のポピュレーション(平衡状態でのポピュレーション)を記載している。さらに、下部のマトリックスにおいては、1=>2のように、どの構造からどの構造への遷移であるかの記載の下に、その遷移の活性化エネルギー(kcal/mol)と遷移の速度定数Kobsと刻み幅Δtをかけた実際の遷移の量が記載されている。
[0013]
(発明の詳細な説明)
 以下、本発明の方法をより詳細に説明する。
 なお、本発明においてRNA上でいう「長さ」とは、配列長であって、配列された塩基の数を意味する。長さの最小単位は塩基1個であり、「長さ」は「塩基数」と読み替えてよい。よって、本明細書では、〔長さが100個〕などと記載し、例えば〔長さrだけ移動させる〕は〔塩基数rだけ移動させる〕を意味する。
[0014]
 上述した通り、本発明は、RNAの二次構造の動態を予測する方法を提供し、当該方法は、以下のステップ(A)~(E)を含む:
 (A)RNAの塩基配列から長さL1を有する区間を選択し、該区間において、平衡状態における2種以上の特定二次構造を決定し、かつ各特定二次構造のギッブス自由エネルギーおよび平衡状態での存在確率を算出するステップ;
 (B)前記特定二次構造同士の間の相互変換における遷移状態の二次構造を決定するとともに、各遷移状態の二次構造のギッブス自由エネルギーを求めるステップ;
 (C)各遷移状態の二次構造およびそれらのギッブス自由エネルギーから、各特定二次構造同士の間の相互変換の経路の活性化エネルギーを算出するステップ;
 (D)各活性化エネルギーを基にして、平衡状態に達するまでの時間的区間内の各時点における各特定二次構造の存在確率を決定するための連立微分方程式を形成するステップ;および、
 (E)各特定二次構造が、前記時間的区間内の各時点において何程の確率で存在し得るかを決定するステップ。
[0015]
 或る区間内のRNA配列は、二次構造を形成するために十分な時間を与えられたときに、その取り得る各二次構造の存在確率が時間的に実質的に変化しない状態に達し、そして該状態においては、各二次構造が、マクスウェル・ボルツマンの分布則に従う存在確率で分布すると考えられる。本明細書でいう「平衡状態」とは、二次構造の存在確率の分布が、実質的にマクスウェル・ボルツマンの分布則に従うものである状態をいう。ここで「実質的にマクスウェル・ボルツマンの分布則に従う」とは、各二次構造の存在確率が、例えば、マクスウェル・ボルツマンの分布則に従う存在確率に対して10%未満の差異、好ましくは1%未満の差異であることを意味する。
[0016]
 本発明の方法は、上記のような平衡状態にある時点よりも前(すなわち、「非平衡状態」)の任意の時点において、2種以上の二次構造が各々いかなる確率で存在し得るかを決定するための方法である。該任意の時点での各二次構造の存在確率を決定することにより、非平衡状態にある時間的区間内において、各二次構造の存在確率がどのように推移していくかを決定すること、すなわち、RNAの二次構造の動態を予測することが可能となる。従って、本明細書において、「RNAの二次構造の動態」とは、RNAが取り得る各二次構造の各時点における存在確率の経時的な変化の態様を意味している。なお、本発明の方法によって、非平衡状態にある時間的区間内でのRNAの二次構造の経時的な変化が予測され得るが、該時間的区間内の特定の時点を指定したときには、当該特定の時点での各二次構造の存在確率が予測されたものと考えることができる。従って、別の観点からは、本発明により、平衡状態に達するまでの時間的区間内の或る時点におけるRNAの二次構造を予測するための方法が提供される。以下で詳細に説明するように、本発明の方法では、非平衡状態にある時間的区間内での存在確率を解析される2種以上の二次構造は、平衡状態において存在し得ると予測された二次構造である。
[0017]
 本発明の方法による予測の対象とすべきRNAとしては、mRNA(本明細書では、特に矛盾のない限りはその前駆RNAをも包含する意味で用いられる)の他、miRNA(本明細書では、特に矛盾のない限りはその前駆RNAをも包含する意味で用いられる)、rRNA、tRNA、ncRNAなどが挙げられる。ncRNAは、タンパク質に翻訳されないRNAであって、機能性RNAなどとも呼ばれ、近年一段と注目されている。
 これらのなかでも、mRNAは、従来より医薬品開発において重要な創薬ターゲットの1つとして位置付けられており、アンチセンス医薬品やsiRNAの効率は、mRNAの二次構造に左右されるという知見もある。また、そうした二次構造に関する知見は、miRNAなどの既知または未知の分子生物学現象に関与する構造の解明に寄与し得ると考えられる一方で、一部の特定のmRNAを除いては、生体内におけるその高次構造がこれまで全く明らかにされていないことから、RNA中の有意味な二次構造の存在を予測、議論すべき対象としてきわめて重要である。
 一方、miRNAはステム・ループ構造を有する前駆体(pre-miRNA)から酵素dicerで切断されて生成し、mRNAを切断分解することなく物理的にその翻訳を抑制するので、この前駆体RNAの構造を安定化することにより、成熟miRNAの生成を抑制して結果的に標的mRNAの翻訳を促進することができる。このような観点から、pre-miRNAも、好ましい予測対象の1つであるといえる。
 また、ncRNAは、何らかの機能を生体内で果たしていると考えられている。まちがいなくそのような機能は、それらRNAが何らかの高次構造をとることで生体内での役割を担っていると考えられている。さらに、こうしたncRNAの中には、rRNAやtRNAのように(おそらく)一意の構造をとるであろうものも含まれる場合があるが、ほとんどのものはその構造的挙動がまったく未知である。このような観点から、ncRNAもまた、本発明にとって好ましい予測対象の1つであるといえる。
[0018]
 上記ステップ(A)では、先ず、RNAの塩基配列から、予測の対象とすべき長さL1の区間が選択される。区間の長さL1は、生物学的な条件をできるだけ模倣できるように設定することが好ましい。例えば、mRNAのような長いRNAが翻訳される過程においては、mRNA上の高次構造とタンパク質をリボソームが一本鎖に直しながら進行する。すなわち、リボソームが通り抜けた直後のmRNAは、いかような構造をもとり得る自由な状態にある。そうした自由な状態である部位は、次々と何らかの高次構造をとっていく。この状態は、現在の構造予測技術(例えば、本発明のプログラムにおける手段P1や、mfoldなどの公知の構造予測ソフトウェア)が想定している状況に近似していると考えられる。特に、mfoldにおいては、実験的に得られたパラメータを使用して構造予測を行っているが、そのパラメータは比較的短いRNAについての実験で得られたものである。つまりmfoldが正しい結果を与えるのは比較的短い(経験的に200塩基以下である)配列に対してだけである。そこで、リボソームの進行速度と、そうした現在の構造予測技術(後述するように、ステップ(A)での特定二次構造の決定は、従来公知の構造予測技術を用いて行われる)の特性とを考慮して、好ましい予測結果が得られるように、区間の長さL1を設定すればよい。
 また、区間の長さL1は、生物学的な機能を有し得る局所的な二次構造(例、ステム・ループ、シュードノット構造、凹形構造(後述)など)を含み得る長さであるべきである。
 以上のことから、具体的な区間の長さL1としては、例えば、塩基数50~400個、好ましくは塩基数100~300個が挙げられる。
[0019]
 次いで、上記のようにして選択された長さL1の区間において、後のステップでの解析に供されるべき平衡状態における2種以上の二次構造が決定される。また、決定された各二次構造のギッブス自由エネルギーおよび平衡状態での存在確率が算出される。ここで決定された平衡状態における二次構造を、本明細書では「特定二次構造」と呼ぶ。
 特定二次構造には、最も安定な(すなわち、ギッブス自由エネルギーが最も低い)二次構造が含まれることが好ましい。また、特定二次構造の種類は2種以上であれば特に制限されないが、典型的には2~20種である。
 非平衡状態において殆ど存在しないと考えられる二次構造をこの段階で除外しておくことが好ましい。そのような除外されるべき二次構造としては、例えば、ギッブス自由エネルギーがあまりに大きい(すなわち、平衡状態における存在確率があまりに低い)二次構造、および、実験的データなどから、非平衡状態において実質的に全く存在し得ないと予想される二次構造などが挙げられる。
 特定二次構造の選択のために、例えば、特定二次構造の個数を予め設定しておき、平衡状態での存在確率が高いもの(つまりギッブス自由エネルギーが低いもの)から優先的に選択すること、平衡状態での存在確率が所定の値以上の特定二次構造のみを選択すること、および予測ソフトウェアの望まれない特性のために算出された特定二次構造を排除することなどの基準を与えることもできる。特定二次構造の選択は、当業者によって適宜為され得る。
[0020]
 本発明の方法においては、長さL1の区間からの二次構造の抽出は、例えば、mfoldなどの従来公知のRNA構造予測ソフトウェアを利用して行えばよい。また、抽出された各二次構造のギッブス自由エネルギーは、例えば、mfoldなどの構造予測ソフトウェアによって算出された値を用いることができる。あるいは、抽出された各二次構造のギッブス自由エネルギーを、以下の参考文献:
   1. Mathews et al, J. Mol. Biol. (1999) 288, 911-940.
   2. Xia T et al, Biochemistry (1998) 37, 14719-14735.
   3. Giese MR et al, Biochemistry (1998) 37, 1094-1100.
   4. Draper DE et al, Methods in Enzymology (1995) 259, 281-242.
   5. Serra MJ et al, Methods in Enzymology (1995) 259, 242-261.
に記載の方法およびパラメータをそのまま使用して算出してもよい。このエネルギー算出方法は、小さな差異はあるものの広く使われており、好ましい結果を与え得ると考えられる。各二次構造に対して算出されたギッブス自由エネルギーを用いてマクスウェル・ボルツマンの分布則を適用することによって、各二次構造の平衡状態での存在確率を算出することができる。
 また、mfoldなどの構造予測方法自体の条件は、目的に適合するよう適宜に変更してよい。例えば、通常37℃の温度での構造予測が行われるが、好熱菌(高い温度の環境で生育する菌)などに見られるRNAの構造解析には、37℃よりはるかに高い温度を指定して構造予測をする。また、実験的データなどが得られている場合などには、mfoldなどの構造予測結果に影響を与えるその他のパラメータを変更してよい(ステム・ループのループの最低長さや、強制的にある構造を組ませるように指定するなど)。
 さらに、そのRNAのある特殊条件での構造解析をするため、推測に沿って、mfoldなどの構造予測結果に影響を与えるパラメータを変更してよい。例えば、ある領域にタンパク質またはアンチセンスなどが結合した条件での構造解析を行うのであれば、その領域を含むある一定領域を該RNA中での構造形成に参加しないものとして、該領域をどことも相互作用しないと指定して構造予測方法を用いてもよい。
[0021]
 次いで、上記ステップ(B)では、ステップ(A)で決定された複数の特定二次構造同士の間の相互変換における遷移状態の二次構造が決定され、かつ各遷移状態の二次構造のギッブス自由エネルギーが求められる。
 さらに、上記ステップ(C)では、ステップ(B)で得られた各遷移状態の二次構造およびそれらのギッブス自由エネルギーから、各特定二次構造同士の間の相互変換の経路の活性化エネルギーが算出される。
[0022]
 ここで、ある2種の特定二次構造について、それらの間の相互変換における遷移状態の二次構造を、それらの特定二次構造の最大公約数的共通構造としている。すなわち、ある特定二次構造AからBへの(およびその逆の)遷移について、最低限ほどかねばならない部分二次構造すべてをシングルストランドにした二次構造を、それらの間の相互変換における遷移状態の二次構造とする。各遷移状態の二次構造のギッブス自由エネルギーの算出は、上記したのと同様の方法で行うことができる。
[0023]
 上述した通り、本発明の方法では、非平衡状態においても、ステップ(A)において平衡状態で存在し得ると決定された二次構造のみが存在すると想定されている。
 上記ステップ(D)では、ステップ(A)で決定された各特定二次構造が、非平衡状態にある時間的区間内の各時点においていかなる確率で存在し得るかを決定するための連立微分方程式が形成される。該連立微分方程式の形成は、ステップ(C)で得られた各活性化エネルギーを用いて行われる。ここで、活性化エネルギーは、基準となる状態(そのギッブス自由エネルギーをΔGとする)から遷移状態(そのギッブス自由エネルギーをΔG TSとする)へのエネルギー変化ΔΔG=ΔG TS-ΔG(正の値を取る)と定義される。該方法をより詳しく説明すると以下の通りである。
[0024]
 いま、ステップ(A)で決定された特定二次構造がN個であるとし、それらを二次構造n(n=1、2、・・・、N)とする。また、二次構造nのギッブス自由エネルギーをG 、時間tでの存在確率をP n(t)とする。平衡状態となる時点はt=∞であると考えることができるため、平衡状態での二次構造nの存在確率をP n(∞)と表す。ここで、以下の式(1)が成り立つ。
[0025]
[数1]


[0026]
 上記式(1)中、Rは気体定数(R=8.31447 J K -1mol -1)、Tは絶対温度(K)である。
 さらに、二次構造nと二次構造mとの間の相互変換における遷移状態n,m の二次構造のギッブス自由エネルギーをG n,mとする。二次構造nから二次構造mへの変換の速度定数(k n→m)(単位:秒 -1)を示すアイリングの絶対反応速度式は以下の式(2)の通りである。但し、nとmが同じ値の場合は、速度定数はゼロとする。
[0027]
[数2]


[0028]
 上記式(2)中、Rは気体定数(R=8.31447 J K -1mol -1)、Tは絶対温度(K)、Kはアイリング定数(K=1.00)、k Bはボルツマン定数(k B=1.3806504x10 -23 J K -1)、hはプランク定数(h=6.62621x10 -34 J s)である。ここで式(2)中の(G n,m-G n)が、二次構造nから二次構造mへの相互変換において遷移状態n,m を通る経路の活性化エネルギーである。
 すると、二次構造nの存在確率の時間変化(dP n(t)/dt)は以下の式(3)のように表される。
[0029]
[数3]


[0030]
 上記式(3)においてn=1、2、・・・、Nとすることで得られる連立微分方程式を、適当な初期条件のもとで解くことによって、平衡状態に達するまでの時間的区間内の各時点における各特定二次構造の存在確率を決定することができる。
[0031]
 現実の解析においては、通常、ステップ(D)で得られる連立微分方程式は数値的に解かれる。すなわち、該連立微分方程式を、存在確率の初期値を与える時点t=t 0から計算の終了時点t=T foldまで適当な時間刻み幅Δt毎に数値的に解いていくことで、任意の時点tにおける各特定二次構造の存在確率P n(t)を得ることができる。このとき、1次近似を仮定すれば、存在確率P n(t)の計算式は、以下の式(4)の通りとなる。
[0032]
[数4]


[0033]
 次いで、上記ステップ(E)では、各特定二次構造が前記時間的区間内の各時点において何程の確率で存在し得るかを決定するために、ステップ(D)で得られた連立微分方程式が解かれる。連立微分方程式の初期値は、生物学的な条件をできるだけ模倣できるように設定することが好ましい。本発明の方法のステップ(E)では、該連立微分方程式は、以下に限定されないが、例えば、上記[1]のステップ(E)に示した(i)~(iv)のうちのいずれか1つの解法によって解かれ得る。
[0034]
 上記(i)の解法は、〔平衡状態における各特定二次構造の存在確率は、形成される水素結合ペアの個数に比例しており、また水素結合ペアの個数が多い特定二次構造ほどその形成に時間がかかるはずであるため、非平衡状態の或る時点t 0においては、各特定二次構造は平衡状態におけるそれらの存在確率の逆数に比例した確率で存在し得る〕との考えに基づいている。そのようにして時点t 0における各特定二次構造の存在確率を初期条件として与え、例えば上記式(4)をコンピュータを用いて数値的に解くことによって、各特定二次構造が非平衡状態にある時間的区間内の各時点において何程の確率で存在し得るかを決定することができる。
[0035]
 上記(ii)の解法は、非平衡状態の或る時点t 0において、全ての特定二次構造が等しい確率で存在し得るとする考えに基づいている。そのようにして時点t 0における各特定二次構造の存在確率を初期条件として与え、例えば上記式(4)をコンピュータを用いて数値的に解くことによって、各特定二次構造が非平衡状態にある時間的区間内の各時点において何程の確率で存在し得るかを決定することができる。
[0036]
 上記(iii)の解法は、時点t 0における各特定二次構造の存在確率を、各々特定の数値として与えるのではなく、各特定二次構造毎に正規分布として与えるものである。或る特定二次構造のための正規分布の平均値は、上記(i)の解法における初期時点t 0でのその存在確率と同様にして得られる。一方、或る特定二次構造のための正規分布の標準偏差は、上記のようにして得られた平均値と、上記(ii)の解法と同様にして求めた存在確率との差の絶対値により与えられる。各特定二次構造毎に、各々の正規分布に従う確率で選択された1個の乱数を与え、さらに、与えられた乱数を全ての特定二次構造について総和することで得られる値で、各特定二次構造のために与えられた乱数を除することによって、各特定二次構造の時点t 0における存在確率が得られる。そのようにして得られた時点t 0における各特定二次構造の存在確率を初期条件として、例えば上記式(4)をコンピュータを用いて数値的に解くことで、「該初期値を用いた場合の」、各特定二次構造の非平衡状態にある時間的区間内での任意の時点での存在確率が求まる。上記の〔正規分布に従った乱数を与えて、該乱数に基づいて時点t 0における各特定二次構造の存在確率を決定し、それを初期条件として連立微分方程式を解く〕という操作を合計でN回繰り返す。各操作毎の、任意の時点tにおける或る特定二次構造の存在確率を算術平均する(すなわち、該存在確率の総和をNで除する)ことで得られる値が、当該特定二次構造の時点tにおける存在確率となる。上記操作の繰り返し回数Nは、大きければ大きいほど乱数の選択に左右されない結果が得られるため好ましいが、実際の解析では、通常1~100の範囲から選択される適当な回数とすればよい。
[0037]
 上記(iv)の解法は、時点t 0における各特定二次構造の存在確率を、各々特定の数値として与えるのではなく、各特定二次構造毎に正規分布として与えるものである。或る特定二次構造のための正規分布の平均値は、上記(ii)の解法における初期時点t 0でのその存在確率と同様にして得られる。一方、或る特定二次構造のための正規分布の標準偏差は、上記のようにして得られた平均値と、上記(i)の解法と同様にして求めた存在確率との差の絶対値により与えられる。各特定二次構造毎に、各々の正規分布に従う確率で選択された1個の乱数を与え、さらに、与えられた乱数を全ての特定二次構造について総和することで得られる値で、各特定二次構造のために与えられた乱数を除することによって、各特定二次構造の時点t 0における存在確率が得られる。そのようにして得られた時点t 0における各特定二次構造の存在確率を初期条件として、例えば上記式(4)をコンピュータを用いて数値的に解くことで、「該初期値を用いた場合の」、各特定二次構造の非平衡状態にある時間的区間内での任意の時点での存在確率が求まる。
 上記の〔正規分布に従った乱数を与えて、該乱数に基づいて時点t 0における各特定二次構造の存在確率を決定し、それを初期条件として連立微分方程式を解く〕という操作を合計でN回繰り返す。各操作毎の、任意の時点tにおける或る特定二次構造の存在確率を算術平均する(すなわち、該存在確率の総和をNで除する)ことで得られる値が、当該特定二次構造の時点tにおける存在確率となる。上記操作の繰り返し回数Nは、大きければ大きいほど乱数の選択に左右されない結果が得られるため好ましいが、実際の解析では、通常1~100の範囲から選択される適当な回数とすればよい。
[0038]
 上記(i)~(iv)の解法において、コンピュータを用いた現実の数値解析では、通常、存在確率の初期値を与える時点t 0をt=0として計算すればよい。また、上記時間刻み幅Δtは、1ピコ秒~1ミリ秒であり得、例えば1ミリ秒(0.001秒)である。
 また、計算の終了時点t=T foldについては、生物学的な条件を考慮に入れて設定することが好ましい。あるいは、mRNAのある時点での構造を見たい場合には、任意の数値を入れることも間違いではない。そこで、T foldの値としては、例えば、0~30秒であり得る。なお、後述する実施例3で示すように、各特定二次構造の存在確率の時間的変遷を解析したとき、おおまかには以下の4つのケースが観察される:
(c1)時点t=T foldまでに、各特定二次構造の存在確率が平衡状態での存在確率に収束するケース;
(c2)各特定二次構造の存在確率が平衡状態での存在確率に向かいつつも、時点t=T foldまでには収束に達しないケース;
(c3)初期時点t=0の直後に、各特定二次構造の存在確率が平衡状態での存在確率に収束するケース;および、
(c4)時点t=T foldまでに、各特定二次構造の存在確率が、平衡状態での存在確率に向かって全く進まないケース。
 上記(c2)のケースは、本発明の方法で想定される時間枠において、もっともドラスティックに各構造の存在確率が時間の影響を受けるケースと言える。
 上記(c3)のケースは、各特定二次構造間の活性化エネルギーが極めて低く、二次構造間の変換が非常に速い(すなわち、速度定数が非常に大きい)ものである。このような場合、理論的に正しくは刻み幅を細かくすることで解決されるが、後述する実施例で行っているように、二次構造間の相互変換を数値的に解かずに、マクスウェル・ボルツマンの分布則によって存在比率を与えてもよい。
 上記(c4)のケースは、各特定二次構造間の活性化エネルギーが極めて高く、二次構造間の変換が非常に遅い(すなわち、速度定数が非常に小さい)ものである。これは、事実上平衡状態に達しない系であり、各特定二次構造の存在確率は初期条件に大きく左右される。この場合、生体内でそれらの特定二次構造が相互変換し、やがて存在確率の分布が平衡状態に達するとする考え自体が正しくなく、むしろ、ある初期時点での二次構造が単純にそのまま持続すると考えるのが妥当であり得る。ここで、自然界の現象として、RNAの二次構造の初期値が毎回同じである保証はないため、その部分の二次構造の存在確率の分布は、毎回異なるものであり得る。二次構造の存在確率の分布が毎回異なるものであるRNA部分の構造が生物学的に意義がある可能性は低いと考えられる。このような考え方に基づいて、本発明の方法に従って、RNAの塩基配列上で生物学的機能を有する部分と有しない部分とを区別することが可能となり得る。
 現実の解析では、単一の小区間であっても上記(c1)~(c4)のケースが混在する場合もある。
[0039]
 本発明の方法ではまた、上記特許文献2で開示された概念に従って、上記長さL1を有する区間を該文献の方法における小区間と考えて該区間において上記ステップ(A)~(E)を行い、続いて該区間をRNA塩基配列上で数塩基ずらすことで得られる別の区間においても同様のステップを行う、という操作を繰り返した後、全てのそのような区間において得られた結果を集積する手法も用いることができる。すなわち、上述した本発明の方法は、一つの実施形態では、上記[2]に示した(I)~(VI)のステップをさらに有する。
[0040]
 以下に、上記[2]に示した実施形態を詳細に説明する。
 なお、以下では、配列中の塩基の位置(配列の順番)を特定し得るように1つ1つの塩基に付与された番号(本明細書では、5’側から3’側に向かって数える)を「塩基アドレス」と呼んで説明する。塩基配列上の位置xは塩基アドレスであり、各小区間の先頭(5’側)からどの位置にあるかを示す相対的な塩基アドレスを用いてもよいが、対象区間Sの先頭(5’側)からどの位置にあるかを示す塩基アドレスを用いても、両者を併用してもよい。
[0041]
 上記(I)のステップでは、二次構造の動態を解析すべき対象区間Sが選択される。本実施形態においても、該長さL1は、上述したものと同様の大きさであってよい。
 対象区間Sの長さLは、長さL1より大きい限り限定されず、用途に応じて任意に決定してよく、例えばRNAの全長としてもよい。例えば、所望の局所的な二次構造の存在をタンパク質への翻訳を抑制する目的で予測する場合、翻訳抑制への寄与が大きいと予測される位置、例えば、5’非翻訳領域(5’UTR)、5’UTRとコード領域(CDS)とにまたがる領域、好ましくは翻訳開始点およびその近傍などが選択され得る。一方、所望の局所的な二次構造の存在をタンパク質への翻訳を増大させる目的で予測する場合、mRNAの安定性への寄与が大きいと予測される位置、例えば、3’非翻訳領域(3’UTR)などが選択され得る。あるいは、フレームシフトなどの生命現象を引き起こすシュードノット構造(特殊なステム・ループ構造)を見つける目的ならばコード領域が主に選択され得る。
[0042]
 上記(II)のステップでは、対象区間S内の一端(すなわち、最も5’側)である位置x=x 1と、位置x=x 1+L1とを両端とする小区間S(x 1)を、上記のステップ(A)で選択される長さL1を有する区間であると想定して、該小区間S(x 1)において上記のステップ(A)~(E)を行うことにより、該小区間S(x 1)において決定された各特定二次構造が、非平衡状態にある時間的区間内の各時点で何程の確率で存在し得るかを決定する。上記のステップ(A)~(E)の実行自体は、上述した通りに行えばよい。なお、(II)のステップおよび以下で説明する(III)のステップにおいては、特定二次構造の個数が1つである小区間が存在してもよい。特定二次構造が1つである小区間においては、二次構造の動態を解析すべき前記時間的区間内の任意の時点において、当該1つの特定二次構造が100%の存在確率で存在するとすればよく、従って、上記のステップ(B)~(E)を行わなくてよい。
[0043]
 上記(III)のステップでは、先ず、上記位置x 1を、所定の微小長さrだけ他端側(すなわち、5’→3’方向)へ移動させた位置をx 2 =x 1+rとし、位置x 2と位置x 2+L1とを両端とする小区間をS(x 2)とする。
 微小長さrは、小区間同士が互いに重なり合うように小区間の長さL1よりも短い値であればよいが、数多く重なり合わせることによって、後に詳述する蓋然性の値自体の信頼性を高めることができる。微小長さrは、塩基1個に近い小さい値が好ましいが、演算処理すべき小区間の数がそれに応じて増加するので、用いる演算装置の処理能力や、必要な信頼性の程度、計算する対象であるRNAのおかれた生物学的環境などに応じて適当な数値を選択すればよい。これらの点からは、微小長さrの塩基数は1~10個、特に、1~3個が好ましい範囲である。なお、微小長さrは、どの程度の精度で蓋然性を算出するかに影響する因子であり、分解能(resolution)と考えることもできる。
[0044]
 上記(II)のステップでの操作と同様にして、小区間S(x 2)においても上記のステップ(A)~(E)を行うことにより、該小区間S(x 2)において決定された各特定二次構造が、非平衡状態にある時間的区間内の各時点で何程の確率で存在し得るかを決定する。
 このように、小区間を微小長さrだけ他端側へ移動させて、そのときの小区間でステップ(A)~(E)を行うという操作を、小区間の一端が対象区間Sの反対端に到達したときの最後の小区間S(x M)まで繰り返し行う。このような繰り返し操作によって、各小区間において決定される各特定二次構造が、非平衡状態にある時間的区間内の各時点で何程の確率で存在し得るかが決定される。x iの添え字iは、位置x iが位置x 1を(i-1)回微小長さrだけ(すなわち、長さ{r×(i-1)}だけ)他端側へ移動させた位置であることを示している。また、小区間S(x i)は、位置x iと位置x i+L1とを両端とする小区間である。
[0045]
 生体内の条件を適切に算出するために、上記(II)および(III)のステップにおいて、上記ステップ(E)における解法の選択、初期条件を与える時点t 0、計算の終了時点T foldを各小区間毎にそれぞれ独立に与えてもよいし、或いは、簡単のためにそれらを全ての小区間で統一させてもよい。
[0046]
 また、上記(II)および(III)のステップでは、対象区間Sの一端から小区間の位置を順次移動させて順番に解析しているように表現しているが、その演算や操作は、必ずしも端から順番通りに行なうことだけを意味するものではない。即ち、対象区間Sから、互いに一部がオーバラップした小区間S(x 1)~S(x M)を順不同に選び出し、上記のステップ(A)~(E)を行ってもよい。すなわち、上記(II)および(III)のステップは、結果として、対象区間Sにそのような互いに重なり合った小区間が設定されてそれぞれについて解析された状態を得ることを意味する。演算やプログラムの進行においては、端から順番の解析が好ましい手順である。
[0047]
 ステップ(IV)では、先ず、上記(II)~(III)のステップでの各小区間において、各小区間毎に決定された各特定二次構造が、任意に選択された時間tにおいて何程の確率で存在し得るかが決定される。その決定のためには、上記(II)~(III)のステップで得られた結果を参照すればよい。
[0048]
 続いて各小区間における各特定二次構造毎に、所定の特徴パラメータ値を有する任意のモチーフM K が存在するか否かが調べられる。ここで、特徴パラメータ値は、モチーフの位置および構造を表すものであり、従って、本明細書におけるモチーフは、構造だけでなく、RNA上での位置をも含む概念であることに留意されたい。求めるべきモチーフの塩基数は、長さL1以下である限り特に制限はなく、後述のモチーフを規定する各特徴パラメータ値のとり得る範囲の組み合わせによってその下限と上限が定まる。例えば、モチーフが、ステムの長さがx1~x2塩基、ループの長さがy1~y2塩基、ステム・ループ内の許容ミスマッチ数がi個のステム・ループの場合、モチーフの長さは[x1×2+y1]~[x2×2+y2+i]塩基となる。好ましいモチーフの塩基数としては、例えば、10~100塩基程度が挙げられ、20~60塩基程度がより好ましい塩基数として挙げられる。特定二次構造におけるモチーフの存在または非存在の決定は、該特定二次構造が、該モチーフの特徴パラメータ値を満たす局所的な二次構造を含むか否かを単純にパターンマッチングすることによって行うことができる。また、明らかなように、当該モチーフの全配列(あるいは、当該モチーフを構成するために必要な全配列;以下同様であり、以後は特に断らない。)を含まない小区間については、小区間内で決定されたいずれの特定二次構造にも当該モチーフは存在し得ないため、パターンマッチングは、当該モチーフの全配列を含む小区間においてのみ行えばよい。
[0049]
 ここで、モチーフとはすなわち局所的な二次構造であり、その二次構造の特徴をパラメータ値で規定し得るものであればよい。また、モチーフを規定するパラメータ値が、特徴パラメータ値である。
 具体的なモチーフとしては、ステム・ループ構造、シュードノット(pseudoknot)構造(偽結び目構造)、凹形構造などが挙げられる。
 凹形構造とは、本発明者が重要な局所的な二次構造として着目したものであって、互いに接近した2つのステム・ループ構造の間に存在し、他の塩基配列と相互作用していない間隙部分を有する構造である。その間隙部分自体は平坦ではあるが、その両側にステムに関与する塩基が存在する点から、モチーフとしている。凹形構造は、グローイン(groin)構造と呼ぶこともできる。
[0050]
 ステム・ループ構造を規定する特徴パラメータ値としては、RNA上での位置、ステムの長さ、ループの長さ、ステム・ループ内の許容ミスマッチの数が挙げられる。
 ステム・ループ構造は、必ず一対以上の塩基対(ステム部分)および折り返し部位(ループ部分)を持つので、必然的にパラメータの最小値は決定される(すなわち、ステムの長さ1以上、ループの長さ3以上、ミスマッチの数0以上)。しかし計算のスピードを上げるため、パラメータの最小値を上げても良い。
 ステムの長さ(塩基対形成に寄与する片方の鎖の塩基数を意味する。G-Uのwobblingを含む)は、特に上限値は定めなくてもよいが、塩基数4~50、より好ましくは5~20である。
 ループの長さは、塩基数3~50、より好ましくは4~10である。
 ステム・ループ内の許容ミスマッチの数は、塩基数0~30、より好ましくは0~10である。
[0051]
 シュードノット構造を規定する特徴パラメータ値としては、ステム・ループを規定するパラメータに加え、主たるステム・ループ(より5’端側に存在するステム・ループ)のループ部分のうち3’端側の配列と相互作用する連続したある塩基数と、その部分に相補的な配列の主たるステム・ループから見た距離(塩基数)があげられる。
 シュードノット構造は、必ず一対以上の塩基対(ステム部分)と折り返し部位(ループ部分)、および一対以上の相互作用部位(ループ部分と相互作用する部位)を持つので、必然的に上記パラメータの最小値は決定される(すなわち、ステムの長さ1以上、ループの長さ3以上、ミスマッチの数0以上、もう一つの連続した塩基数1以上、距離2以上)。しかし計算のスピードを上げるため、パラメータの最小値を上げても良い。
 ループ部分のうち3’端側の配列と相互作用する連続したある塩基数は3~50、より好ましくは、5~20である。
 その部分に相補的な配列の主たるステム・ループの3’端からの距離(塩基数)は4~50、より好ましくは5~20である。
[0052]
 凹形構造では、中心となる間隙部分の両側にそれぞれステムa、bが存在する。両側のステムa、bは、それぞれ別個に塩基対をなしているが、凹形構造を規定するには、ステムa、bのそれぞれの塩基対のうちの内側の塩基(すなわち、それぞれの間隙部分側の塩基)の数を指定するだけで充分である。よって、本発明では、凹形構造を規定する好ましい特徴パラメータ値として、中心となる間隙部分の塩基数と、その両側にある2つのステム・ループ構造のそれぞれのステムの内側部分の塩基数の和、および該ステムの内側部分に許容されるミスマッチの数を採用する。
 凹形構造は、必ず間隙部分をはさんで形成されている2個以上の塩基対に参加する塩基(間隙部分をはさむ2本のステムで、間隙部分側の塩基)を持つので、必然的に上記パラメータの最小値は決定される(すなわち、間隙部分0以上、ステムの内側部分の塩基数2以上、ミスマッチの数0以上)。しかし計算のスピードを上げるため、パラメータの最小値を上げても良い。
 間隙部分の塩基数は0~50、より好ましくは0~10である。
 間隙部分の両側にある2つのステムのそれぞれの内側部分の塩基数の和は、4~50、より好ましくは4~30である。
 ミスマッチは0~10、より好ましくは0~5である。
[0053]
 なお、モチーフの「位置」そのものは、点(塩基1つ分の位置)ではなく、モチーフを直線状に展開したときの塩基長をもった区間である。演算などの取り扱い上、モチーフの位置を1点の塩基アドレスで表した方が便利な場合には、その位置が占める配列区間の先頭などの1点の塩基の塩基アドレスと、モチーフの塩基長とによって、モチーフの位置を表現してもよい。
[0054]
 上記のようにして、各小区間毎に、時間tでの各特定二次構造の存在確率と、各特定二次構造内にモチーフM が存在するか否かとが決定され、続いて、これらを用いて各小区間内において、時間tにおいてモチーフM k が何程の蓋然性で存在するかが決定される。第1の蓋然性(VS)は、例えば以下のようにして算出することができる。
 一般に、或る小区間S(x i)(i=1,2,3…M)について、時間tにおける特定二次構造j(j=1,2,3…N)の存在確率をP i,j(t)とし、或る小区間iについて得られた特定二次構造全ての存在確率の和を、
[0055]
[数5]


[0056]
とする。
 さらに、そのようにして求められたP i,total(t)を用いて、
[0057]
[数6]


[0058]
を計算する。このようにして求められるVS i(M k)(t)が、時間tにおける小区間iでのモチーフM k についての第1の蓋然性となる。また上記の場合、特定二次構造jにM k が発見された場合であっても、そのM k が特定二次構造jにおいて最も5’側または最も3’側の構造であることを理由に、発見されなかったものとして扱うこともできる。
 ところで、これが平衡状態にある場合には、以下の式:
[0059]
[数7]


[0060]
が成り立つ。ここでのΔG i,jとは、或る小区間iについて、特定二次構造jのギッブス自由エネルギーであり、ΔG i,minは或る小区間iについて、最も安定な特定二次構造のギッブス自由エネルギーである。
 上記式(7)は、各特定二次構造の平衡状態での存在確率を求めるに際して、マクスウェル・ボルツマンの分布則を適用したものであり、これと上記式(5)および(6)を用いて求められる平衡状態での第1の蓋然性は、上記特許文献2に開示された方法で求められる第1の蓋然性に対応するものである。本発明者は、上記特許文献2で開示した発明の方法では、この方針を採らなかったが、それは、小区間内の特定二次構造を全て発見できていないだろうことを意識した上で、簡便性を優先したためである。本発明の方法においては、時間の概念を入れて厳密性を追及するため、上記式(7)を利用することが好ましい。
 ここで例えば、或る小区間S(x i)において、上記(II)または(III)のステップで4つの特定二次構造s1、s2、s3およびs4が得られており、時間tにおけるそれらの存在確率をそれぞれ、P s1(t)、P s2(t)、P s3(t)およびP s4(t)とし、かつ、得られた特定二次構造全ての存在確率の和をP total(t)とする。さらに、特定二次構造s1、s2およびs3はモチーフM k1 を有するが特定二次構造s4はモチーフM k1 を有さず、かつ特定二次構造s1およびs4はモチーフM k2 を有するが特定二次構造s2およびs3はモチーフM k2 を有さないものとする。この場合、時間tにおいて、モチーフM k1 が存在する蓋然性は;
  {P s1(t)+P s2(t)+P s3(t)}/P total(t)、
モチーフM k2が存在する蓋然性は;
  {P s1(t)+P s4(t)}/P total(t)
となる。
 決定された蓋然性の値が、小区間S(x i)での時間tにおける上記各モチーフに関する第1の蓋然性である。上記の例の場合、
 VS i(M k1)(t)={P s1(t)+P s2(t)+P s3(t)}/P total(t)、
 VS i(M k2)(t)={P s1(t)+P s4(t)}/P total(t)
となる。
[0061]
 次に、上記(V)のステップでは、上記(IV)のステップにおいて各小区間内の各特定二次構造におけるモチーフM k の存在および非存在を調べた結果に基づいて、各小区間を、(a)当該小区間においてモチーフM k を有する特定二次構造が少なくとも1つ存在する、または(b)当該小区間においてモチーフM k を有する特定二次構造が存在しない、のいずれであるかについて評価する。当該評価は、前記の小区間S(x 1)からS(x M)の全ての小区間に対して行うべきであるが、実際には、当該モチーフの全配列を含まない小区間にはモチーフM k を有する特定二次構造が存在し得ないため、当該モチーフの全配列を含む小区間のみ考慮すればよい。モチーフM k に対して,どれ程の小区間が上記(a)に評価されるかは様々であって、1つの小区間のみの場合から、モチーフM k の全配列を含む全ての小区間の場合まで、様々である。ここで、上記(a)に評価される小区間の頻度を求め、これをモチーフM k に関する第2の蓋然性とする。
 この第2の蓋然性は、上記(a)に評価される小区間の個数で表してもよいし、当該個数を、モチーフM k の全配列を含む全ての小区間の個数で割って100を乗じた数値で表してもよい。
 より好ましくは、第2の蓋然性は以下の通りに算出される。すなわち、注目しているモチーフM k に関して、該モチーフの全配列が含まれる全ての小区間S(x i)(i=1,2,3…M)の各々に対して重み係数W(i,M k)を与え、下記式:
[0062]
[数8]


[0063]
(式中、”{“記号の表記は、小区間S(x i)について、いずれかの特定二次構造にモチーフM が発見された場合はW(i,M k)を、いずれの特定二次構造にもモチーフM k が発見されない場合は0を加えることを意味する。)を計算する。ここで、式(8)におけるW(M k)は:
[0064]
[数9]


[0065]
で与えられる。重み係数W(i,M k)は、通常全てのiに対して1とすればよいが、変数としてもよい。W(i,M k)を定数1とした場合には、W(M k)は、モチーフM k の全配列が含まれる小区間の総数に等しい。さらに、VD i(M k)を、モチーフM k の全配列が含まれる全ての小区間にわたって総和することで、
[0066]
[数10]


[0067]
を計算する。このようにして得られるVD(M k)を、第2の蓋然性とすることができる。VD i(M k)の算出は一見冗長であるが、後述する第3の蓋然性の算出を容易にし得るため、好ましい。
 第1の蓋然性VSが、各小区間内で評価した場合のモチーフの存在の蓋然性であるのに対し、第2の蓋然性VDは、各小区間を互いに関連させることで得られたモチーフの存在の蓋然性である。
[0068]
 上記(VI)のステップでは、第1の蓋然性と第2の蓋然性とに基づいて、時間tにおけるモチーフM k の存在を決定する。例えば、第1の蓋然性と第2の蓋然性とが共に高いものが、実際に存在する可能性の高いモチーフであると決定することができる。あるいは、第1の蓋然性と第2の蓋然性とが共に高い位置がどこにあるかをより効果的に示すために、上記特許文献2に記載したのと同様にして、これら2つの蓋然性を1つのグラフ内で関連づけて視覚的に表示してもよい。視覚的な表示の方法については、本発明のプログラムにおいて後述する。さらに代替的には、第1の蓋然性と第2の蓋然性とから新たな値(これを「第3の蓋然性VI」とする)を算出し、第3の蓋然性に基づいてモチーフの存在を決定することもできる。第3の蓋然性としては、例えば、第1の蓋然性と第2の蓋然性とを掛け合わせた後に総和をとり、その結果を100で割ることで得られる値を用いることができ、そのようにして算出された第3の蓋然性の値が大きい場合に、モチーフが存在する蓋然性が高いと決定され得る。その場合、第3の蓋然性を、例えば、上記したVS i(M k)(t)およびVD i(M k)を用いて、:
[0069]
[数11]


[0070]
を計算することによって求めることができる。
[0071]
 次に、本発明の方法を実施するための装置の構成を具体的に示す。
 図1は当該装置の構成の一例を模式的に示すブロック図である。上記[9]および図1に示すとおり、当該装置は、第一演算部(Y1)、第二演算部(Y2)、第三演算部(Y3)、第四演算部(Y4)および第五演算部(Y5)を少なくとも有して構成される。
[0072]
 第一演算部(Y1)は、上記ステップ(A)を行うための演算部である。第一演算部(Y1)には、予測の対象とすべき長さL1のRNA塩基配列、温度など、予測に必要な初期値が入力される。RNA塩基配列は、オペレータによるデータ名の入力に基づいて、自体の記憶装置(HDD、CD-ROMなど)や外部データベースにアクセスしデータを取り込むように構成されていてもよい。これらの初期値については、本発明の方法において説明したとおりである。第一演算部(Y1)は、入力されたこれらの初期値に基づいて、予測の対象とする区間において、平衡状態における複数の特定二次構造ならびにそれらのギッブス自由エネルギーおよび存在確率を算出する演算を行う。その後の解析に利用すべき二次構造の選択を可能とするために、第一演算部(Y1)は、初期値として二次構造の選択の基準の入力を可能とするか、あるいは、予測の対象とされた区間での特定二次構造の候補ならびにそれらのギッブス自由エネルギーおよび/または存在確率を出力し、オペレータによる特定二次構造の選択の入力を可能とするものであることが好ましい。
 第二演算部(Y2)は、上記ステップ(B)を行うための演算部である。第二演算部(Y2)では、第一演算部(Y1)で算出された特定二次構造に対して、それらの間での相互変換における遷移状態の二次構造が決定され、さらに各遷移状態の二次構造のギッブス自由エネルギーが求められる。
 第三演算部(Y3)は、上記ステップ(C)を行うための演算部である。第三演算部(Y3)では、第二演算部(Y2)で求められた各遷移状態の二次構造およびそれらのギッブス自由エネルギーから、各特定二次構造同士の間の相互変換の経路の活性化エネルギーが算出される。
 第四演算部(Y4)は、上記ステップ(D)を行うための演算部である。第四演算部(Y4)では、第三演算部(Y3)で算出された各活性化エネルギーを基にして、平衡状態に達するまでの時間的区間内の各時点における、第一演算部(Y1)で決定された各特定二次構造の存在確率を決定するための連立微分方程式が形成される。
 第五演算部(Y5)は、上記ステップ(E)を行うための演算部である。第五演算部(Y5)には、第四演算部(Y4)で得られた連立微分方程式を数値的に解くために、初期条件、終了時点T fold、時間刻み幅Δtなど、解析に必要な値が適宜入力される。本発明の方法において説明したように、本発明では、連立微分方程式は上記(i)~(iv)のうちのいずれか1つの解法によって解かれ得る。従って、第五演算部(Y5)は、これら4つの演算プロセスを組み込んでおり、かつオペレータによるこれら4つの解法のうちのいずれを選択するかの入力を可能とするものであることが好ましい。上記(iii)または(iv)の解法が選択された場合には、本発明の方法によって説明したような、操作の繰り返し回数Nが入力される。
 上述した第一から第五演算部でなされる演算のプロセスは、本発明の方法において説明したとおりである。
[0073]
 図1に示すように、第一から第五演算部は、一体のまたはそれぞれ独立したコンピュータとして構成するのが装置としての好ましい態様である。その場合、前記コンピュータに付帯する入力装置および出力装置によって入出力がなされ得るが、入力装置をさらにクライアント用のコンピュータとしてもよい。あるいは、当該装置は、予測結果の出力のための出力部(Y6)をさらに有してもよい。出力部(Y6)は、液晶ディスプレイモニターなどの各種画像表示装置、プリンターなどの各種印刷出力装置などであってよく、それぞれの出力装置には専用のコンピュータが付帯していてもよい。
[0074]
 本発明の装置は、上述した(I)~(VI)のステップをさらに有する実施形態で本発明の方法を行うために、上記[10]に示すように、さらに上記(y1)~(y6)の手段を有してもよい。その場合の装置の構成の一例を図2に示す。
[0075]
 (y1)の手段は、上記(I)のステップを行うための手段である。(y1)の手段には、解析の対象とすべきRNAの一次構造データが入力され、かつ該一次構造データから選択される長さLの対象区間Sが入力される。一次構造データは、必ずしもRNA全体のデータである必要はなく、予測の対象とすべき対象区間Sの一次構造データのみが入力されてもよい。また、必要なRNAの一次構造データは、自体の記憶装置に格納しておくだけでなく、オペレータの命令に応じて、CD-ROMや、外部データベース(外部コンピュータ)へのアクセスによって供給されるように構成してもよい。そのような一次構造データの供給源としては、公知の遺伝子データベース(例えば、米国のNational Center for Biotechnology Information (NCBI) が提供するGenBank)などが挙げられる。
 (y2)の手段は、上記(II)のステップを行うための手段である。(y2)の手段は、(y1)の手段で決定された対象区間S内の位置xの一端を位置x=x 1として、位置x 1と位置x 1+L1とを両端とする小区間をS(x 1)とし、上記の手段(Y1)で選択される区間が前記小区間S(x 1)であるとして、上記の手段(Y1)~(Y5)を用いて、該小区間S(x 1)での前記各特定二次構造が、非平衡状態にある時間的区間内の各時点で何程の確率で存在し得るかを決定する。
 (y3)の手段は、上記(III)のステップを行うための手段である。(y3)の手段は、前記位置x 1を、所定の微小長さrだけ他端側へ移動させた位置をx 2 =x 1+rとし、位置x 2と位置x 2+L1とを両端とする小区間をS(x 2)とし、前記小区間S(x 1)の代わりに該小区間S(x 2)を用いる以外は前記(y2)の手段と同様にして、該小区間S(x 2)での前記各特定二次構造が、非平衡状態にある時間的区間内の各時点で何程の確率で存在し得るかを決定する、という操作を、前記位置x=x 1 、x 2 に続いて、小区間の一端が対象区間の他端に達するまで、位置x=x 3 、...、x M について繰り返し行い、それによって、最初の小区間S(x 1)から最後の小区間S(x M)までの各小区間毎の前記各特定二次構造が、非平衡状態にある時間的区間内の各時点で何程の確率で存在し得るかを決定する。
 (y4)の手段は、上記(IV)のステップを行うための手段である。(y4)の手段には、存在を解析すべきモチーフの特徴パラメータ値が入力される。モチーフの特徴パラメータ値については、本発明の方法で説明したとおりである。(y4)の手段は、該特徴パラメータ値を有するモチーフM k が、任意の時間tにおいて、各小区間毎に決定された各特定二次構造中に存在するか否かを調べ、その結果と、時間tでの各特定二次構造の存在確率とに基づいて、各小区間において、時間tにおいて該モチーフが何程の蓋然性で存在するか(すなわち、上述した「第1の蓋然性」)を決定する。
 (y5)の手段は、上記(V)のステップを行うための手段である。(y5)の手段は、前記(y4)の手段を用いて調べた、各小区間における各特定二次構造における前記モチーフM の存在および非存在の結果に基づいて、各小区間を、(a)当該小区間においてモチーフM k を有する特定二次構造が少なくとも1つ存在する、または(b)当該小区間においてモチーフM k を有する特定二次構造が存在しない、のいずれであるかについて評価する。そして、上記(a)に評価される小区間の頻度を求めることによって、モチーフM k に関する第2の蓋然性を算出する。その詳細については、本発明の方法において上述したとおりである。
 (y6)の手段は、上記(VI)のステップを行うための手段であり、前記(y4)および(y5)の手段で求められた、第1の蓋然性と、第2の蓋然性とに基づいて、前記時間tにおける前記モチーフM k の存在を決定するための出力を行う。(VI)のステップを行うに際して、複数のモチーフに対して蓋然性を決定することが有用であり得る。その場合、複数のモチーフについての蓋然性の算出結果を同時に視覚化する機能を(y6)の手段が有することが有利である。そのため、(y6)の手段は、例えば、各モチーフに対応する各位置毎に第1の蓋然性と第2の蓋然性とを共に出力する、これら2つの蓋然性を1つのグラフ内で関連づけて視覚的に表示する、または、各モチーフに対応する各位置毎に第1の蓋然性と第2の蓋然性とから第3の蓋然性を算出してこれを出力する、などの機能を有し得る。これらの機能については、本発明のプログラムにおいて後述する。
[0076]
 次に、本発明の方法を実施するためのプログラムの構成を具体的に示す。
 当該プログラムは、上記した本発明の装置をコンピュータで構成した場合に、該コンピュータで本発明の方法を実行させるように構成されたプログラムである。当該プログラムは、1台のコンピュータ内で実行される態様だけではなく、手元のコンピュータ内にメインのプログラムが置かれ、他のコンピュータ内の演算プログラムへ初期値が送られ、演算結果がメインのプログラムに戻るというように、複数のコンピュータ内で実行される態様であってもよい。
[0077]
 当該プログラムは、図1にそのフローを示すとおり、基本的には、上述した本発明の方法の流れと全く同様であって、コンピュータを、少なくとも、上記(A)のステップを行う手段(P1)、上記(B)のステップを行う手段(P2)、上記(C)のステップを行う手段(P3)、上記(D)のステップを行う手段(P4)、および上記(E)のステップを行う手段(P5)として機能させるように構成される。
[0078]
 手段(P1)では、予測の対象とすべき長さL1のRNA塩基配列、温度など、予測に必要な初期値が入力される。RNA塩基配列は、オペレータによるデータ名の入力に基づいて、自体の記憶装置(HDD、CD-ROMなど)や外部データベースにアクセスしデータを取り込むプログラム構成であってもよい。これらの初期値については、本発明の方法において説明したとおりである。手段(P1)は、入力されたこれらの初期値に基づいて、予測の対象とする区間において、平衡状態における複数の特定二次構造ならびにそれらのギッブス自由エネルギーおよび存在確率を算出する演算を行う。その後の解析に利用すべき二次構造の選択を可能とするために、手段(P1)は、初期値として二次構造の選択の基準の入力を可能とするか、あるいは、予測の対象とされた区間での特定二次構造の候補ならびにそれらのギッブス自由エネルギーおよび/または存在確率を出力し、オペレータによる特定二次構造の選択の入力を可能とするものであることが好ましい。
 手段(P2)では、手段(P1)を用いて算出された特定二次構造に対して、それらの間での相互変換における遷移状態の二次構造が決定され、さらに各遷移状態の二次構造のギッブス自由エネルギーが求める。
 手段(P3)は、手段(P2)を用いて求められた各遷移状態の二次構造およびそれらのギッブス自由エネルギーから、各特定二次構造同士の間の相互変換の経路の活性化エネルギーを算出する。
 手段(P4)は、手段(P3)を用いて算出された各活性化エネルギーを基にして、平衡状態に達するまでの時間的区間内の各時点における、手段(P1)を用いて決定された各特定二次構造の存在確率を決定するための連立微分方程式を形成する。
 手段(P5)では、手段(P4)を用いて得られた連立微分方程式を数値的に解くために、初期条件、終了時点T fold、時間刻み幅Δtなど、解析に必要な値が適宜入力される。本発明の方法において説明したように、本発明では、連立微分方程式は上記(i)~(iv)のうちのいずれか1つの解法によって解かれ得る。従って、手段(P5)は、これら4つの演算プロセスを実装しており、かつオペレータによるこれら4つの解法のうちのいずれを選択するかの入力を可能とするものであることが好ましい。上記(iii)または(iv)の解法が選択された場合には、本発明の方法によって説明したような、操作の繰り返し回数Nが入力される。
 上述した手段(P1)から(P5)が行う演算のプロセスは、本発明の方法において説明したとおりである。
[0079]
 当該プログラムは、コンピュータをさらに、上記手段(P5)を用いて得られた各特定二次構造の存在確率の時間変化を出力するための手段(P6)として機能させるように構成されることが好ましい。手段(P6)は、各特定二次構造毎に、各時点における存在確率を数値として出力するものであってもよく、あるいは、横軸を時間とし、縦軸には、各時点毎に、各特定二次構造の存在確率を各々異なる色を有する存在確率に比例した長さの棒として、連続的かつ重ならないように一定の順序で並べる(その場合、全ての棒について合計した長さは常に一定である)などして、二次構造の挙動の視覚的な理解を可能とするように工夫してもよい。後者に類似した方法の一例は、図8~11に示されている。
[0080]
 本発明のプログラムは、上記[6]に示すように、コンピュータをさらに、上述した(I)~(VI)のステップをそれぞれ行うための(p1)~(p6)の手段として機能させるように構成することもできる。そのような構成によって、上述した(I)~(VI)のステップをさらに有する実施形態で本発明の方法を行うことが可能となる。その場合のプログラムの構成の一例を図2に示す。
[0081]
 (p1)の手段は、上記(I)のステップをコンピュータに行わせるための手段である。(p1)の手段では、解析の対象とすべきRNAの一次構造データが入力され、かつ該一次構造データから選択される長さLの対象区間Sが入力される。一次構造データは、必ずしもRNA全体のデータである必要はなく、予測の対象とすべき対象区間Sの一次構造データのみが入力されてもよい。また、必要なRNAの一次構造データは、コンピュータの記憶装置に格納しておいたものだけでなく、オペレータの命令に応じて、CD-ROMや、外部データベース(外部コンピュータ)へのアクセスによって供給されるように構成してもよい。そのような一次構造データの供給源としては、公知の遺伝子データベース(例えば、米国のNational Center for Biotechnology Information (NCBI) が提供するGenBank)などが挙げられる。
 (p2)の手段は、上記(II)のステップをコンピュータに行わせるための手段である。(p2)の手段は、(p1)の手段を用いて選択された対象区間S内の位置xの一端を位置x=x 1として、位置x 1と位置x 1+L1とを両端とする小区間をS(x 1)とし、上記の手段(P1)で選択される区間が前記小区間S(x 1)であるとして、上記の手段(P1)~(P5)を用いて、該小区間S(x 1)での前記各特定二次構造が、非平衡状態にある時間的区間内の各時点で何程の確率で存在し得るかを決定する。
 (p3)の手段は、上記(III)のステップをコンピュータに行わせるための手段である。(p3)の手段は、前記位置x 1を、所定の微小長さrだけ他端側へ移動させた位置をx 2 =x 1+rとし、位置x 2と位置x 2+L1とを両端とする小区間をS(x 2)とし、前記小区間S(x 1)の代わりに該小区間S(x 2)を用いる以外は前記(p2)の手段と同様にして、該小区間S(x 2)での前記各特定二次構造が、非平衡状態にある時間的区間内の各時点で何程の確率で存在し得るかを決定する、という操作を、前記位置x=x 1 、x 2 に続いて、小区間の一端が対象区間の他端に達するまで、位置x=x 3 、...、x M について繰り返し行い、それによって、最初の小区間S(x 1)から最後の小区間S(x M)までの各小区間毎の前記各特定二次構造が、非平衡状態にある時間的区間内の各時点で何程の確率で存在し得るかを決定する。
 (p4)の手段は、上記(IV)のステップをコンピュータに行わせるための手段である。(p4)の手段では、存在を解析すべきモチーフの特徴パラメータ値が入力される。モチーフの特徴パラメータ値については、本発明の方法で説明したとおりである。(p4)の手段は、該特徴パラメータ値を有するモチーフM が、任意の時間tにおいて、各小区間毎に決定された各特定二次構造中に存在するか否かを調べ、その結果と、時間tでの各特定二次構造の存在確率とに基づいて、各小区間において、時間tにおいて該モチーフが何程の蓋然性で存在するか(すなわち、上述した「第1の蓋然性」)を決定する。
 (p5)の手段は、上記(V)のステップをコンピュータに行わせるための手段である。(p5)の手段は、前記(p4)の手段を用いて調べた、各小区間における各特定二次構造における前記モチーフM k の存在および非存在の結果に基づいて、各小区間を、(a)当該小区間においてモチーフM k を有する特定二次構造が少なくとも1つ存在する、または(b)当該小区間においてモチーフM k を有する特定二次構造が存在しない、のいずれであるかについて評価する。そして、上記(a)に評価される小区間の頻度を求めることによって、モチーフM k に関する第2の蓋然性を算出する。その詳細については、本発明の方法において上述したとおりである。
 (p6)の手段は、上記(VI)のステップをコンピュータに行わせるための手段であり、前記(p4)および(p5)の手段で求められた、第1の蓋然性と、第2の蓋然性とに基づいて、前記時間tにおける前記モチーフM の存在を決定するための出力をコンピュータに行わせる。(VI)のステップを行うに際して、複数のモチーフに対して蓋然性を決定することが有用であり得る。その場合、複数のモチーフについての蓋然性の算出結果を同時に視覚化する機能を(p6)の手段が有することが有利である。そのため、(p6)の手段は、コンピュータを、例えば、各モチーフに対応する各位置毎に第1の蓋然性と第2の蓋然性とを共に出力するように機能させる、これら2つの蓋然性を1つのグラフ内で関連づけて視覚的に表示するように機能させる、または、各モチーフに対応する各位置毎に第1の蓋然性と第2の蓋然性とから第3の蓋然性を算出してこれを出力するように機能させる、などのプログラム構成であり得る。
 第1、第2の蓋然性の出力の仕方は、利用者が両方の値を考慮し易いように表示され得るものであればよい。例えば、第1、第2の蓋然性のそれぞれの値を、対象区間Sに対応させて単純に表形式で並列的に出力し、それぞれの値について、規定値よりも高い値にマーキングを施して示すという態様や、第1、第2の蓋然性のそれぞれに基準を設け、両方の基準を満たす位置をピックアップして示すという態様であってもよい。
 第3の蓋然性の出力の仕方は、第1、第2の蓋然性の出力をもって成しても良いし、利用者が第3の蓋然性の値を考慮し易いように表示され得るものであればよい。例えば、対象区間Sに対応させて単純に第3の蓋然性の値を順番にリスト出力する態様や、あるいは、対象区間S内の各位置を横軸、各位置に対応する第3の蓋然性を縦軸としたグラフで表示する態様であってもよい。
 また、第1の蓋然性と第2の蓋然性とを1つのグラフ内で関連づけて視覚的に表示させる方法としては、上記特許文献2に記載されたものと同様の方法が挙げられる。該方法は、対象区間Sが一方の軸に対応したグラフ中に、第1、第2の蓋然性を同時に描き、二次元、三次元のグラフィックなデザインにて出力するものである。
 そのような表示例を次に示す。なお、以下の説明は、「特定の時点tにおける」第1、第2の蓋然性を、グラフ中に表示する例であることに留意されたい。適切な特定の時点を選択して該特定の時点tでのモチーフの存在または非存在を決定することにより、生体内でRNAが実際に取っている可能性が高い二次構造であって、RNAの機能の発現に重要であると考えられる当該二次構造を予測することが可能となる。計算上は、通常、そのような特定の時点tをT foldとした上で、時点t=T foldにおける第1、第2の蓋然性をグラフ表示すればよい。
[0082]
 本例は、対象区間Sが一方の軸に対応したグラフ中に、棒グラフの形式で第1の蓋然性を表し、同時に、第2の蓋然性をも示す例である。
 例えば、(p4)の手段の演算において、ある小区間内のある位置に第1の蓋然性が与えられる場合、グラフ内の対象区間S上の対応する位置に棒グラフの棒(グラフ棒)を表示し、そのグラフ棒の長さで個々の第1の蓋然性の大きさを表す。
 小区間S(x 1)~S(x M)においてそれぞれ得られた1以上の位置を、グラフの対象区間Sに対応させて、第1の蓋然性の値をグラフ棒として描いていくと、位置によっては、グラフ棒同士の重なり合いが多数生じることになる。
 そこで、グラフ棒同士の重なり合いの程度が高いほど、すなわち、第2の蓋然性が高いほど、その重なり合った部分の色彩および/または模様を変えて表示する。
[0083]
 色彩の変化は、色相(単色光の波長に相当するもの)、彩度(あざやかさ、すなわち白みを帯びていない度合)、明度(明るさ、すなわち光の強弱)から選ばれる任意の要素を変化させればよい。例えば、第2の蓋然性がより高いグラフ棒の色をより濃く示したり、グラフ棒の色を青色から赤色へ変化させて示す、モニター画面では最重要部分を点滅させるなど、見る者が識別し易いように種々の工夫を凝らして表示してよい。表示装置が白黒の場合には、ハッチングやあみがけなどの種々の模様が識別には有用となる場合もある。
 このように表示することによって、第1の蓋然性については、グラフ棒の長さで直感でき、第2の蓋然性については、グラフ棒の色彩や模様の違いで直感できる。例えば、グラフ棒がより長くより濃い位置に、実際にモチーフが存在する蓋然性が高いと予想できる。
 このような予測結果を表示するための手段は、演算手段とは別の独立した表示専用のプログラムとしてもよい。
[0084]
 表示におけるオプションとしては、例えば、第2の蓋然性の値(グラフ棒同士の多重度)の閾値を変更できるようにし、その閾値以上の多重度を持っていないモチーフは表示しないようにしてもよい。それによって、色の薄い構造は表示されなくなる。
 また、第1の蓋然性の閾値を変更できるようにし、その閾値よりも低い第1の蓋然性を示す構造については表示しないようにしてもよい。第1の蓋然性を各モチーフのΔGを用いた関数f(ΔG)で表した場合には、関数f(ΔG)の値がより小さい方が、第1の蓋然性が高いことを意味する場合がある。その場合には、閾値よりも高いf(ΔG)の構造を表示させないようにする。これらの操作によって、存在する確からしさが低いモチーフを考慮に入れないようにすることができる。
 上記のような閾値による表示変更は、第3の蓋然性について行ってもよい。その操作によって、総合的な観点から存在する確からしさが低いモチーフを考慮に入れないようにすることができる。
[0085]
 モチーフの予測結果をスペクトルとして表示する方法は、本発明の有用性が顕著となる特徴的な表示方法である。しかし、グラフ棒が長くかつ充分に重なって示された部分が、同じモチーフの重なり合いに由来するものなのか、あるいは、互いに異なるモチーフの抽出結果がたまたま重なって表示されているのかがわからない場合がある。
 予測の信頼性を向上させるには、グラフ棒が多く重なって示された部分が、同じモチーフが重なり合ったことに由来することが明示される必要がある。以下に、その明示のための方法を二つ提案する。
[0086]
 1つの方法は、モチーフがステム・ループ構造や凹形構造などステムを持つ構造の場合に有用な方法であって、グラフに表示される全ての位置にあるモチーフについて塩基の状態を調べ、次の2種類の塩基(i)、(ii)を、互いに識別可能なように(例えば、色彩などを変えて)グラフに重ねて表示するという方法である。
 (i)常に塩基対(ステム)を形成しない塩基(フリー塩基と呼ぶ)。
 (ii)あるモチーフ中においては塩基対を形成しない塩基であるが、他のモチーフ中においては塩基対形成に参加している塩基(偽フリー塩基と呼ぶ)。
 例えば、モチーフが最も単純なステム・ループ構造であって、同じ構造のモチーフだけが重なり合って第2の蓋然性が高く表示された場合、その位置の中央には、フリー塩基(ループ部分に対応する)だけが表示され、偽フリー塩基は表示されない可能性が高いはずである(注意すべき例外は後述する)。このように、フリー塩基と偽フリー塩基とを参照すれば、第2の蓋然性の高さが、同じモチーフの重なり合いに由来するかどうかを判断する指標になる。
[0087]
 もう1つの方法は、グラフに表示された位置の配列を、もう一度構造予測にかけるなどしてそのモチーフ自体の安定性値ΔG値を出し、それをグラフに重ねて表示する方法である。凹形構造についてはモチーフ自体が安定性値を持たないので、間隙部分をはさむステム・ループ構造二本の安定性値ΔGの和をもってモチーフの安定性値としてもよい。
 第2の蓋然性の高さが、同じモチーフの重なり合いに由来するならば、その重なり合っている構造のΔG値はすべて同一となり、逆に棒グラフが複数の構造に由来していれば重なり合っている構造のΔG値はそれぞれ互いに異なっている可能性が非常に高い。
 また、この方法を用いて算出された各モチーフM k のΔGについて閾値を与え、画面上で不安定な構造を表示させないようにすることもできる。
[0088]
 次の二点に注意が必要である。
 一点は、対象区間Sの両端では、第2の蓋然性が低くなる(重なり合いの度合いが少なくなる)ことである。これはただ単に、端部では、小区間の重なり合う数が少なくなることに起因する。そのため、対象区間Sの端部に近い二次構造を議論するときには、第2の蓋然性の低さが、小区間の重なり合う数が少ないことによるのか、その構造の第2の蓋然性が低いためなのかを詳細に調べればよい。或いは、上述したように、第1の蓋然性の算出にあたって、最初からこうした構造を除いても良い。
 他の一点は、フリー塩基と偽フリー塩基とを参照する場合に、偽フリー塩基が確認されたからといって、その位置のモチーフが全て存在の蓋然性が低いとは限らないことである。例えば、一つの配列に対して十分な確からしさで存在するモチーフM 1が予測されていたとする。ところが、フレームの関係で該モチーフM 1を形成する配列の全長が含まれないフレームにおいては、該モチーフM 1を取れず、かわりに該配列が参加する形で別のモチーフM 2が予測される場合がある。こうした場合も最後の出力の結果としては、違う構造の重なりとして出力され、この二つの構造について互いに偽フリー塩基が出力される。しかしながら、この場合はモチーフM 1の蓋然性がより高いことは自明である。このようなM 2については上述の第一あるいは第2の蓋然性に対する閾値を設定することで除くことができる場合が多いが、より厳密には詳細を調べればよい。
[0089]
 また、RNA上の各位置に対して、その位置を含む全ての小区間にわたって、各小区間を代表するギッブス自由エネルギーを総和し、それをさらに、その位置を含む小区間の総数で割ることで得られた値を、その位置でのRNAの一般的な構造化の度合いとして表示してもよい。ここで、小区間を代表するギッブス自由エネルギーとは、その小区間において求められた各特定二次構造のギッブス自由エネルギーに、各特定二次構造の時間tでの存在確率をかけ、該小区間での全ての特定二次構造にわたって総和したものである。
 つまり、時間tでの或る小区間iについて、時間tにおける特定二次構造jの存在確率をP i,j(t)とし、その構造j(j=1,2,・・・,N)のギッブス自由エネルギーをΔG i,jとした場合、ある小区間iを代表するギッブス自由エネルギーΔG i(t)は、
[0090]
[数12]


[0091]
となる。さらに、或る位置xに関して、位置xを含む全ての小区間i(i=1,2,・・・,M)にわたって、各小区間を代表するギッブス自由エネルギーを総和し、それをさらに、その位置を含む小区間の総数Mで割ればよい。
[0092]
[数13]


[0093]
 こうして求められたΔG(x)(t)をスペクトル中に折れ線グラフとして表示してもよい。さらに、モチーフの存在する位置のΔG(x)(t)について閾値を与え、それを満たさない位置にかかるモチーフは画面上で表示させないようにしてもよい。
[0094]
 本プログラムの実施に際して、種々のパラメータの入力は、研究者が理解しやすい形で為されることが望ましい。さらに、その入力メソッドが、研究者の計算に対する理解を助けるように配置されていることが望ましい。そこで、たとえば計算の概要を示すようなアイコンと共に、それらの計算ステップで必要とされるパラメータを直観的に入力できるパネルが用意されていても良い。そのようなパネルの一例を図3に示す。
[0095]
 この計算は科学計算であるが、本プログラムはそうした科学計算に慣れていない実験系研究者をも対象にしている。そのため、計算の進捗にあたっては、上述したパラメータ設定のアイコンと対応するアイコンを表示し、計算がどのように行われているかを示して、研究者の計算に対する理解を助けてもよい。各ステップの進捗がバーで示されると同時に、終了したステップにはチェックマークが入り、計算の進行を体感できるようにしても良い。その一例を図4に示す。
[0096]
 また、本計算で行われる連立微分方程式の数値計算は、新たに示す「数値解析をイメージさせるアイコン」と共に進捗バーを表示することで、計算の意味と進捗を体感できるようにしても良い。その一例を図5に示す。

実施例

[0097]
実施例1: ADM mRNAの構造
 アドレノメデュリン(ADM)の5’UTRには特徴的なステム・ループ構造があり、その構造がADMの翻訳を制御しているといわれている(Brenet et al., Oncogene 25(49): 6510-9 (2006))。その構造は実験的に見つかっており、そうした有意な構造をコンピュータを用いて予測することは非常に意味がある。
[0098]
 そこで、本実施例1では、ADM mRNA配列(GenBank Accession No. NM_001124)について、3塩基(r=3)ずつずらしながら200塩基(L1=200)の小区間で区切って構造予測を行い、各小区間について一つまたは複数の構造予測結果(すなわち、特定二次構造)を得た。
 この構造予測結果から、妥当性が高いステム・ループを抽出した(抽出法の詳細については、上記Brenet et al., Oncogene 25(49): 6510-9 (2006)を参照)。そのとき、ステム・ループの評価には、各小区間での構造予測結果の存在確率(ポピュレーション)が必要になる。この「各小区間での構造予測結果の存在確率(ポピュレーション)」を以下の3つの計算方法で取得した。
 第1の計算方法は、S. Nakamura, J. Biochem., Apr 22 (2009) [Epub ahead of print], PMID: 19386779(http://jb.oxfordjournals.org/cgi/reprint/mvp064?ijkey=UAuE0vg4VWCEPzL&keytype=refから入手可能)に記載されている既存の方法である。すなわち、構造予測に用いたソフトウェアGCG-mfoldで出力されるエネルギー値を用いて、それぞれの構造群の相互変換が平衡状態にあるとしてポピュレーションを計算した。
 第2の計算方法では、用いるエネルギー値の算出を、新たにRNA構造予測システムGeno Poemics TM(GP)内部に搭載されているエネルギー計算ルーチンを用いて(計算手法は、以下の論文: Mathews et al, J. Mol. Biol. (1999) 288, 911-940.; Xia T et al, Biochemistry (1998) 37, 14719-14735.; Giese MR et al, Biochemistry (1998) 37, 1094-1100.; Draper DE et al, Methods in Enzymology (1995) 259, 281-242.;および、Serra MJ et al, Methods in Enzymology (1995) 259, 242-261.を参考にした)、構造群の相互変換が平衡状態にあるとしてポピュレーションを計算した。
 第3の計算方法では、用いるエネルギー値の算出を、新たにGP内部に搭載されているエネルギー計算ルーチンを用いて、構造群の相互変換が平衡状態に達していない(非平衡状態)としてポピュレーションを計算した。
 この非平衡状態の計算では、各構造のポピュレーションを計算するにあたり、初期値から微分方程式を数値的に解いた。初期値としては、各小区間において得られた構造予測結果が等しい比率で存在すると仮定することで得られる存在確率を用いた。さらに、時間刻み幅Δt=0.001秒、計算の終了時点T fold=3秒、温度37℃に設定の上、反応速度kを算出するために、アイリングの絶対反応速度式を使用した(プランク定数h(6.62621x10 -34 J s)、ボルツマン定数k B(1.3806504x10 -23 J K -1)、気体定数R(8.31447 J K -1 mol -1)、アイリング定数K(1.00))。
 微分方程式を数値的に解くにあたっては一次近似を採用した。さらに、ある構造間、例えば構造Aと構造Bとの間でk A→B×Δtとk B→A×Δtのどちらか一方が1.0以上であった場合には、その構造間での変換はΔtの時間内に既に平衡に達しているとした(その構造間では数値的に解かずにマクスウェル・ボルツマン分布によるポピュレーション配分を行った)。ここで、上記の関係が繋がって複数にわたる場合には、その複数の繋がった構造間で平衡に達していると考える(たとえばAとB、BとCに上記の関係があった場合には、A、BおよびCの3つの構造間で平衡に達していると考える。一方、AとB、CとDのそれぞれに上記の関係があった場合には、A-B間とC-D間のそれぞれの間が独立して平衡状態に達していると考える)。
[0099]
 上記3つの方法でのそれぞれの計算結果について、GPを用いてそれぞれスペクトル絵を作成した(第1から第3の計算方法が、それぞれスペクトル1~3に対応している)。簡単のため、VI=90%の閾値を適用した結果を図6に示す(すなわち、第3の蓋然性が90%以上のモチーフのみ表示している)。非平衡状態での計算から得られたスペクトル3については、時点t=T foldにおける結果を示している。ここで、上記Brenet et al., Oncogene 25(49): 6510-9 (2006)では、40-80付近のステム・ループ構造が重要であると記載されている。スペクトルを見ると、スペクトル1~3のいずれにおいてもそのステム・ループ(該図中では、楕円で囲まれている)が抽出されていることが分かる。しかし、スペクトル3については、出力されているステム・ループ数がそもそも少ない。これは、スペクトル3はスペクトル1、2に比べてノイズが少ないことを意味している。
[0100]
 以上より、本発明に従った非平衡状態での計算によって、mRNAの5’UTRを解析するにあたって、従来の方法と比べた改善が得られることが分かった。また、本実施例では注目すべきステム・ループ以外をノイズと称したが、このように厳しい条件下でも残ってくるピークは、何らかの未知の生体内機能を担っている有意味なものである可能性が考えられる。
[0101]
実施例2: CD83 mRNAの構造
 CD83のCDSには特徴的なステム・ループ構造があり、その構造がHuRと相互作用することでCD83の翻訳が制御されているといわれている(Prechtel et al., J. Biol. Chem. 281(16): 10912-25 (2006))。その構造は実験的に見つかっており、そうした有意な構造をコンピュータを用いて予測することは非常に意味がある。
[0102]
 そこで、本実施例2では、CD83 mRNA配列(GenBank Accession No. NM_004233)についても同様に、3塩基(r=3)ずつずらしながら200塩基(L1=200)の小区間で区切って構造予測を行い、各小区間について一つまたは複数の構造予測結果(すなわち、特定二次構造)を得た。この構造予測結果から、妥当性が高いステム・ループを抽出した(抽出法の詳細については上記Brenet et al., Oncogene 25(49): 6510-9 (2006)を参照)。そのとき、ステム・ループの評価には、各小区間での構造予測結果の存在確率(ポピュレーション)が必要になる。この「各小区間での構造予測結果の存在確率(ポピュレーション)」を次の3つの計算方法で取得した。
 第1の計算方法は、S. Nakamura, J. Biochem., Apr 22 (2009) [Epub ahead of print], PMID: 19386779(http://jb.oxfordjournals.org/cgi/reprint/mvp064?ijkey=UAuE0vg4VWCEPzL&keytype=refから入手可能)に記載されている既存の方法である。すなわち、構造予測に用いたソフトウェアGCG-mfoldで出力されるエネルギー値を用いて、それぞれの構造群の相互変換が平衡状態にあるとしてポピュレーションを計算した。
 第2の計算方法では、用いるエネルギー値の算出を、新たにRNA構造予測システムGeno Poemics TM(GP)内部に搭載されているエネルギー計算ルーチンを用いて(計算手法は、以下の論文: Mathews et al, J. Mol. Biol. (1999) 288, 911-940.; Xia T et al, Biochemistry (1998) 37, 14719-14735.; Giese MR et al, Biochemistry (1998) 37, 1094-1100.; Draper DE et al, Methods in Enzymology (1995) 259, 281-242.;および、Serra MJ et al, Methods in Enzymology (1995) 259, 242-261.を参考にした)、構造群の相互変換が平衡状態にあるとしてポピュレーションを計算した。
 第3の計算方法では、用いるエネルギー値の算出を、新たにGP内部に搭載されているエネルギー計算ルーチンを用いて、構造群の相互変換が平衡状態に達していない(非平衡状態)としてポピュレーションを計算した。この非平衡状態での計算は、上記実施例1の方法と全く同じ方法を採用した。
[0103]
 上記3つの計算方法でのそれぞれの計算結果について、GPを用いてそれぞれスペクトル絵を作成した(第1から第3の計算方法が、それぞれスペクトル4~6に対応している)。簡単のため、VI=90%の閾値を適用した結果を図7に示す(すなわち、第3の蓋然性が90%以上のモチーフのみ表示している)。非平衡状態での計算から得られたスペクトル6については、時点t=T foldにおける結果を示している。ここで、Prechtel AT et al., J Biol Chem. 2006; 281(16): 10912-25.では、540-580付近のステム・ループ構造がHuRとの相互作用に重要であると記載されている。また同文献では、440-510付近のステム・ループはHuRとの相互作用において重要ではないと記載されている。ここで、ノイズの減小と同時に、これら二つのステム・ループの差別化が目で見えれば、本発明に従った非平衡状態での計算の優位性が示される。
 スペクトルを見ると、どのスペクトルでも540-580付近のステム・ループ(該図中では、楕円で囲まれている)が抽出されていることが分かる。しかし、スペクトル4では440-510付近のステム・ループ(該図中では、×印が付されている)が出力されているのに対して、スペクトル5、6では、明確に440-510付近のステム・ループの出力がない。一方、ノイズについて見ると、スペクトル4とスペクトル5ではピークの数が等しく、ノイズの点では同程度である。しかし、スペクトル6においては、明確にピーク数が減少して4本となっている。これらのことから、本発明に従った非平衡状態での計算によって、重要な部分が差別化できると共に、ノイズの減少も達成できることが示された。
[0104]
 以上より、本発明に従った非平衡状態での計算によって、mRNAのCDSを解析するにあたって、従来の方法と比べた改善が得られることが分かった。また、本実施例では注目すべきステム・ループ以外をノイズと称したが、このように厳しい条件下でも残ってくるピークは、何らかの未知の生体内機能を担っている有意味なものである可能性が考えられる。
[0105]
実施例3: ADM mRNAの構造解析からのmRNAの挙動についての検討
 実施例1で示したように、ADM mRNAのモチーフについてスペクトルを得ることができた。本実施例3では、本発明に従った計算手法が、mRNAの構造的挙動を解析するにあたって、これまでにない視点を提供できることを示す。
[0106]
 本手法は、mRNAの各小区間について予測される安定な二次構造の存在確率を、連立微分方程式を解いていくことで算出する。すなわち、これは存在確率の時間変化を追っていることに他ならない。実施例1および2のように特定の時点でのスペクトル絵を生成するのではなく、本実施例では、個別の小区間毎の特定二次構造の存在確率の時間変化としてプロットした(その例を、図8~11に示す)。
[0107]
 このプロットの様子は下記の(1)~(4)のケースに大別された。計算にあたっては、T foldを3.0秒としている。
(1)T fold = 3.0秒以内で、ほぼ平衡状態に達する系(例:図8)。
(2)平衡状態に向かって進みつつも、T fold = 3.0秒では平衡に達していない系(例:図9)。図9に示した2つの例では、プロット線は、収束地点(平衡状態)に向かって矢印の方に進んではいるが、収束までにはあと数十秒~数分かかると思われる。
(3)0.001秒以内(最小刻み幅Δt内)に平衡状態に達してしまう系(例:図10)。構造間の活性化エネルギーが極めて低く、構造間の変換がとても速い(速度定数kが大きい)。図10に示した例では、活性化エネルギーは2.57-5.72 kcal/molである。
(4)全く平衡状態に向かって進んで行かない系(例:図11)。構造間の活性化エネルギーが高く、構造間の変換がとても遅い(速度定数kが小さい)。図11に示した例では、活性化エネルギーは23.8-51.18 kcal/molである。事実上平衡状態に達しない系であり、ポピュレーションは初期値に大きく引っ張られる。この系のポピュレーションを考えるに、各特定二次構造が平衡状態にあると仮定することは本質的に正しくないと考えられる。
 実際には、単一の小区間であってもこの4種類が混在している場合もある。例えば、ある小区間について6種の特定二次構造が予想されたとして、そのうち5種についてはそれぞれの構造間の変換が速くそれらの間で平衡に達するが、残りの1種がその他との構造変換を事実上起こさず全く平衡に向かって進まないという系もある。
[0108]
 ここで上記(4)のケースについて検討する。
 上述したように、この系を平衡状態にあると仮定してポピュレーション計算することは本質的に正しくない。系が平衡状態にあると仮定してポピュレーション計算してよいかどうかをチェックするためには、二次構造のエネルギー量だけでなく、二次構造間の遷移状態の二次構造を決定し、各相互変換における活性化エネルギーを評価する必要がある。それは、本発明の手法を用いることで達成することができる。
 この系でのポピュレーションは初期値(各構造の最初のポピュレーション)に大きく左右される。ここで、自然界の現象としてmRNAの挙動を見た場合に、毎回同じ初期値をとっている保証はない(というよりも同じ初期値をとっている理由がない)。とすると、毎回構造群が同じポピュレーションをとるわけではないということになる。さらに「毎回同じポピュレーションになるわけではない」部位の構造が生物学的に意味を持つ可能性は低いと考えられる。従って、このような構造は基本的に無視されるべきだろう。
 実は、実施例1および2において、本発明の方法が有意にノイズが少ないのは、半ば自動的必然的に「毎回同じポピュレーションになるわけではない部位の構造を重視しない」という機能が働いているからである。つまり、実施例1および2においては、予測構造のポピュレーションに対して同一の初期値を与えている。平衡状態であれば予測構造のうち最安定なものが極めて優勢なポピュレーションをとるのであろうが(その場合、その最安定構造に含まれるモチーフは妥当性が高いものとして出力される)、平衡に達しない状態では、どの構造も優勢にならないのである。結果として、妥当な閾値(例えば、VI=90%)を設定すると、この部位の構造は完全に無視される。もちろん、平衡状態に達しない構造の中に共通モチーフがあって、それが有意なポピュレーションをとるのであれば、その共通モチーフは出力される。
 GPシステムの考え方は、構造が一つに定まらないものを、複数の構造がどれくらいのポピュレーションで存在するかというプロファイルを定量的に算出し、そこから共通モチーフを抽出する(そしてそのモチーフの妥当性を算出する)ことにある。本実施例の手法は、そのプロファイル自体が定まるかどうかをチェックする方法として用いることができる。
[0109]
 振り返るに、あるRNAの安定構造を論じるにあたって、従来、その安定構造のギッブス自由エネルギー量だけを評価することが多かった。そこに、複数の安定構造が同時に存在し得るという事実を正面から捉え、複数の安定構造のポピュレーションのプロファイルから、共通モチーフの妥当性を定量的に算出および画像化するのを可能にしたのがGPシステムである。さらに、各安定構造間の遷移状態の構造を決定し、各相互変換の活性化エネルギーを評価することで、複数の安定構造のポピュレーションプロファイルが定まるかどうかを算出する(定まるのであればより確からしい値を算出する)のが本発明の方法である。
 つまり、mRNAの構造を調べるということについては、その安定性だけを論じる時代、安定性だけから妥当性まで算出する時代、そして安定構造だけで論じるのではなく安定構造間の遷移状態の構造をも加味して妥当性を算出する時代、と順を追って進んできたとも言える。本発明の方法は、前記最後の、安定構造間の遷移状態の構造をも加味した妥当性の算出を実現できる現時点で唯一無二のものである。

産業上の利用可能性

[0110]
 本発明の方法により、従来の方法では扱うことができなかった、RNAの二次構造がマクスウェル・ボルツマンの分布則に従う存在確率の分布を示すよりも前の時点における該RNAの二次構造の動態を予測することが可能となる。それにより、生体内でRNAが実際に取っている二次構造をより反映する形で、RNAの二次構造を予測することができる。本発明はまた、該方法をコンピュータに実行させるためのプログラム(該プログラムを記録した記録媒体を含む)、および該方法を実施するための装置を提供する。
 本出願は日本で出願された特願2009-262258(出願日:2009年11月17日)を基礎としており、その内容は本明細書に全て包含されるものである。

請求の範囲

[請求項1]
 RNAの二次構造の動態を予測する方法であって、
 (A)RNAの塩基配列から長さL1を有する区間を選択し、該区間において、平衡状態における2種以上の特定二次構造を決定し、かつ、各特定二次構造のギッブス自由エネルギーと平衡状態での存在確率とを算出するステップと、
 (B)前記特定二次構造同士の間の相互変換における遷移状態の二次構造を決定するとともに、各遷移状態の二次構造のギッブス自由エネルギーを求めるステップと、
 (C)各遷移状態の二次構造およびそれらのギッブス自由エネルギーから、各特定二次構造同士の間の相互変換の経路の活性化エネルギーを算出するステップと、
 (D)各活性化エネルギーを基にして、平衡状態に達するまでの時間的区間内の各時点における各特定二次構造の存在確率を決定するための連立微分方程式を形成するステップと、
 (E)下記(i)~(iv)のうちのいずれか1つの解法によって、各特定二次構造が、前記時間的区間内の各時点において何程の確率で存在し得るかを決定するステップと、
を有することを特徴とする、RNAの二次構造の動態を予測する方法。
  (i)前記平衡状態における各特定二次構造の存在確率の逆数を、前記時間的区間の或る時点t 0における各特定二次構造の存在比率とし、そこから導かれる該時点t 0での各特定二次構造の存在確率を初期値として前記連立微分方程式を解き、求める存在確率とする解法。
  (ii)前記時間的区間の或る時点t 0において各特定二次構造の存在比率が全て等しいものであるとして、各特定二次構造の該時点t 0での存在確率を決定し、これを初期値として、前記連立微分方程式を解き、求める存在確率とする解法。
  (iii)前記時間的区間の或る時点t 0における各特定二次構造の存在確率を、それぞれ正規分布として与えるものとし、それぞれの正規分布は、前記平衡状態における各特定二次構造の存在確率の逆数を、前記時間的区間の或る時点t 0における各特定二次構造の存在比率とすることで導かれる存在確率を平均値とし、かつ、該平均値と、各特定二次構造が等しい比率で存在した場合の存在確率との差の絶対値を標準偏差とするものとし、
 各特定二次構造毎に前記正規分布に従った1個の乱数を与えて対応する存在比率を求め、そこから導かれる各存在確率を初期値として前記連立微分方程式を解き、該初期値を用いた場合の存在確率を求める、という操作をN回繰り返した上で、各時点毎に、該N回の操作毎の存在確率の総和をNで除して平均し、求める存在確率とする解法。
  (iv)前記時間的区間の或る時点t 0における各特定二次構造の存在確率を、それぞれ正規分布として与えるものとし、それぞれの正規分布は、前記平衡状態における各特定二次構造が等しい比率で存在した場合の存在確率を平均値とし、かつ、該平均値と、前記平衡状態における各特定二次構造の存在確率の逆数を、前記時間的区間の或る時点t 0における各特定二次構造の存在比率とすることで導かれる存在確率との差の絶対値を標準偏差とするものとし、
 各特定二次構造毎に前記正規分布に従った1個の乱数を与えて対応する存在比率を求め、そこから導かれる各存在確率を初期値とし、前記連立微分方程式を解き、該初期値を用いた場合の存在確率を求める、という操作をN回繰り返した上で、各時点毎に、該N回の操作毎の存在確率の総和をNで除して平均し、求める存在確率とする解法。
[請求項2]
 以下の(I)~(VI)のステップをさらに有するものである、請求項1に記載の方法。
 (I)前記長さL1よりも大きい長さLを有する対象区間Sを与えるステップ。
 (II)前記対象区間S内の位置xの一端を位置x=x 1として、位置x 1と位置x 1+L1とを両端とする小区間をS(x 1)とし、
 上記のステップ(A)で選択される区間が前記小区間S(x 1)であるとして、上記のステップ(A)~(E)を行って、該小区間S(x 1)での前記各特定二次構造が、前記時間的区間内の各時点で何程の確率で存在し得るかを決定するステップ。
 (III)前記位置x 1を、所定の微小長さrだけ他端側へ移動させた位置をx 2 =x 1+rとし、位置x 2と位置x 2+L1とを両端とする小区間をS(x 2)とし、前記小区間S(x 1)の代わりに該小区間S(x 2)を用いる以外は前記(II)のステップと同様にして、該小区間S(x 2)での前記各特定二次構造が、前記時間的区間内の各時点で何程の確率で存在し得るかを決定する、という操作を、前記位置x=x 1 、x 2に続いて、小区間の一端が対象区間Sの他端に達するまで、位置x=x 3 、...、x Mについて繰り返し行い、それによって、最初の小区間S(x 1)から最後の小区間S(x M)までの各小区間毎の前記各特定二次構造が、前記時間的区間内の各時点で何程の確率で存在し得るかを決定するステップ。
 (IV)前記(II)~(III)のステップの結果に基づいて、前記の小区間S(x 1)からS(x M)のうちの任意の小区間S(x i)において、任意の時間tにおける前記各特定二次構造の存在確率を決定し、
 前記各特定二次構造において、位置および構造を表す所定の特徴パラメータ値を有する任意のモチーフM k が存在するか否かを調べ、その結果と上記で得られた前記時間tでの各特定二次構造の存在確率とに基づいて、該小区間S(x i)内において、前記時間tにおいて該モチーフM k が何程の蓋然性で存在するかを決定し、これを、小区間S(X i)における該モチーフM k に関する第1の蓋然性とする、
 という操作を、前記の小区間S(x 1)からS(x M)の全てに対して行なうステップ。
 (V)前記(IV)のステップにおいて各小区間内の各特定二次構造における前記モチーフM k の存在および非存在を調べた結果に基づいて、各小区間において、モチーフM k を有する特定二次構造が少なくとも1つ存在するか否かを決定する、
 という操作を、前記の小区間S(x 1)からS(x M)の全てに対して行い、そのような特定二次構造が少なくとも1つ存在する小区間の頻度を、モチーフM k に関する第2の蓋然性とするステップ。
 (VI)前記(IV)および(V)のステップで求められた、第1の蓋然性と、第2の蓋然性とに基づいて、前記時間tにおける前記モチーフM k の存在を決定するステップ。
[請求項3]
 前記モチーフM k が、ステム・ループ構造を含むものである、請求項2に記載の方法。
[請求項4]
 前記モチーフM k が、2つのステム・ループ構造の間にあって他の塩基配列と相互作用していない間隙部分と、その両側にある前記ステム・ループ構造のそれぞれの内側のステム部分とを有してなる凹形構造を含むものである、請求項2または3に記載の方法。
[請求項5]
 RNAの二次構造の動態を予測するためのコンピュータプログラムであって、
 コンピュータを、
 (P1)RNAの塩基配列から長さL1を有する区間を選択し、該区間において、平衡状態における2種以上の特定二次構造を決定し、かつ各特定二次構造のギッブス自由エネルギーおよび平衡状態での存在確率を算出する手段、
 (P2)前記特定二次構造同士の間の相互変換における遷移状態の二次構造を決定するとともに、各遷移状態の二次構造のギッブス自由エネルギーを求める手段、
 (P3)各遷移状態の二次構造およびそれらのギッブス自由エネルギーから、各特定二次構造同士の間の相互変換の経路の活性化エネルギーを算出する手段、
 (P4)各活性化エネルギーを基にして、平衡状態に達するまでの時間的区間内の各時点における各特定二次構造の存在確率を決定するための連立微分方程式を形成する手段、
 (P5)下記(i)~(iv)のうちのいずれか1つの解法によって、各特定二次構造が、前記時間的区間内の各時点において何程の確率で存在し得るかを決定する手段、
として機能させるための前記コンピュータプログラム。
  (i)前記平衡状態における各特定二次構造の存在確率の逆数を、前記時間的区間の或る時点t 0における各特定二次構造の存在比率とし、そこから導かれる該時点t 0での各特定二次構造の存在確率を初期値として前記連立微分方程式を解き、求める存在確率とする解法。
  (ii)前記時間的区間の或る時点t 0において各特定二次構造の存在比率が全て等しいものであるとして、各特定二次構造の該時点t 0での存在確率を決定し、これを初期値として、前記連立微分方程式を解き、求める存在確率とする解法。
  (iii)前記時間的区間の或る時点t 0における各特定二次構造の存在確率を、それぞれ正規分布として与えるものとし、それぞれの正規分布は、前記平衡状態における各特定二次構造の存在確率の逆数を、前記時間的区間の或る時点t 0における各特定二次構造の存在比率とすることで導かれる存在確率を平均値とし、かつ、該平均値と、各特定二次構造が等しい比率で存在した場合の存在確率との差の絶対値を標準偏差とするものとし、
 各特定二次構造毎に前記正規分布に従った1個の乱数を与えて対応する存在比率を求め、そこから導かれる各存在確率を初期値として前記連立微分方程式を解き、該初期値を用いた場合の存在確率を求める、という操作をN回繰り返した上で、各時点毎に、該N回の操作毎の存在確率の総和をNで除して平均し、求める存在確率とする解法。
  (iv)前記時間的区間の或る時点t 0における各特定二次構造の存在確率を、それぞれ正規分布として与えるものとし、それぞれの正規分布は、前記平衡状態における各特定二次構造が等しい比率で存在した場合の存在確率を平均値とし、かつ、該平均値と、前記平衡状態における各特定二次構造の存在確率の逆数を、前記時間的区間の或る時点t 0における各特定二次構造の存在比率とすることで導かれる存在確率との差の絶対値を標準偏差とするものとし、
 各特定二次構造毎に前記正規分布に従った1個の乱数を与えて対応する存在比率を求め、そこから導かれる各存在確率を初期値とし、前記連立微分方程式を解き、該初期値を用いた場合の存在確率を求める、という操作をN回繰り返した上で、各時点毎に、該N回の操作毎の存在確率の総和をNで除して平均し、求める存在確率とする解法。
[請求項6]
 コンピュータを、さらに以下の(p1)~(p6)の手段として機能させるためのものである、請求項5に記載のコンピュータプログラム。
 (p1)前記長さL1よりも大きい長さLを有する対象区間Sを入力する手段。
 (p2)前記対象区間S内の位置xの一端を位置x=x 1として、位置x 1と位置x 1+L1とを両端とする小区間をS(x 1)とし、
 上記の手段(P1)で選択される区間が前記小区間S(x 1)であるとして、上記の手段(P1)~(P5)を用いて、該小区間S(x 1)での前記各特定二次構造が、前記時間的区間内の各時点で何程の確率で存在し得るかを決定する手段。
 (p3)前記位置x 1を、所定の微小長さrだけ他端側へ移動させた位置をx 2 =x 1+rとし、位置x 2と位置x 2+L1とを両端とする小区間をS(x 2)とし、前記小区間S(x 1)の代わりに該小区間S(x 2)を用いる以外は前記(p2)の手段と同様にして、該小区間S(x 2)での前記各特定二次構造が、前記時間的区間内の各時点で何程の確率で存在し得るかを決定する、という操作を、前記位置x=x 1 、x 2 に続いて、小区間の一端が対象区間Sの他端に達するまで、位置x=x 3 、...、x M について繰り返し行い、それによって、最初の小区間S(x 1)から最後の小区間S(x M)までの各小区間毎の前記各特定二次構造が、前記時間的区間内の各時点で何程の確率で存在し得るかを決定する手段。
 (p4)前記(p2)~(p3)の手段を用いて得られた結果に基づいて、前記の小区間S(x 1)からS(x M)のうちの任意の小区間S(x i)において、任意の時間tにおける前記各特定二次構造の存在確率を決定し、
 前記各特定二次構造において、位置および構造を表す所定の特徴パラメータ値を有する任意のモチーフM k が存在するか否かを調べ、その結果と上記で得られた前記時間tでの各特定二次構造の存在確率とに基づいて、該小区間S(x i)内において、前記時間tにおいて該モチーフM k が何程の蓋然性で存在するかを決定し、これを、小区間S(X i)における該モチーフM k に関する第1の蓋然性とする、
 という操作を、前記の小区間S(x 1)からS(x M)の全てに対して行なう手段。
 (p5)前記(p4)の手段を用いて各小区間内の各特定二次構造における前記モチーフM k の存在および非存在を調べた結果に基づいて、各小区間において、モチーフM k を有する特定二次構造が少なくとも1つ存在するか否かを決定する、
 という操作を、前記の小区間S(x 1)からS(x M)の全てに対して行い、そのような特定二次構造が少なくとも1つ存在する小区間の頻度を、モチーフM k に関する第2の蓋然性とする手段。
 (p6)前記(p4)および(p5)の手段で求められた、第1の蓋然性と、第2の蓋然性とに基づいて、前記時間tにおける前記モチーフM k の存在を決定するための出力を行う手段。
[請求項7]
 前記モチーフM k が、ステム・ループ構造を含むものである、請求項6に記載のコンピュータプログラム。
[請求項8]
 前記モチーフM k が、2つのステム・ループ構造の間にあって他の塩基配列と相互作用していない間隙部分と、その両側にある前記ステム・ループ構造のそれぞれの内側のステム部分とを有してなる凹形構造を含むものである、請求項6または7に記載のコンピュータプログラム。
[請求項9]
 RNAの二次構造の動態を予測するための装置であって、
 (Y1)RNAの塩基配列から長さL1を有する区間を選択し、該区間において、平衡状態における2種以上の特定二次構造およびそれらの存在確率を算出する第一演算部と、
 (Y2)前記特定二次構造同士の間の相互変換における遷移状態の二次構造を決定するとともに、各遷移状態の二次構造のギッブス自由エネルギーを求める第二演算部と、
 (Y3)各遷移状態の二次構造およびそれらのギッブス自由エネルギーから、各特定二次構造同士の間の相互変換の経路の活性化エネルギーを算出する第三演算部と、
 (Y4)各活性化エネルギーを基にして、平衡状態に達するまでの時間的区間内の各時点における各特定二次構造の存在確率を決定するための連立微分方程式を形成する第四演算部と、
 (Y5)下記(i)~(iv)のうちのいずれか1つの解法によって、各特定二次構造が、前記時間的区間内の各時点において何程の確率で存在し得るかを決定する第五演算部と、
を有する、前記装置。
  (i)前記平衡状態における各特定二次構造の存在確率の逆数を、前記時間的区間の或る時点t 0における各特定二次構造の存在比率とし、そこから導かれる該時点t 0での各特定二次構造の存在確率を初期値として前記連立微分方程式を解き、求める存在確率とする解法。
  (ii)前記時間的区間の或る時点t 0において各特定二次構造の存在比率が全て等しいものであるとして、各特定二次構造の該時点t 0での存在確率を決定し、これを初期値として、前記連立微分方程式を解き、求める存在確率とする解法。
  (iii)前記時間的区間の或る時点t 0における各特定二次構造の存在確率を、それぞれ正規分布として与えるものとし、それぞれの正規分布は、前記平衡状態における各特定二次構造の存在確率の逆数を、前記時間的区間の或る時点t 0における各特定二次構造の存在比率とすることで導かれる存在確率を平均値とし、かつ、該平均値と、各特定二次構造が等しい比率で存在した場合の存在確率との差の絶対値を標準偏差とするものとし、
 各特定二次構造毎に前記正規分布に従った1個の乱数を与えて対応する存在比率を求め、そこから導かれる各存在確率を初期値として前記連立微分方程式を解き、該初期値を用いた場合の存在確率を求める、という操作をN回繰り返した上で、各時点毎に、該N回の操作毎の存在確率の総和をNで除して平均し、求める存在確率とする解法。
  (iv)前記時間的区間の或る時点t 0における各特定二次構造の存在確率を、それぞれ正規分布として与えるものとし、それぞれの正規分布は、前記平衡状態における各特定二次構造が等しい比率で存在した場合の存在確率を平均値とし、かつ、該平均値と、前記平衡状態における各特定二次構造の存在確率の逆数を、前記時間的区間の或る時点t 0における各特定二次構造の存在比率とすることで導かれる存在確率との差の絶対値を標準偏差とするものとし、
 各特定二次構造毎に前記正規分布に従った1個の乱数を与えて対応する存在比率を求め、そこから導かれる各存在確率を初期値とし、前記連立微分方程式を解き、該初期値を用いた場合の存在確率を求める、という操作をN回繰り返した上で、各時点毎に、該N回の操作毎の存在確率の総和をNで除して平均し、求める存在確率とする解法。
[請求項10]
 さらに以下の(y1)~(y6)の手段を有するものである、請求項9に記載の装置。
 (y1)前記長さL1よりも大きい長さLを有する対象区間Sを入力する手段。
 (y2)前記対象区間S内の位置xの一端を位置x=x 1として、位置x 1と位置x 1+L1とを両端とする小区間をS(x 1)とし、
 上記の手段(Y1)で選択される区間が前記小区間S(x 1)であるとして、上記の手段(Y1)~(Y5)を用いて、該小区間S(x 1)での前記各特定二次構造が、前記時間的区間内の各時点で何程の確率で存在し得るかを決定する手段。
 (y3)前記位置x 1を、所定の微小長さrだけ他端側へ移動させた位置をx 2 =x 1+rとし、位置x 2と位置x 2+L1とを両端とする小区間をS(x 2)とし、前記小区間S(x 1)の代わりに該小区間S(x 2)を用いる以外は前記(y2)の手段と同様にして、該小区間S(x 2)での前記各特定二次構造が、前記時間的区間内の各時点で何程の確率で存在し得るかを決定する、という操作を、前記位置x=x 1 、x 2 に続いて、小区間の一端が対象区間Sの他端に達するまで、位置x=x 3 、...、x M について繰り返し行い、それによって、最初の小区間S(x 1)から最後の小区間S(x M)までの各小区間毎の前記各特定二次構造が、前記時間的区間内の各時点で何程の確率で存在し得るかを決定する手段。
 (y4)前記(y2)~(y3)の手段によって決定された結果に基づいて、前記の小区間S(x 1)からS(x M)のうちの任意の小区間S(x i)において、任意の時間tにおける前記各特定二次構造の存在確率を決定し、
 前記各特定二次構造において、位置および構造を表す所定の特徴パラメータ値を有する任意のモチーフM k が存在するか否かを調べ、その結果と上記で得られた前記時間tでの各特定二次構造の存在確率とに基づいて、該小区間S(x i)内において、前記時間tにおいて該モチーフM k が何程の蓋然性で存在するかを決定し、これを、小区間S(X i)における該モチーフM k に関する第1の蓋然性とする、
 という操作を、前記の小区間S(x 1)からS(x M)の全てに対して行なう手段。
 (y5)前記(y4)の手段を用いて各小区間内の各特定二次構造における前記モチーフM の存在および非存在を調べた結果に基づいて、各小区間において、モチーフM を有する特定二次構造が少なくとも1つ存在するか否かを決定する、
 という操作を、前記の小区間S(x 1)からS(x M)の全てに対して行い、そのような特定二次構造が少なくとも1つ存在する小区間の頻度を、モチーフM k に関する第2の蓋然性とする手段。
 (y6)前記(y4)および(y5)の手段で求められた、第1の蓋然性と、第2の蓋然性とに基づいて、前記時間tにおける前記モチーフM k の存在を決定するための出力を行う手段。
[請求項11]
 前記モチーフM k が、ステム・ループ構造を含むものである、請求項10に記載の装置。
[請求項12]
 前記モチーフM k が、2つのステム・ループ構造の間にあって他の塩基配列と相互作用していない間隙部分と、その両側にある前記ステム・ループ構造のそれぞれの内側のステム部分とを有してなる凹形構造を含むものである、請求項10または11に記載の装置。

図面

[ 図 1]

[ 図 2]

[ 図 3]

[ 図 4]

[ 図 5]

[ 図 6]

[ 図 7]

[ 図 8]

[ 図 9]

[ 図 10]

[ 図 11]