遺伝子相互作用推定方法、遺伝子相互作用推定プログラム、及び遺伝子相互作用推定装置、並びに、結合サイト推定方法、結合サイト推定プログラム、及び結合サイト推定装置
【課題】 ＤＮＡマイクロアレイを用いながらも、高精度に遺伝子相互作用の推定を行うことができる遺伝子相互作用推定方法を提供する。
【解決手段】 遺伝子相互作用推定方法は、発現データの相関を利用した第１のフィルタリング処理、結合サイトを推定する全く新規の手法である第２のフィルタリング処理、及び実験操作を加えた遺伝子発現データを利用した第３のフィルタリング処理による３つの独立した推定方法を融合させて行う。第２のフィルタリング処理としては、遺伝子ａと、当該遺伝子ａと相関係数の絶対値が大きい上位所定個数の遺伝子群Ｂの各遺伝子ｂとのペアについて、これらのウィンドウに含まれる近接した開始位置をもつ所定長ｋの文字列同士の類似度の最大値をその位置の類似度と定義したときの各位置に対する類似度の合計であるウィンドウ類似度を算出する。
【解決手段】 遺伝子相互作用推定方法は、発現データの相関を利用した第１のフィルタリング処理、結合サイトを推定する全く新規の手法である第２のフィルタリング処理、及び実験操作を加えた遺伝子発現データを利用した第３のフィルタリング処理による３つの独立した推定方法を融合させて行う。第２のフィルタリング処理としては、遺伝子ａと、当該遺伝子ａと相関係数の絶対値が大きい上位所定個数の遺伝子群Ｂの各遺伝子ｂとのペアについて、これらのウィンドウに含まれる近接した開始位置をもつ所定長ｋの文字列同士の類似度の最大値をその位置の類似度と定義したときの各位置に対する類似度の合計であるウィンドウ類似度を算出する。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、遺伝子が因子の制御を受けるとき、その因子の制御を受ける他の遺伝子に与える影響を推定する遺伝子相互作用推定方法、遺伝子相互作用推定プログラム、及び遺伝子相互作用推定装置に関する。
【背景技術】
【０００２】
近年のゲノム科学の進展にともない、ある環境における多数の遺伝子の発現強度を定量的に測定する技術であるＤＮＡマイクロアレイ技術により、多数の遺伝子発現データを同時に取得することが可能となっている。また、遺伝子の機能を解明する研究も、個々の遺伝子について解析する段階から、遺伝子の相互作用によって細胞内の活動がどのように行われるかを知る段階に移りつつある。
【０００３】
ここで、遺伝子は、例えば図２７に示すように、上流部分に制御領域を有し、そこに制御因子が結合することにより、遺伝子の発現が促進されたり抑制されたりする。また、遺伝子の発現によって因子（タンパク質）が生成されると、その因子が他の遺伝子の発現に影響を与えるが、各遺伝子が、自身を含めた他の遺伝子にいかなる影響を与えるかを推定する問題を遺伝子相互作用推定問題という。
【０００４】
かかる遺伝子相互作用推定の手法としては、発現データのみを利用したモデル推定によるものが提案されている。すなわち、ＤＮＡマイクロアレイによる発現データのみを利用し、発現データを最もよく説明可能なモデルを推定するものである。
【０００５】
まず、この種の遺伝子相互作用推定の手法としては、ブール代数を用いたブーリアンネットワークが提案されている（例えば、非特許文献１等参照。）。ブーリアンネットワークは、同期して動作する論理回路のようなモデルである。例えば図２８に示すモデルの場合には、１個のノード（頂点）が１個の遺伝子に対応し、各頂点は、“０（発現していない）”又は“１（発現している）”の状態をとる。また、各頂点の状態は、単位時刻毎に同期して変化し、遺伝子発現の制御規則は、ブール関数の形式で表現される。かかる手法においては、発現の時系列データを与え、データとモデルとの不一致を最小にするようなブーリアンネットワークを推定する問題を解くことになる。しかしながら、かかる手法においては、定量的な発現強度を扱うことができないという問題がある。
【０００６】
また、この種の遺伝子相互作用推定の手法としては、ベイジアンネットワークを用いる手法が提案されている（例えば、非特許文献２等参照。）。確率変数を頂点で表し、因果関係や相関関係といった依存する関係を有する変数間にリンクを張ったグラフ構造による確率モデルをグラフィカルモデルというが、このうち、例えば図２９に示すように、特にリンクの向きが因果関係の方向であり、このリンクを辿ったパスが循環しないものをベイジアンネットワークという。ベイジアンネットワークは、確率変数間の定性的な依存関係をグラフ構造によって表し、その確率変数間に定義される条件付き確率によって確率変数間の定量的な依存関係を表す。かかる手法においては、各遺伝子の発現強度を確率変数とし、その間の依存関係をベイジアンネットワークとして推定することになる。しかしながら、かかる手法においては、同じデータを説明可能な複数のモデルが存在してしまうという問題がある。
【０００７】
さらに、この種の遺伝子相互作用推定の手法としては、グラフィカルガウシアンモデリングを用いる手法が提案されている（例えば、非特許文献３等参照。）。相関係数によるクラスタリングの問題点は、直接の相互作用によって相関がみられる場合と、いくつかの他の遺伝子の発現を介した相互作用によって相関がみられる場合や、ある共通の遺伝子の発現に影響を受ける２つの遺伝子にみられる擬似相関の場合とを区別できないことにある。ここで、ｎ個の遺伝子がある環境での発現パターンは、ｎ次元の正規分布からのサンプルであると仮定すると、偏相関係数が“０”となるか否かによって擬似相関を排除することができる。これがグラフィカルガウシアンモデリングを用いる手法である。しかしながら、かかる手法においては、発現の挙動が非常に類似した遺伝子が存在した場合には、相関係数行列が正則でなくなり、偏相関係数を算出することができなくなる。そのため、かかる手法においては、前処理としてクラスタリングを行う必要があり、遺伝子グループ単位でしか依存関係を推定することができないという問題がある。
【０００８】
さらにまた、この種の遺伝子相互作用推定の手法としては、動的システムモデルを用いる手法が提案されている（例えば、非特許文献４等参照。）。この手法は、微分方程式等によるシステムモデルを仮定し、そのパラメータ推定を発現データとの一致度を目的関数とする最適化問題として定式化し、遺伝的アルゴリズムによって解くものである。しかしながら、この手法においては、計算コストが高く、数千もの遺伝子ネットワークを一度に解くことができず、分割や階層化等の処理を行う必要があるという問題がある。また、かかる手法においては、ノイズをも説明するモデルを出力してしまう可能性があり、いかにしてノイズを排除するかという問題もある。
【０００９】
このように、ＤＮＡマイクロアレイは、遺伝子相互作用の解析には有力な手段ではあるものの、上述した手法は、いずれも発現データのみを利用することから、ＤＮＡマイクロアレイによる発現データに含まれる多くのノイズの問題、実験コストの観点から数千もの遺伝子数に比べて数百程度の少数の実験データしか得られないという問題、転写因子と転写因子結合サイトとの関係を得ることができないという問題があり、推定精度が高くない。
【００１０】
この問題を解決するために、ＤＮＡマイクロアレイによる発現データと遺伝子の上流部分にある制御領域に含まれるパターンの解析とを組み合わせることにより、転写因子とこれによって制御される遺伝子群の関係を導き出す手法が提案されている（例えば、非特許文献５及び非特許文献６並びに特許文献１等参照。）。
【００１１】
具体的には、非特許文献５には、系列の発現データを用いて、転写因子の発現パターンと共通のモチーフを制御領域に有する遺伝子の発現パターンの合成との相関をみることにより、複数の転写因子の組み合わせによる影響を除外し、転写因子と結合サイトとの関係を推定する手法が開示されている。また、非特許文献６には、転写因子の発現パターンと各遺伝子の発現パターンとの相関によって遺伝子をクラスタリングし、その後モチーフをみつける手法が開示されている。さらに、特許文献１には、遺伝子発現量データと遺伝子配列データとに基づいてコンピュータによって遺伝子間の制御関係を推定する手法が開示されている。すなわち、この手法においては、遺伝子発現量の時系列データから、任意に早発現遺伝子と遅発現遺伝子とを選択し、その２遺伝子に関する相互相関関数を求め、その無相関検定を行って帰無仮説が棄却された場合には、その２遺伝子を記録する。そして、この手法は、早発現遺伝子と遅発現遺伝子とを選択し尽くすまで、この手順を繰り返し行い、その後、転写因子の上流配列に当該転写因子のコンセンサス制御配列の数が多い遺伝子と、先に記録された早発現遺伝子と遅発現遺伝子とをそれぞれ照合し、一致した場合に、その２遺伝子を記録して、これら２遺伝子には制御関係があると推定するものである。
【００１２】
【非特許文献１】阿久津達也著、「遺伝子ネットワークの推定アルゴリズム」、数理科学１９９９年６月号、株式会社サイエンス社、１９９９年、Ｎｏ．４３２、ｐ．４０−４６
【非特許文献２】N. Friedman, M. Linial, I. Nachman and D. Pe’er著, 「Using Bayesian networksto analyzeexpression data」, Journal of Computational Biology,vol. 7, no. 3/4, 2000年, p. 601-620
【非特許文献３】H. Toh, K. Horimoto著, 「Inference of a genetic network by a conbined approach of cluster analysis and graphicalGaussian modeling」, Bioinformatics, vol. 18, no. 2,2002年, p. 287-297
【非特許文献４】岡本正宏著、「Ｓ−ｓｙｓｔｅｍによる遺伝子の相互作用推定」、ゲノム情報生物学、株式会社中山書店、２０００年、ｐ．１６５−１８８
【非特許文献５】K. Birnbaum, P. N. Benfeyand D. E. Shasha著, 「cisElement/Transcription Factor Analysis (cis/TF): AMethod for Discovering Transcription Factor/cisElement Relationships」, Genome Res. 2001 September,vol. 11, no. 9, 2001年, p. 1567-1573
【非特許文献６】Z. Zhu, Y. Pilpel and G. M. Church著, 「Computational Identification ofTranscription Factor Binding Sites via a Transcription-factor-centricClustering (TFCC) Algorithm」, J. Mol. Biol. vol. 318,2002年, p. 7181
【特許文献１】特開２００３−１４１１２３号公報
【発明の開示】
【発明が解決しようとする課題】
【００１３】
ところで、各遺伝子の発現レベルは、様々な要因によって変化するが、転写因子の発現強度の変化は、ノイズ等のバックグラウンドの変化に対して極めて小さい。これに対して、上述した非特許文献５及び非特許文献６並びに特許文献１に記載された従来の手法は、いずれも発現データ間の相関を求めることが基本であることから、ノイズの影響を大きく受けるという問題があった。
【００１４】
また、多くの転写因子は、常に発現しており、その後の作用によって制御されていることが知られているが、従来の手法は、かかる転写因子に対しては適用することができないという問題もあった。
【００１５】
さらに、特開２０００−３４２２５７号公報には、エンハンサー又はプロモーターの内部構造を探究することにより、遺伝子の転写制御をより明確にすることができる手法が開示されている。すなわち、この手法においては、遺伝子の制御因子の結合部位構造を推定するにあたって、遺伝子をコードしているコード領域の上流域又は下流域にある転写因子となるタンパク質が結合するエンハンサー又はプロモーター領域内の制御因子の結合部位の制御構造を推定したい遺伝子を設定し、前記エンハンサー又はプロモーター領域において前記遺伝子に関与する制御因子、若しくは仮説的に導入する制御因子の結合部位について、その遺伝子座及びその他の遺伝子の発現の要因をパラメータとする計算モデルを構築し、この構築した計算モデルによって遺伝子の転写量を算出し、パラメータ探索アルゴリズムを使用して、実験で得られている前記設定された遺伝子の発現が得られるように計算モデルのパラメータを探索し、前記エンハンサー又はプロモーターのマイクロストラクチャーを推定している。しかしながら、この手法は、遺伝子の制御関係そのものを推定するものではない。
【００１６】
本発明は、このような実情に鑑みてなされたものであり、ＤＮＡマイクロアレイを用いながらも、高精度に遺伝子相互作用の推定を行うことができる遺伝子相互作用推定方法、遺伝子相互作用推定プログラム、及び遺伝子相互作用推定装置を提供することを目的とする。
【課題を解決するための手段】
【００１７】
本願出願人は、遺伝子が因子の制御を受けるとき、その因子の制御を受ける他の遺伝子に与える影響を推定する遺伝子相互作用推定問題について鋭意研究を重ねた結果、発現データの相関を利用した第１のフィルタリング処理、結合サイトを推定する全く新規の手法である第２のフィルタリング処理、及び実験操作を加えた遺伝子発現データを利用した第３のフィルタリング処理による３つの独立した推定方法を融合させる、という発想をなすに至り、さらに、第２のフィルタリング処理として行うべき具体的処理として極めて斬新な手法を考案することにより、本発明をなすに至った。
【００１８】
すなわち、上述した目的を達成する本発明にかかる遺伝子相互作用推定方法は、遺伝子ａが因子Ｆの制御を受けるとき、前記因子Ｆの制御を受ける他の遺伝子ｂに与える影響を推定する遺伝子相互作用推定方法において、ウィンドウ類似度算出手段により、前記遺伝子ａと、当該遺伝子ａと相関係数の絶対値が大きい上位所定個数の遺伝子群Ｂの各遺伝子ｂとのペアについて、これら２つの遺伝子ａ，ｂの制御領域からそれぞれ所定長ｌの連続領域をウィンドウとして切り出し、これらのウィンドウに含まれる近接した開始位置をもつ所定長ｋの文字列同士の類似度の最大値をその位置の類似度と定義したときの各位置に対する類似度の合計であるウィンドウ類似度を算出するウィンドウ類似度算出工程を備えることを特徴としている。
【００１９】
このような本発明にかかる遺伝子相互作用推定方法においては、ウィンドウ類似度という概念を取り入れることにより、遺伝子ａ，ｂの制御領域全体ではなく局所的な類似度を発見することが可能となる。
【００２０】
ここで、前記ウィンドウ類似度算出工程では、統計的に特異的に出現頻度が高い文字列については、それに対する類似度を０として前記ウィンドウ類似度が算出される。すなわち、本願出願人は、制御領域におけるオリゴヌクレオチドの出現頻度の解析を行い、統計的に出現頻度が特異的に高いオリゴヌクレオチドは結合文字列になる可能性が低いことを新たに発見したが、本発明にかかる遺伝子相互作用推定方法においては、この事実を領域の類似度算出に反映させることにより、結合サイト以外で類似領域が発見される可能性を低減することができる。
【００２１】
また、本発明にかかる遺伝子相互作用推定方法は、領域探索手段により、前記遺伝子ａの制御領域のうち、最大値ＭａｘａＢ［ｉ］が位置ｉにおける前記ウィンドウ類似度の分布において統計的に特異的に大きく、且つ、当該遺伝子ａの制御領域の位置ｉの変化に対する前記遺伝子群Ｂの各遺伝子ｂの前記ウィンドウ類似度の最大値Ｍａｘａｂ［ｉ］をみたときに極大値を与える位置であるピーク位置を含むような領域Ｒを探索する領域探索工程を備えることを特徴としている。すなわち、本発明にかかる遺伝子相互作用推定方法においては、遺伝子ａ，ｂに共通の結合サイトが存在するとした場合には、このピーク位置からのウィンドウ内に存在する可能性が高いことから、かかるピーク位置を含むような領域Ｒを探索することにより、結合サイトの候補を絞り込むことが可能となる。
【００２２】
さらに、本発明にかかる遺伝子相互作用推定方法は、遺伝子群抽出手段により、前記領域探索工程にて求められた前記領域Ｒにピーク位置が存在し、且つ、類似度が所定値よりも高い遺伝子群ＢＲを求めて抽出する遺伝子群抽出工程を備えることを特徴としている。これにより、本発明にかかる遺伝子相互作用推定方法においては、遺伝子群Ｂのうち、制御領域内に所定値以上の類似度を有する遺伝子群ＢＲを抽出することができる。
【００２３】
このように、本発明にかかる遺伝子相互作用推定方法においては、第２のフィルタリング処理としてこれらの工程を行い、特定のパターンを仮定せずに、指定した遺伝子ａと類似した制御領域を有する遺伝子を網羅的に探索する。このようにして求められた類似度が高い領域は、共通の因子が結合するサイトであるものと推定することができる。
【００２４】
さらにまた、本発明にかかる遺伝子相互作用推定方法は、前記ウィンドウ類似度算出工程に先だって、遺伝子群選択手段により、前記遺伝子ａと相関係数の絶対値が大きい上位所定個数の前記遺伝子群Ｂを選択する遺伝子群選択工程を備えることを特徴としている。すなわち、本発明にかかる遺伝子相互作用推定方法においては、共通の制御を受ける遺伝子のフィルタリングにのみ相関係数を用いる。このように、本発明にかかる遺伝子相互作用推定方法においては、計算コストが高い制御領域の解析処理、すなわち、第２のフィルタリング処理に先だって、第１のフィルタリング処理として発現データの相関を利用した処理を行うことにより、計算コストが高い第２のフィルタリング処理の計算時間を大幅に削減することが可能となる。
【００２５】
また、本発明にかかる遺伝子相互作用推定方法は、第２の遺伝子群抽出手段により、前記遺伝子群抽出工程にて抽出された前記遺伝子群ＢＲのうち、特定の実験操作において大きく発現強度が変化しているもののみを選択する第２の遺伝子群抽出工程を備えることを特徴としている。具体的には、前記第２の遺伝子群抽出工程では、前記実験操作として遺伝子破壊を行うことによって得られているＤＮＡマイクロアレイ破壊株データにおいて前記遺伝子ａに強く影響を与えた因子群が求められ、求められた各因子Ｆについて、当該因子Ｆの破壊株における発現強度の変化比に基づいて、前記遺伝子群ＢＲのうち、前記因子Ｆの影響を強く受けた遺伝子群が遺伝子群ＢＲＦとして抽出される。このように、本発明にかかる遺伝子相互作用推定方法においては、第３のフィルタリング処理として、第２のフィルタリング処理によって抽出された遺伝子群ＢＲのうち、特定の実験操作において大きく発現強度が変化しているもののみを選択することにより、実験操作によって生成される因子又はそれに強く影響を受ける未知の因子を転写因子として推定することができる。
【００２６】
さらに、本発明にかかる遺伝子相互作用推定方法は、アラインメント手段により、前記遺伝子ａ及び前記遺伝子群ＢＲＦの各遺伝子について、マッチした位置からのウィンドウを取り出し、マルチプルアラインメントを行うアラインメント工程を備えることを特徴としている。これにより、本発明にかかる遺伝子相互作用推定方法においては、共通パターンの存在が確認できた場合には、推定が成功したものと判断することが可能となる。
【００２７】
また、上述した目的を達成する本発明にかかる遺伝子相互作用推定プログラムは、遺伝子ａが因子Ｆの制御を受けるとき、前記因子Ｆの制御を受ける他の遺伝子ｂに与える影響を推定するコンピュータ実行可能な遺伝子相互作用推定プログラムにおいて、前記コンピュータを、前記遺伝子ａと、当該遺伝子ａと相関係数の絶対値が大きい上位所定個数の遺伝子群Ｂの各遺伝子ｂとのペアについて、これら２つの遺伝子ａ，ｂの制御領域からそれぞれ所定長ｌの連続領域をウィンドウとして切り出し、これらのウィンドウに含まれる近接した開始位置をもつ所定長ｋの文字列同士の類似度の最大値をその位置の類似度と定義したときの各位置に対する類似度の合計であるウィンドウ類似度を算出するウィンドウ類似度算出手段として機能させることを特徴としている。
【００２８】
さらに、上述した目的を達成する本発明にかかる遺伝子相互作用推定装置は、遺伝子ａが因子Ｆの制御を受けるとき、前記因子Ｆの制御を受ける他の遺伝子ｂに与える影響を推定する遺伝子相互作用推定装置において、前記遺伝子ａと、当該遺伝子ａと相関係数の絶対値が大きい上位所定個数の遺伝子群Ｂの各遺伝子ｂとのペアについて、これら２つの遺伝子ａ，ｂの制御領域からそれぞれ所定長ｌの連続領域をウィンドウとして切り出し、これらのウィンドウに含まれる近接した開始位置をもつ所定長ｋの文字列同士の類似度の最大値をその位置の類似度と定義したときの各位置に対する類似度の合計であるウィンドウ類似度を算出するウィンドウ類似度算出手段を備えることを特徴としている。
【００２９】
このような本発明にかかる遺伝子相互作用推定プログラム及び遺伝子相互作用推定装置においては、ウィンドウ類似度という概念を取り入れることにより、遺伝子ａ，ｂの制御領域全体ではなく局所的な類似度を発見することが可能となる。
【発明の効果】
【００３０】
本発明によれば、個々の文字列が制御領域にどの程度共通に含まれるのかを調べるのではなく、ある程度の長さの部分領域（ウィンドウ）を文字列の集合体として捉え、ウィンドウが全体として類似しているか否かをも網羅的に計算することにより、文字列の組み合わせによる発現制御にも対応することが可能となる。そして、本発明によれば、ＤＮＡマイクロアレイデータの影響及び結合パターンに含まれるあいまいさの影響を低減することができ、ＤＮＡマイクロアレイを用いながらも、高精度に遺伝子相互作用の推定を行うことができる。
【発明を実施するための最良の形態】
【００３１】
以下、本発明を適用した具体的な実施の形態について図面を参照しながら詳細に説明する。
【００３２】
この実施の形態は、遺伝子ａが因子Ｆの制御（正の制御又は負の制御）を受けるとき、因子Ｆの制御を受ける他の遺伝子ｂにいかなる影響を与えるかを推定する遺伝子相互作用推定方法である。特に、この遺伝子相互作用推定方法は、入力されたデータに対して３つの独立したフィルタリング処理を施すことにより、ＤＮＡマイクロアレイを用いながらも、高精度に遺伝子相互作用の推定を行うことができるものである。
【００３３】
なお、以下では、被推定生物に実験操作を加えた遺伝子発現データとして、枯菌草バチルス・スブチリス（Bacillus subtillis）の遺伝子発現データを用いた実際のデータを用いた説明を適宜行うものとする。
【００３４】
この遺伝子相互作用推定方法を実行する遺伝子相互作用推定装置は、一般的には、プロセッサによって所定の遺伝子相互作用推定プログラムを実行するコンピュータを用いて実現される。図１に、本発明を適用した遺伝子相互作用推定方法における一連の処理を示す。
【００３５】
まず、この遺伝子相互作用推定方法においては、ＤＮＡマイクロアレイ破壊株データ、ＤＮＡ塩基配列データ、及び相互作用の推定対象となる遺伝子ａのデータが遺伝子相互作用推定装置に入力される。なお、ＤＮＡマイクロアレイデータとは、ＤＮＡマイクロアレイ技術によって得られた遺伝子ａの発現データであり、ＤＮＡマイクロアレイ破壊株データとは、ＤＮＡマイクロアレイデータの１種であり、破壊株と野生株との発現量の比を観測したデータである。ここで、破壊株とは、実験的手法によって標的となる遺伝子ａを破壊した株のことであり、野生株とは、遺伝子ａを破壊していない株のことである。また、株とは、遺伝的形質が同じ生物の別の個体である。さらに、ＤＮＡは、４つの塩基文字（ａ，ｔ，ｇ，ｃ）から構成される１つの塩基配列とその相補配列とが互いに結合し、二重らせん構造をとっている物質であるが、ＤＮＡ塩基配列データとは、かかる配列の片側の塩基配列を表したデータである。
【００３６】
この遺伝子相互作用推定方法においては、かかる入力に対して、ＤＮＡマイクロアレイによる発現データのバックグラウンドノイズの影響を低減するために、発現データの相関を利用した第１のフィルタリング処理を行う。すなわち、この遺伝子相互作用推定方法においては、ステップＳ１において、遺伝子相互作用推定装置における遺伝子群選択手段により、遺伝子ａと相関係数の絶対値が大きい上位所定個数の遺伝子群Ｂを選択する。具体的には、この遺伝子相互作用推定方法においては、遺伝子ａと相関係数の絶対値が大きい上位３００の遺伝子群Ｂを選択する。例えば、遺伝子ａが“ａｈｐＣ”である場合には、図２に示すように、相関係数の絶対値が大きい順に、“ａｈｐＣ”、“ａｈｐＦ”、“ｙｕｒＶ”、“ｙｕｒＷ”、“ｙｕｒＵ”、・・・といったように、上位３００の遺伝子が選択され、これらが遺伝子群Ｂとされる。
【００３７】
すなわち、この遺伝子相互作用推定方法においては、共通の制御を受ける遺伝子のフィルタリングにのみ相関係数を用いる。さらに換言すれば、この遺伝子相互作用推定方法においては、共通の制御を受ける遺伝子の候補を絞り込むために、相関係数の値そのものには拘泥せずに、ある値以上の相関係数を有するものは全て共通の制御を受ける遺伝子の候補として残すように処理を行う。
【００３８】
続いて、この遺伝子相互作用推定方法においては、第１のフィルタリング処理によって抽出された遺伝子群Ｂを対象として第２のフィルタリング処理を施し、制御領域内の類似部分領域の探索による結合サイト（binding site）の推定を行う。
【００３９】
まず、この遺伝子相互作用推定方法においては、第２のフィルタリング処理の一環として、図１中ステップＳ２において、遺伝子相互作用推定装置におけるウィンドウ類似度算出手段により、遺伝子ａと遺伝子群Ｂの各遺伝子ｂとのペアについて、制御領域のウィンドウ類似度を算出する。
【００４０】
ここで、ウィンドウ類似度とは、例えば図３に示すように、遺伝子の制御領域全体ではなく、類似部分の多少のずれを許容しつつ局所的な類似度を発見するために新たに考案した概念である。このとき、この遺伝子相互作用推定方法においては、期待値よりも出現頻度が特異的に高い文字列は結合サイトにはならず、かかる文字列が比較すべき２つの遺伝子の制御領域に共通に含まれていたとしても共通の因子が結合するわけではない、という本願出願人の発見に基づいて、結合サイトになりそうな文字列の類似性のみを考慮し、遺伝子の制御に関係していない文字列が類似していても、それらは無視する。
【００４１】
具体的には、ウィンドウ類似度は、例えば図４に示すように、比較する２つの遺伝子ａ，ｂの制御領域からそれぞれ長さｌ（例えば３０程度）の連続領域をウィンドウとして切り出し、これらのウィンドウに含まれる近接した開始位置をもつ長さｋ（例えば６程度）の文字列（オリゴヌクレオチド）同士の類似度の最大値をその位置の類似度と定義したときの各位置に対する類似度の合計値として算出される。なお、オリゴヌクレオチド同士の類似度は、一致する文字数に関する単調増加関数として定義する。このとき、前後に数個（例えば３程度）ずれた範囲での類似度を全て算出し、その最大値をその位置の類似度と定義する。これにより、この遺伝子相互作用推定方法は、モチーフアラインメントを近似的に行っていることになり、結合文字列に頻繁にみられる、間に可変長の任意文字列を含むパターンにも対応することが可能となる。また、統計的に特異的に出現頻度が高いオリゴヌクレオチドについては、結合サイトにはならないとの発見に基づいて、それに対する類似度を“０”とし、計算対象から除外する。実際には、例えば図４に示すように、遺伝子ａの制御領域のｉ文字目を開始位置とするウィンドウをＷａ［ｉ］とすると、このウィンドウＷａ［ｉ］に対する遺伝子ｂのウィンドウＷｂ［ｊ］の類似度ｗｓｉｍａ，ｂ（ｉ，ｊ）は、文字列の長さを“６”とすると、
ｗｓｉｍａ，ｂ（ｉ，ｊ）：＝Σｋ＝０，ｌ−６ｓｉｍａ，ｂ（ｉ＋ｋ，ｊ＋ｋ）
として算出される。
【００４２】
ここで、１つのウィンドウのペアに対する類似度の計算量は、ウィンドウの長さｌに比例すると考えられ、ウィンドウのペアの総数は、遺伝子ａ，ｂの制御領域全体の長さをｎとすると、高々ｎ２個である。したがって、全てのウィンドウのペアに対する計算量は、計算量がｌである１つのウィンドウのペアがｎ２個存在するため、Ｏ（ｎ２・ｌ）となるはずである。これに対して、この遺伝子相互作用推定方法においては、遺伝子ａ，ｂの制御領域全体の長さをｎとすると、全ての位置ｉ，ｊに対するｓｉｍａ，ｂ（ｉ，ｊ）は、予め計算量Ｏ（ｎ２）で算出することができ、さらに、全てのウィンドウのペアに対する類似度は、ｗｓｉｍａ，ｂ（ｉ＋１，ｊ＋１）＝ｗｓｉｍａ，ｂ（ｉ，ｊ）−ｓｉｍａ，ｂ（ｉ，ｊ）＋ｓｉｍａ，ｂ（ｉ＋ｌ−５，ｊ＋ｌ−５）の関係を用いると計算量Ｏ（ｎ２）で算出することができる。したがって、遺伝子相互作用推定方法においては、全てのウィンドウのペアに対する類似度は、ウィンドウの長さｌには依存せずに計算量Ｏ（ｎ２）で算出することができ、非常に少ない計算量で済ませることができる。
【００４３】
このように、この遺伝子相互作用推定方法においては、ウィンドウ類似度という概念を取り入れることにより、遺伝子ａ，ｂの制御領域全体ではなく局所的な類似度を発見する。また、この遺伝子相互作用推定方法においては、ウィンドウ類似度を算出する際に、遺伝子ａの制御領域における各位置ｉに対する遺伝子群Ｂの類似度の最大値ＭａｘａＢ［ｉ］と平均値ＡｖｇａＢ［ｉ］、及び遺伝子群Ｂの各遺伝子ｂの類似度の最大値Ｍａｘａｂ［ｉ］を算出する。なお、位置ｉとは、制御領域を文字列として取り出したとき、その先頭からｉ番目の文字の場所のことを意味する。また、位置ｉからの長さｌのウィンドウとは、ｉ番目からｉ＋ｌ−１番目までの部分文字列のことを意味する。
【００４４】
続いて、この遺伝子相互作用推定方法においては、第２のフィルタリング処理の一環として、図１中ステップＳ３において、遺伝子相互作用推定装置における領域探索手段により、遺伝子ａの制御領域のうち、最大値ＭａｘａＢ［ｉ］が位置ｉにおけるウィンドウ類似度の分布において統計的に特異的に大きく、且つ、ピーク位置を含むような領域Ｒを探索する。より具体的には、この遺伝子相互作用推定方法においては、遺伝子ａの制御領域のうち、最大値ＭａｘａＢ［ｉ］と平均値ＡｖｇａＢ［ｉ］との差分ＭａｘａＢ［ｉ］−ＡｖｇａＢ［ｉ］が所定値よりも大きく、且つ、ピーク位置を含むような領域Ｒを探索する。なお、ピーク位置とは、例えば図５中矢印で示すように、遺伝子ａの位置を前後にずらして当該遺伝子ａの制御領域の位置ｉの変化に対する遺伝子群Ｂの各遺伝子ｂのウィンドウ類似度の最大値Ｍａｘａｂ［ｉ］をみたときに極大値を与える位置のことである。このとき、この遺伝子相互作用推定方法においては、所定の閾値以上の類似度を有するピーク位置のみを対象とするために、最大値ＭａｘａＢ［ｉ］と平均値ＡｖｇａＢ［ｉ］との差分ＭａｘａＢ［ｉ］−ＡｖｇａＢ［ｉ］が所定値よりも大きい連続領域に含まれるピーク位置を対象とする。そして、この遺伝子相互作用推定方法においては、遺伝子ａ，ｂに共通の結合サイトが存在するとした場合には、このピーク位置からのウィンドウ内に存在する可能性が高いことから、かかるピーク位置を含むような領域Ｒを探索する。
【００４５】
具体的には、この遺伝子相互作用推定方法においては、例えば図６に示すように、遺伝子ａの制御領域の位置ｉの変化に対する類似度の最大値ＭａｘａＢ［ｉ］と平均値ＡｖｇａＢ［ｉ］との差分をプロットし、同図中塗りつぶしプロットで示されるピーク位置を求める。そして、この遺伝子相互作用推定方法においては、類似度が１８０００以上のピーク位置を対象とし、最大値ＭａｘａＢ［ｉ］と平均値ＡｖｇａＢ［ｉ］との差分ＭａｘａＢ［ｉ］−ＡｖｇａＢ［ｉ］の値が１５０００以上の値を連続的に有してピーク位置を含むような領域を探索するものとすると、例えば図７中斜線部に示すように、位置ｉ＝４１１８０５６〜４１１８０９６を領域Ｒとして求める。
【００４６】
なお、領域Ｒを探索するために、最大値ＭａｘａＢ［ｉ］ではなく差分ＭａｘａＢ［ｉ］−ＡｖｇａＢ［ｉ］を用いるのは、以下の理由による。上述したように、統計的に特異的に出現頻度が高いオリゴヌクレオチドについては、それに対する類似度を“０”とし、計算対象から除外しているため、各ウィンドウ毎に類似度の算出に用いられる文字列の総数は異なることになる。したがって、高い類似度の文字列がウィンドウ内に少数しか存在しない場合と、低い類似度の文字列がウィンドウ内に多数存在する場合とにおいて、同程度のウィンドウ類似度が算出される可能性がある。遺伝子相互作用推定方法においては、高い類似度を有するものを選ぶのが目的であるため、前者のみを選別する必要がある。ここで、前者と後者とでは、後者の方が平均値が高くなるため、最大値ＭａｘａＢ［ｉ］と平均値ＡｖｇａＢ［ｉ］との差分をみることにより、これら前者と後者とを簡便に区別することができる。領域Ｒを探索するために、最大値ＭａｘａＢ［ｉ］ではなく差分ＭａｘａＢ［ｉ］−ＡｖｇａＢ［ｉ］を用いるのは、このような理由によるものである。ただし、遺伝子相互作用推定方法においては、最大値ＭａｘａＢ［ｉ］が位置ｉにおけるウィンドウ類似度の分布において統計的に特異的に大きいと判断したものを選ぶことができるのであれば、差分ＭａｘａＢ［ｉ］−ＡｖｇａＢ［ｉ］を用いる方法以外であっても適用することができる。
【００４７】
遺伝子相互作用推定方法においては、このようにして領域Ｒを求めると、第２のフィルタリング処理の一環として、図１中ステップＳ４において、遺伝子相互作用推定装置における遺伝子群抽出手段により、求めた領域Ｒにピーク位置が存在し、且つ、類似度が所定値よりも高い遺伝子群ＢＲを求める。すなわち、この遺伝子相互作用推定方法においては、このステップＳ４を行うことにより、第１のフィルタリング処理によって抽出された遺伝子群Ｂのうち、制御領域内に所定値以上の類似度を有する遺伝子群ＢＲを抽出することができる。
【００４８】
このように、遺伝子相互作用推定方法においては、第２のフィルタリング処理として、特定のパターンを仮定せずに、指定した遺伝子ａと類似した制御領域を有する遺伝子を網羅的に探索する。このとき、この遺伝子相互作用推定方法においては、類似度を算出する際に、部分文字列の出現頻度に関する統計的特異性に基づいて、結合サイトに含まれる可能性が低い部分を除外する。このようにして求められた類似度が高い領域は、共通の因子が結合するサイトであるものと推定することができる。
【００４９】
また、遺伝子相互作用推定方法においては、図１中ステップＳ４において、遺伝子群ＢＲを求めると、被推定生物に実験操作を加えた遺伝子発現データを利用して転写因子を推定する処理を行う。すなわち、この遺伝子相互作用推定方法においては、実験操作を加えない個体と実験操作を加えた個体との遺伝子発現量の変化を比較することにより、実験操作の直接的な影響を知ることができる。具体的には、遺伝子相互作用推定方法においては、遺伝子相互作用推定装置における第２の遺伝子群抽出手段により、第３のフィルタリング処理として、第２のフィルタリング処理によって抽出された遺伝子群ＢＲのうち、特定の実験操作において大きく発現強度が変化しているもののみを選択する。ここでは、実験操作として遺伝子破壊を行うことによって得られているＤＮＡマイクロアレイ破壊株データにおいて遺伝子ａに強く影響を与えた因子群を求める。そして、この遺伝子相互作用推定方法においては、求めた各因子Ｆについて、遺伝子群ＢＲのうち、因子Ｆの影響を強く受けた遺伝子群を抽出し、それをＢＲＦとして求める。なお、遺伝子群ＢＲＦは、因子Ｆによる野生株の発現強度と比較した破壊株の発現強度の変化が大きく、且つ、その順位も上位であるものである。なお、発現強度の変化比は、ｌｏｇ（（破壊株の発現強度）／（野生株の発現強度））の絶対値で表される。
【００５０】
例えば、図８に示す遺伝子群ＢＲが抽出されているものとする。この場合、遺伝子相互作用推定方法においては、図９に示すように、遺伝子ａとしての“ａｈｐＣ”の発現強度が変化した破壊株データに対する破壊遺伝子が作り出す因子を、当該遺伝子ａに強く影響を与えた因子群として求める。そして、この遺伝子相互作用推定方法においては、例えば図１０に示すように、求めた因子Ｆとしての“ＰｅｒＲ”を作り出す遺伝子“ｐｅｒＲ”の破壊株における発現強度の変化比を求める。なお、同図中括弧内の数字は、遺伝子内での発現強度の変化の大きさの順位を示している。これにより、この遺伝子相互作用推定方法においては、遺伝子群ＢＲのうち、遺伝子ａとしての“ａｈｐＣ”に強く影響を与えた因子Ｆの１つである“ＰｅｒＲ”の影響を強く受けた遺伝子群ＢＲＦとして、図８及び図１０中太字で示すように、“ｋａｔＡ”、“ｙｆｍＪ”、“ｍｒｇＡ”、“ｈｅｍＡ”、“ｐｏｎＡ”、“ｙｋｖＷ”、及び“ｙｄｊＬ”を抽出することができる。
【００５１】
このように、遺伝子相互作用推定方法においては、第３のフィルタリング処理として、第２のフィルタリング処理によって抽出された遺伝子群ＢＲのうち、特定の実験操作において大きく発現強度が変化しているもののみを選択することにより、実験操作によって生成される因子又はそれに強く影響を受ける未知の因子を転写因子として推定することができる。
【００５２】
そして、この遺伝子相互作用推定方法においては、図１中ステップＳ５において、遺伝子相互作用推定装置におけるアラインメント手段により、遺伝子ａ及び遺伝子群ＢＲＦの各遺伝子について、マッチした位置からのウィンドウを取り出し、例えば図１１に示すように、マルチプルアラインメント（multiple alignment）を行う。そして、この遺伝子相互作用推定方法においては、共通パターンの存在が確認できた場合には、推定が成功したものと判断する。これにより、この遺伝子相互作用推定方法においては、例えば図１２に示すように、遺伝子ａを制御する因子Ｆ、制御方向（正又は負）、因子Ｆによって制御される遺伝子ａ以外の遺伝子群ＢＲＦ、及び遺伝子ａ及び遺伝子群ＢＲＦの各遺伝子の結合サイトを示す情報を出力し、一連の処理を終了する。
【００５３】
このように、遺伝子相互作用推定方法においては、３つのフィルタリング処理を行うことによって得られる複数の結果を複合的に利用し、遺伝子ａが因子Ｆの制御（正の制御又は負の制御）を受けるとき、因子Ｆの制御を受ける他の遺伝子ｂにいかなる影響を与えるかを推定する。このとき、この遺伝子相互作用推定方法においては、相関係数等の信頼性が低いデータについては大まかな選択を行うためにのみ用い、実験操作データ等の信頼性が高いデータについては値そのものを選択基準として用いることにより、結果の信頼性を高めることができる。
【００５４】
なお、上述した説明においては、図２及び図７乃至図１２等を用いて、遺伝子ａとして“ａｈｐＣ”が因子Ｆとしての“ＰｅｒＲ”の制御を受ける場合について具体的な結果を示しているが、以下、このような遺伝子相互作用推定方法を具体的に実行した他の実験結果についても参考までに示す。
【００５５】
まず、遺伝子ａとして、“ｐｕｒＡ”を用いた場合の具体的な実験結果について説明する。
【００５６】
まず、実験においては、図１中ステップＳ１における処理と同様に、遺伝子ａと相関係数の絶対値が大きい上位３００の遺伝子群Ｂを選択した。これにより、例えば図１３に示すように、相関係数の絶対値が大きい順に、遺伝子“ｐｕｒＡ”、“ｙｕｍＤ”、“ｇｌｙＡ”、“ｙｋｂＡ”、“ｐｕｒＱ”、・・・が遺伝子群Ｂとして選択された。
【００５７】
続いて、実験においては、図１中ステップＳ２における処理と同様に、遺伝子ａと遺伝子群Ｂの各遺伝子ｂとのペアについて、制御領域のウィンドウ類似度を算出し、ステップＳ３における処理と同様に、遺伝子ａの制御領域のうち、最大値ＭａｘａＢ［ｉ］と平均値ＡｖｇａＢ［ｉ］との差分ＭａｘａＢ［ｉ］−ＡｖｇａＢ［ｉ］が所定値よりも大きく、且つ、ピーク位置を含むような領域Ｒを探索した。ここでは、類似度が１８０００以上のピーク位置を対象とし、最大値ＭａｘａＢ［ｉ］と平均値ＡｖｇａＢ［ｉ］との差分ＭａｘａＢ［ｉ］−ＡｖｇａＢ［ｉ］の値が１５０００以上の値を連続的に有してピーク位置を含むような領域を探索するように条件を設定したため、図１４中斜線部に示すように、位置ｉ＝４１５６０１９〜４１５６０２９が領域Ｒとして求められた。
【００５８】
続いて、実験においては、図１中ステップＳ４における処理と同様に、求めた領域Ｒにピーク位置が存在し、且つ、類似度が所定値よりも高い遺伝子群ＢＲを抽出した。これにより、図１５に示す遺伝子群ＢＲが抽出された。そして、実験においては、ＤＮＡマイクロアレイ破壊株データにおいて遺伝子ａに強く影響を与えた因子群を求め、求めた各因子Ｆについて、遺伝子群ＢＲのうち、因子Ｆの影響を強く受けた遺伝子群を抽出し、それをＢＲＦとして求めた。これにより、図１６に示すように、遺伝子ａとしての“ｐｕｒＡ”の発現強度が変化した破壊株が、当該遺伝子“ｐｕｒＡ”に強く影響を与えた因子群として求められ、遺伝子群ＢＲのうち、当該遺伝子ｐｕｒＡ”に強く影響を与えた因子Ｆの１つである“ＰｕｒＲ”の影響を強く受けた遺伝子群ＢＲＦとして、図１５及び図１７中太字で示すように、“ｙｄｅＱ”及び“ｇｌｙＡ”が抽出された。
【００５９】
そして、実験においては、図１中ステップＳ５における処理と同様に、遺伝子ａ及び遺伝子群ＢＲＦの各遺伝子について、マッチした位置からのウィンドウを取り出し、マルチプルアラインメントを行い、図１８に示すような結果が得られた。
【００６０】
これにより、遺伝子ａを制御する因子Ｆ、制御方向（正又は負）、因子Ｆによって制御される遺伝子ａ以外の遺伝子群ＢＲＦ、及び遺伝子ａ及び遺伝子群ＢＲＦの各遺伝子の結合サイトを示す情報として、図１９に示すような結果が出力として求められた。
【００６１】
つぎに、遺伝子ａとして、“ｐｈｏＤ”を用いた場合の具体的な実験結果について説明する。
【００６２】
まず、実験においては、図１中ステップＳ１における処理と同様に、遺伝子ａと相関係数の絶対値が大きい上位３００の遺伝子群Ｂを選択した。これにより、例えば図２０に示すように、相関係数の絶対値が大きい順に、遺伝子“ｐｈｏＤ”、“ｙｑｇＧ”、“ｔｕａＤ”、“ｔｕａＢ”、“ｙｑｇＫ”、・・・が遺伝子群Ｂとして選択された。
【００６３】
続いて、実験においては、図１中ステップＳ２における処理と同様に、遺伝子ａと遺伝子群Ｂの各遺伝子ｂとのペアについて、制御領域のウィンドウ類似度を算出し、ステップＳ３における処理と同様に、遺伝子ａの制御領域のうち、最大値ＭａｘａＢ［ｉ］と平均値ＡｖｇａＢ［ｉ］との差分ＭａｘａＢ［ｉ］−ＡｖｇａＢ［ｉ］が所定値よりも大きく、且つ、ピーク位置を含むような領域Ｒを探索した。ここでは、類似度が１８０００以上のピーク位置を対象とし、最大値ＭａｘａＢ［ｉ］と平均値ＡｖｇａＢ［ｉ］との差分ＭａｘａＢ［ｉ］−ＡｖｇａＢ［ｉ］の値が１５０００以上の値を連続的に有してピーク位置を含むような領域を探索するように条件を設定したため、図２１中斜線部に示すように、位置ｉ＝２８３４９３〜２８３５０２が領域Ｒとして求められた。
【００６４】
続いて、実験においては、図１中ステップＳ４における処理と同様に、求めた領域Ｒにピーク位置が存在し、且つ、類似度が所定値よりも高い遺伝子群ＢＲを抽出した。これにより、図２２に示す遺伝子群ＢＲが抽出された。そして、実験においては、ＤＮＡマイクロアレイ破壊株データにおいて遺伝子ａに強く影響を与えた因子群を求め、求めた各因子Ｆについて、遺伝子群ＢＲのうち、因子Ｆの影響を強く受けた遺伝子群を抽出し、それをＢＲＦとして求めた。これにより、図２３に示すように、遺伝子ａとしての“ｐｈｏＤ”の発現強度が変化した破壊株が、当該遺伝子“ｐｈｏＤ”に強く影響を与えた因子群として求められ、遺伝子群ＢＲのうち、遺伝子“ｐｈｏＤ”に強く影響を与えた因子Ｆの１つである“ＰｈｏＰ”の影響を強く受けた遺伝子群ＢＲＦとして、図２２及び図２４中太字で示すように、“ｐｈｏＡ”及び“ｐｈｏＢ”が抽出された。
【００６５】
そして、実験においては、図１中ステップＳ５における処理と同様に、遺伝子ａ及び遺伝子群ＢＲＦの各遺伝子について、マッチした位置からのウィンドウを取り出し、マルチプルアラインメントを行い、図２５に示すような結果が得られた。
【００６６】
これにより、遺伝子ａを制御する因子Ｆ、制御方向（正又は負）、因子Ｆによって制御される遺伝子ａ以外の遺伝子群ＢＲＦ、及び遺伝子ａ及び遺伝子群ＢＲＦの各遺伝子の結合サイトを示す情報として、図２６に示すような結果が出力として求められた。
【００６７】
以上説明したように、本発明の実施の形態として示した遺伝子相互作用推定方法においては、発現データの相関を利用した第１のフィルタリング処理、結合サイトを推定する全く新規の手法である第２のフィルタリング処理、及び実験操作を加えた遺伝子発現データを利用した第３のフィルタリング処理による３つの独立した推定方法を融合させることにより、ＤＮＡマイクロアレイデータの影響及び結合パターンに含まれるあいまいさの影響を低減することができ、ＤＮＡマイクロアレイを用いながらも、高精度に遺伝子相互作用の推定を行うことができる。
【００６８】
特に、この遺伝子相互作用推定方法においては、計算コストが高い制御領域の解析処理、すなわち、第２のフィルタリング処理に先だって、第１のフィルタリング処理として発現データの相関を利用した処理を行うことにより、計算コストが高い第２のフィルタリング処理の計算時間を大幅に削減することが可能となる。なお、本願出願人は、この遺伝子相互作用推定方法を実現するプログラムを、インテル社製ＣＰＵ“Ｐｅｎｔｉｕｍ（登録商標）４、クロック３ＧＨｚ”の５倍程度の速度を有するＡＭＤ社製ＣＰＵ“Ｏｐｔｅｒｏｎ（登録商標）２．４ＧＨｚ”を４並列構成としたコンピュータを用いて実行した場合には、４０００種類の遺伝子全ての計算に約１０日程度要することを確認している。これに対して、本願出願人は、第１のフィルタリング処理を行わない場合には、この日数の１０倍以上の時間が必要となることも確認している。
【００６９】
また、この遺伝子相互作用推定方法においては、第２のフィルタリング処理として、プロモーター領域での発見で頻繁に行われるように、個々のオリゴヌクレオチドが制御領域にどの程度共通に含まれるのかを調べるのではなく、ある程度の長さの部分領域（ウィンドウ）をオリゴヌクレオチドの集合体として捉え、ウィンドウが全体として類似しているか否かをも網羅的に計算する。これにより、この遺伝子相互作用推定方法においては、オリゴヌクレオチドの組み合わせによる発現制御にも対応することが可能となる。さらに、本願出願人は、制御領域におけるオリゴヌクレオチドの出現頻度の解析を行い、統計的に出現頻度が特異的に高いオリゴヌクレオチドは結合文字列になる可能性が低いことを新たに発見したが、この遺伝子相互作用推定方法においては、この事実を領域の類似度算出に反映させることにより、結合サイト以外で類似領域が発見される可能性を低減することができる。さらにまた、この遺伝子相互作用推定方法においては、結合サイトはピーク位置に存在すると予想できることから、結合サイトの候補を絞り込むことも可能となる。
【００７０】
また、この遺伝子相互作用推定方法においては、第３のフィルタリング処理については、単一のＤＮＡマイクロアレイデータのみに依存し、ノイズの影響を受けやすいことから、最終段階に適用するが、かかる第３のフィルタリング処理を行うことにより、ＤＮＡマイクロアレイによる発現データのバックグラウンドノイズに対して十分に大きな発現強度の変化を与えることができ、また、各遺伝子に対する因子の直接的影響を導き出すことができる。さらに、この遺伝子相互作用推定方法においては、第３のフィルタリング処理を行うにあたって、値のみならず順位も考慮することにより、相対評価を行っていることと等価な状況を作り出すことができる。
【００７１】
なお、本発明は、上述した実施の形態に限定されるものではなく、その趣旨を逸脱しない範囲で適宜変更が可能であることはいうまでもない。
【図面の簡単な説明】
【００７２】
【図１】本発明の実施の形態として示す遺伝子相互作用推定方法における一連の処理をフローチャートである。
【図２】遺伝子ａと相関係数の絶対値が大きい上位所定個数の遺伝子群Ｂの具体例を説明する図である。
【図３】ウィンドウ類似度について説明するための図である。
【図４】ウィンドウ類似度の具体的な算出方法について説明するための図である。
【図５】ピーク位置について説明するための図である。
【図６】ピーク位置を具体的に求める様子について説明するための図である。
【図７】領域Ｒを具体的に求める様子について説明するための図である。
【図８】抽出された遺伝子群ＢＲの具体例を説明する図である。
【図９】遺伝子ａとしての“ａｈｐＣ”の発現強度が変化した破壊株の具体例を説明する図である。
【図１０】因子Ｆとしての遺伝子“ｐｅｒＲ”の破壊株における発現強度の変化比の具体例を説明する図である。
【図１１】マルチプルアラインメントの具体的な方法について説明するための図である。
【図１２】本発明の実施の形態として示す遺伝子相互作用推定方法における一連の処理を経て出力される情報の具体例を説明する図である。
【図１３】遺伝子ａとしての“ｐｕｒＡ”と相関係数の絶対値が大きい上位所定個数の遺伝子群Ｂの具体例を説明する図である。
【図１４】遺伝子ａとしての“ｐｕｒＡ”と図１３に示す遺伝子群Ｂとに基づいて求められた領域Ｒについて説明するための図である。
【図１５】図１４に示す領域Ｒから抽出された遺伝子群ＢＲの具体例を説明する図である。
【図１６】遺伝子ａとしての“ｐｕｒＡ”の発現強度が変化した破壊株の具体例を説明する図である。
【図１７】因子Ｆとしての遺伝子“ｐｕｒＲ”の破壊株における発現強度の変化比の具体例を説明する図である。
【図１８】遺伝子ａとしての“ｐｕｒＡ”及び遺伝子群ＢＲＦの各遺伝子について、マッチした位置からのウィンドウを取り出し、マルチプルアラインメントを行った結果について説明するための図である。
【図１９】図１３乃至図１８のデータに基づいて出力された実験結果としての情報の具体例を説明する図である。
【図２０】遺伝子ａとしての“ｐｈｏＤ”と相関係数の絶対値が大きい上位所定個数の遺伝子群Ｂの具体例を説明する図である。
【図２１】遺伝子ａとしての“ｐｈｏＤ”と図２０に示す遺伝子群Ｂとに基づいて求められた領域Ｒについて説明するための図である。
【図２２】図２１に示す領域Ｒから抽出された遺伝子群ＢＲの具体例を説明する図である。
【図２３】遺伝子ａとしての“ｐｈｏＤ”の発現強度が変化した破壊株の具体例を説明する図である。
【図２４】因子Ｆとしての遺伝子“ｐｈｏＰ”の破壊株における発現強度の変化比の具体例を説明する図である。
【図２５】遺伝子ａとしての“ｐｈｏＤ”及び遺伝子群ＢＲＦの各遺伝子について、マッチした位置からのウィンドウを取り出し、マルチプルアラインメントを行った結果について説明するための図である。
【図２６】図２０乃至図２５のデータに基づいて出力された実験結果としての情報の具体例を説明する図である。
【図２７】遺伝子の構造について説明する図である。
【図２８】ブーリアンネットワークの具体的なモデルについて説明する図である。
【図２９】ベイジアンネットワークの具体的なモデルについて説明する図である。
【符号の説明】
【００７３】
ａ，ｂ 遺伝子
ＡｖｇａＢ［ｉ］ 類似度の平均値
Ｂ，ＢＲ，ＢＲＦ 遺伝子群
Ｆ 因子
ＭａｘａＢ［ｉ］ 類似度の最大値
Ｒ 領域
Ｗａ［ｉ］，Ｗｂ［ｊ］ ウィンドウ
ｗｓｉｍａ，ｂ（ｉ，ｊ） ウィンドウ類似度
【技術分野】
【０００１】
本発明は、遺伝子が因子の制御を受けるとき、その因子の制御を受ける他の遺伝子に与える影響を推定する遺伝子相互作用推定方法、遺伝子相互作用推定プログラム、及び遺伝子相互作用推定装置に関する。
【背景技術】
【０００２】
近年のゲノム科学の進展にともない、ある環境における多数の遺伝子の発現強度を定量的に測定する技術であるＤＮＡマイクロアレイ技術により、多数の遺伝子発現データを同時に取得することが可能となっている。また、遺伝子の機能を解明する研究も、個々の遺伝子について解析する段階から、遺伝子の相互作用によって細胞内の活動がどのように行われるかを知る段階に移りつつある。
【０００３】
ここで、遺伝子は、例えば図２７に示すように、上流部分に制御領域を有し、そこに制御因子が結合することにより、遺伝子の発現が促進されたり抑制されたりする。また、遺伝子の発現によって因子（タンパク質）が生成されると、その因子が他の遺伝子の発現に影響を与えるが、各遺伝子が、自身を含めた他の遺伝子にいかなる影響を与えるかを推定する問題を遺伝子相互作用推定問題という。
【０００４】
かかる遺伝子相互作用推定の手法としては、発現データのみを利用したモデル推定によるものが提案されている。すなわち、ＤＮＡマイクロアレイによる発現データのみを利用し、発現データを最もよく説明可能なモデルを推定するものである。
【０００５】
まず、この種の遺伝子相互作用推定の手法としては、ブール代数を用いたブーリアンネットワークが提案されている（例えば、非特許文献１等参照。）。ブーリアンネットワークは、同期して動作する論理回路のようなモデルである。例えば図２８に示すモデルの場合には、１個のノード（頂点）が１個の遺伝子に対応し、各頂点は、“０（発現していない）”又は“１（発現している）”の状態をとる。また、各頂点の状態は、単位時刻毎に同期して変化し、遺伝子発現の制御規則は、ブール関数の形式で表現される。かかる手法においては、発現の時系列データを与え、データとモデルとの不一致を最小にするようなブーリアンネットワークを推定する問題を解くことになる。しかしながら、かかる手法においては、定量的な発現強度を扱うことができないという問題がある。
【０００６】
また、この種の遺伝子相互作用推定の手法としては、ベイジアンネットワークを用いる手法が提案されている（例えば、非特許文献２等参照。）。確率変数を頂点で表し、因果関係や相関関係といった依存する関係を有する変数間にリンクを張ったグラフ構造による確率モデルをグラフィカルモデルというが、このうち、例えば図２９に示すように、特にリンクの向きが因果関係の方向であり、このリンクを辿ったパスが循環しないものをベイジアンネットワークという。ベイジアンネットワークは、確率変数間の定性的な依存関係をグラフ構造によって表し、その確率変数間に定義される条件付き確率によって確率変数間の定量的な依存関係を表す。かかる手法においては、各遺伝子の発現強度を確率変数とし、その間の依存関係をベイジアンネットワークとして推定することになる。しかしながら、かかる手法においては、同じデータを説明可能な複数のモデルが存在してしまうという問題がある。
【０００７】
さらに、この種の遺伝子相互作用推定の手法としては、グラフィカルガウシアンモデリングを用いる手法が提案されている（例えば、非特許文献３等参照。）。相関係数によるクラスタリングの問題点は、直接の相互作用によって相関がみられる場合と、いくつかの他の遺伝子の発現を介した相互作用によって相関がみられる場合や、ある共通の遺伝子の発現に影響を受ける２つの遺伝子にみられる擬似相関の場合とを区別できないことにある。ここで、ｎ個の遺伝子がある環境での発現パターンは、ｎ次元の正規分布からのサンプルであると仮定すると、偏相関係数が“０”となるか否かによって擬似相関を排除することができる。これがグラフィカルガウシアンモデリングを用いる手法である。しかしながら、かかる手法においては、発現の挙動が非常に類似した遺伝子が存在した場合には、相関係数行列が正則でなくなり、偏相関係数を算出することができなくなる。そのため、かかる手法においては、前処理としてクラスタリングを行う必要があり、遺伝子グループ単位でしか依存関係を推定することができないという問題がある。
【０００８】
さらにまた、この種の遺伝子相互作用推定の手法としては、動的システムモデルを用いる手法が提案されている（例えば、非特許文献４等参照。）。この手法は、微分方程式等によるシステムモデルを仮定し、そのパラメータ推定を発現データとの一致度を目的関数とする最適化問題として定式化し、遺伝的アルゴリズムによって解くものである。しかしながら、この手法においては、計算コストが高く、数千もの遺伝子ネットワークを一度に解くことができず、分割や階層化等の処理を行う必要があるという問題がある。また、かかる手法においては、ノイズをも説明するモデルを出力してしまう可能性があり、いかにしてノイズを排除するかという問題もある。
【０００９】
このように、ＤＮＡマイクロアレイは、遺伝子相互作用の解析には有力な手段ではあるものの、上述した手法は、いずれも発現データのみを利用することから、ＤＮＡマイクロアレイによる発現データに含まれる多くのノイズの問題、実験コストの観点から数千もの遺伝子数に比べて数百程度の少数の実験データしか得られないという問題、転写因子と転写因子結合サイトとの関係を得ることができないという問題があり、推定精度が高くない。
【００１０】
この問題を解決するために、ＤＮＡマイクロアレイによる発現データと遺伝子の上流部分にある制御領域に含まれるパターンの解析とを組み合わせることにより、転写因子とこれによって制御される遺伝子群の関係を導き出す手法が提案されている（例えば、非特許文献５及び非特許文献６並びに特許文献１等参照。）。
【００１１】
具体的には、非特許文献５には、系列の発現データを用いて、転写因子の発現パターンと共通のモチーフを制御領域に有する遺伝子の発現パターンの合成との相関をみることにより、複数の転写因子の組み合わせによる影響を除外し、転写因子と結合サイトとの関係を推定する手法が開示されている。また、非特許文献６には、転写因子の発現パターンと各遺伝子の発現パターンとの相関によって遺伝子をクラスタリングし、その後モチーフをみつける手法が開示されている。さらに、特許文献１には、遺伝子発現量データと遺伝子配列データとに基づいてコンピュータによって遺伝子間の制御関係を推定する手法が開示されている。すなわち、この手法においては、遺伝子発現量の時系列データから、任意に早発現遺伝子と遅発現遺伝子とを選択し、その２遺伝子に関する相互相関関数を求め、その無相関検定を行って帰無仮説が棄却された場合には、その２遺伝子を記録する。そして、この手法は、早発現遺伝子と遅発現遺伝子とを選択し尽くすまで、この手順を繰り返し行い、その後、転写因子の上流配列に当該転写因子のコンセンサス制御配列の数が多い遺伝子と、先に記録された早発現遺伝子と遅発現遺伝子とをそれぞれ照合し、一致した場合に、その２遺伝子を記録して、これら２遺伝子には制御関係があると推定するものである。
【００１２】
【非特許文献１】阿久津達也著、「遺伝子ネットワークの推定アルゴリズム」、数理科学１９９９年６月号、株式会社サイエンス社、１９９９年、Ｎｏ．４３２、ｐ．４０−４６
【非特許文献２】N. Friedman, M. Linial, I. Nachman and D. Pe’er著, 「Using Bayesian networksto analyzeexpression data」, Journal of Computational Biology,vol. 7, no. 3/4, 2000年, p. 601-620
【非特許文献３】H. Toh, K. Horimoto著, 「Inference of a genetic network by a conbined approach of cluster analysis and graphicalGaussian modeling」, Bioinformatics, vol. 18, no. 2,2002年, p. 287-297
【非特許文献４】岡本正宏著、「Ｓ−ｓｙｓｔｅｍによる遺伝子の相互作用推定」、ゲノム情報生物学、株式会社中山書店、２０００年、ｐ．１６５−１８８
【非特許文献５】K. Birnbaum, P. N. Benfeyand D. E. Shasha著, 「cisElement/Transcription Factor Analysis (cis/TF): AMethod for Discovering Transcription Factor/cisElement Relationships」, Genome Res. 2001 September,vol. 11, no. 9, 2001年, p. 1567-1573
【非特許文献６】Z. Zhu, Y. Pilpel and G. M. Church著, 「Computational Identification ofTranscription Factor Binding Sites via a Transcription-factor-centricClustering (TFCC) Algorithm」, J. Mol. Biol. vol. 318,2002年, p. 7181
【特許文献１】特開２００３−１４１１２３号公報
【発明の開示】
【発明が解決しようとする課題】
【００１３】
ところで、各遺伝子の発現レベルは、様々な要因によって変化するが、転写因子の発現強度の変化は、ノイズ等のバックグラウンドの変化に対して極めて小さい。これに対して、上述した非特許文献５及び非特許文献６並びに特許文献１に記載された従来の手法は、いずれも発現データ間の相関を求めることが基本であることから、ノイズの影響を大きく受けるという問題があった。
【００１４】
また、多くの転写因子は、常に発現しており、その後の作用によって制御されていることが知られているが、従来の手法は、かかる転写因子に対しては適用することができないという問題もあった。
【００１５】
さらに、特開２０００−３４２２５７号公報には、エンハンサー又はプロモーターの内部構造を探究することにより、遺伝子の転写制御をより明確にすることができる手法が開示されている。すなわち、この手法においては、遺伝子の制御因子の結合部位構造を推定するにあたって、遺伝子をコードしているコード領域の上流域又は下流域にある転写因子となるタンパク質が結合するエンハンサー又はプロモーター領域内の制御因子の結合部位の制御構造を推定したい遺伝子を設定し、前記エンハンサー又はプロモーター領域において前記遺伝子に関与する制御因子、若しくは仮説的に導入する制御因子の結合部位について、その遺伝子座及びその他の遺伝子の発現の要因をパラメータとする計算モデルを構築し、この構築した計算モデルによって遺伝子の転写量を算出し、パラメータ探索アルゴリズムを使用して、実験で得られている前記設定された遺伝子の発現が得られるように計算モデルのパラメータを探索し、前記エンハンサー又はプロモーターのマイクロストラクチャーを推定している。しかしながら、この手法は、遺伝子の制御関係そのものを推定するものではない。
【００１６】
本発明は、このような実情に鑑みてなされたものであり、ＤＮＡマイクロアレイを用いながらも、高精度に遺伝子相互作用の推定を行うことができる遺伝子相互作用推定方法、遺伝子相互作用推定プログラム、及び遺伝子相互作用推定装置を提供することを目的とする。
【課題を解決するための手段】
【００１７】
本願出願人は、遺伝子が因子の制御を受けるとき、その因子の制御を受ける他の遺伝子に与える影響を推定する遺伝子相互作用推定問題について鋭意研究を重ねた結果、発現データの相関を利用した第１のフィルタリング処理、結合サイトを推定する全く新規の手法である第２のフィルタリング処理、及び実験操作を加えた遺伝子発現データを利用した第３のフィルタリング処理による３つの独立した推定方法を融合させる、という発想をなすに至り、さらに、第２のフィルタリング処理として行うべき具体的処理として極めて斬新な手法を考案することにより、本発明をなすに至った。
【００１８】
すなわち、上述した目的を達成する本発明にかかる遺伝子相互作用推定方法は、遺伝子ａが因子Ｆの制御を受けるとき、前記因子Ｆの制御を受ける他の遺伝子ｂに与える影響を推定する遺伝子相互作用推定方法において、ウィンドウ類似度算出手段により、前記遺伝子ａと、当該遺伝子ａと相関係数の絶対値が大きい上位所定個数の遺伝子群Ｂの各遺伝子ｂとのペアについて、これら２つの遺伝子ａ，ｂの制御領域からそれぞれ所定長ｌの連続領域をウィンドウとして切り出し、これらのウィンドウに含まれる近接した開始位置をもつ所定長ｋの文字列同士の類似度の最大値をその位置の類似度と定義したときの各位置に対する類似度の合計であるウィンドウ類似度を算出するウィンドウ類似度算出工程を備えることを特徴としている。
【００１９】
このような本発明にかかる遺伝子相互作用推定方法においては、ウィンドウ類似度という概念を取り入れることにより、遺伝子ａ，ｂの制御領域全体ではなく局所的な類似度を発見することが可能となる。
【００２０】
ここで、前記ウィンドウ類似度算出工程では、統計的に特異的に出現頻度が高い文字列については、それに対する類似度を０として前記ウィンドウ類似度が算出される。すなわち、本願出願人は、制御領域におけるオリゴヌクレオチドの出現頻度の解析を行い、統計的に出現頻度が特異的に高いオリゴヌクレオチドは結合文字列になる可能性が低いことを新たに発見したが、本発明にかかる遺伝子相互作用推定方法においては、この事実を領域の類似度算出に反映させることにより、結合サイト以外で類似領域が発見される可能性を低減することができる。
【００２１】
また、本発明にかかる遺伝子相互作用推定方法は、領域探索手段により、前記遺伝子ａの制御領域のうち、最大値ＭａｘａＢ［ｉ］が位置ｉにおける前記ウィンドウ類似度の分布において統計的に特異的に大きく、且つ、当該遺伝子ａの制御領域の位置ｉの変化に対する前記遺伝子群Ｂの各遺伝子ｂの前記ウィンドウ類似度の最大値Ｍａｘａｂ［ｉ］をみたときに極大値を与える位置であるピーク位置を含むような領域Ｒを探索する領域探索工程を備えることを特徴としている。すなわち、本発明にかかる遺伝子相互作用推定方法においては、遺伝子ａ，ｂに共通の結合サイトが存在するとした場合には、このピーク位置からのウィンドウ内に存在する可能性が高いことから、かかるピーク位置を含むような領域Ｒを探索することにより、結合サイトの候補を絞り込むことが可能となる。
【００２２】
さらに、本発明にかかる遺伝子相互作用推定方法は、遺伝子群抽出手段により、前記領域探索工程にて求められた前記領域Ｒにピーク位置が存在し、且つ、類似度が所定値よりも高い遺伝子群ＢＲを求めて抽出する遺伝子群抽出工程を備えることを特徴としている。これにより、本発明にかかる遺伝子相互作用推定方法においては、遺伝子群Ｂのうち、制御領域内に所定値以上の類似度を有する遺伝子群ＢＲを抽出することができる。
【００２３】
このように、本発明にかかる遺伝子相互作用推定方法においては、第２のフィルタリング処理としてこれらの工程を行い、特定のパターンを仮定せずに、指定した遺伝子ａと類似した制御領域を有する遺伝子を網羅的に探索する。このようにして求められた類似度が高い領域は、共通の因子が結合するサイトであるものと推定することができる。
【００２４】
さらにまた、本発明にかかる遺伝子相互作用推定方法は、前記ウィンドウ類似度算出工程に先だって、遺伝子群選択手段により、前記遺伝子ａと相関係数の絶対値が大きい上位所定個数の前記遺伝子群Ｂを選択する遺伝子群選択工程を備えることを特徴としている。すなわち、本発明にかかる遺伝子相互作用推定方法においては、共通の制御を受ける遺伝子のフィルタリングにのみ相関係数を用いる。このように、本発明にかかる遺伝子相互作用推定方法においては、計算コストが高い制御領域の解析処理、すなわち、第２のフィルタリング処理に先だって、第１のフィルタリング処理として発現データの相関を利用した処理を行うことにより、計算コストが高い第２のフィルタリング処理の計算時間を大幅に削減することが可能となる。
【００２５】
また、本発明にかかる遺伝子相互作用推定方法は、第２の遺伝子群抽出手段により、前記遺伝子群抽出工程にて抽出された前記遺伝子群ＢＲのうち、特定の実験操作において大きく発現強度が変化しているもののみを選択する第２の遺伝子群抽出工程を備えることを特徴としている。具体的には、前記第２の遺伝子群抽出工程では、前記実験操作として遺伝子破壊を行うことによって得られているＤＮＡマイクロアレイ破壊株データにおいて前記遺伝子ａに強く影響を与えた因子群が求められ、求められた各因子Ｆについて、当該因子Ｆの破壊株における発現強度の変化比に基づいて、前記遺伝子群ＢＲのうち、前記因子Ｆの影響を強く受けた遺伝子群が遺伝子群ＢＲＦとして抽出される。このように、本発明にかかる遺伝子相互作用推定方法においては、第３のフィルタリング処理として、第２のフィルタリング処理によって抽出された遺伝子群ＢＲのうち、特定の実験操作において大きく発現強度が変化しているもののみを選択することにより、実験操作によって生成される因子又はそれに強く影響を受ける未知の因子を転写因子として推定することができる。
【００２６】
さらに、本発明にかかる遺伝子相互作用推定方法は、アラインメント手段により、前記遺伝子ａ及び前記遺伝子群ＢＲＦの各遺伝子について、マッチした位置からのウィンドウを取り出し、マルチプルアラインメントを行うアラインメント工程を備えることを特徴としている。これにより、本発明にかかる遺伝子相互作用推定方法においては、共通パターンの存在が確認できた場合には、推定が成功したものと判断することが可能となる。
【００２７】
また、上述した目的を達成する本発明にかかる遺伝子相互作用推定プログラムは、遺伝子ａが因子Ｆの制御を受けるとき、前記因子Ｆの制御を受ける他の遺伝子ｂに与える影響を推定するコンピュータ実行可能な遺伝子相互作用推定プログラムにおいて、前記コンピュータを、前記遺伝子ａと、当該遺伝子ａと相関係数の絶対値が大きい上位所定個数の遺伝子群Ｂの各遺伝子ｂとのペアについて、これら２つの遺伝子ａ，ｂの制御領域からそれぞれ所定長ｌの連続領域をウィンドウとして切り出し、これらのウィンドウに含まれる近接した開始位置をもつ所定長ｋの文字列同士の類似度の最大値をその位置の類似度と定義したときの各位置に対する類似度の合計であるウィンドウ類似度を算出するウィンドウ類似度算出手段として機能させることを特徴としている。
【００２８】
さらに、上述した目的を達成する本発明にかかる遺伝子相互作用推定装置は、遺伝子ａが因子Ｆの制御を受けるとき、前記因子Ｆの制御を受ける他の遺伝子ｂに与える影響を推定する遺伝子相互作用推定装置において、前記遺伝子ａと、当該遺伝子ａと相関係数の絶対値が大きい上位所定個数の遺伝子群Ｂの各遺伝子ｂとのペアについて、これら２つの遺伝子ａ，ｂの制御領域からそれぞれ所定長ｌの連続領域をウィンドウとして切り出し、これらのウィンドウに含まれる近接した開始位置をもつ所定長ｋの文字列同士の類似度の最大値をその位置の類似度と定義したときの各位置に対する類似度の合計であるウィンドウ類似度を算出するウィンドウ類似度算出手段を備えることを特徴としている。
【００２９】
このような本発明にかかる遺伝子相互作用推定プログラム及び遺伝子相互作用推定装置においては、ウィンドウ類似度という概念を取り入れることにより、遺伝子ａ，ｂの制御領域全体ではなく局所的な類似度を発見することが可能となる。
【発明の効果】
【００３０】
本発明によれば、個々の文字列が制御領域にどの程度共通に含まれるのかを調べるのではなく、ある程度の長さの部分領域（ウィンドウ）を文字列の集合体として捉え、ウィンドウが全体として類似しているか否かをも網羅的に計算することにより、文字列の組み合わせによる発現制御にも対応することが可能となる。そして、本発明によれば、ＤＮＡマイクロアレイデータの影響及び結合パターンに含まれるあいまいさの影響を低減することができ、ＤＮＡマイクロアレイを用いながらも、高精度に遺伝子相互作用の推定を行うことができる。
【発明を実施するための最良の形態】
【００３１】
以下、本発明を適用した具体的な実施の形態について図面を参照しながら詳細に説明する。
【００３２】
この実施の形態は、遺伝子ａが因子Ｆの制御（正の制御又は負の制御）を受けるとき、因子Ｆの制御を受ける他の遺伝子ｂにいかなる影響を与えるかを推定する遺伝子相互作用推定方法である。特に、この遺伝子相互作用推定方法は、入力されたデータに対して３つの独立したフィルタリング処理を施すことにより、ＤＮＡマイクロアレイを用いながらも、高精度に遺伝子相互作用の推定を行うことができるものである。
【００３３】
なお、以下では、被推定生物に実験操作を加えた遺伝子発現データとして、枯菌草バチルス・スブチリス（Bacillus subtillis）の遺伝子発現データを用いた実際のデータを用いた説明を適宜行うものとする。
【００３４】
この遺伝子相互作用推定方法を実行する遺伝子相互作用推定装置は、一般的には、プロセッサによって所定の遺伝子相互作用推定プログラムを実行するコンピュータを用いて実現される。図１に、本発明を適用した遺伝子相互作用推定方法における一連の処理を示す。
【００３５】
まず、この遺伝子相互作用推定方法においては、ＤＮＡマイクロアレイ破壊株データ、ＤＮＡ塩基配列データ、及び相互作用の推定対象となる遺伝子ａのデータが遺伝子相互作用推定装置に入力される。なお、ＤＮＡマイクロアレイデータとは、ＤＮＡマイクロアレイ技術によって得られた遺伝子ａの発現データであり、ＤＮＡマイクロアレイ破壊株データとは、ＤＮＡマイクロアレイデータの１種であり、破壊株と野生株との発現量の比を観測したデータである。ここで、破壊株とは、実験的手法によって標的となる遺伝子ａを破壊した株のことであり、野生株とは、遺伝子ａを破壊していない株のことである。また、株とは、遺伝的形質が同じ生物の別の個体である。さらに、ＤＮＡは、４つの塩基文字（ａ，ｔ，ｇ，ｃ）から構成される１つの塩基配列とその相補配列とが互いに結合し、二重らせん構造をとっている物質であるが、ＤＮＡ塩基配列データとは、かかる配列の片側の塩基配列を表したデータである。
【００３６】
この遺伝子相互作用推定方法においては、かかる入力に対して、ＤＮＡマイクロアレイによる発現データのバックグラウンドノイズの影響を低減するために、発現データの相関を利用した第１のフィルタリング処理を行う。すなわち、この遺伝子相互作用推定方法においては、ステップＳ１において、遺伝子相互作用推定装置における遺伝子群選択手段により、遺伝子ａと相関係数の絶対値が大きい上位所定個数の遺伝子群Ｂを選択する。具体的には、この遺伝子相互作用推定方法においては、遺伝子ａと相関係数の絶対値が大きい上位３００の遺伝子群Ｂを選択する。例えば、遺伝子ａが“ａｈｐＣ”である場合には、図２に示すように、相関係数の絶対値が大きい順に、“ａｈｐＣ”、“ａｈｐＦ”、“ｙｕｒＶ”、“ｙｕｒＷ”、“ｙｕｒＵ”、・・・といったように、上位３００の遺伝子が選択され、これらが遺伝子群Ｂとされる。
【００３７】
すなわち、この遺伝子相互作用推定方法においては、共通の制御を受ける遺伝子のフィルタリングにのみ相関係数を用いる。さらに換言すれば、この遺伝子相互作用推定方法においては、共通の制御を受ける遺伝子の候補を絞り込むために、相関係数の値そのものには拘泥せずに、ある値以上の相関係数を有するものは全て共通の制御を受ける遺伝子の候補として残すように処理を行う。
【００３８】
続いて、この遺伝子相互作用推定方法においては、第１のフィルタリング処理によって抽出された遺伝子群Ｂを対象として第２のフィルタリング処理を施し、制御領域内の類似部分領域の探索による結合サイト（binding site）の推定を行う。
【００３９】
まず、この遺伝子相互作用推定方法においては、第２のフィルタリング処理の一環として、図１中ステップＳ２において、遺伝子相互作用推定装置におけるウィンドウ類似度算出手段により、遺伝子ａと遺伝子群Ｂの各遺伝子ｂとのペアについて、制御領域のウィンドウ類似度を算出する。
【００４０】
ここで、ウィンドウ類似度とは、例えば図３に示すように、遺伝子の制御領域全体ではなく、類似部分の多少のずれを許容しつつ局所的な類似度を発見するために新たに考案した概念である。このとき、この遺伝子相互作用推定方法においては、期待値よりも出現頻度が特異的に高い文字列は結合サイトにはならず、かかる文字列が比較すべき２つの遺伝子の制御領域に共通に含まれていたとしても共通の因子が結合するわけではない、という本願出願人の発見に基づいて、結合サイトになりそうな文字列の類似性のみを考慮し、遺伝子の制御に関係していない文字列が類似していても、それらは無視する。
【００４１】
具体的には、ウィンドウ類似度は、例えば図４に示すように、比較する２つの遺伝子ａ，ｂの制御領域からそれぞれ長さｌ（例えば３０程度）の連続領域をウィンドウとして切り出し、これらのウィンドウに含まれる近接した開始位置をもつ長さｋ（例えば６程度）の文字列（オリゴヌクレオチド）同士の類似度の最大値をその位置の類似度と定義したときの各位置に対する類似度の合計値として算出される。なお、オリゴヌクレオチド同士の類似度は、一致する文字数に関する単調増加関数として定義する。このとき、前後に数個（例えば３程度）ずれた範囲での類似度を全て算出し、その最大値をその位置の類似度と定義する。これにより、この遺伝子相互作用推定方法は、モチーフアラインメントを近似的に行っていることになり、結合文字列に頻繁にみられる、間に可変長の任意文字列を含むパターンにも対応することが可能となる。また、統計的に特異的に出現頻度が高いオリゴヌクレオチドについては、結合サイトにはならないとの発見に基づいて、それに対する類似度を“０”とし、計算対象から除外する。実際には、例えば図４に示すように、遺伝子ａの制御領域のｉ文字目を開始位置とするウィンドウをＷａ［ｉ］とすると、このウィンドウＷａ［ｉ］に対する遺伝子ｂのウィンドウＷｂ［ｊ］の類似度ｗｓｉｍａ，ｂ（ｉ，ｊ）は、文字列の長さを“６”とすると、
ｗｓｉｍａ，ｂ（ｉ，ｊ）：＝Σｋ＝０，ｌ−６ｓｉｍａ，ｂ（ｉ＋ｋ，ｊ＋ｋ）
として算出される。
【００４２】
ここで、１つのウィンドウのペアに対する類似度の計算量は、ウィンドウの長さｌに比例すると考えられ、ウィンドウのペアの総数は、遺伝子ａ，ｂの制御領域全体の長さをｎとすると、高々ｎ２個である。したがって、全てのウィンドウのペアに対する計算量は、計算量がｌである１つのウィンドウのペアがｎ２個存在するため、Ｏ（ｎ２・ｌ）となるはずである。これに対して、この遺伝子相互作用推定方法においては、遺伝子ａ，ｂの制御領域全体の長さをｎとすると、全ての位置ｉ，ｊに対するｓｉｍａ，ｂ（ｉ，ｊ）は、予め計算量Ｏ（ｎ２）で算出することができ、さらに、全てのウィンドウのペアに対する類似度は、ｗｓｉｍａ，ｂ（ｉ＋１，ｊ＋１）＝ｗｓｉｍａ，ｂ（ｉ，ｊ）−ｓｉｍａ，ｂ（ｉ，ｊ）＋ｓｉｍａ，ｂ（ｉ＋ｌ−５，ｊ＋ｌ−５）の関係を用いると計算量Ｏ（ｎ２）で算出することができる。したがって、遺伝子相互作用推定方法においては、全てのウィンドウのペアに対する類似度は、ウィンドウの長さｌには依存せずに計算量Ｏ（ｎ２）で算出することができ、非常に少ない計算量で済ませることができる。
【００４３】
このように、この遺伝子相互作用推定方法においては、ウィンドウ類似度という概念を取り入れることにより、遺伝子ａ，ｂの制御領域全体ではなく局所的な類似度を発見する。また、この遺伝子相互作用推定方法においては、ウィンドウ類似度を算出する際に、遺伝子ａの制御領域における各位置ｉに対する遺伝子群Ｂの類似度の最大値ＭａｘａＢ［ｉ］と平均値ＡｖｇａＢ［ｉ］、及び遺伝子群Ｂの各遺伝子ｂの類似度の最大値Ｍａｘａｂ［ｉ］を算出する。なお、位置ｉとは、制御領域を文字列として取り出したとき、その先頭からｉ番目の文字の場所のことを意味する。また、位置ｉからの長さｌのウィンドウとは、ｉ番目からｉ＋ｌ−１番目までの部分文字列のことを意味する。
【００４４】
続いて、この遺伝子相互作用推定方法においては、第２のフィルタリング処理の一環として、図１中ステップＳ３において、遺伝子相互作用推定装置における領域探索手段により、遺伝子ａの制御領域のうち、最大値ＭａｘａＢ［ｉ］が位置ｉにおけるウィンドウ類似度の分布において統計的に特異的に大きく、且つ、ピーク位置を含むような領域Ｒを探索する。より具体的には、この遺伝子相互作用推定方法においては、遺伝子ａの制御領域のうち、最大値ＭａｘａＢ［ｉ］と平均値ＡｖｇａＢ［ｉ］との差分ＭａｘａＢ［ｉ］−ＡｖｇａＢ［ｉ］が所定値よりも大きく、且つ、ピーク位置を含むような領域Ｒを探索する。なお、ピーク位置とは、例えば図５中矢印で示すように、遺伝子ａの位置を前後にずらして当該遺伝子ａの制御領域の位置ｉの変化に対する遺伝子群Ｂの各遺伝子ｂのウィンドウ類似度の最大値Ｍａｘａｂ［ｉ］をみたときに極大値を与える位置のことである。このとき、この遺伝子相互作用推定方法においては、所定の閾値以上の類似度を有するピーク位置のみを対象とするために、最大値ＭａｘａＢ［ｉ］と平均値ＡｖｇａＢ［ｉ］との差分ＭａｘａＢ［ｉ］−ＡｖｇａＢ［ｉ］が所定値よりも大きい連続領域に含まれるピーク位置を対象とする。そして、この遺伝子相互作用推定方法においては、遺伝子ａ，ｂに共通の結合サイトが存在するとした場合には、このピーク位置からのウィンドウ内に存在する可能性が高いことから、かかるピーク位置を含むような領域Ｒを探索する。
【００４５】
具体的には、この遺伝子相互作用推定方法においては、例えば図６に示すように、遺伝子ａの制御領域の位置ｉの変化に対する類似度の最大値ＭａｘａＢ［ｉ］と平均値ＡｖｇａＢ［ｉ］との差分をプロットし、同図中塗りつぶしプロットで示されるピーク位置を求める。そして、この遺伝子相互作用推定方法においては、類似度が１８０００以上のピーク位置を対象とし、最大値ＭａｘａＢ［ｉ］と平均値ＡｖｇａＢ［ｉ］との差分ＭａｘａＢ［ｉ］−ＡｖｇａＢ［ｉ］の値が１５０００以上の値を連続的に有してピーク位置を含むような領域を探索するものとすると、例えば図７中斜線部に示すように、位置ｉ＝４１１８０５６〜４１１８０９６を領域Ｒとして求める。
【００４６】
なお、領域Ｒを探索するために、最大値ＭａｘａＢ［ｉ］ではなく差分ＭａｘａＢ［ｉ］−ＡｖｇａＢ［ｉ］を用いるのは、以下の理由による。上述したように、統計的に特異的に出現頻度が高いオリゴヌクレオチドについては、それに対する類似度を“０”とし、計算対象から除外しているため、各ウィンドウ毎に類似度の算出に用いられる文字列の総数は異なることになる。したがって、高い類似度の文字列がウィンドウ内に少数しか存在しない場合と、低い類似度の文字列がウィンドウ内に多数存在する場合とにおいて、同程度のウィンドウ類似度が算出される可能性がある。遺伝子相互作用推定方法においては、高い類似度を有するものを選ぶのが目的であるため、前者のみを選別する必要がある。ここで、前者と後者とでは、後者の方が平均値が高くなるため、最大値ＭａｘａＢ［ｉ］と平均値ＡｖｇａＢ［ｉ］との差分をみることにより、これら前者と後者とを簡便に区別することができる。領域Ｒを探索するために、最大値ＭａｘａＢ［ｉ］ではなく差分ＭａｘａＢ［ｉ］−ＡｖｇａＢ［ｉ］を用いるのは、このような理由によるものである。ただし、遺伝子相互作用推定方法においては、最大値ＭａｘａＢ［ｉ］が位置ｉにおけるウィンドウ類似度の分布において統計的に特異的に大きいと判断したものを選ぶことができるのであれば、差分ＭａｘａＢ［ｉ］−ＡｖｇａＢ［ｉ］を用いる方法以外であっても適用することができる。
【００４７】
遺伝子相互作用推定方法においては、このようにして領域Ｒを求めると、第２のフィルタリング処理の一環として、図１中ステップＳ４において、遺伝子相互作用推定装置における遺伝子群抽出手段により、求めた領域Ｒにピーク位置が存在し、且つ、類似度が所定値よりも高い遺伝子群ＢＲを求める。すなわち、この遺伝子相互作用推定方法においては、このステップＳ４を行うことにより、第１のフィルタリング処理によって抽出された遺伝子群Ｂのうち、制御領域内に所定値以上の類似度を有する遺伝子群ＢＲを抽出することができる。
【００４８】
このように、遺伝子相互作用推定方法においては、第２のフィルタリング処理として、特定のパターンを仮定せずに、指定した遺伝子ａと類似した制御領域を有する遺伝子を網羅的に探索する。このとき、この遺伝子相互作用推定方法においては、類似度を算出する際に、部分文字列の出現頻度に関する統計的特異性に基づいて、結合サイトに含まれる可能性が低い部分を除外する。このようにして求められた類似度が高い領域は、共通の因子が結合するサイトであるものと推定することができる。
【００４９】
また、遺伝子相互作用推定方法においては、図１中ステップＳ４において、遺伝子群ＢＲを求めると、被推定生物に実験操作を加えた遺伝子発現データを利用して転写因子を推定する処理を行う。すなわち、この遺伝子相互作用推定方法においては、実験操作を加えない個体と実験操作を加えた個体との遺伝子発現量の変化を比較することにより、実験操作の直接的な影響を知ることができる。具体的には、遺伝子相互作用推定方法においては、遺伝子相互作用推定装置における第２の遺伝子群抽出手段により、第３のフィルタリング処理として、第２のフィルタリング処理によって抽出された遺伝子群ＢＲのうち、特定の実験操作において大きく発現強度が変化しているもののみを選択する。ここでは、実験操作として遺伝子破壊を行うことによって得られているＤＮＡマイクロアレイ破壊株データにおいて遺伝子ａに強く影響を与えた因子群を求める。そして、この遺伝子相互作用推定方法においては、求めた各因子Ｆについて、遺伝子群ＢＲのうち、因子Ｆの影響を強く受けた遺伝子群を抽出し、それをＢＲＦとして求める。なお、遺伝子群ＢＲＦは、因子Ｆによる野生株の発現強度と比較した破壊株の発現強度の変化が大きく、且つ、その順位も上位であるものである。なお、発現強度の変化比は、ｌｏｇ（（破壊株の発現強度）／（野生株の発現強度））の絶対値で表される。
【００５０】
例えば、図８に示す遺伝子群ＢＲが抽出されているものとする。この場合、遺伝子相互作用推定方法においては、図９に示すように、遺伝子ａとしての“ａｈｐＣ”の発現強度が変化した破壊株データに対する破壊遺伝子が作り出す因子を、当該遺伝子ａに強く影響を与えた因子群として求める。そして、この遺伝子相互作用推定方法においては、例えば図１０に示すように、求めた因子Ｆとしての“ＰｅｒＲ”を作り出す遺伝子“ｐｅｒＲ”の破壊株における発現強度の変化比を求める。なお、同図中括弧内の数字は、遺伝子内での発現強度の変化の大きさの順位を示している。これにより、この遺伝子相互作用推定方法においては、遺伝子群ＢＲのうち、遺伝子ａとしての“ａｈｐＣ”に強く影響を与えた因子Ｆの１つである“ＰｅｒＲ”の影響を強く受けた遺伝子群ＢＲＦとして、図８及び図１０中太字で示すように、“ｋａｔＡ”、“ｙｆｍＪ”、“ｍｒｇＡ”、“ｈｅｍＡ”、“ｐｏｎＡ”、“ｙｋｖＷ”、及び“ｙｄｊＬ”を抽出することができる。
【００５１】
このように、遺伝子相互作用推定方法においては、第３のフィルタリング処理として、第２のフィルタリング処理によって抽出された遺伝子群ＢＲのうち、特定の実験操作において大きく発現強度が変化しているもののみを選択することにより、実験操作によって生成される因子又はそれに強く影響を受ける未知の因子を転写因子として推定することができる。
【００５２】
そして、この遺伝子相互作用推定方法においては、図１中ステップＳ５において、遺伝子相互作用推定装置におけるアラインメント手段により、遺伝子ａ及び遺伝子群ＢＲＦの各遺伝子について、マッチした位置からのウィンドウを取り出し、例えば図１１に示すように、マルチプルアラインメント（multiple alignment）を行う。そして、この遺伝子相互作用推定方法においては、共通パターンの存在が確認できた場合には、推定が成功したものと判断する。これにより、この遺伝子相互作用推定方法においては、例えば図１２に示すように、遺伝子ａを制御する因子Ｆ、制御方向（正又は負）、因子Ｆによって制御される遺伝子ａ以外の遺伝子群ＢＲＦ、及び遺伝子ａ及び遺伝子群ＢＲＦの各遺伝子の結合サイトを示す情報を出力し、一連の処理を終了する。
【００５３】
このように、遺伝子相互作用推定方法においては、３つのフィルタリング処理を行うことによって得られる複数の結果を複合的に利用し、遺伝子ａが因子Ｆの制御（正の制御又は負の制御）を受けるとき、因子Ｆの制御を受ける他の遺伝子ｂにいかなる影響を与えるかを推定する。このとき、この遺伝子相互作用推定方法においては、相関係数等の信頼性が低いデータについては大まかな選択を行うためにのみ用い、実験操作データ等の信頼性が高いデータについては値そのものを選択基準として用いることにより、結果の信頼性を高めることができる。
【００５４】
なお、上述した説明においては、図２及び図７乃至図１２等を用いて、遺伝子ａとして“ａｈｐＣ”が因子Ｆとしての“ＰｅｒＲ”の制御を受ける場合について具体的な結果を示しているが、以下、このような遺伝子相互作用推定方法を具体的に実行した他の実験結果についても参考までに示す。
【００５５】
まず、遺伝子ａとして、“ｐｕｒＡ”を用いた場合の具体的な実験結果について説明する。
【００５６】
まず、実験においては、図１中ステップＳ１における処理と同様に、遺伝子ａと相関係数の絶対値が大きい上位３００の遺伝子群Ｂを選択した。これにより、例えば図１３に示すように、相関係数の絶対値が大きい順に、遺伝子“ｐｕｒＡ”、“ｙｕｍＤ”、“ｇｌｙＡ”、“ｙｋｂＡ”、“ｐｕｒＱ”、・・・が遺伝子群Ｂとして選択された。
【００５７】
続いて、実験においては、図１中ステップＳ２における処理と同様に、遺伝子ａと遺伝子群Ｂの各遺伝子ｂとのペアについて、制御領域のウィンドウ類似度を算出し、ステップＳ３における処理と同様に、遺伝子ａの制御領域のうち、最大値ＭａｘａＢ［ｉ］と平均値ＡｖｇａＢ［ｉ］との差分ＭａｘａＢ［ｉ］−ＡｖｇａＢ［ｉ］が所定値よりも大きく、且つ、ピーク位置を含むような領域Ｒを探索した。ここでは、類似度が１８０００以上のピーク位置を対象とし、最大値ＭａｘａＢ［ｉ］と平均値ＡｖｇａＢ［ｉ］との差分ＭａｘａＢ［ｉ］−ＡｖｇａＢ［ｉ］の値が１５０００以上の値を連続的に有してピーク位置を含むような領域を探索するように条件を設定したため、図１４中斜線部に示すように、位置ｉ＝４１５６０１９〜４１５６０２９が領域Ｒとして求められた。
【００５８】
続いて、実験においては、図１中ステップＳ４における処理と同様に、求めた領域Ｒにピーク位置が存在し、且つ、類似度が所定値よりも高い遺伝子群ＢＲを抽出した。これにより、図１５に示す遺伝子群ＢＲが抽出された。そして、実験においては、ＤＮＡマイクロアレイ破壊株データにおいて遺伝子ａに強く影響を与えた因子群を求め、求めた各因子Ｆについて、遺伝子群ＢＲのうち、因子Ｆの影響を強く受けた遺伝子群を抽出し、それをＢＲＦとして求めた。これにより、図１６に示すように、遺伝子ａとしての“ｐｕｒＡ”の発現強度が変化した破壊株が、当該遺伝子“ｐｕｒＡ”に強く影響を与えた因子群として求められ、遺伝子群ＢＲのうち、当該遺伝子ｐｕｒＡ”に強く影響を与えた因子Ｆの１つである“ＰｕｒＲ”の影響を強く受けた遺伝子群ＢＲＦとして、図１５及び図１７中太字で示すように、“ｙｄｅＱ”及び“ｇｌｙＡ”が抽出された。
【００５９】
そして、実験においては、図１中ステップＳ５における処理と同様に、遺伝子ａ及び遺伝子群ＢＲＦの各遺伝子について、マッチした位置からのウィンドウを取り出し、マルチプルアラインメントを行い、図１８に示すような結果が得られた。
【００６０】
これにより、遺伝子ａを制御する因子Ｆ、制御方向（正又は負）、因子Ｆによって制御される遺伝子ａ以外の遺伝子群ＢＲＦ、及び遺伝子ａ及び遺伝子群ＢＲＦの各遺伝子の結合サイトを示す情報として、図１９に示すような結果が出力として求められた。
【００６１】
つぎに、遺伝子ａとして、“ｐｈｏＤ”を用いた場合の具体的な実験結果について説明する。
【００６２】
まず、実験においては、図１中ステップＳ１における処理と同様に、遺伝子ａと相関係数の絶対値が大きい上位３００の遺伝子群Ｂを選択した。これにより、例えば図２０に示すように、相関係数の絶対値が大きい順に、遺伝子“ｐｈｏＤ”、“ｙｑｇＧ”、“ｔｕａＤ”、“ｔｕａＢ”、“ｙｑｇＫ”、・・・が遺伝子群Ｂとして選択された。
【００６３】
続いて、実験においては、図１中ステップＳ２における処理と同様に、遺伝子ａと遺伝子群Ｂの各遺伝子ｂとのペアについて、制御領域のウィンドウ類似度を算出し、ステップＳ３における処理と同様に、遺伝子ａの制御領域のうち、最大値ＭａｘａＢ［ｉ］と平均値ＡｖｇａＢ［ｉ］との差分ＭａｘａＢ［ｉ］−ＡｖｇａＢ［ｉ］が所定値よりも大きく、且つ、ピーク位置を含むような領域Ｒを探索した。ここでは、類似度が１８０００以上のピーク位置を対象とし、最大値ＭａｘａＢ［ｉ］と平均値ＡｖｇａＢ［ｉ］との差分ＭａｘａＢ［ｉ］−ＡｖｇａＢ［ｉ］の値が１５０００以上の値を連続的に有してピーク位置を含むような領域を探索するように条件を設定したため、図２１中斜線部に示すように、位置ｉ＝２８３４９３〜２８３５０２が領域Ｒとして求められた。
【００６４】
続いて、実験においては、図１中ステップＳ４における処理と同様に、求めた領域Ｒにピーク位置が存在し、且つ、類似度が所定値よりも高い遺伝子群ＢＲを抽出した。これにより、図２２に示す遺伝子群ＢＲが抽出された。そして、実験においては、ＤＮＡマイクロアレイ破壊株データにおいて遺伝子ａに強く影響を与えた因子群を求め、求めた各因子Ｆについて、遺伝子群ＢＲのうち、因子Ｆの影響を強く受けた遺伝子群を抽出し、それをＢＲＦとして求めた。これにより、図２３に示すように、遺伝子ａとしての“ｐｈｏＤ”の発現強度が変化した破壊株が、当該遺伝子“ｐｈｏＤ”に強く影響を与えた因子群として求められ、遺伝子群ＢＲのうち、遺伝子“ｐｈｏＤ”に強く影響を与えた因子Ｆの１つである“ＰｈｏＰ”の影響を強く受けた遺伝子群ＢＲＦとして、図２２及び図２４中太字で示すように、“ｐｈｏＡ”及び“ｐｈｏＢ”が抽出された。
【００６５】
そして、実験においては、図１中ステップＳ５における処理と同様に、遺伝子ａ及び遺伝子群ＢＲＦの各遺伝子について、マッチした位置からのウィンドウを取り出し、マルチプルアラインメントを行い、図２５に示すような結果が得られた。
【００６６】
これにより、遺伝子ａを制御する因子Ｆ、制御方向（正又は負）、因子Ｆによって制御される遺伝子ａ以外の遺伝子群ＢＲＦ、及び遺伝子ａ及び遺伝子群ＢＲＦの各遺伝子の結合サイトを示す情報として、図２６に示すような結果が出力として求められた。
【００６７】
以上説明したように、本発明の実施の形態として示した遺伝子相互作用推定方法においては、発現データの相関を利用した第１のフィルタリング処理、結合サイトを推定する全く新規の手法である第２のフィルタリング処理、及び実験操作を加えた遺伝子発現データを利用した第３のフィルタリング処理による３つの独立した推定方法を融合させることにより、ＤＮＡマイクロアレイデータの影響及び結合パターンに含まれるあいまいさの影響を低減することができ、ＤＮＡマイクロアレイを用いながらも、高精度に遺伝子相互作用の推定を行うことができる。
【００６８】
特に、この遺伝子相互作用推定方法においては、計算コストが高い制御領域の解析処理、すなわち、第２のフィルタリング処理に先だって、第１のフィルタリング処理として発現データの相関を利用した処理を行うことにより、計算コストが高い第２のフィルタリング処理の計算時間を大幅に削減することが可能となる。なお、本願出願人は、この遺伝子相互作用推定方法を実現するプログラムを、インテル社製ＣＰＵ“Ｐｅｎｔｉｕｍ（登録商標）４、クロック３ＧＨｚ”の５倍程度の速度を有するＡＭＤ社製ＣＰＵ“Ｏｐｔｅｒｏｎ（登録商標）２．４ＧＨｚ”を４並列構成としたコンピュータを用いて実行した場合には、４０００種類の遺伝子全ての計算に約１０日程度要することを確認している。これに対して、本願出願人は、第１のフィルタリング処理を行わない場合には、この日数の１０倍以上の時間が必要となることも確認している。
【００６９】
また、この遺伝子相互作用推定方法においては、第２のフィルタリング処理として、プロモーター領域での発見で頻繁に行われるように、個々のオリゴヌクレオチドが制御領域にどの程度共通に含まれるのかを調べるのではなく、ある程度の長さの部分領域（ウィンドウ）をオリゴヌクレオチドの集合体として捉え、ウィンドウが全体として類似しているか否かをも網羅的に計算する。これにより、この遺伝子相互作用推定方法においては、オリゴヌクレオチドの組み合わせによる発現制御にも対応することが可能となる。さらに、本願出願人は、制御領域におけるオリゴヌクレオチドの出現頻度の解析を行い、統計的に出現頻度が特異的に高いオリゴヌクレオチドは結合文字列になる可能性が低いことを新たに発見したが、この遺伝子相互作用推定方法においては、この事実を領域の類似度算出に反映させることにより、結合サイト以外で類似領域が発見される可能性を低減することができる。さらにまた、この遺伝子相互作用推定方法においては、結合サイトはピーク位置に存在すると予想できることから、結合サイトの候補を絞り込むことも可能となる。
【００７０】
また、この遺伝子相互作用推定方法においては、第３のフィルタリング処理については、単一のＤＮＡマイクロアレイデータのみに依存し、ノイズの影響を受けやすいことから、最終段階に適用するが、かかる第３のフィルタリング処理を行うことにより、ＤＮＡマイクロアレイによる発現データのバックグラウンドノイズに対して十分に大きな発現強度の変化を与えることができ、また、各遺伝子に対する因子の直接的影響を導き出すことができる。さらに、この遺伝子相互作用推定方法においては、第３のフィルタリング処理を行うにあたって、値のみならず順位も考慮することにより、相対評価を行っていることと等価な状況を作り出すことができる。
【００７１】
なお、本発明は、上述した実施の形態に限定されるものではなく、その趣旨を逸脱しない範囲で適宜変更が可能であることはいうまでもない。
【図面の簡単な説明】
【００７２】
【図１】本発明の実施の形態として示す遺伝子相互作用推定方法における一連の処理をフローチャートである。
【図２】遺伝子ａと相関係数の絶対値が大きい上位所定個数の遺伝子群Ｂの具体例を説明する図である。
【図３】ウィンドウ類似度について説明するための図である。
【図４】ウィンドウ類似度の具体的な算出方法について説明するための図である。
【図５】ピーク位置について説明するための図である。
【図６】ピーク位置を具体的に求める様子について説明するための図である。
【図７】領域Ｒを具体的に求める様子について説明するための図である。
【図８】抽出された遺伝子群ＢＲの具体例を説明する図である。
【図９】遺伝子ａとしての“ａｈｐＣ”の発現強度が変化した破壊株の具体例を説明する図である。
【図１０】因子Ｆとしての遺伝子“ｐｅｒＲ”の破壊株における発現強度の変化比の具体例を説明する図である。
【図１１】マルチプルアラインメントの具体的な方法について説明するための図である。
【図１２】本発明の実施の形態として示す遺伝子相互作用推定方法における一連の処理を経て出力される情報の具体例を説明する図である。
【図１３】遺伝子ａとしての“ｐｕｒＡ”と相関係数の絶対値が大きい上位所定個数の遺伝子群Ｂの具体例を説明する図である。
【図１４】遺伝子ａとしての“ｐｕｒＡ”と図１３に示す遺伝子群Ｂとに基づいて求められた領域Ｒについて説明するための図である。
【図１５】図１４に示す領域Ｒから抽出された遺伝子群ＢＲの具体例を説明する図である。
【図１６】遺伝子ａとしての“ｐｕｒＡ”の発現強度が変化した破壊株の具体例を説明する図である。
【図１７】因子Ｆとしての遺伝子“ｐｕｒＲ”の破壊株における発現強度の変化比の具体例を説明する図である。
【図１８】遺伝子ａとしての“ｐｕｒＡ”及び遺伝子群ＢＲＦの各遺伝子について、マッチした位置からのウィンドウを取り出し、マルチプルアラインメントを行った結果について説明するための図である。
【図１９】図１３乃至図１８のデータに基づいて出力された実験結果としての情報の具体例を説明する図である。
【図２０】遺伝子ａとしての“ｐｈｏＤ”と相関係数の絶対値が大きい上位所定個数の遺伝子群Ｂの具体例を説明する図である。
【図２１】遺伝子ａとしての“ｐｈｏＤ”と図２０に示す遺伝子群Ｂとに基づいて求められた領域Ｒについて説明するための図である。
【図２２】図２１に示す領域Ｒから抽出された遺伝子群ＢＲの具体例を説明する図である。
【図２３】遺伝子ａとしての“ｐｈｏＤ”の発現強度が変化した破壊株の具体例を説明する図である。
【図２４】因子Ｆとしての遺伝子“ｐｈｏＰ”の破壊株における発現強度の変化比の具体例を説明する図である。
【図２５】遺伝子ａとしての“ｐｈｏＤ”及び遺伝子群ＢＲＦの各遺伝子について、マッチした位置からのウィンドウを取り出し、マルチプルアラインメントを行った結果について説明するための図である。
【図２６】図２０乃至図２５のデータに基づいて出力された実験結果としての情報の具体例を説明する図である。
【図２７】遺伝子の構造について説明する図である。
【図２８】ブーリアンネットワークの具体的なモデルについて説明する図である。
【図２９】ベイジアンネットワークの具体的なモデルについて説明する図である。
【符号の説明】
【００７３】
ａ，ｂ 遺伝子
ＡｖｇａＢ［ｉ］ 類似度の平均値
Ｂ，ＢＲ，ＢＲＦ 遺伝子群
Ｆ 因子
ＭａｘａＢ［ｉ］ 類似度の最大値
Ｒ 領域
Ｗａ［ｉ］，Ｗｂ［ｊ］ ウィンドウ
ｗｓｉｍａ，ｂ（ｉ，ｊ） ウィンドウ類似度
【特許請求の範囲】
【請求項１】
遺伝子ａが因子Ｆの制御を受けるとき、前記因子Ｆの制御を受ける他の遺伝子ｂに与える影響を推定する遺伝子相互作用推定方法において、
ウィンドウ類似度算出手段により、前記遺伝子ａと、当該遺伝子ａと相関係数の絶対値が大きい上位所定個数の遺伝子群Ｂの各遺伝子ｂとのペアについて、これら２つの遺伝子ａ，ｂの制御領域からそれぞれ所定長ｌの連続領域をウィンドウとして切り出し、これらのウィンドウに含まれる近接した開始位置をもつ所定長ｋの文字列同士の類似度の最大値をその位置の類似度と定義したときの各位置に対する類似度の合計であるウィンドウ類似度を算出するウィンドウ類似度算出工程を備えること
を特徴とする遺伝子相互作用推定方法。
【請求項２】
前記ウィンドウ類似度算出工程では、統計的に特異的に出現頻度が高い文字列については、それに対する類似度を０として前記ウィンドウ類似度が算出されること
を特徴とする請求項１記載の遺伝子相互作用推定方法。
【請求項３】
前記ウィンドウ類似度算出工程では、前記遺伝子ａ，ｂの制御領域全体の長さｎに対して、全てのウィンドウのペアに対する類似度が前記ウィンドウの長さｌに依存しない計算量Ｏ（ｎ２）で算出されること
を特徴とする請求項１記載の遺伝子相互作用推定方法。
【請求項４】
前記ウィンドウ類似度算出工程では、前記遺伝子ａの制御領域を文字列として取り出したときにその先頭からｉ番目の文字の場所である当該制御領域における各位置ｉに対する前記遺伝子群Ｂの類似度の最大値ＭａｘａＢ［ｉ］と平均値ＡｖｇａＢ［ｉ］とが算出されること
を特徴とする請求項１記載の遺伝子相互作用推定方法。
【請求項５】
領域探索手段により、前記遺伝子ａの制御領域のうち、最大値ＭａｘａＢ［ｉ］が位置ｉにおける前記ウィンドウ類似度の分布において統計的に特異的に大きく、且つ、当該遺伝子ａの制御領域の位置ｉの変化に対する前記遺伝子群Ｂの各遺伝子ｂの前記ウィンドウ類似度の最大値Ｍａｘａｂ［ｉ］をみたときに極大値を与える位置であるピーク位置を含むような領域Ｒを探索する領域探索工程を備えること
を特徴とする請求項４記載の遺伝子相互作用推定方法。
【請求項６】
前記領域探索工程では、最大値ＭａｘａＢ［ｉ］が位置ｉにおける前記ウィンドウ類似度の分布において統計的に特異的に大きい連続領域に含まれるピーク位置を対象として、前記領域Ｒが探索されること
を特徴とする請求項５記載の遺伝子相互作用推定方法。
【請求項７】
前記領域探索工程では、最大値ＭａｘａＢ［ｉ］と平均値ＡｖｇａＢ［ｉ］との差分ＭａｘａＢ［ｉ］−ＡｖｇａＢ［ｉ］が所定値よりも大きい場合に、当該最大値ＭａｘａＢ［ｉ］が位置ｉにおける前記ウィンドウ類似度の分布において統計的に特異的に大きいと判断すること
を特徴とする請求項５又は請求項６記載の遺伝子相互作用推定方法。
【請求項８】
遺伝子群抽出手段により、前記領域探索工程にて求められた前記領域Ｒにピーク位置が存在し、且つ、類似度が所定値よりも高い遺伝子群ＢＲを求めて抽出する遺伝子群抽出工程を備えること
を特徴とする請求項５記載の遺伝子相互作用推定方法。
【請求項９】
前記ウィンドウ類似度算出工程に先だって、遺伝子群選択手段により、前記遺伝子ａと相関係数の絶対値が大きい上位所定個数の前記遺伝子群Ｂを選択する遺伝子群選択工程を備えること
を特徴とする請求項８記載の遺伝子相互作用推定方法。
【請求項１０】
第２の遺伝子群抽出手段により、前記遺伝子群抽出工程にて抽出された前記遺伝子群ＢＲのうち、特定の実験操作において大きく発現強度が変化しているもののみを選択する第２の遺伝子群抽出工程を備えること
を特徴とする請求項９記載の遺伝子相互作用推定方法。
【請求項１１】
前記第２の遺伝子群抽出工程では、前記実験操作として遺伝子破壊を行うことによって得られているＤＮＡマイクロアレイ破壊株データにおいて前記遺伝子ａに強く影響を与えた因子群が求められ、求められた各因子Ｆについて、当該因子Ｆの破壊株における発現強度の変化比に基づいて、前記遺伝子群ＢＲのうち、前記因子Ｆの影響を強く受けた遺伝子群が遺伝子群ＢＲＦとして抽出されること
を特徴とする請求項１０記載の遺伝子相互作用推定方法。
【請求項１２】
アラインメント手段により、前記遺伝子ａ及び前記遺伝子群ＢＲＦの各遺伝子について、マッチした位置からのウィンドウを取り出し、マルチプルアラインメントを行うアラインメント工程を備えること
を特徴とする請求項１１記載の遺伝子相互作用推定方法。
【請求項１３】
遺伝子ａが因子Ｆの制御を受けるとき、前記因子Ｆの制御を受ける他の遺伝子ｂに与える影響を推定するコンピュータ実行可能な遺伝子相互作用推定プログラムにおいて、
前記コンピュータを、
前記遺伝子ａと、当該遺伝子ａと相関係数の絶対値が大きい上位所定個数の遺伝子群Ｂの各遺伝子ｂとのペアについて、これら２つの遺伝子ａ，ｂの制御領域からそれぞれ所定長ｌの連続領域をウィンドウとして切り出し、これらのウィンドウに含まれる近接した開始位置をもつ所定長ｋの文字列同士の類似度の最大値をその位置の類似度と定義したときの各位置に対する類似度の合計であるウィンドウ類似度を算出するウィンドウ類似度算出手段として機能させること
を特徴とする遺伝子相互作用推定プログラム。
【請求項１４】
遺伝子ａが因子Ｆの制御を受けるとき、前記因子Ｆの制御を受ける他の遺伝子ｂに与える影響を推定する遺伝子相互作用推定装置において、
前記遺伝子ａと、当該遺伝子ａと相関係数の絶対値が大きい上位所定個数の遺伝子群Ｂの各遺伝子ｂとのペアについて、これら２つの遺伝子ａ，ｂの制御領域からそれぞれ所定長ｌの連続領域をウィンドウとして切り出し、これらのウィンドウに含まれる近接した開始位置をもつ所定長ｋの文字列同士の類似度の最大値をその位置の類似度と定義したときの各位置に対する類似度の合計であるウィンドウ類似度を算出するウィンドウ類似度算出手段を備えること
を特徴とする遺伝子相互作用推定装置。
【請求項１】
遺伝子ａが因子Ｆの制御を受けるとき、前記因子Ｆの制御を受ける他の遺伝子ｂに与える影響を推定する遺伝子相互作用推定方法において、
ウィンドウ類似度算出手段により、前記遺伝子ａと、当該遺伝子ａと相関係数の絶対値が大きい上位所定個数の遺伝子群Ｂの各遺伝子ｂとのペアについて、これら２つの遺伝子ａ，ｂの制御領域からそれぞれ所定長ｌの連続領域をウィンドウとして切り出し、これらのウィンドウに含まれる近接した開始位置をもつ所定長ｋの文字列同士の類似度の最大値をその位置の類似度と定義したときの各位置に対する類似度の合計であるウィンドウ類似度を算出するウィンドウ類似度算出工程を備えること
を特徴とする遺伝子相互作用推定方法。
【請求項２】
前記ウィンドウ類似度算出工程では、統計的に特異的に出現頻度が高い文字列については、それに対する類似度を０として前記ウィンドウ類似度が算出されること
を特徴とする請求項１記載の遺伝子相互作用推定方法。
【請求項３】
前記ウィンドウ類似度算出工程では、前記遺伝子ａ，ｂの制御領域全体の長さｎに対して、全てのウィンドウのペアに対する類似度が前記ウィンドウの長さｌに依存しない計算量Ｏ（ｎ２）で算出されること
を特徴とする請求項１記載の遺伝子相互作用推定方法。
【請求項４】
前記ウィンドウ類似度算出工程では、前記遺伝子ａの制御領域を文字列として取り出したときにその先頭からｉ番目の文字の場所である当該制御領域における各位置ｉに対する前記遺伝子群Ｂの類似度の最大値ＭａｘａＢ［ｉ］と平均値ＡｖｇａＢ［ｉ］とが算出されること
を特徴とする請求項１記載の遺伝子相互作用推定方法。
【請求項５】
領域探索手段により、前記遺伝子ａの制御領域のうち、最大値ＭａｘａＢ［ｉ］が位置ｉにおける前記ウィンドウ類似度の分布において統計的に特異的に大きく、且つ、当該遺伝子ａの制御領域の位置ｉの変化に対する前記遺伝子群Ｂの各遺伝子ｂの前記ウィンドウ類似度の最大値Ｍａｘａｂ［ｉ］をみたときに極大値を与える位置であるピーク位置を含むような領域Ｒを探索する領域探索工程を備えること
を特徴とする請求項４記載の遺伝子相互作用推定方法。
【請求項６】
前記領域探索工程では、最大値ＭａｘａＢ［ｉ］が位置ｉにおける前記ウィンドウ類似度の分布において統計的に特異的に大きい連続領域に含まれるピーク位置を対象として、前記領域Ｒが探索されること
を特徴とする請求項５記載の遺伝子相互作用推定方法。
【請求項７】
前記領域探索工程では、最大値ＭａｘａＢ［ｉ］と平均値ＡｖｇａＢ［ｉ］との差分ＭａｘａＢ［ｉ］−ＡｖｇａＢ［ｉ］が所定値よりも大きい場合に、当該最大値ＭａｘａＢ［ｉ］が位置ｉにおける前記ウィンドウ類似度の分布において統計的に特異的に大きいと判断すること
を特徴とする請求項５又は請求項６記載の遺伝子相互作用推定方法。
【請求項８】
遺伝子群抽出手段により、前記領域探索工程にて求められた前記領域Ｒにピーク位置が存在し、且つ、類似度が所定値よりも高い遺伝子群ＢＲを求めて抽出する遺伝子群抽出工程を備えること
を特徴とする請求項５記載の遺伝子相互作用推定方法。
【請求項９】
前記ウィンドウ類似度算出工程に先だって、遺伝子群選択手段により、前記遺伝子ａと相関係数の絶対値が大きい上位所定個数の前記遺伝子群Ｂを選択する遺伝子群選択工程を備えること
を特徴とする請求項８記載の遺伝子相互作用推定方法。
【請求項１０】
第２の遺伝子群抽出手段により、前記遺伝子群抽出工程にて抽出された前記遺伝子群ＢＲのうち、特定の実験操作において大きく発現強度が変化しているもののみを選択する第２の遺伝子群抽出工程を備えること
を特徴とする請求項９記載の遺伝子相互作用推定方法。
【請求項１１】
前記第２の遺伝子群抽出工程では、前記実験操作として遺伝子破壊を行うことによって得られているＤＮＡマイクロアレイ破壊株データにおいて前記遺伝子ａに強く影響を与えた因子群が求められ、求められた各因子Ｆについて、当該因子Ｆの破壊株における発現強度の変化比に基づいて、前記遺伝子群ＢＲのうち、前記因子Ｆの影響を強く受けた遺伝子群が遺伝子群ＢＲＦとして抽出されること
を特徴とする請求項１０記載の遺伝子相互作用推定方法。
【請求項１２】
アラインメント手段により、前記遺伝子ａ及び前記遺伝子群ＢＲＦの各遺伝子について、マッチした位置からのウィンドウを取り出し、マルチプルアラインメントを行うアラインメント工程を備えること
を特徴とする請求項１１記載の遺伝子相互作用推定方法。
【請求項１３】
遺伝子ａが因子Ｆの制御を受けるとき、前記因子Ｆの制御を受ける他の遺伝子ｂに与える影響を推定するコンピュータ実行可能な遺伝子相互作用推定プログラムにおいて、
前記コンピュータを、
前記遺伝子ａと、当該遺伝子ａと相関係数の絶対値が大きい上位所定個数の遺伝子群Ｂの各遺伝子ｂとのペアについて、これら２つの遺伝子ａ，ｂの制御領域からそれぞれ所定長ｌの連続領域をウィンドウとして切り出し、これらのウィンドウに含まれる近接した開始位置をもつ所定長ｋの文字列同士の類似度の最大値をその位置の類似度と定義したときの各位置に対する類似度の合計であるウィンドウ類似度を算出するウィンドウ類似度算出手段として機能させること
を特徴とする遺伝子相互作用推定プログラム。
【請求項１４】
遺伝子ａが因子Ｆの制御を受けるとき、前記因子Ｆの制御を受ける他の遺伝子ｂに与える影響を推定する遺伝子相互作用推定装置において、
前記遺伝子ａと、当該遺伝子ａと相関係数の絶対値が大きい上位所定個数の遺伝子群Ｂの各遺伝子ｂとのペアについて、これら２つの遺伝子ａ，ｂの制御領域からそれぞれ所定長ｌの連続領域をウィンドウとして切り出し、これらのウィンドウに含まれる近接した開始位置をもつ所定長ｋの文字列同士の類似度の最大値をその位置の類似度と定義したときの各位置に対する類似度の合計であるウィンドウ類似度を算出するウィンドウ類似度算出手段を備えること
を特徴とする遺伝子相互作用推定装置。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図１６】
【図１７】
【図１８】
【図１９】
【図２０】
【図２１】
【図２２】
【図２３】
【図２４】
【図２５】
【図２６】
【図２７】
【図２８】
【図２９】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図１６】
【図１７】
【図１８】
【図１９】
【図２０】
【図２１】
【図２２】
【図２３】
【図２４】
【図２５】
【図２６】
【図２７】
【図２８】
【図２９】
【公開番号】特開２００８−２７１５１（Ｐ２００８−２７１５１Ａ）
【公開日】平成２０年２月７日（２００８．２．７）
【国際特許分類】
【出願番号】特願２００６−１９８４８５（Ｐ２００６−１９８４８５）
【出願日】平成１８年７月２０日（２００６．７．２０）
【出願人】（３０４０２４４３０）国立大学法人北陸先端科学技術大学院大学 (169)
【公開日】平成２０年２月７日（２００８．２．７）
【国際特許分類】
【出願日】平成１８年７月２０日（２００６．７．２０）
【出願人】（３０４０２４４３０）国立大学法人北陸先端科学技術大学院大学 (169)
[ Back to top ]