JP4459597B2 - 遺伝子発現作用推定方法、遺伝子発現作用推定装置及び遺伝子発現作用推定プログラム - Google Patents

Description

本発明は、複数の遺伝子間において、個々の遺伝子が他の遺伝子から時間的に影響を受けて変化する、その変化作用量を推測する方法に関するものである。

従来、こうした複数の遺伝子群において、個々の遺伝子が他の遺伝子から影響を受ける様子、つまり遺伝子相互作用を推定する方法として、閾値検定、統計的手法を用いたＢａｙｅｓｉａｎネットワーク、二項関係を用いたブーリアンネットワーク、更には微分方程式を用いる方法等、各種の技術がある。
なかでも同値類のフィードバック制御が同定可能な微分方程式モデルを使用した推定方法がよく知られている。

なお、出願人は、出願時点で先行技術文献を知得していない。

従来の方法は、推定のために各遺伝子の作用係数に初期値を与えて推定を開始するが、その初期値の与え方に関して適切な方法を示していないので、最適方法を推定するまでに非常に時間がかかるという課題がある。
また計算量が膨大になり、または繰返し推定の中で数値の発散を防ぐために制限を加える必要がある、そのため遺伝子数を制限する必要がある、などの課題もある。

この発明は上記の課題を解決するためになされたもので、適切な係数の初期値が容易に得られ、かつ計算量を抑えた遺伝子発現相互作用推定方法を得る。

この発明に係る遺伝子発現作用推定方法は、遺伝子モデルを作成し、操作されていない遺伝子群からなる野生株と特定遺伝子を破壊した破壊株の時間列データを基に遺伝子発現作用を推定する方法において、
遺伝子モデルに対し、野生株と特定の破壊株との時間列データを同時に満足する発現係数を推定する手順を野生株に含まれる推定対象個別遺伝子の数だけ繰返すループとして持つ個別遺伝子への発現係数推定ステップと、
特定の破壊株との推定を終えると、次の破壊株を取上げて個別遺伝子への発現係数推定ステップを繰返すことを、所定数の破壊株について順次行う、繰り返しループとして持つ遺伝子発現係数行列推定ステップとを備えて、
上記２つのステップを経て、遺伝子発現係数行列の推定初期値を得るようにした。

野生株と破壊株とのタイムコース値の組合わせから遺伝子間相互作用行列を初期値推定するようにしたので、推定収束までの時間が短縮され、かつ安定で妥当な推定ができる効果がある。

実施の形態１．
遺伝子の相互作用、または他の遺伝子から受ける影響度を推定するために、最小二乗法による準最適化手法を適用して初期値を得て、またその初期値に対して線形微分方程式モデルによる準最適化を適用して、制約を少なくして推定所要時間を短縮する方法を説明する。
本実施の形態において、線形微分方程式モデルとして次の（１）式を設定する。

ここでｘは遺伝子発現量ベクトル、Ａは各遺伝子間の相互作用を示す係数行列、ｗは不確定要素によるバイアス係数ベクトルである。不確定要素に関しては、時間の多項式関数とする場合もあるが、ここではバイアスモデルとして記述する。また上記（１）式から、最小二乗法では遺伝子発現量ベクトルを（１）式の両辺を積分した次の（２）式とする。

ｘ₀ は遺伝子発現量初期値（初期値ベクトル）、ｔは時間を示す。
また最小二乗法は、遺伝子発現量ベクトル（遺伝子ベクトル）ｘの一つの遺伝子に対し、時間で変化するタイムコースデータを基に次の（３）式で各遺伝子Ｙが変化すると設定する。そうしてその時間変化において、他の遺伝子からの相互作用を示す行列Ｂと、Ｗの推定を行う。
Ｙ＝Ｙ₀ ＋Ｂ・Ｘ＋Ｗ・ｔ （３）
ここでＹは対象遺伝子群中における一つの遺伝子のタイムコースデータであり、Ｘは個別の遺伝子Ｙの集合体としての遺伝子群で、遺伝子発現量ｘの時間積分値である。またＹ₀ はその初期値である。
こうして設定したモデルに基づく推定方法をフローにしたものを図１と図２に示す。またこの動作を展開する計算回路構成の例を図３に示す。図３において、推定を行うための回路構成として、実際のデータを読込む入力装置５１、本装置の解析結果である、遺伝子間の影響度をパスウェイとして得て、このパスウェイを表示する表示装置５２と、この入出力間で本処理を行うハードウェア構成の遺伝子発現作用推定装置からなる。この遺伝子発現作用推定装置は、中央処理装置を有する計算機で実現できるものである。遺伝子発現作用推定装置の詳細構成として、以下の構成要素を持っている。入力側から順に、タイムコース・データ入力処理部４１、これで得たデータを補間する補間処理部４２、補間結果から正規化する正規化処理部４３、モデルの係数初期値を推定する初期値推定部４４、得られた初期値に基づいて相互作用係数を推定する相互作用係数推定処理部４５、更に最適化処理をする係数最適化処理部４６、入力のタイムコース・データを記憶するタイムコース・データ記憶部４８、モデルの数式を表現する遺伝子発現モデルを記憶する解析手法記憶部４７、得られた遺伝子間の影響を図式化するパスウェイ描画処理部４９で、それぞれ対応プログラムを持ったハードウェアである。あるいは、上記各部は、プログラムのみにより実現されてもよいし、ハードウェアのみにより実現されてもよいし、プログラムとハードウェアとの組み合わせにより実現されてもよい。

