(57)【要約】【課題】本発明の実施形態は、遺伝子ネットワークを含む、複雑な生物学的情報の解析に対する新規の推論法の適用を含む。新規方法は、ベイズ推論法の修正、および発現遺伝子間の因果関係を決定するための、およびある実施形態においては、調節遺伝子の上流エフェクターを決定するための当該方法の適用を含む。【解決手段】ベイズ法の追加的な修正は、発現遺伝子間の因果関係を推論するための経時的データの使用および遺伝子破壊データの使用を含む。他の実施形態は、より正確に発現遺伝子間のネットワーク情報を提供するためのブートストラップ法の使用およびエッジ効果の決定を含む。遺伝子ネットワークについての情報は記憶装置に格納され得るとともに出力装置に転送され得るか、または遠隔地に転送され得る。【選択図】なし【特許請求の範囲】
【請求項１】
遺伝子ネットワークを構築する方法であって、
１つまたは複数の種類の生物学的データを値の表現に変換するステップと、
前記１つまたは複数の種類の生物学的データからの前記値の表現の少なくとも１つをベイズ計算モデルにおけるベイズ事前確率として使用して前記遺伝子ネットワークを構築するステップと、
を含んでなる方法。
【請求項２】
前記１つまたは複数の種類の生物学的データが遺伝子発現データを含んでなる、請求項１に記載の方法。
【請求項３】
前記遺伝子発現データが転写産物の検出により得られる、請求項２に記載の方法。
【請求項４】
少なくとも２種類の遺伝子発現データを含んでなる、請求項１に記載の方法。
【請求項５】
前記少なくとも２種類の遺伝子発現データの少なくとも１種類が経時的遺伝子発現データである、請求項４に記載の方法。
【請求項６】
前記少なくとも２種類の遺伝子発現データの少なくとも１種類が遺伝子ノックダウン発現データである、請求項４に記載の方法。
【請求項７】
前記２種類の遺伝子発現データが経時的遺伝子発現データおよび遺伝子ノックダウン発現データである、請求項４に記載の方法。
【請求項８】
前記値の表現が行列表現である、請求項１に記載の方法。
【請求項９】
前記値が離散値または連続値である、請求項１に記載の方法。
【請求項１０】
前記変換するステップが、
遺伝子の第１の集合について前記経時的遺伝子発現データから遺伝子発現行列を作成するステップであって、前記データが前記第１の集合内の各遺伝子の発現の時間経過に基づく発現結果を含み、前記第１の集合内の前記遺伝子の互いに対する各時点の変動の平均効果および測度を定量化するステップ、
前記第１の集合内の前記遺伝子間のネットワーク関係を生成するステップ、
第１のベイズ事前確率として前記経時的発現データに基づくベイズ計算モデルを提供するステップであって、前記ベイズモデルがＢＮＲＣdynamic基準を最小化するステップを含んでなるステップ、
遺伝子の第２の集合の前記遺伝子ノックダウン発現データから遺伝子発現行列を作成するステップであって、前記データが遺伝子の第３の集合内の各遺伝子の破壊に基づく発現結果を含み、前記第２の集合内の前記遺伝子の各々の発現に対する前記第３の集合内の前記遺伝子の各々の破壊の変動の平均効果および測度を定量化するステップ、
前記第２の集合内の遺伝子と前記第３の集合内の遺伝子との間のネットワーク関係を生成するステップ、および
前記第２の集合内の遺伝子と前記第３の集合内の遺伝子との間の前記ネットワーク関係の行列表現を第２のベイズ事前確率として提供するステップ、
を含んでなる、請求項７に記載の方法。
【請求項１１】
前記変動の測度が分散である、請求項１０に記載の方法。
【請求項１２】
ＢＮＲＣdynamic基準を最小化する前記ステップが、多項式基底、フーリエ級数、ウェーブレット基底、回帰スプライン基底およびＢスプラインからなる群より選択される非線形曲線フィッティング法を使用するステップを含んでなる、請求項１０に記載の方法。
【請求項１３】
前記非線形曲線フィッティング法がノンパラメトリック法である、請求項１２に記載の方法。
【請求項１４】
ＢＮＲＣdynamic基準を最小化するための前記ノンパラメトリック法が不均一誤差分散を使用するステップを含む、請求項１３に記載の方法。
【請求項１５】
ベイズ計算モデルにおけるエッジの信頼性がブートストラップ法を使用して決定される、請求項１０に記載の方法。
【請求項１６】
前記ブートストラップ法が、
ａ）第１の集合についての経時的遺伝子発現データから何回も復元無作為抽出することによりブートストラップ遺伝子発現行列を提供するステップ、
ｂ）前記ブートストラップ遺伝子発現行列に基づき前記遺伝子ネットワークを推定するステップ、
ｃ）ステップａ）およびｂ）をＢ回反復し、それによりＢ個の遺伝子ネットワークを生成するステップ、および
ｄ）エッジの信頼性を前記Ｂ個の遺伝子ネットワークから計算するステップ、
を含んでなる、請求項１５に記載の方法。
【請求項１７】
遺伝子ノックダウン発現データ行列を前記第１および第２のベイズ事前確率と組み合わせて前記遺伝子ネットワークを構築するステップをさらに含んでなる、請求項１０に記載の方法。
【請求項１８】
コンピュータシステムと連動した使用向けのコンピュータプログラム製品であって、前記コンピュータプログラム製品がコンピュータにより読取可能な記憶メディアおよびそのなかに埋め込まれるコンピュータプログラム機構を含んでなり、前記コンピュータプログラム機構が遺伝子ネットワークを構築するための構築モジュールを含んでなり、
（ａ）１つまたは複数の種類の生物学的データをそれぞれ値の表現に変換するための命令、
（ｂ）ベイズ計算モデルにおけるベイズ事前確率として各値の表現を使用して前記遺伝子ネットワークを構築するための命令
を含んでなる、コンピュータプログラム製品。
【請求項１９】
コンピュータ読取可能メモリであって、記憶メディアであるとともに、
プロセッサにより実行可能なそのなかに埋め込まれる命令を含むコンピュータプログラムであって、前記プロセッサが前記命令を実行するとき、
１つまたは複数の種類の生物学的データをそれぞれ値の表現に変換するステップ、
各値の表現をベイズ計算モデルにおけるベイズ事前確率として使用して前記遺伝子ネットワークを構築するステップ、
を含む複数のステップを実施するコンピュータプログラム、
を含んでなる、コンピュータ読取可能メモリ。
【請求項２０】
請求項１に記載の方法により構築される遺伝子ネットワークモデルを含んでなるデータベース。
【請求項２１】
薬剤により影響を受ける遺伝子ネットワークを同定する方法であって、
薬剤への曝露により生成される遺伝子の第１の集合についての経時的遺伝子発現データを提供するステップ、
遺伝子の第２の集合についての遺伝子ノックダウン発現データを提供するステップ、
前記薬剤により影響を受ける遺伝子ネットワークを経時的遺伝子発現データおよび遺伝子ノックダウン発現データの組み合わせに基づき同定するステップであって、前記経時的発現データおよび遺伝子ノックダウン発現データの少なくとも一方がベイズ計算モデルにおけるベイズプライアとして使用されるステップ、
を含んでなる、方法。
【請求項２２】
遺伝子ネットワークを含む系内の標的遺伝子を同定する方法であって、請求項１に記載のステップを含んでなり、請求項１に記載の方法により構築される前記遺伝子ネットワークの親遺伝子が前記標的遺伝子であるよう選択される、方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、遺伝子ネットワークモデルを構築するとともに遺伝子間の関係を決定するためのシステムおよび方法に関する。
【０００２】
本出願は、参照により全体として本明細書に援用される、２００５年１０月１２日付けで出願された米国仮特許出願第６０／７２５，７０１号明細書の利益を主張する。
【背景技術】
【０００３】
関連技術の説明
生命科学、医薬、創薬ならびに新薬開発および製薬業界における最近の研究および展開の最も重要な側面の１つは、大量の生データを解釈するとともにかかるデータに基づく結論を引き出すための方法および装置を開発する必要性である。バイオインフォマティクスはシステム生物学の理解に大きく貢献してきたとともに生物系の構成要素間の複雑な関係のさらなる理解さえも生み出すことを約束している。特に、発現遺伝子を迅速に検出するとともに遺伝子の発現を定量化するための新規方法の出現に伴い、バイオインフォマティクスは、たとえ特定の遺伝子が生物体の生物学上果たす正確な役割について確実な知識がなくとも、潜在的な治療標的を予測するために使用され得る。
【０００４】
遺伝系のシミュレーションはシステム生物学の中心的論題である。シミュレーションは生物学的知識に基づき得るため、ネットワーク推定方法は従来未知であった関係を予測または推論することにより生物学的シミュレーションを支援できる。
【０００５】
特に、マイクロアレイ技術の開発により様々な生物体からの多数の遺伝子の発現の研究が可能となっている。多量の生データは生物体由来の数多くの遺伝子から得られ得るとともに、遺伝子発現は、突然変異、疾患または薬物のいずれかによる介入によって研究され得る。特定の遺伝子の発現が特定の疾患または特定の介入に対する応答において増加するという所見は、ひいては効力、安定性、または他の特性を増加させるべく修飾され得るという所見につながり得るとともに、該修飾化合物はアッセイで再試験され得る。この手法は、候補により影響を受ける作用または細胞経路の機序についての情報はほとんどまたは全くもたらさない。
【０００６】
創薬の第２の手法としては、既知の分子標的、典型的にはクローン化遺伝子配列または単離酵素もしくはタンパク質に対する特異的効果に関する多数の化合物の試験が挙げられる。例えば、ハイスループットアッセイが展開され得、ここでは多数の化合物が、特異的プロモーターからの転写レベルまたは同定されたタンパク質の結合を変化させる能力について試験され得る。ハイスループットスクリーニングの使用は薬物候補を同定するための強力な方法論ではあるものの、これもまた細胞または生物体レベルでの化合物の効果についての情報、特に影響を受ける実際の細胞経路に関する情報をほとんどまたは全く提供しない。
【０００７】
事実上、候補薬物の有効性および毒性の経路分析が新薬開発プロセスの相当な割合を費やす（例えば、ドヴィータ（ＤｅＶｉｔａ）ら「Ｃａｎｃｅｒ：Ｐｒｉｎｃｉｐｌｅｓ  ＆  Ｐｒａｃｔｉｃｅ  ｏｆ  Ｏｎｃｏｌｏｇｙ」第５版，１９９７年，Ｌｉｐｐｉｎｃｏｔｔ−Ｒａｖｅｎ  Ｐｕｂｌｉｓｈｅｒｓ，フィラデルフィアのなかの、オリフ（Ｏｌｉｆｆ）ら，１９９７年，「Ｍｏｌｅｃｕｌａｒ  Ｔａｒｇｅｔｓ  ｆｏｒ  Ｄｒｕｇ  Ｄｅｖｅｌｏｐｍｅｎｔ」を参照）。従って、この分析を改良する方法には相当な現在価値がある。
【０００８】
過去には、対象の生体系内部に現れる経路および機序については（薬物作用の経路を含め）、既知の入力に対する系の応答を単に観測することによりある程度一部の情報を収集することが可能とされてきた。観測される応答は典型的には、遺伝子発現（すなわち、ｍＲＮＡ存在量）および／またはタンパク質存在量であった。入力は、遺伝子突然変異（遺伝子欠失など）を含む実験的摂動、薬物処置、および成長環境条件の変化である。
【０００９】
しかしながら、系の詳細を単純に観測された入力−出力関係から推論しようという試みは、通常、絶望的な課題であった。たとえ経路仮説が利用可能であったとしても、特異的実験が適切または有効な試験または経路仮説の確証を提供するかどうかを判定することは容易ではなかった。かつ、かかる実験によってさえも、経路仮説の観点においてそれらの結果をいかに解釈するべきかについて、必ずしも認知されているとは限らなかった。
【００１０】
多くの努力および精巧な測定にもかかわらず、タンパク質およびｍＲＮＡ濃度などの単純な観測からは、生体系の経路の再構築において具体的な進展はほとんどなく、それらの時間依存相互作用についてはさらに少なかった（マカダムス（ＭｃＡｄａｍｓ）ら，１９９５年，「Ｃｉｒｃｕｉｔ  ｓｉｍｕｌａｔｉｏｎ  ｏｆ  ｇｅｎｅｔｉｃ  ｎｅｔｗｏｒｋｓ」，Ｓｃｉｅｎｃｅ  ２６９：６５０−６５６頁；レイニッツ（Ｒｅｉｎｉｔｚ）ら，１９９５年，「Ｍｅｃｈａｎｉｓｍ  ｏｆ  ｅｖｅ  ｓｔｒｉｐｅ  ｆｏｒｍａｔｉｏｎ」，Ｍｅｃｈａｎｉｓｍｓ  ｏｆ  Ｄｅｖｅｌｏｐｍｅｎｔ  ４９：１３３−１５８頁）。
【００１１】
この問題に対する一手法が他の研究領域からモデリングツールをもたらし、この問題に影響を与えた。例えば、電気工学界に身近なブール表現およびネットワークモデルが生体系に適用された（ミクレキー（Ｍｉｋｕｌｅｃｋｙ），１９９０年，「Ｍｏｄｅｌｉｎｇ  ｉｎｔｅｓｔｉｎａｌ  ａｂｓｏｒｐｔｉｏｎ  ａｎｄ  ｏｔｈｅｒ  ｎｕｔｒｉｔｉｏｎ−ｒｅｌａｔｅｄ  ｐｒｏｃｅｓｓｅｓ  ｕｓｉｎｇ  ＰＳＰＩＣＥ  ａｎｄ  ＳＴＥＬＬＡ」，Ｊ．ｏｆ  Ｐｅｄ．Ｇａｓｔｒｏｅｎｔｅｒｏｌｏｇｙ  ａｎｄ  Ｎｕｔｒｉｔｉｏｎ  １１：７−２０頁；マカダムス（ＭｃＡｄａｍｓ）ら，１９９５年，「Ｃｉｒｃｕｉｔ  ｓｉｍｕｌａｔｉｏｎ  ｏｆ  ｇｅｎｅｔｉｃ  ｎｅｔｗｏｒｋｓ」，Ｓｃｉｅｎｃｅ  ２６９：６５０−６５６頁）。一適用は、発生の間、特に逐次的な生物発生の間の遺伝子転写の制御に対するものであった（ユウ（Ｙｕｈ）ら，１９９８年，「Ｇｅｎｏｍｉｃ  Ｃｉｓ−ｒｅｇｕｌａｔｏｒｙ  ｌｏｇｉｃ：Ｅｘｐｅｒｉｍｅｎｔａｌ  ａｎｄ  ｃｏｍｐｕｔａｔｉｏｎａｌ  ａｎａｌｙｓｉｓ  ｏｆ  ａ  ｓｅａ  ｕｒｃｈｉｎ  ｇｅｎｅ」，Ｓｃｉｅｎｃｅ  ２７９：１８９６−１９０２頁）。
【発明の開示】
【発明が解決しようとする課題】
【００１２】
生物体内の生物学的経路のモデルを開発および試験するなかで認められる困難により、生物測定技術において近年なされた大幅な進歩の有効な使用が妨げられてきた。従来のバイオインフォマティクス技術は何千もの分子シグナルおよびそれらの同時調節パターンの同時研究を可能としてきた一方、細胞内の分子シグナル間の因果関係は解明できないという弱点があった。プロセスを駆動および調節するのは典型的には、個々のシグナルの独立した作用よりむしろ、細胞内で作動する全てのシグナル間の因果関係の組み合わせであることから、これは重大な障害である。
【００１３】
このように、遺伝子間の因果関係を決定するための方法を開発するため多くの努力が費やされており、それらの遺伝子は生物学的現象の中心であるとともに、それらの遺伝子の発現は研究中の生物学的プロセスにおいては周辺的である。かかる周辺的な遺伝子の発現は生物学的または病態生理学的条件のマーカーとしては有用であり得るが、かかる遺伝子が生理学的または病態生理学的条件の中心でない場合、かかる遺伝子に基づく新薬開発は労力に値しないこともある。対照的に、プロセスの中心として同定される遺伝子については、薬物または他の介入の開発が、遺伝子の変質した発現に関連する病態向け処置の開発に不可欠であり得る。
【００１４】
発現遺伝子間の関係を決定するために、ますます数学的方法が用いられるようになっている。しかしながら、遺伝子調節ネットワークを遺伝子発現データから正確に引き出すことは困難であり得る。微分方程式モデル（チェン（Ｃｈｅｎ）ら，１９９９年；ドホーン（ｄｅ  Ｈｏｏｎ）ら，２００３年）、状態空間モデル（ランゲル（Ｒａｎｇｅｌ）ら，２００４年），ブールネットワークモデル（アクツ（Ａｋｕｔｓｕ）ら，１９９８年；シュムレビッチ（Ｓｈｍｕｌｅｖｉｃｈ）ら，２００２年）およびベイジアンネットワークモデル（フリードマン（Ｆｒｉｅｄｍａｎ）ら，２０００年；イモト（Ｉｍｏｔｏ）ら，２００２年）などのいくつかの数学的方法が遺伝子調節ネットワークを推論するために使用されている。参照によって全体的に本明細書に援用される、２００２年９月２６日付けで出願された米国特許出願第１０／２５９，７２３号明細書、２００３年１１月２５日付けで出願された米国特許出願第１０／７２２，０３３号明細書、および２００３年１１月１８日付けで出願された米国特許出願第１０／７１６，３３０号明細書もまた参照されたい。
【課題を解決するための手段】
【００１５】
発明の概要
特定の実施形態において、本発明は、ノンパラメトリック回帰によるベイジアンネットワークモデルにおける経時的発現データの使用を含む。動的ベイジアンネットワークモデルにより経時的データから推定されるベイジアンネットワークモデルが遺伝子ノックダウン発現データと組み合わされるとともに、ノックダウンマイクロアレイから調節関係が推定されることで化学的化合物の系に対する影響を反映する正確な遺伝子ネットワークを構築し得る。本発明は、薬剤により影響を受ける遺伝子ネットワークを同定し、遺伝子ネットワークを包含する系内の標的遺伝子を同定するため、または生データを関係者から受け取るサービスを提供するとともに本発明に従い構築される遺伝子ネットワークモデルに基づき関係者向けの標的遺伝子を同定するために適用され得る。
【発明を実施するための最良の形態】
【００１６】
発明の詳細な説明
本発明は、異なる種類の生物学的情報を組み合わせて正確な遺伝子ネットワークを構築するための計算戦略を利用する。本発明のある実施形態において、本方法は、ドラッガブル遺伝子ネットワーク、すなわち化学的化合物により最も強く影響を受けるものを発見するために使用される。本方法を説明するため、我々は２種類のマイクロアレイデータを使用する：１つは化学的化合物による処置後の時間に伴う転写産物存在量の応答を測定することにより得られる遺伝子発現データである。もう１つは遺伝子ノックダウン発現データであり、ここでは各マイクロアレイについて１個の遺伝子がノックダウンされている。図１は我々の戦略の概念図を提供する。
【００１７】
まず我々は、ＧTとして示される遺伝子間の動的関係を、経時的データに基づき動的ベイジアンネットワークを使用して推定する。次に、遺伝子ノックダウン発現データにおいて、ノックダウンされた遺伝子の情報が分かるため、ノックダウンされた遺伝子とその被調節因子との間の可能性のある調節関係が導き出され得る。この情報をＲとする。最後に、遺伝子ネットワークＧKが、ＧT及びＲと共に、ＸKとして示される遺伝子ノックダウンデータにより、マルチソースの生物学的情報に基づくベイジアンネットワークを使用して推定される。遺伝子ネットワークをマルチソースの生物学的情報に基づき推定するための鍵となる着想はＧTおよびＲをＧKのベイズ事前確率として使用することである。本発明においては、グラフの事前確率を拡張することにより連続値として表される事前情報を使用する。遺伝子ネットワークの推定後、我々は、推定された遺伝子ネットワークから生物学的に妥当性のある情報を引き出すための、Ｇ．ＮＥＴの拡張版であるｉＮＥＴと呼ばれる遺伝子ネットワーク解析ツールを開発した。ｉＮｅｔツールは、アノテーション付きの遺伝子ネットワークの可視化を伴う遺伝子間の様々なパス探索向けのコンピューティング環境を提供する。
【００１８】
本発明の方法はドラッガブル遺伝子ネットワークを、組織特異的ネットワークのサブネットワークである有向グラフとして推定できる。この方法においては、エッジ方向が極めて重要な情報であるとともに化合物関連遺伝子の選択が必要となる。我々の方法はまた、様々な種類の生物学的データを使用することもでき、これにより推定されたネットワークの確度が高まるとともに、単に影響を受けた遺伝子の同定のみならず、ネットワークとしてのそれらの依存性の解明も可能となる。
【００１９】
遺伝子ネットワークの構築方法
本発明の一実施形態において、我々はベイジアンネットワークおよび動的ベイジアンネットワークを使用して、それぞれ遺伝子ノックダウンおよび経時的マイクロアレイデータから遺伝子ネットワークを推定する。本節においては、これらの２つのネットワークモデルを簡単に説明するとともに、次に我々がどのようにマルチソースの生物学的情報を組み合わせてより正確な遺伝子ネットワークを推定するのかについて明らかにする。
【００２０】
ｐ個の確率変数の集合である観測データＸ、Ｘ＝｛Ｘ1，．．．，Ｘp）があるとともに有向グラフＧとして示されるｐ個の確率変数間の依存性が未知であり、及びそれをＸから推定したいと仮定する。マイクロアレイデータに基づく遺伝子ネットワーク推定において、遺伝子は特異的ＲＮＡ種の存在量を表す確率変数と見なされるとともに、Ｘはマイクロアレイデータである。ベイズ手法から、観測データを条件とするグラフの事後確率を最大化することにより最適グラフが選択される。ベイズ定理により、グラフの事後確率は次式
【数１】
2009521015.tif 000002
（式中、ｐ（Ｇ）はグラフの事前確率であり、ｐ（Ｘ｜Ｇ）はＧを条件とするデータＸの尤度であり、及びｐ（Ｘ）は正規化定数であるとともにＧの選択に依存しない）で表すことができる。従って、ｐ（Ｇ｜Ｘ）に基づくグラフ選択のためにはｐ（Ｇ）を設定するとともにｐ（Ｘ｜Ｇ）を計算する必要がある。
【００２１】
グラフｐ（Ｇ）の事前確率によりマイクロアレイデータ以外の生物学的データを使用して遺伝子ネットワークを推定することが可能であるとともに尤度ｐ（Ｘ｜Ｇ）がベイジアンネットワークおよび動的ベイジアンネットワークにより、それぞれ遺伝子ノックダウンおよび経時的マイクロアレイデータから計算され得る。当業者が認識するであろうとおり、本発明は遺伝子ノックダウンデータおよび経時的マイクロアレイデータ以外の生物学的データに広範に適用され得る。次節においては、我々のｐ（Ｇ｜Ｘ）の構築方法について明らかにする。
【００２２】
ベイジアンネットワーク
ベイジアンネットワークは確率変数における因果関係を表すグラフィカルモデルである。ベイジアンネットワークにおいては、接続されたノード間のマルコフ関係をエンコードする有向非循環グラフを使用する。確率変数の集合Ｘ＝｛Ｘ1，．．．，Ｘp）があるとともに有向非循環グラフＧKを表すことによりＸにおける因果関係があると仮定する。ひいてはベイジアンネットワークにより、次式の条件付き確率の積
【数２】
2009521015.tif 000003
（式中Ｐａjは、ＧKにおけるＸjの直接の親に対応する確率変数の集合である）によって同時確率を計算することが可能となる。遺伝子ネットワーク推定においては、遺伝子は、グラフ中でノードとして示される、特異的ＲＮＡ種の存在量を表す確率変数と見なされるとともに、遺伝子間の相互作用はノード間の直接的なエッジにより表される。
【００２３】
ＸKがＮ×ｐの遺伝子ノックダウンデータ行列であり、その第（ｉ，ｊ）番目の要素ｘj｜Ｄiは第Ｄi番目の遺伝子がノックダウンされる時の第ｊ番目の遺伝子の発現データに対応するものとする（ｊ＝１，．．．，ｐおよびｉ＝１，．．．，Ｎ）。
【００２４】
ここで、第ｉ番目のノックダウンマイクロアレイがノックダウンする第Ｄi番目の遺伝子により測定されると仮定する。マイクロアレイデータは連続変数をとるため、次式の密度
【数３】
2009521015.tif 000004
（式中、Θ＝（θ’1，．．．，θ’p）はパラメータベクトルであり、Ｐａj｜Djは第ｉ番目のノックダウンマイクロアレイにより測定されるＰａjの発現値ベクトルである）を使用することにより分解（１）が表される。従って、グラフＧKの構築は条件付き確率ｆj（ｊ＝１，．．．，ｐ）をモデル化することに等しく、これは本質的に回帰問題と同じである。
【００２５】
【数４】
2009521015.tif 000005
の構築について、次式
【数５】
2009521015.tif 000006
（式中、
【数６】
2009521015.tif 000007
はＰａj｜Dj，  εj｜Djの第ｋ番目の要素であり、ｉ＝１，．．．Ｎについて〜ｉ．ｉ．ｄ．Ｎ（０、σ2）、およびｍjk（ｋ＝ｌ，．．．，｜Ｐａj｜）は
【数７】
2009521015.tif 000008
のときＢスプラインにより構築される平滑化関数である）の形のＢスプラインを伴うノンパラメトリック回帰モデルを仮定する。
【００２６】
ここで、
【数８】
2009521015.tif 000009
は、それぞれパラメータおよびＢスプラインである。
【００２７】
従って、尤度ｐ（ＸK｜ＧK）は次式
【数９】
2009521015.tif 000010
（式中、ｐ（Θ｜λ，ＧK）はハイパーパラメータλにより特定されるパラメータΘの事前分布である）により得られる。高次元積分はラプラス近似による解析法を用いて漸近的に近似され得るとともにイモト（Ｉｍｏｔｏ）らが、次式
【数１０】
2009521015.tif 000011
の形のＢＮＲＣと命名されたグラフの選択基準を定義した
（式中、
【数１１】
2009521015.tif 000012
であり、ｒはΘの次元であるとともに、Θはｌλ（Θ｜ＸK）の最頻値である）。ネットワーク構造が学習されることにより欲張り山登り法によりＢＮＲＣ（ＧK）が減少する。欲張り山登り法により得られる解は最適として保証され得ないことに留意すべきである。より良い解を見つけるためには、欲張り法を反復するとともにＧKとして最良のものを選択する。尤度ｐ（ＸK｜ＧK）が数種のネットワーク構造についてほぼ同じ値となることはかなり頻繁に起こり、様々な種類の生物学的情報に基づく有効なｐ（ＧK）の構築が鍵となる技法である。ｐ（ＧK）の構築方法は「遺伝子ネットワーク推定のためのマルチソース生物学的情報の組み合わせ」と題する節において明らかにする。
【００２８】
動的ベイジアンネットワーク
動的ベイジアンネットワークは経時的データに基づく確率変数における依存性を表す。Ｘ（ｔ）＝｛Ｘ1（ｔ），．．．，Ｘp（ｔ）｝を時刻ｔ（ｔ＝１，．．．，Ｔ）におけるｐ個の確率変数の集合とする。動的ベイジアンネットワークにおいて、ｐ個のノードを含む有向グラフは完全な２通りのグラフとして書き換えられ、これによりＸ（ｔ）からＸ（ｔ＋１）（ｔ＝１，．．．，Ｔ−１）までの直接的なエッジが可能となる。
【００２９】
従って上記に定義される２通りのグラフを推定することによりｐ個の確率変数間の因果関係の有向グラフＧTが構築される。ＧT構造のもと、ひいては次式の分解
【数１２】
2009521015.tif 000013
（式中、Ｐａj（ｔ）は、ＧTにおけるＸjの直接の親に対応する時刻ｔの確率変数の集合である）が得られる。
【００３０】
ＸTはＴ×ｐの経時的データ行列であり、その第（ｔ、ｊ）番目の要素ｘj（ｔ）が時刻ｔにおける第ｊ番目の遺伝子の発現データに対応するものとする（ｊ＝１、．．．，ｐおよびｔ＝１，．．．，Ｔ）。ベイジアンネットワークにおいて記載されるとおり、（３）における分解が次式の密度
【数１３】
2009521015.tif 000014
（式中、Ξ＝（ξ’1，．．．，ξ’p）’はパラメータベクトルであり、ｐａj（ｔ）は時刻ｔにおいて測定されるＸjの直接の親の発現値ベクトルである）を使用することにより適用される。
【００３１】
ここで、
【数１４】
2009521015.tif 000015
と定める。ベイジアンネットワークと同じ方法でＢスプラインを用いるノンパラメトリック回帰を使用してｆDBNを構築できる。それゆえ、（２）においてｆBNをｆDBNに置き換えることによりキム（Ｋｉｍ）らは、ＢＮＲＣdynamicと命名された動的ベイジアンネットワークについてのグラフ選択基準を提示し、適用に成功した。
【００３２】
遺伝子ネットワーク推定のためのマルチソース生物学的情報の組み合わせ
イモト（Ｉｍｏｔｏ）らは、より正確な遺伝子ネットワークを推定することを目的として生物学的知識を発現データと組み合わせるための一般的枠組みを提示した。イモト（Ｉｍｏｔｏ）らにおいては、生物学的知識は、例えば既知または未知といった２進値として表されるとともに、ｐ（Ｇ）を構築するために使用される。しかしながら、現実には、実際の生物学的知識における確かさは様々である。バーナード（Ｂｅｒｎａｒｄ）およびハルトミンク（Ｈａｒｔｅｍｉｎｋ）は、ｐ値の集合（連続情報）である結合部位データを使用してｐ（Ｇ）を構築した。本発明の方法において、我々は連続的および離散的事前情報を含むマルチソース情報を使用してｐ（Ｇ）を構築する。
【００３３】
Ｚkが第ｋ番目の事前情報を表す行列であり、式中、第（ｉ，ｊ）番目の要素
【数１５】
2009521015.tif 000016
が「遺伝子ｉ→遺伝子ｊ」の情報を表すものとする。例えば、（１）Ｚkとして事前ネットワークＧpriorを使用する場合、
【数１６】
2009521015.tif 000017
はｅ（ｉ，ｊ）∈Ｇpriorならば１を、または
【数１７】
2009521015.tif 000018
ならば０をとる。ここで、ｅ（ｉ，ｊ）は遺伝子ｉから遺伝子ｊへの直接的なエッジを示す。
【００３４】
（２）Ｚkとして遺伝子ノックダウンデータを使用することにより、
【数１８】
2009521015.tif 000019
はノックダウンする遺伝子ｉにより遺伝子ｊがどのように変化するかについて示す値を表す。遺伝子ｉノックダウンデータについての遺伝子ｊの対数比の絶対値を
【数１９】
2009521015.tif 000020
として使用できる。Ｇの隣接行列Ｅ＝（ｅij）l≦i,j≦pを使用して（式中、ｅ（ｉ，ｊ）∈Ｇのとき（ｅij）＝１またはその他のときに０）、次式の確率関数
【数２０】
2009521015.tif 000021
（式中、πIJ＝Ｐｒ（ｅij＝１）である）を有するｅijに対してベルヌーイ分布を仮定する。
【００３５】
πijの構築については、πij＝｛１＋ｅｘｐ（−ηij）｝-1のときの線形予測量
【数２１】
2009521015.tif 000022
（式中、ωkおよびｃk（ｋ＝１，．．．，Ｋ）は、それぞれ、重みパラメータおよび基準パラメータである）を用いるロジスティックモデルを使用する。次に事前情報Ｚk（ｋ＝１，．．．，Ｋ）に基づくグラフの事前確率を次式
【数２２】
2009521015.tif 000023
により定義する。
【００３６】
当業者が認識するであろうとおり、本発明の方法は１つまたは複数の種類のデータを事前確率として適合させることができる。グラフのこの事前確率はエッジｅ（ｉ，ｊ）（ｉ，ｊ＝１、．．．，ｐ）が互いに独立していることを仮定する。現実には、ｅijの間にはｐ（ｅij＝１）＜ｐ（ｅij＝１｜ｅki＝１）などのいくつかの依存性があるが、かかる情報をｐ（Ｇ）に加えてかかる情報の質による早熟性を検討する。
【００３７】
本方法を実施するための例示的コンピュータシステム
図９および以下の考察は、本明細書に記載される方法を実行するコンピュータプログラムに好適なコンピューティング環境についての簡潔で一般的な説明を提供することを意図している。遺伝子ネットワークを構築するための本方法は、プログラムモジュールに体系化されるコンピュータにより実行可能な命令において実行される。プログラムモジュールとしては、タスクを処理するとともにデータタイプを実行して上述される技法を実行する、ルーチン、プログラム、オブジェクト、コンポーネント、およびデータ構造が挙げられる。
【００３８】
図９はデスクトップコンピュータの典型的な構成を示すが、本発明は、マルチプロセッサシステム、マイクロプロセッサベースの、またはプログラム可能な家庭用電化製品、ミニコンピュータ、メインフレームコンピュータなどを含む、他のコンピュータシステム構成で実行されてもよい。本発明はまた分散コンピューティング環境で使用されてもよく、ここではタスクは演算装置により並列に実行されることでパフォーマンスが高まる。例えば、大きな集合の非線形モデルの有効性測定に関連するタスクは、複数のコンピュータ上、単一のコンピュータにおける複数のプロセッサ上、またはその双方で同時に実施され得る。分散コンピューティング環境においては、プログラムモジュールはローカルおよび遠隔の双方の記憶装置内に位置する。
【００３９】
図９に示されるコンピュータシステムは本発明の実行に好適であるとともに、処理ユニット２９２１、システムメモリ２９２２、および、システムメモリの処理ユニット２９２１への接続を含む、様々なシステムコンポーネントを相互接続するシステムバス２９２３を伴うコンピュータ２９２０を含む。システムバスは、メモリバスまたはメモリコントローラ、周辺機器用バス、およびバスアーキテクチャを使用するローカルバスを含む、数種のバス構造のうち任意のものを備えてもよい。システムメモリは、読取専用メモリ（ＲＯＭ）２９２４およびランダムアクセスメモリ（ＲＡＭ）２９２５を含む。
【００４０】
不揮発性システム２９２６（例えば、ＢＩＯＳ）はＲＯＭ２９２４内に格納されるとともに基本ルーチンを含んで起動中などにパーソナルコンピュータ２９２０内の要素間の情報を転送することができる。パーソナルコンピュータ２９２０は、ハードディスクドライブ２９２７、例えば、リムーバブルディスク２９２９から読み出す、またはそこへ書き込むための磁気ディスクドライブ２９２８、および、例えば、ＣＤ−ＲＯＭディスク２９３１を読み取る、または他の光メディアから読み出す、またはそこへ書き込むための光ディスクドライブ２９３０をさらに含むことができる。ハードディスクドライブ２９２７、磁気ディスクドライブ２９２８、および光ディスクドライブ２９３０は、それぞれ、ハードディスクドライブインタフェース２９３２、磁気ディスクドライブインタフェース２９３３、および光学ドライブインタフェース２９３４によりシステムバス２９２３に接続される。ドライブおよびそれらの関連コンピュータ読取可能メディアは、データ、データ構造、コンピュータにより実行可能な命令（動的リンクライブラリおよび実行可能ファイルなどのプログラムコードを含む）などのパーソナルコンピュータ２９２０用の不揮発性ストレージを提供する。
【００４１】
上記のコンピュータ読取可能メディアの記述はハードディスク、リムーバブル磁気ディスク、およびＣＤを参照しているが、これはまた、磁気カセット、フラッシュメモリカード、デジタルビデオディスクなどの、コンピュータにより読取可能な他の種類のメディアも含み得る。
【００４２】
オペレーティングシステム２９３５、１つまたは複数のアプリケーションプログラム２９３６、他のプログラムモジュール２９３７、およびプログラムデータ２９３８を含む、数多くのプログラムモジュールがドライブおよびＲＡＭ２９２５内に格納される。ユーザはコマンドおよび情報をパーソナルコンピュータ２９２０にキーボード２９４０およびマウス２９４２などのポインティングデバイスを通じて入力してもよい。他の入力装置（図示せず）としては、マイクロフォン、ジョイスティック、ゲームパッド、衛星放送受信アンテナ、スキャナなどが挙げられ得る。
【００４３】
これらおよび他の入力装置はしばしばシステムバスに連結されるシリアルポートインタフェース２９４６を通じて処理ユニット２９２１に接続されるが、パラレルポート、ゲームポート、またはユニバーサルシリアルバス（ＵＳＢ）などの、他のインタフェースにより接続されてもよい。モニタ２９４７または他の種類のディスプレイ装置もまたディスプレイコントローラまたはビデオアダプタ２９４８などのインタフェースを介してシステムバス２９２３に接続される。モニタに加え、パーソナルコンピュータは典型的には、スピーカおよびプリンタなどの他の周辺出力装置（図示せず）を含む。
【００４４】
上記のコンピュータシステムは単に例として提供されるものである。本発明は多種多様な他の構成で実行され得る。さらに、遺伝子の関連性の定量化に関するデータを収集および解析するうえでは多種多様な手法が可能である。例えば、必要に応じて異なるコンピュータシステム上でデータが収集され、非線形モデルが構築され、モデルの有効性が測定され、および結果が提示され得る。加えて、様々なソフトウェア態様がハードウェアで、またその逆も然り、実装され得る。
【実施例】
【００４５】
以下の実施例は、例示としてのみ、および限定としてではなく提供される。当業者は、変更または修正されても本質的に同様の結果を与え得るであろう様々な主要ではないパラメータを容易に認識するであろう。
【００４６】
ドラッガブル遺伝子ネットワークを生成するためのヒト内皮細胞への適用
本発明の動作方法を実証するため、我々はヒト内皮細胞からの発現データを解析し、ヒト内皮細胞転写産物の抗高脂血症薬フェノフィブラートによる処置に対する応答を明らかにする新規の経時的データを生成した。ヒト内皮細胞における２７０個の遺伝子のノックダウン実験から新規データも生成した。フェノフィブラート関連遺伝子ネットワークは、本発明の方法によりフェノフィブラートの経時的データおよび２７０個の遺伝子のノックダウン発現データに基づき推定される。推定された遺伝子ネットワークは、フェノフィブラートにより活性化されることで知られるＰＰＡＲ−αに関連する遺伝子調節関係を明らかにする。我々のコンピュータ解析は、遺伝子ノックダウンおよび薬物投与の経時的マイクロアレイに基づくこの計算戦略が、ドラッガブル遺伝子の発見において、新規の、および改良されたツールを提供することを示唆する。
【００４７】
実施例１：フェノフィブラートの経時的データ
我々は２５μＭのフェノフィブラートに対するヒト内皮細胞遺伝子の時間応答を測定した。２０，４６９個のプローブの発現レベルがＣｏｄｅＬｉｎｋ（商標）Ｈｕｍａｎ  Ｕｎｉｓｅｔ  Ｉ  ２０Ｋにより６時点（０、２、４、６、８および１８時間）で測定された。ここで時刻０は本観測の開始点およびフェノフィブラートに対する曝露の直前を意味する。さらに、我々はこの経時的データを重複データとして測定することにより実験の質を確認した。
【００４８】
我々のフェノフィブラートの経時的データは重複データであるとともに６時点を含むことから、経時的データセットの作成には２6＝６４通りの可能な組み合わせがある。６４個のデータセットにおける親−子関係に対しては同じ回帰関数をフィットさせるべきである。この制約のもと、ノンパラメトリック回帰モデルを６４個のデータセットの接続されたデータに対してフィットさせることを検討する。すなわち、遺伝子ｉ→遺伝子ｊを検討するならば、モデルｘj(c)（ｔ）＝ｍj（ｘi(c)（ｔ−ｌ））＋εj（ｔ）（式中、ｘj(c)（ｔ）は第ｃ番目のデータセット（ｃ＝１，．．．，６４）の時刻ｔにおける遺伝子ｊの発現データである）をフィットさせることとなる。ベイジアンネットワークにおいて、推定されたエッジの信頼性はブートストラップ法を使用して測定され得る。経時的データについては、ブロックリサンプリングなどブートストラップ法のいくつかの修正が提示されているが、これらの方法を短い時間経過により生成される少数のデータ点に対して適用することは困難である。
【００４９】
しかしながら、上記の経時的モデリングを使用することにより、我々はブートストラップに基づく方法を次のとおり定義できる：Ｄ＝｛Ｄ（１），．．．，Ｄ（６４））を全遺伝子の組み合わせ経時的データとする。無作為にＤ（ｃ）を復元再抽出するとともにブートストラップ標本Ｄ*＝｛Ｄ*（１），．．．，Ｄ*（６４）｝を定義する。次にＤ*に基づき遺伝子ネットワークを再推定する。ブートストラップ反復を１０００回繰り返すと
【数２３】
2009521015.tif 000024
（式中
【数２４】
2009521015.tif 000025
は第Ｂ番目のブートストラップ標本に基づき推定されたグラフである）が得られる。ひいてはエッジの推定された信頼性が
【数２５】
2009521015.tif 000026
のときの第１の事前情報Ｚ1の行列表現として使用され得る。
【００５０】
実施例２：ｓｉＲＮＡによる遺伝子ノックダウンデータ
例示的遺伝子ネットワークの構築について、我々はｓｉＲＮＡを使用することにより２７０個の遺伝子のノックダウンデータを新規に作成した。２０，４６９個のプローブを各ノックダウンマイクロアレイについてｓｉＲＮＡトランスフェクションの２４時間後にＣｏｄｅＬｉｎｋ（商標）Ｈｕｍａｎ  Ｕｎｉｓｅｔ  Ｉ  ２０Ｋにより測定した。ノックダウン遺伝子は主に転写因子およびシグナル伝達分子であった。
【００５１】
【数２６】
2009521015.tif 000027
を第ｉ番目のノックダウンマイクロアレイの元の強度ベクトルとする。各マイクロアレイの正規化された発現値について、発現値ベクトルｖ＝（ｖ1，．．．，ｖp）’の中央値を対照データとして計算した（式中、ｖj＝中央値）。
【００５２】
【数２７】
2009521015.tif 000028
ローエス（ｌｏｅｓｓ）正規化法をＭＡ変換データに適用し、逆変換を正規化された
【数２８】
2009521015.tif 000029
に対し適用することにより正規化された強度ｘj|D_iが得られた。正規化された
【数２９】
2009521015.tif 000030
を対数比と呼ぶ。
【００５３】
２７０個の遺伝子のノックダウンマイクロアレイデータにおいて、各マイクロアレイについてどの遺伝子がノックダウンされているかは既知である。従って、遺伝子Ｄiをノックダウンした場合に、発現レベルを有意に変化させる遺伝子が遺伝子Ｄiの直接的な被調節因子と考えられ得る。我々はこの情報を次のとおり補正対数比を計算することにより測定した：対数比の変動は標本および対照の強度のそれらの和に依存する。正規化されたＭＡ変換データから、次に条件付き分散ｓj＝Ｖａｒ［ｌｏｇ（ｘj|D_i／ｖj）｜ｌｏｇ（ｘj|D_i・ｖj）］が得られ、および対数比はＶａｒ_ｚi(2)j＝１を満たす補正されたｚi(2)j＝ｌｏｇ（ｘj|D_i／ｖj）であり得る。
【００５４】
実施例３：遺伝子ネットワークモデルを生成するためのフェノフィブラートの経時的データ、ｓｉＲＮＡによる遺伝子ノックダウンデータ、およびノックダウンデータ行列の組み合わせ
【００５５】
【表１】
2009521015.tif 000031
【００５６】
フェノフィブラートの経時的データおよび２７０個の遺伝子のノックダウンデータからフェノフィブラート関連遺伝子ネットワークを推定するため、まずフェノフィブラートと関連する可能性のある遺伝子の集合を次のとおり定義した：第１に、その分散の補正された対数比｜ｌｏｇ（ｘj|D_i／ｖj）／ｓj｜が１．５より大きい遺伝子の集合を各時点から抽出した。次にＧＯ  Ｔｅｒｍ  Ｆｉｎｄｅｒを使用して有意な選択遺伝子クラスタを特定した。表１は１８時間において有意な遺伝子クラスタを示す。第１の列は発現値がどのように変化したかを示す、すなわち、
【数３０】
2009521015.tif 000032
はそれぞれ、「過剰発現した」および「抑制された」ことを意味する。
【００５７】
【数３１】
2009521015.tif 000033
を伴うクラスタのＧＯアノテーションは主に細胞周期に関連し、これらのクラスタ内の遺伝子は偏在的に発現するとともにこれは共通の生物学的機能である。他方、
【数３２】
2009521015.tif 000034
を伴うクラスタのＧＯアノテーションは主に脂質代謝に関連する。生物学においては、フェノフィブラートが曝露のおよそ１２時間後に作用することが報告されている。遺伝子選択についての我々の第１の解析は、フェノフィブラートが脂質代謝に関連する遺伝子に影響を及ぼすことを示唆するとともにこれは生物学的事実と一致する。我々はまた、８時間時点マイクロアレイからの遺伝子にも焦点を当てた。８時間時点マイクロアレイからの選択された遺伝子中には特異的な機能を伴うクラスタは所見され得なかったが、脂質代謝に関連する遺伝子が一部存在した。従って我々は８および１８時間時点マイクロアレイからの遺伝子を使用した。最後に、２６７個のノックダウン遺伝子（３個の遺伝子は我々のチップ上に発見されなかった）を上記の選択された遺伝子に追加し、合計１１９２個の遺伝子がフェノフィブラートに関連する可能性のある遺伝子として定義されたとともに次のネットワーク解析に使用された。
【００５８】
推定された動的ネットワークおよびノックダウン遺伝子情報をそれぞれ第１および第２の事前情報Ｚ1およびＺ2の行列表現に変換することにより、Ｚ1、Ｚ2およびノックダウンデータ行列ＸKに基づき遺伝子ネットワーク
【数３３】
2009521015.tif 000035
を推定した。推定された遺伝子ネットワークから生物学的情報を引き出すため、我々はまず脂質代謝関連遺伝子に焦点を当てたが、これはこの機能に関連するクラスタが１８時間マイクロアレイにおいて有意に変化したためである。推定された遺伝子ネットワークにおいては、４２個の脂質代謝関連遺伝子があったが、それらのなかでＰＰＡＲ−α（ヒト（Ｈｏｍｏ  ｓａｐｉｅｎｓ）ペルオキシソーム増殖剤応答性受容体α）が唯一の転写因子である。実際に、ＰＰＡＲ−αはフェノフィブラートの標的として知られている。従って、次に我々はＰＰＡＲ−αの下流ノードに焦点を当てた。
【００５９】
図２はＰＰＡＲ−αの下流ノード（４９１個の遺伝子）のコンピュータレンダリングを提供する。ここでは、ＰＰＡＲ−αの下流４段にある遺伝子がＰＰＡＲ−αの候補被調節因子であると考える。ＰＰＡＲ−αの候補被調節因子のなかには、２１個の脂質代謝関連遺伝子およびＰＰＡＲ−αに関連することが以前実験的に同定された１１個の分子がある。実際に、ＰＰＡＲ−αはフェノフィブラートにより活性化されることが知られており、これは我々のネットワークの確度を支持するものである。
【００６０】
特に、ＰＰＡＲ−αを有する１つのサブネットワークをルートノードとして図３に示す。ＰＰＡＲ−αを標的とするフェノフィブラートの薬効の１つは、ＬＤＬコレステロールの低減である。ＬＤＬＲおよびＶＬＤＬＲは主にコレステロールの輸送に寄与するとともに、我々の推定ネットワークにおいてそれらはＰＰＡＲ−αの子、つまりＰＰＡＲ−αの候補被調節因子である。ＬＤＬＲに関しては、ＰＰＡＲ−αとの関係性が報告されている。そのうえ、コレステロール代謝に関連する数個の遺伝子が我々のネットワークにおいてはＰＰＡＲ−αの子である。我々はＳＴＡＴ５ＢおよびＧＬＳがＰＰＡＲ−αの子であることを我々のネットワークから引き出すことができたが、これについてはそれらのＰＰＡＲ−αとの調節関係が報告されている。従って、既知のＰＰＡＲ−α被調節因子に関わる多くの直接的および間接的な関係が内皮細胞内でフェノフィブラート処置により誘発されることを我々のネットワークが示しているのは、驚くことではない。ＰＰＡＲ−αの上流ノードにおいては、ヘテロ二量体を形成するＰＰＡＲ−αおよびＲＸＲ−αが親を共有する。本発明の方法を使用して、我々はフェノフィブラート関連遺伝子ネットワークを生成するとともに、それにより事前の生物学的知識なしにＰＰＡＲ−αがフェノフィブラート調節の鍵分子の１つであると推定することができた。
【００６１】
薬理ゲノミクスの観点から、ドラッガブル遺伝子ネットワークを知ることは非常に重要である。我々の遺伝子ネットワークは、化学的化合物の作用機序を予測し、より有効な薬物標的を発見し、および所与の薬物に対する曝露による生じる副作用を予測するうえでの可能性を有する。本発明は化学的化合物に関連する遺伝子ネットワークを発見するための計算法を提供する。我々は本目的上、遺伝子ノックダウンマイクロアレイデータおよび経時的応答マイクロアレイデータを使用するとともに観測データから得られる複数の情報を組み合わせることによりベイズ統計の枠組みのもと正確な遺伝子ネットワークを推定する。我々はヒト内皮細胞における遺伝子ネットワーク推論の実例を使用して本方法の全体プロセスを説明した。
【００６２】
フェノフィブラートの経時的データおよびヒト内皮細胞内の遺伝子ノックダウンからのデータを使用して、我々は、ＰＰＡＲ−αの既知のアゴニストである薬物フェノフィブラートに関連する遺伝子ネットワークを推定することに成功した。推定された遺伝子ネットワークにおいて、ＰＰＡＲ−αは脂質代謝関連遺伝子を含む多くの直接的および間接的被調節因子を有するとともにこの結果はＰＰＡＲ−αが推定されたフェノフィブラート関連ネットワークのトリガとして働くことを示す。ＰＰＡＲ−αの候補被調節因子には多くの既知の関係性があるとともに我々は推定されたネットワークにおいてＰＰＡＲ−αとＲＸＲ−αとの間の関係を発見できた。ペルオキシソーム増殖剤応答性受容体（ＰＰＡＲ）は内皮細胞および数種の他の細胞型により発現するリガンド活性化転写因子である。これらは天然脂肪酸および合成フィブラートなどのリガンドにより活性化される。活性化されると、これらはレチノイドＸ受容体（ＲＸＲ）とヘテロ二量体化して標的遺伝子の転写を活性化する。これらの遺伝子の多くは、炭水化物およびグルコース代謝を制御するとともに炎症反応を下方調節するタンパク質をコードする。
【００６３】
ヒト内皮細胞アポトーシスの研究への適用
本発明はまた、ヒト内皮細胞（ＥＣ）アポトーシス過程の研究においても適用された。
【００６４】
実施例４：アポトーシスの経時的遺伝子アレイデータ
ＥＣアポトーシス中のトランスクリプトーム変化の動態を理解するため、我々は経時実験を実施した。ＨＵＶＥＣ（１０ドナー由来のプール培養物）が上述される研究において使用されるものと同一のＳＦＤ条件に曝露された。ＲＮＡが培養物から調製された後アポトーシスが誘導された（時刻０）とともに０．５、１．５、３、６、９、１２および２４時間後、ＣｏｄｅＬｉｎｋ  Ｈｕｍａｎ  Ｕｎｉｓｅｔ  ２０Ｋ遺伝子アレイとハイブリダイズされた。この実験は独立して３回繰り返された。ＥＣアポトーシス中の転写産物の調節は図４（ａ）〜（ｄ）に示される散布図で可視化され得る。次にＥＣアポトーシスの時間経過にわたり一貫して調節された２７６個の転写産物の一部をさらなる解析のため選択した。
【００６５】
これを行うため、３回全ての反復実験において３点以上の隣接する時点でＺ≧１．５σだけ調節されなかった全ての転写産物を我々の分析から除外した（方法を参照）。これらの２７６個の転写産物から、ｋ平均クラスタリングが使用され、同様の経時的プロファイルを伴う８群の転写産物が（未分類群に加え）選択された（図５）。これらのプロファイルおよび図５ａ〜ｄに示される散布図から、数多くの転写産物が１．５時間以降に調節され始めるとともに概して変化速度は１２時間後に低下するように思われることが確認され得る。前の研究において同定された個々の転写産物をプロットした時、我々は、概して、増殖因子をコードする転写産物（Ａｎｇ−２およびＩＬ−８など）が最も早期に調節される転写産物の１つであった一方、細胞周期に関連する転写産物（サイクリンＡ２、Ｈ、およびＥ、ＣＤＣ６、ＣＤＣ２８およびキネシン様紡錘体タンパク質など）は後に調節された（図６）ことに注目した。これらの転写産物の多くは２８時間ＳＦＤ後に調節されることが、我々の以前のアフィメトリックス研究により検証されている。これらのデータは、一部の機能クラスがＥＣ培養物のＳＦＤ誘導アポトーシス中、共通の調節パターンに従うことを示唆している。より精緻な分析（実施例５を参照）からは、これらの転写産物の共通の上流調節因子が示唆される。
【００６６】
実施例５：遺伝子ネットワーク生成のためのアポトーシスの経時的遺伝子アレイデータおよび遺伝子破壊実験からのデータの組み合わせ
我々は遺伝子ネットワーク（図８（ｃ））を上述されるアポトーシスの経時的遺伝子アレイデータから生成した。このネットワークは、実施例４のＨＵＶＥＣアポトーシスの経時的遺伝子アレイデータの中央値に基づき動的ベイジアンネットワーク推論法（キム（Ｋｉｍ）ら、２００３年）を使用して生成されたもので、このネットワークは８回の破壊実験からの新規の遺伝子アレイデータを（「ベイズプライア」として）本発明により提供される方法で組み込んだ点で図７（ａ）に示されるものとは異なる。これらの８回の破壊実験において、細胞周期制御に関連する特異的ＲＮＡ（ＣＤＣ４５Ｌ、ＣＣＮＥ１、ＣＥＮＰＡ、ＣＥＮＰＦ、ＣＤＣ２、ＣＤＣ２５Ｃ、ＭＣＭ６およびＣＤＣ６）の存在量はｓｉＲＮＡプールを使用するＨＵＶＥＣにおいて６５％超低減された。遺伝子アレイ解析から、これらのｓｉＲＮＡ処置により多数の転写産物が調節された（おそらく、標的ＲＮＡ下方調節の結果として、または可能性としてある場合においてはまた、予期せぬ標的外のｓｉＲＮＡ効果の結果としても）ことが示された。ｓｉＲＮＡ処置のＨＵＶＥＣトランスクリプトームに対する効果の例が図８（ａ）〜（ｂ）に示される。ベイズプライアとしてのこの破壊情報の組込みにより修正される推定経時的遺伝子ネットワークが図８（ｃ）に示される。
【００６７】
図８（ｃ）に示される修正遺伝子ネットワークに含まれる４８８本のエッジのうち、３３８本は図７（ａ）に示される未修正の遺伝子ネットワークと共有される。修正ネットワークには存在するが未修正ネットワークには存在しない１５０本のエッジのうち９４本は、ｓｉＲＮＡ破壊実験において標的とされた８個のＲＮＡのうちの１個を親として有する。これらの特定のエッジが破壊遺伝子アレイデータから受け継がれた因果情報を有するように思われるのは、それらの全てがｓｉＲＮＡ処置による破壊親の子に対する効果を正確に予測するためである。従って、破壊データのベイズプライアとしての包含は、動的な経時的遺伝子ネットワークの、転写産物間の特異的な有向関係を予測する能力を高める働きをする。当業者はまた、本発明が生物学的データベースアノテーションおよびプロテオミクス実験からのデータの組込みに適用され得るとともに、それにより生成される遺伝子ネットワークの予測力をさらに高め得ることも認識するであろう。
【００６８】
前述の発明は、理解を明確にする目的から例示および例としていくらか詳細に記載されているが、添付の特許請求の範囲の精神または範囲から逸脱することなく本発明に対し特定の変更および修正がなされてもよいことは本発明の教示に照らして当業者には直ちに明らかである。本明細書に引用される全ての刊行物、特許、および特許出願は、あらゆる目的上、全体として参照により本明細書によって援用される。
【００６９】
Akutsu, T., Kuhara, S., Maruyama, O. and Miyano, S. 1998 A System for Identifying Genetic Networks from Gene Expression Patterns Produced by Gene Disruptions and Overexpressions. Genome Inform Ser Workshop Genome Inform 9, 151-160.
Chen, T., He, H. L. and Church, G. M. 1999 Modeling gene expression with differential equations. Pac Symp Biocomput, 29-40.
de Hoon, M. J., Imoto, S., Kobayashi, K., Ogasawara, N. and Miyano, S. 2003 Inferring gene regulatory networks from time-ordered gene expression data of Bacillus subtilis using differential equations. Pac Symp Biocomput, 17-28.
Friedman, N., Linial, M., Nachman, I. and Pe'er, D. 2000 Using Bayesian networks to analyze expression data. J Comput Biol 7, 601-20.
Imoto, S., Goto, T. and Miyano, S. 2002 Estimation of genetic networks and functional structures between genes by using Bayesian networks and nonparametric regression. Pac Symp Biocomput, 175-86.
Imoto, S., Tamada, Y., Araki, H., Yasuda, K., Print, C, Charnock-Jones, D., Sanders D, Savoie, C., Tashiro, K., Kuhara, S. and Miyano, S. 2006 Computational strategy for discovering druggable gene networks from genome-wide RNA expression profiles. Pacific Symposium on Biocomputing 11, 559-571.
Johnson, N. A., Sengupta, S., Saidi, S. A., Lessan, K., Charnock-Jones, S. D., Scott, L., Stephens, R., Freeman, T. C., Tom, B. D., Harris, M., Denyer, G., Sundaram, M., Sasisekharan, R., Smith, S. K. and Print, C. G. 2004 Endothelial cells preparing to die by apoptosis initiate a program of transcriptome and glycome regulation. Faseb J 18, 188-90.
Rangel., C, Angus, J., Ghahramani, Z., Lioumi, M., Sotheran, E., Gaiba, A., Wild, D. L. and Falciani, F. 2004 Modeling T-cell activation using gene expression profiling and statespace models. Bioinformatics 20, 1361-72.
Schoenfeld, J., Lessan, K., Johnson, N. A., Charnock-Jones, D. S., Evans, A., Vourvouhaki, E., Scott, L., Stephens, R., Freeman, T. C., Saidi, S. A., Tom, B., Weston, G. C., Rogers, P., Smith, S. K. and Print, C. G. 2004 Bioinformatic analysis of primary endothelial cell gene array data illustrated by the analysis of transcriptome changes in endothelial cells exposed to VEGF-A and PlGF. Angiogenesis 7, 143-56.
Shmulevich, I., Dougherty, E. R., Kim, S. and Zhang, W. 2002 Probabilistic Boolean Networks: a rule-based uncertainty model for gene regulatory networks. Bioinformatics 18, 261-74.
T. Akutsu et al., Pac. Symp. Biocomput., 4:17-28, 1999.
K. Basso et al., Nat. Genet., 37:382-390, 2005.
A. Bernard and A.J. Hartemink, Pac. Symp. Biocomput., 10:459-470, 2005
A. Cabrero et al., Curr. Drug Targets Inflamm. Allergy, 1:243-248, 2002.
T. Chen et al., Pac. Symp. Biocomput., 4:29-40, 1999.
D. di Bernardo et al., Nat. Genet., 37:382-390, 2005.
N. Friedman et al., J. Comp. Biol., 7:601-620, 2000.
K. Goya et al., Arterioscler. Thromb. Vasc. Biol., 24:658-663, 2004.
A.J. Hartemink et al., Pac. Symp. Biocomput., 7:437-449, 2002.
K. Hayashida et al., Biochem. Biophys. Res. Commun., 323:1116-1123, 2004.
D. Heckerman et al., Machine Learning, 20:197-243, 1995.
S. Imoto et al., Pac. Symp. Biocomput., 7:175-186, 2002.
S. Imoto et al., J. Bioinform. Comp. Biol., 2:77-98, 2004.
S. Imoto et al., J. Bioinform. Comp. Biol., 1:459-474, 2003.
K.K. Islam et al., Biochim. Biophys. Acta., 1734:259-268, 2005.
S. Kersten et al., FASEB J., 15:1971-1978, 2001.
S. Kim et al., Biosystems, 75:57-65, 2004.
T.I. Lee et al., Science, 298:799-804, 2002.
M.J. Marton et al., Nat. Med., 4:1293-1301, 1998.
C.J. Savoie et al., DNA Res., 10:19-25, 2003.
J.M. Shipley and D.J. Waxman, Mol. Pharmacol., 64:355-364, 2003.
Y. Tamada et al., Genome Informatics, 16:182-191, 2005.
E.P. van Someren et al., Pharmacogenomics, 3:507-525, 2002.
【図面の簡単な説明】
【００７０】
【図１】本発明の遺伝子ネットワーク構築方法における様々なステップを示す概略図である。
【図２】本発明の方法に従い構築される遺伝子ネットワーク内のＰＰＡＲ−αの下流ノードのコンピュータ生成による可視化である。
【図３】ＰＰＡＲ−αに関連するサブネットワークを示す。
【図４】ＥＣアポトーシス中のトランスクリプトームの経時変化を示す。散布図は点により表される各転写産物を示し、各ｘ軸上に時刻０における存在量が示されるとともにｙ軸上に他の時点における存在量が示される；図４（ａ）；１．５時間ＳＦＤ、図４（ｂ）：６時間ＳＦＤ、図４（ｃ）：９時間ＳＦＤおよび図４（ｄ）：２４時間ＳＦＤ。調節されない当該転写産物がほぼ斜め方向に残留する。培養物がＳＦＤ条件に２４時間曝露された時、時刻０の健常培養物と比較して上方調節されたと思われる転写産物が白色で強調表示されている。これらの転写産物の時間に伴う漸進的な上方調節が確認され得る。
【図５】ｋ平均クラスタリングの適用を示し、これにより異なる時間的調節パターンに従う８群の転写産物（６８５個の高度に調節された転写産物の最終候補リストから）が明らかとなった。
【図６】転写産物存在量のＳＦＤ依存調節の動態を示す。値は、３通りの経時的データの正規化された転写産物存在量の中央値の倍数変化（０時間に対する）を表す。負値は下方調節を表す。正値は上方調節を表す。
【図７】動的な経時的遺伝子調節ネットワークの例を示し、（ａ）は３通りのアポトーシスの経時的データの中央値から生成される動的ベイジアン遺伝子ネットワークを表すグラフを示す。点は転写産物（「ノード」）を表すとともに点間の線は潜在的な因果相互作用（「エッジ」）を表す。（ｂ）は、ネットワーク内の親ＲＮＡ転写産物（ＣＤＫＮ１Ｃ、黒）およびその推定上の子（Ｃ１ＱＴＮＦ５、緑；ＡＫＲ１Ｃ３、青；ＭＬＦ１、橙およびＳＵＶ３９Ｈ１、赤）についてのアポトーシスの経時的遺伝子アレイ転写産物存在量データを示すグラフを示す。
【図８】ｓｉＲＮＡ破壊実験からの事前情報を含めることにより修正される動的な経時的遺伝子調節ネットワークの例を示す。（ａ）は、非トランスフェクトＨＵＶＥＣにおける転写産物存在量（ｘ軸）をルシフェラーゼに対する対照ｓｉＲＮＡをトランスフェクトされるモックＨＵＶＥＣにおける転写産物存在量（ｙ軸）と比較する散布図である−転写産物存在量の調節はほとんど起こらなかったと思われる。（ｂ）は、ルシフェラーゼに対するｓｉＲＮＡをトランスフェクトされるＨＵＶＥＣにおける転写産物存在量（ｘ軸）をＮＦκＢｐ１０５に対するｓｉＲＮＡをトランスフェクトされるＨＵＶＥＣにおける転写産物存在量（ｙ軸）と比較する散布図である−ＮＦκＢｐ１０５（円で囲まれている）は４倍超下方調節されたとともに多数の追加的転写産物における存在量もまた調節された。（ｃ）は、図８（ａ）〜（ｂ）に示されるものと同様の、８回のｓｉＲＮＡ破壊実験からの遺伝子アレイデータをベイズプライアとして組み込むことにより修正された、アポトーシスの生成された経時的遺伝子ネットワークを表すグラフを示す。
【図９】本発明の実装に動作環境として好適なコンピュータシステムを図示するブロック図である。
【図１０】ベイジアンネットワークおよび動的ベイジアンネットワークを使用して遺伝子ネットワークモデルを構築するための本発明の方法を示すフローチャートである。
【図１】
2009521015.tif 000036
【図２】
2009521015.tif 000037
【図３】
2009521015.tif 000038
【図４】
2009521015.tif 000039
【図５】
2009521015.tif 000040
【図６】
2009521015.tif 000041
【図７】
2009521015.tif 000042
【図８】
2009521015.tif 000043
【図９】
2009521015.tif 000044
【図１０】
2009521015.tif 000045