この計算機による推定方法を図１と図２に基づいて説明する。
先ずタイムコースデータである実験値により着目する個別の遺伝子Ｙの時間変化データを読み込む。図４は実験等で得られる観測値であり、図において、左端の列は説明のためのインデックスであり、図では１から２７まで記されている。またこの図では説明を簡単にするために遺伝子を４個（Ｇ０１，Ｇ０２，Ｇ０３，Ｇ０４）の例を示している。インデックス５ないし７は操作されていない自然の個別遺伝子からなる遺伝子群を野生株として、その野生株の時間変化データを示している。そして左端の列が時刻、この例では３つの時刻（１．１，１．２，１．３）を表し、次の列（１．１１，１．２１，１．３１）は、遺伝子Ｇ０１の各左側時刻におけるタイムコース値で、次の１．１２ないし１．３２は、遺伝子Ｇ０２の同じくサンプル時刻のタイムコース値、同様に１．１３ないし１．３３はＧ０３のタイムコース値、１．１４ないし１．３４はＧ０４のタイムコース値である。更にインデックス１０ないし１２は、インデックス８に示すように、Ｇ０１遺伝子を破壊した破壊株の３つの時刻（１．１，１．２，１．３）におけるタイムコース値で、２．１１，２．１２，２．１３，２．１４の各下への列が、それぞれ遺伝子Ｇ０１，Ｇ０２，Ｇ０３，Ｇ０４の時刻変化値である。以下、インデックス１５ないし１７は遺伝子Ｇ０２を破壊した破壊株、２０ないし２２はＧ０３を破壊した破壊株、２５ないし２７はＧ０４を破壊した破壊株のタイムコース値である。これをタイムコース・データ入力処理部４１がＳ１１で行う。

Ｓ１２において、計算機はタイムコースデータを計算機処理のために絶対値が大きいデータで正規化し、補間処理部４２が行うＳ１３で、対数値を基に図４の経過時間の中間値をスプライン（平滑化）補間処理する。正規化処理部４３が行うＳ１４で、補間後の対数値データを元の物理量に戻す。
Ｓ１５では、サンプル時刻毎に各遺伝子のそれまでのデータを積分して行列Ｘと、またベクトルＹを生成する。Ｓ１５ではこの他に、生成された行列Ｘにおいて、推定しない相互作用の要素はゼロに固定し、またＸ行列、Ｙベクトルの有効桁を計算の精度と安定化のために、調整する。

本実施の形態で特徴的なループを説明する。Ｓ１６では、Ｓ２１で囲まれる大きなループを形成している。この大きなループは、野生株の時間変化と、１つの破壊株（図４では、例えばインデックス８に注記したＧ０１を破壊した破壊株で、そのタイムコース値はインデックス１０ないし１２に示される）をループの１つとし、これを野生株と次の破壊株（Ｇ０２を破壊した破壊株）の組み、というように全組合わせによるＢとＷを推定して行くことを表している。
Ｓ１７とＳ２０で囲まれる小さなループは、例えばＧ０１を破壊した破壊株と野生株のタイムコース値からＢとＷを得ることを表している。つまりＳ１７で、先ず最初の野生株の遺伝子Ｇ０１を取り上げて、図４のインデックス５ないし７のタイムコース値と１０ないし１２のタイムコース値とから、Ｇ０１に注目して、Ｓ１８により最小二乗法により暫定的に（３）式の行列ＢとベクトルＷ内の遺伝子Ｇ０１に関する相互作用の係数を推定する。こうして図４の例では、インデックス５ないし７のタイムコース値と１０ないし１２のタイムコース値による推定では、遺伝子数が４つあるので４つのＢとＷが推定される。つまりこの小さなループでは、図４のインデックス１０ないし１２を用いた推定で、一回目のループではＧ０１に着目した推定をし、次の小さなループの２回目ではＧ０２に着目した推定を行い、三回目の小さなループにおいてはＧ０３に着目して推定し、最後に四回目の小さなループではｇ０４に着目した推定を行う。その結果、Ｓ１９で、得られたＢとＷ内の計数を用いて、（１）式の行列Ａとｗの各値を暫定的に推定する。

次いでＳ１６に戻って、タイムコースデータ中のその他の破壊株、図４の例ではインデックス１３，１８，２３に記された他の破壊株について、Ｓ１７ないしＳ１９のループを繰り返す。即ち大きなループでは、順次他の破壊株と野生株とを用いた推定に移ることを表している。この大きなループを含む二重ループの処理を行うのが初期値推定部４４である。
こうして順次破壊株を取り出して、それぞれ破壊株と野生株と同時に満足するＢとＷの推定を行う。その結果、全破壊株が全遺伝子に与える相互作用が推定できる。このとき遺伝子Ｙは、ループ中に含まれる遺伝子の発現量の時系列データを格納したベクトルとなる。またＸは遺伝子Ｙベクトルのタグ時刻における発現量の積分値を格納している。
Ｓ２１で全破壊株が及ぼす遺伝子のＢとＷの係数の推定が終わると、Ｓ２３で相互作用係数推定処理部４５が、破壊株の影響がない係数について、残りのＡ，ｗの係数を野生株のタイムコースデータを入力して最小二乗法により推定する。こうして係数行列Ａとベクトルｗが暫定的に準最適化用の初期値として確定する。

Ｓ２４以降は、係数最適化処理部４６が最終的な係数を確定するための処理を行う。Ｓ２４では、Ｓ２３で暫定的に確定したｘのモデルについて、Ａとｗの各要素の係数を微小変動させた、２個以上の幾つかのモデル、例えば５個のモデルを線形微分方程式モデルとして生成する。
Ｓ２５では、確定までの演算時間を制限するために、最大試行回数を設定し、Ｓ２６では、Ｓ２４で生成した５個のモデルのうち、最初のモデルを取り出して、Ｓ２７とＳ２８の処理をする。Ｓ２７では、モデルの微分方程式を数値積分して、各個体の各遺伝子のタイムコース値を求める。
Ｓ２８で、上記積分で求めたタイムコース値と、対応するＳ１１での入力実験値とを比較して、推定相対誤差を求める。
Ｓ２９からＳ２６に戻って、Ｓ２４で生成した５個全てのモデルについて、Ｓ２７とＳ２８の処理を繰り返し、各モデルの推定相対誤差を求める。

Ｓ３０で所定の誤差以下のモデルがあれば、そのモデルのＡ推定行列とｗベクトルであるとして推定を終える。
Ｓ３０でどのモデルも所定の誤差以下でなければ、Ｓ３１で５個のモデルのうち、最も誤差の少ないモデルを個体として選択し、再び５個のモデルを生成する。そしてＳ３２でその５個のモデルについてそれぞれ行列Ａとベクトルｗの異なる係数を微小変動させる。Ｓ３３では、最大試行回数内で所定の回数Ｓ３０で誤差が閾値以下にならなかった場合に、Ｓ３２の微小変動を小さくする。
Ｓ２５とＳ３４とで囲まれる試行ループにより、Ｓ３０で所定の誤差以下のモデルが得られる可能性が高く、誤差が少ないモデルが定まれば、そのモデルが表す行列Ａとベクトルｗが最終推定結果となる。

図５は、このように推定してＳ３０で誤差が少ないとして得た最終のモデルを推定値で表し、図４のタイムコース実験値とを比較した図である。
また図６は、推定の結果得られる、遺伝子間の相互作用、つまり他の遺伝子への影響関係を図式化したものである。この図はパスウェイ描画処理部４９で得られる。図１と図２による推定の結果、Ｂが定まり、遺伝子Ｃ３は、遺伝子Ｂ６からプラスの相互作用を受け、Ｂ７からマイナスの作用を受け、Ｃ２と組になってＣ１にプラス作用を与えることが推定される。
こうして遺伝子を微分方程式でモデル化し、野生株と破壊株とのタイムコース値の組合わせから微分方程式モデルの遺伝子間相互作用行列を初期値推定するようにしたので、推定収束までのシミュレーション回数が短縮され、かつ安定で妥当な推定ができる。
なお、図２と図３の各ステップを手順として持つプログラムを用意して、記録媒体またはオンラインロードにより他の計算機にインストールし、そのプログラムを計算機の中央処理装置により実行させて上記解析を行わせるようにしてもよい。

この発明の実施の形態１における行列とベクトル推定処理フローを示す図である。 実施の形態１における行列とベクトル推定処理フローの続きを示す図である。 実施の形態１における行列とベクトル推定処理を行う回路構成を示す図である。 実施の形態１における実験で得られるタイムコース値の例を示す図である。 実験値のタイムコース値と実施の形態１で得られた推定値に基づいたタイムコース値との比較を示す図である。 実施の形態１での推定結果による遺伝子間の相互作用を示すパスウェイ図である。

符号の説明

Ｓ１６ 特定破壊株を取上げるステップ、Ｓ１７ 推定対象の個別遺伝子を選択するステップ、Ｓ１８ 最小二乗法による行列Ｂ、ベクトルＷの係数推定ステップ、Ｓ１９ 行列Ａ、ベクトルｗ内の要素への登録ステップ、Ｓ２０ 特定破壊株の中で未推定の遺伝子があるかを調べるステップ、Ｓ２１ 未推定の破壊株があるかを調べるステップ、Ｓ１７，Ｓ２０ 特定破壊株の中で個別遺伝子の係数を順次遺伝子毎に推定する小さなループ、Ｓ１６，Ｓ２０ 破壊株毎に順次推定する大きなループ、Ｓ２２ 全破壊株について推定が終わったかを確認するステップ、Ｓ２３ 行列Ａとｗの未推定係数を推定するステップ、Ｓ２４ 推定対象の複数の線形微分方程式モデル作成ステップ、Ｓ２６ Ｓ２４で作成した一つの線形微分方程式モデルを選択するステップ、Ｓ２７ モデルの積分により各遺伝子の時間列値を得るステップ、Ｓ２８ 時間積分値と入力遺伝子時間列値との誤差算出ステップ、Ｓ２９ Ｓ２４で作成した線形微分方程式モデルを順次取上げて行き、全モデルが終わったかを確認するステップ、Ｓ３０ Ｓ２４で作成したモデルの中で所定の誤差以下のモデルがあるかを調べるステップ、Ｓ３１ 誤差最小のモデルに対してＳ２４で作成した数のモデル数を用意するステップ、Ｓ３２ 行列Ａ、ベクトルｗの要素の値をモデル毎に違えて変化させるステップ、Ｓ３３ 変動率を小さくするステップ、Ｓ３４ 設定した最大試行回数を超えていないかを確認するステップ、４１ タイムコース・データ入力処理部、４２ 補間処理部、４３ 正規化処理部、４４ 初期値推定部、４５ 相互作用係数推定処理部、４６ 係数最適化処理部、４７ 解析手法記憶部、４８ タイムコース・データ記憶部、４９ パスウェイ描画処理部。

Claims (7)

1. 遺伝子群において個々の遺伝子の発現が他の遺伝子から受ける作用を数１１で表わされる線形微分方程式によりモデル化した遺伝子モデルのパラメータである行列Ａ及びベクトルｗを算出する遺伝子発現作用推定装置において、
   上記遺伝子群からなる野生株について、上記遺伝子群に含まれる複数の遺伝子の発現量を時系列で観測した時間変化データを読み込み、上記複数の遺伝子のうちいずれかの遺伝子を破壊した破壊株について、上記複数の遺伝子の発現量を時系列で観測した時間変化データを読み込む入力装置と、
   上記入力装置が読み込んだ上記野生株と上記破壊株とについての時間変化データに基づいて、中央処理装置が、数１２を満たす行列Ｂ及びベクトルＷを最小二乗法により算出し、算出した行列Ｂの要素のうち、上記破壊株で破壊された遺伝子からの作用を表わす要素を、上記行列Ａの対応する要素の推定初期値とする相互作用係数初期値推定部とを備えたことを特徴とする遺伝子発現作用推定装置。
   

   

   ただし、ｘは上記複数の遺伝子の発現量を要素とするベクトル、Ａは上記複数の遺伝子間の作用を表わす係数を要素とする行列、ｗは不確定要素を表わすバイアス係数を要素とするベクトルを表わす。
   

   

   ただし、Ｙは上記複数の遺伝子のうち一の遺伝子の時系列の発現量を要素とするベクトル、Ｙ _０ はＹの初期値、Ｂは上記一の遺伝子が他の遺伝子から受ける作用を表わす係数を要素とする行列、Ｘは上記複数の遺伝子の発現量を時間積分した積分値を要素とする行列、Ｗはベクトル、ｔは時間を表わす。
2. 上記遺伝子発現作用推定装置は、更に、相互作用係数推定部を備え、
   上記相互作用係数推定部は、中央処理装置が、上記行列Ａの要素のうち、上記破壊株で破壊された遺伝子からの作用を表わす要素を、上記相互作用係数初期値推定部が推定した推定初期値に固定し、上記入力装置が読み込んだ上記野生株についての時間変化データに基づいて、行列Ａ及びベクトルｗを最小二乗法により算出することを特徴とする請求項１に記載の遺伝子発現作用推定装置。
3. 上記遺伝子発現作用推定装置は、更に、相互作用係数準最適化処理部を備え、
   上記相互作用係数準最適化処理部は、中央処理装置が、上記相互作用係数推定部が推定した行列Ａ及びベクトルｗを初期値とし、上記行列Ａ及びベクトルｗの要素を変動させて、誤差を算出し、算出した誤差が所定の誤差以下である場合に、上記行列Ａ及びベクトルｗを決定することを特徴とする請求項１または請求項２に記載の遺伝子発現作用推定装置。
4. 計算機が実行することにより、上記計算機が請求項１乃至請求項３のいずれかに記載の遺伝子発現作用推定装置として機能することを特徴とする遺伝子発現作用推定プログラム。
5. 遺伝子発現作用推定装置が、遺伝子群において個々の遺伝子の発現が他の遺伝子から受ける作用を数１３で表わされる線形微分方程式によりモデル化した遺伝子モデルのパラメータである行列Ａ及びベクトルｗを算出する遺伝子発現作用推定方法において、
   入力装置が、上記遺伝子群からなる野生株について、上記遺伝子群に含まれる複数の遺伝子の発現量を時系列で観測した時間変化データを読み込み、
   入力装置が、上記複数の遺伝子のうちいずれかの遺伝子を破壊した破壊株について、上記複数の遺伝子の発現量を時系列で観測した時間変化データを読み込み、
   中央処理装置が、読み込んだ上記野生株と上記破壊株とについての時間変化データに基づいて、数１４を満たす行列Ｂ及びベクトルＷを最小二乗法により算出し、
   中央処理装置が、算出した行列Ｂの要素のうち、上記破壊株で破壊された遺伝子からの作用を表わす要素を、上記行列Ａの対応する要素の推定初期値とすることを特徴とする遺伝子発現作用推定方法。
   

   

   ただし、ｘは上記複数の遺伝子の発現量を要素とするベクトル、Ａは上記複数の遺伝子間の作用を表わす係数を要素とする行列、ｗは不確定要素を表わすバイアス係数を要素とするベクトルを表わす。
   

   

   ただし、Ｙは上記複数の遺伝子のうち一の遺伝子の時系列の発現量を要素とするベクトル、Ｙ _０ はＹの初期値、Ｂは上記一の遺伝子が他の遺伝子から受ける作用を表わす係数を要素とする行列、Ｘは上記複数の遺伝子の発現量を時間積分した積分値を要素とする行列、Ｗはベクトル、ｔは時間を表わす。
6. 上記遺伝子発現作用推定方法は、更に、
   中央処理装置が、上記行列Ａの要素のうち、上記破壊株で破壊された遺伝子からの作用を表わす要素を、推定した推定初期値に固定し、読み込んだ上記野生株についての時間変化データに基づいて、行列Ａ及びベクトルｗを最小二乗法により算出することを特徴とする請求項５に記載の遺伝子発現作用推定方法。
7. 上記遺伝子発現作用推定方法は、更に、
   中央処理装置が、推定した行列Ａ及びベクトルｗを初期値とし、上記行列Ａ及びベクトルｗの要素を変動させて、誤差を算出し、算出した誤差が所定の誤差以下である場合に、上記行列Ａ及びベクトルｗの要素を決定することを特徴とする請求項５または請求項６に記載の遺伝子発現作用推定方法。

Images