(57)【要約】  遺伝子発現測定から遺伝子ネットワークを正確に推定することは、生物情報科学分野における大きい課題である。探索空間を生物学的に有意味な部分空間まで縮小して、生物学的に関連する情報により制約された推定モデルを用いて、線形時間で部分空間内における最適解を発見する一般的な方法を提示する。イースト菌および枯草菌データへ適用して本方法の有効性を示した。また、遺伝子ネットワーク関係のデータおよび結果を与えて保存すべく適合されたシステムおよび記憶媒体を提供する。
2007519988.tif 000002
【特許請求の範囲】
【請求項１】
遺伝子ネットワークを推定する方法であって、
（ａ）探索空間の定義を含む、有機体の可能な遺伝子ネットワークの推定モデルを提供するステップと、
（ｂ）前記探索空間の生物学的に関連する部分空間を選択するステップと、
（ｃ）小規模な遺伝子ネットワークを最適に計算するアルゴリズムを繰り返し適用することにより、前記選択された部分空間の最適解を計算するステップと、を含む方法。
【請求項２】
前記推定モデルがベイジアンネットワーク推定モデルである、請求項１に記載の方法。
【請求項３】
前記生物学的に関連する部分空間が、前記有機体の代謝経路に関係する遺伝子群を含む、請求項１に記載の方法。
【請求項４】
前記アルゴリズムが、
（ａ）全てのｇ∈Ｇに対してＦ（ｇ，φ）＝ｓ（ｇ，φ）を計算するステップと、
（ｂ）全てのＡ⊆Ｇ、Ａ≠φ、および全てのｇ∈Ｇに対してＦ（ｇ，Ａ）を、
【数１】
2007519988.tif 000029
として計算するステップと、
（ｃ）Ｍ（φ）＝φに設定するステップと、
（ｄ）全てのＡ⊆Ｇ、Ａ≠φに対して、
（ｉ）ｇ＊＝ａｒｇ  ｍｉｎｇ⊆Ａ（Ｆ（ｇ，Ａ−｛ｇ｝）＋ＱＡ−｛ｇ｝（Ｍ（Ａ−｛ｇ｝）））を計算すると共に、
（ｉｉ）全てのｌ≦ｉ＜｜Ａ｜に対して、Ｍ（Ａ）（ｉ）＝Ｍ（Ａ−｛ｇ＊｝）（ｉ）且つＭ（Ａ）（｜Ａ｜）＝ｇ＊を設定するステップと、
（ｅ）ＱＧ（Ｍ（Ｇ））を返すステップと、を含む、請求項１に記載の方法。
【請求項５】
前記アルゴリズムが、
（ａ）ステップ１およびステップ２でのＦの計算において、全てのｇ∈Ｓｉおよび全てのＡ⊆Ｃｇに対してＦ（ｇ，Ａ）だけを計算するステップと、
（ｂ）ステップ４ａの項Ｆ（ｇ，Ａ−｛ｇ｝）をＦ（ｇ，（Ｃｇ＝Ｓｉ）∪（Ｃｇ∩Ａ））で置換するステップと、に従い変更される、請求項４に記載の方法。
【請求項６】
最適ネットワークＮが定義ｓｃｏｒｅ（Ｎ）＝Σｇ∈Ｇs（ｇ，ＰＮ（ｇ））を有する、請求項１に記載の方法。
【請求項７】
前記アルゴリズムが、
（ａ）Ｇ内の遺伝子群をクラスタリングして、どのクラスタも遺伝子ｃ個より大きくならないようにするステップと、
（ｂ）前記クラスタを降順Ｃ１，・・・，Ｃｎにソートするステップと、
（ｃ）各ｉ∈｛１，・・・，ｎ｝および各ｇ∈Ｃｉに対して、Ｃ１∪．．．∪Ｃｎから親候補を最大ｍ個まで選択するステップと、
（ｄ）定理１．２を用いて最適遺伝子ネットワーク・モデルを計算するステップと、を含む、請求項１に記載の方法。
【請求項８】
前記アルゴリズムが、
（ａ）｜Ｃｉ｜≦ｃであるグループＣｉにおけるＧの遺伝子群をグループ化し、それらを生物学的知識に従いＣ１，・・・，Ｃｎにソートするステップと、
（ｂ）各ｉ∈｛１，・・・，ｎ｝および各遺伝子ｇ∈Ｃｉに対して、Ｃ１∪．．．∪Ｃｉから親候補を最大ｍ個まで選択するステップと、
（ｃ）定理２を用いて最適遺伝子ネットワーク・モデルを計算するステップと、を含む、請求項１に記載の方法。
【請求項９】
前記アルゴリズムが、
（ａ）全てのｇ∈Ｇに対してＦ（ｇ，φ）＝ｓ（ｇ，φ）を計算するステップと、
（ｂ）全てのＡ⊆Ｇ、Ａ≠φ、および全てのｇ∈Ｇに対してＦ（ｇ，Ａ）を、
【数２】
2007519988.tif 000030
として計算するステップと、
（ｃｃ）Ｍ（φ）＝φに設定するステップと、
（ｄ）全てのＡ⊆Ｇ、Ａ≠φに対して、
（ｉ）ｇ＊＝ａｒｇ  ｍｉｎｇ⊆Ａ（Ｆ（ｇ，Ａ−｛ｇ｝）＋ＱＡ−｛ｇ｝（Ｍ（Ａ−｛ｇ｝）））を計算すると共に、
（ｉｉ）全てのｌ≦ｉ＜｜Ａ｜に対して、Ｍ（Ａ）（ｉ）＝Ｍ（Ａ−｛ｇ＊｝）（ｉ）且つＭ（Ａ）（｜Ａ｜）＝ｇ＊を設定するステップと、
（ｅ）ＱＧ（Ｍ（Ｇ））を返すステップと、を含む、請求項１に記載の方法。
【請求項１０】
前記アルゴリズムが、
（ａ）全てのｇ∈Ｇに対してＦｍ（ｇ，φ，１）＝φ、Ｓｍ（ｇ，φ，１）＝ｓ（ｇ，φ）を設定するステップと、
（ｂ）全てのｇ∈Ｇ、全てのＡ⊆Ｇ、Ａ≠Ｎ、および全てのｎ≦ｍに対して、
（ｉ）｛Ｂ⊆Ａ｜Ｂ＝Ａ∨Ｂ＝Ｆｍ（ｇ，Ａ−｛ｈ｝，ｐ）、ｈ∈Ａ，ｐ≦ｍ｝−｛Ｆｍ（ｇ，Ａ，ｐ）｜ｐ＜ｎ｝｝から、ｓ（ｇ，Ｂ＊）が最小化されるようにＢ＊⊆Ａ選択すると共に、
（ｉｉ）Ｆｍ（ｇ，Ａ，ｎ）＝Ｂ＊，Ｓｍ（ｇ，Ａ，ｎ）＝ｓ（ｇ，Ｂ＊）を設定するステップと、
（ｃ）Ｍｍ（φ，１）＝φおよびＤｍ（φ，１）＝φに設定するステップと、
（ｄ）全てのＡ⊆Ｇ、φ、および全てのｎ≦ｍに対して、
（ｉ）三つ組（ｇ，ｐ，ｑ）∈Ａ×ＩＮ≦ｍ×ＩＮ≦ｍを、ｓｃｏｒｅ（ＱＡ−｛ｇ｝（Ｍｍ（Ａ−｛ｇ｝，ｐ），Ｄｍ（Ａ−｛ｇ｝，ｐ）））＋Ｓｍ（ｇ，Ａ−｛ｇ｝，ｑ）が最小化されて、（ｇ，ｐ，ｑ）が、ｒ＜ｎに対してＱＡ（Ｍｍ（Ａ，ｒ），Ｄｍ（Ａ，ｒ））とは異なるネットワークを誘導するように選択し、
（ｉｉ）ｉ＜｜Ａ｜に対してＭｍ（Ａ，ｎ）（ｉ）＝Ｍｍ（Ａ−｛ｇ｝，ｐ）（ｉ）、およびＭｍ（Ａ，ｎ）（｜Ａ｜）＝ｇを設定すると共に、
（ｉｉｉ）ｖをＤｍ（Ａ−｛ｇ｝，ｐ）を表わすものとし、全てのＩ＜｜Ａ｜およびＷ｜Ａ｜＝ｑに対しｗ∈ＩＮ｜Ａ｜をｗｉ＝ｖｉとして設定し、Ｄｍ（Ａ，ｎ）＝ｗを設定するステップと、
（ｅ）全てのｉ≦ｍに対してＱＧ（Ｍｍ（Ｇ，ｉ），Ｄｍ（Ｇ，ｉ））を返すステップと、を含む、請求項１に記載の方法。
【請求項１１】
列挙された遺伝子ネットワークの信頼性が、
（ａ）最も可能性のある遺伝子ネットワーク・モデルＮｉ，ｌ≦ｉ≦ｎを列挙するステップと、
（ｂ）全てのｇ，ｈ∈０Ｇに対して、ネットワークＮｉにおける辺（ｇ，ｈ）の出現を数えるステップと、
（ｃ）少なくともｃ回出現した全ての辺（ｇ，ｈ）を選択するステップと、
（ｄ）｜Ｍ｜＝ｋである選択された辺の集合の全ての部分集合Ｍに対して、ＭＩの全ての辺を含んでいるネットワークを数えるステップと、
（ｅ）少なくともｃ回出現した全てのモチーフＭを返すステップと、を含む、請求項１〜１０のいずれか一項に記載の方法。
【請求項１２】
ＢＲＮＣスコア、ＢＤｅスコア、およびＭＤＬスコアからなるグループから選択されたスコアリング関数を計算するステップを更に含む、請求項１〜１１のいずれか一つに記載の方法。
【請求項１３】
実質的に本明細書に記載されているように遺伝子ネットワークを決定する方法。
【請求項１４】
請求項１〜１１のいずれか一項に記載の方法を用いて得られた結果を含む記憶媒体。
【請求項１５】
実質的に本明細書に記載されている方法を用いて得られた結果を含む記憶媒体。
【請求項１６】
遺伝子ネットワーク関係を決定するシステムであって、
有機体の遺伝子群の定量的発現データを提供する入力装置と、
前記有機体の遺伝子群の定量的発現データを受信すべく適合させた記憶装置と、
前記遺伝子群内のネットワーク関係のベイジアンネットワーク解析を実行し、それにより前記ネットワーク関係を反映するデータ集合を生成すべく適合させたプロセッサと、
前記ネットワーク関係の前記データ集合を表示する出力装置と、を含むシステム。
【請求項１７】
実質的に本明細書に記載されているように遺伝子ネットワーク関係を決定するシステム。
【発明の詳細な説明】
【技術分野】
【０００１】
優先権主張
本出願は、米国仮出願第６０／５２９，１４６号「推定方法および生物学的制約を用いた遺伝子ネットワークの推定」（Ｓａｓｃｈａ  Ｏｔｔ、Ｓａｔｏｒｕ  Ｍｉｙａｎｏ、２００３年１２月１２日出願、代理人整理番号ＧＥＮＮ−０１０１１ＵＳ０）に基づく優先権を主張し、その全文を本明細書に引用する。
【０００２】
本発明の分野
本発明は、生物学的制約を利用する数学的推定方法を用いて生物系における遺伝子ネットワークを推定することに関するものである。特に、本発明は、検索分野が遺伝子システムに対する生物学的制約の解釈により限定される、ベイズ推定を用いる遺伝子ネットワーク推定に関するものである。
【背景技術】
【０００３】
背景
遺伝子発現測定から遺伝子ネットワークを推定することは、システムズ・バイオロジーにおける大きな課題である。遺伝子ネットワークが正しく推定できたならば、細胞プロセスのより良い理解につながり、従って、新薬発見、疾病研究その他の領域に応用できる。発現レベルの測定からの遺伝子ネットワークの推定は焦点となっている研究の一つである。
【０００４】
遺伝子ネットワークのモデル化に頻繁に用いられる方式にベイズ・ネットワークがある。ベイズ・ネットワークにおいて、ネットワークの挙動は、一つの有機体から得られた関連する遺伝子群の全ての遺伝子の発現レベルについての同時確率分布により記述される。同時確率分布は、有向非巡回グラフを用いて条件付き確率分布に分解することができ、これをネットワークと呼ぶことができる。所与の遺伝子発現測定データに対して、ネットワークを選択するために多くのスコア関数が提案されてきた。従前の知識を利用して新規なスコア関数の研究が進行中である。ネットワークは、データが与えられると、ネットワークの尤度に基づくスコア関数を用いてスコア付けされる。スコア関数を選択したならば、遺伝子ネットワークを推定する作業は、最小スコアを有するネットワークを発見することである。
【発明の開示】
【発明が解決しようとする課題】
【０００５】
しかし、ＮＰ困難である遺伝子ネットワーク推定問題および超指数関数的大きさの探索空間に直面すると、遺伝子ネットワークを推定するためにヒューリスティック（発見的）アルゴリズム、例えばグリーディ（貪欲）アルゴリズムまたはシミュレーテッド・アニーリングを適用する者もあった。ヒューリスティック・アルゴリズムは精度を一切保証しないため、使用した計算方法によりネットワーク推定にどの程度偏りが生じるかは不明なままである。
【０００６】
約３０個前後の遺伝子を含む遺伝子ネットワークに対する何らかのデータに関して最適遺伝子ネットワーク・モデルを推定することが可能である。この結果は、遺伝子ｇおよびｇの親集合にスコアを割り当てる任意のスコア関数ｓ：Ｇ×２Ｇに対して成立する。しかし、最適遺伝子ネットワーク・モデルは最尤モデルであるものの、あるデータに関してほぼ同じ尤度を有しながら大幅に異なるモデルが依然として存在し得る。これは特に遺伝子発現測定が一般に当該遺伝子ネットワークに関する部分的な情報しか提供しないためである。更に、遺伝子が１０個しかない遺伝子ネットワークの場合でさえ、約４．１７×１０１８通りの可能なネットワーク・モデルが存在する。従って、たとえ最適遺伝子ネットワーク・モデルであっても一般に対象遺伝子ネットワークには合致しない。
【課題を解決するための手段】
【０００７】
概要
本発明の特定の態様において、制約的な意味で推定されたネットワークの最適性を提供しつつ、遺伝子ネットワークを効率的に推定する方法を提示する。最初に、可能な遺伝子ネットワークの推定モデルを生成する。第二に、探索空間の部分空間のうち生物学的に関連する部分空間を選択する。次いで、小さい遺伝子ネットワークを最適に計算するアルゴリズムを繰り返し適用することにより、選択された部分空間内の最適解を発見する。部分空間の選択は、生物学的知識を適用するか、または、事前知識が利用できない場合、例えばクラスタリング法を用いて行なうことができる。
【０００８】
特定の実施形態において、このような方法により、生物学的には些細だが計算的に効果的な仕方で探索空間を制約することにより生物学的知識を利用することができる。本方法は、制約された探索空間に関する最適性を向上させ、従って、推定のうちヒューリスティック／経験的に行なわれる部分と、最適に行なわれる部分を明確に区別する。
【０００９】
別の実施形態では、ベイジアンネットワークを用いて、遺伝子ネットワークの挙動を全ての遺伝子に対する同時確率分布としてモデル化される。これにより、遺伝子相互作用を極めて一般的にモデリング可能にする。同時確率分布は、他の遺伝子による遺伝子の調整を表わす条件付き確率の積として分解することができる。この分解は、有向非巡回グラフとして表わすことができる。ベイジアンネットワーク・モデルを用いて、生物学的に有り得る遺伝子ネットワークを発見することができる。
【００１０】
遺伝子発現情報は、例えばマイクロアレイ解析を含む、当分野で公知の方法を用いて生成することができる。これらの方法は当分野で公知であるため、ここで詳細に述べる必要はない。しかし、ネットワーク関係の信頼性が遺伝子発現に関して得られた情報の精度および／または信頼性に依存することを理解されたい。従って、いくつかの実施形態において、遺伝子発現データの信頼性を高めるために、反復的な解析および管理された研究が望ましい。
【００１１】
別の態様において、最適および準最適なネットワークを尤度順に列挙する理論的基礎およびアルゴリズムを提供する。更に別の実施形態において、最適遺伝子ネットワーク・モデルを、既存の遺伝子ネットワークに関する利用可能な知識と比較した。更に別の実施形態において、関連する種同士の遺伝子ネットワーク・モデルを比較するために、最適遺伝子ネットワーク・モデルの信頼性が高い部分を抽出する方法を提示し、この方法を枯草菌および大腸菌の利用可能なデータに適用した。
【００１２】
上記の結果、本発明の方法により識別された部分ネットワークが最適ネットワーク・モデル自体よりもはるかに良く知見と一致することが示された。本発明では上位ｎ個の最適遺伝子ネットワーク・モデルに基礎をおくため、ヒューリスティックな方法だけに依存する方法よりも、これらの推定の方がより信頼性の高い結論を導くことができる。
【００１３】
本発明はまた、多くの異なる遺伝子からの遺伝子発現に関するデータを受信すべく適合させた入力装置と、遺伝子ネットワーク関係を推定するプログラムを実行すべく適合させたプロセッサとを少なくとも含む、遺伝子ネットワークの推定解析を実行するシステム含む。追加的な実施形態において、システムは、遺伝子ネットワーク解析結果を遠隔拠点へ表示、印刷、またを送信すべく適合された出力装置を含む。
【００１４】
本発明について、本発明の実施形態の個別の記述および図面を参照しつつ述べる。
【発明を実施するための最良の形態】
【００１５】
詳細な説明
細胞内の遺伝子相互作用の解明は、近年の研究における重要な焦点である（２０）。特定の実施形態において、ベイジアンネットワークを用いて遺伝子ネットワークをモデル化することができる（４，５，６，１１，１２，１４，１６，１９）。ベイジアンネットワークにおいて、ネットワークの挙動は、全ての遺伝子の発現レベルに対する同時確率分布により記述される。同時確率分布は、有向非巡回グラフすなわち「ネットワーク」を用いて条件付き確率分布に分解する必要がある。遺伝子発現測定データが与えられている場合にネットワークを選択するための多くのスコア関数が提案されてきた（３，４，１１）。事前知識を利用する新規スコア関数に関する研究途上にある（１３，２３）。スコア関数を選択したならば、遺伝子ネットワークを推定する作業は、スコアが最小であるネットワークを発見することである。
【００１６】
しかし、ＮＰ困難遺伝子ネットワーク推定問題（１）および超指数大きさ検索空間（１７）に直面して、研究者は遺伝子ネットワークを推定するためにグリーディ・アルゴリズム（１１）等のヒューリスティック・アルゴリズムまたはシミュレーテッド・アニーリング（７）を適用してきた。ヒューリスティック・アルゴリズムは、自身の精度を一切保証しないため、計算方法によりネットワーク推定にどの程度偏りが生じるかは不明なままである。
【００１７】
本発明の特定の態様において、制約的な意味で推定されたネットワークの最適性を提供しつつ、遺伝子ネットワークを効率的に推定する方法を提示する。最初に、可能な遺伝子ネットワークの推定モデルを生成する。第二に、探索空間の生物学的に関連する部分空間を選択する。次いで、小さい遺伝子ネットワークを最適に計算するアルゴリズムを繰り返し適用することにより選択された部分空間内の最適解を発見する（１５）。部分空間の選択は、生物学的知識を適用するか、または、事前知識が利用できない場合、例えばクラスタリング法を用いて行なうことができる。
【００１８】
特定の実施形態において、このような方法により、生物学的には些細だが計算的に効果的な仕方で探索空間を制約することにより生物学的知識を利用することができる。本方法は、制約された探索空間に関する最適性を向上させ、従って、推定のうちヒューリスティック／経験的に行なわれる部分と、最適に行なわれる部分を明確に区別する。
【実施例】
【００１９】
実施例
以下の例は、本発明の態様を例示することを目的としており、限定を意図するものではない。更に、各々の例において、参考文献を括弧で引用し、各々の例の末尾に参考文献の一覧を掲載する。本発明の他の実施形態は、当業者が無用な実験を行なうことなく製作することができる。これらの実施形態は全て本発明の範囲内に含まれる。
【００２０】
例１：生物学的制約を用いた最適遺伝子ネットワークの発見
選択された部分空間内で最適解を発見するために必要な実行時間はＯ（ｎ）であり、ここで、ｎは遺伝子の個数である。従って、いくらでも大きい遺伝子ネットワークの推定が可能である。本発明者らは、この方法を、対応関係ｓ：Ｇ×２Ｇ→ＩＲを有する任意のスコア関数ｓを用いて適用することができる点に注意したが、これはベイジアンネットワーク・フレームワーク内の全てのスコア関数に当てはまり、且つ大部分の他の関数にも当てはまる。
【００２１】
発明者は上述の方法を実行して、イースト菌データおよび枯草菌データに適用した。枯草菌の推定されたネットワークを生物学的知識と比較すると顕著な合致が得られた。
【００２２】
必要な表記法は第二節で解説され、本方法は第三節で詳細に説明されて、ネットワーク推定の結果について第四節で議論する。
【００２３】
１．１．準備
実施例を通じて、遺伝子ネットワークを予測する遺伝子の集合をＧで表記し、遺伝子ｇ∈Ｇおよび親遺伝子の集合Ａ⊆Ｇにスコアを割り当てるスコア関数ｓ：Ｇ×２Ｇ→ＩＲが与えられていると仮定する。所与のネットワークＮに対し、Ｎのスコアをｓｃｏｒｅ（Ｎ）＝ｄｅｆΣｇ∈Ｇｓ（ｇ，ＰＮ（ｇ））として定義し、ここに、ＰＮ（ｇ）はＮにおけるｇの親の集合を表わす。（１５）からいくつかの表記を導入する。
【００２４】
定義１．１：
全てのｇ∈ＧおよびＡ⊆ＧについてＦ：Ｇ×２Ｇ→ＩＲを、Ｆ（ｇ，Ａ）＝ｄｅｆｍｉｎＢ⊆Ａｓ（ｇ，Ｂ）と定義する。選択がＡに制約されている場合、Ｆ（ｇ，Ａ）がｇの親の最適な選択である。πを用いて集合Ａ⊆Ｇの順列π｛１，・・・，｜Ａ｜｝→Ａを表わし、ΠＡでＡのあらゆる順列の集合を表わす。
【００２５】
定義１．２：
Ａ⊆Ｇとする。全てのπ∈ΠＡについてＱＡΠＡ→ＩＲを次式で表わす。
【数１】
2007519988.tif 000003
【００２６】
定義１．３：
全てのＡ⊆Ｇについて、Ｍ：２Ｇ→Ａ⊆ＧΠＡを次式で定義する。
【数２】
2007519988.tif 000004
【００２７】
文献（１５）において、ＱＡ（π）はＭ（Ａ）がトポロジカルソートである全てのネットワーク空間における最適ネットワークであり、ｌはＡにおける少なくとも１個の最適ネットワークのトポロジカルソートである順列であることが明らかにされた。
【００２８】
文献（１５）から以下のアルゴリズムを掲示し、本研究でサブルーチンとして使用する。
【００２９】
アルゴリズム１．１
ステップ１：全てのｇ∈Ｇに対して、Ｆ（ｇ，φ）＝ｓ（ｇ，Ｎ）を計算する。
ステップ２：全てのＡ⊆Ｇ、Ａ≠φおよびｇ∈Ｇに対して、Ｆ（ｇ，Ａ）を
【数３】
2007519988.tif 000005
として計算する。
ステップ３：Ｍ（φ）＝φに設定する。
ステップ４：全てのＡ⊆Ｇ、Ａ≠φに対して、以下の２ステップを実行する。
ステップ４ａ：ｇ＊＝ａｒｇ  ｍｉｎｇ⊆Ａ（Ｆ（ｇ，Ａ−｛ｇ｝）＋ＱＡ−｛ｇ｝（Ｍ（Ａ−｛ｇ｝）））を計算する。
ステップ４ｂ：全てのｌ≦ｉ＜｜Ａ｜に対して、Ｍ（Ａ）（ｉ）＝Ｍ（Ａ−｛ｇ＊｝（ｉ）且つＭ（Ａ）（｜Ａ｜）＝ｇ＊と設定する。
ステップ５：ＱＧ（Ｍ（Ｇ））を返す。
【００３０】
アルゴリズム１．１はステップ１およびステップ２で関数Ｆを計算し、ステップ３およびステップ４でＦを利用して関数Ｍを計算する。ステップ５において、Ｍを用いて最適ネットワークのスコアを計算した。このアルゴリズムは、最高約３０個の遺伝子の遺伝子ネットワークの最適遺伝子ネットワーク・モデルを発見すべく適用されたが、大規模な遺伝子ネットワークの場合には信頼性が下がる。
【００３１】
定理１．１（１５）  アルゴリズム１．１はＯ（ｎ・２ｎ）の動的計画法のステップを用いて最適ネットワークを発見する。
【００３２】
１．２  方法
最初に定理を証明し、次いで大規模遺伝子ネットワークを推定するための２種のアルゴリズムをそこから導く。
【００３３】
１．２．１理論的基礎
以下の導出定理は、本実施例の方法の根拠を提供する。
【００３４】
定理１．２
ｍ、ｃ∈ＩＮとする。各ｇ∈Ｇに対して、Ｃｇ⊆Ｇ−｛ｇ｝を親候補の集合とする。以下の２条件が成立すると仮定する。
１．全てのｇ∈Ｇに対して、｜Ｃｇ｜≦ｍ
２．Ｃｇ、ｇ∈Ｇにより導出されるグラフの強連結成分はｃを超えない。
次いで、選択された親候補に関する最適ネットワークが時間Ｏ（｜Ｇ｜）で発見できる。
【００３５】
証明
Ｌ＝（Ｇ，Ｅ）をＣｇ、ｇ∈Ｇにより導出された遺伝子Ｇの集合上のグラフを表わすものとする。すなわちＥは、ｈ∈Ｇｇが成立する辺（ｈ，ｇ）だけを含む。Ｓ１，・・・，Ｓｎ⊆ＧがＬの強連結成分を表わすとする。ここで、Ｌ’＝（｛Ｓ１，・・・，Ｓｎ｝、Ｅ’）を、次式で定義されるＥとの強連結成分上のグラフを表わすとする。
【数４】
2007519988.tif 000006
【００３６】
従って、強連結成分の定義により、Ｌ’には巡回路が存在しない。一般性を失うことなく、Ｓ１，・・・，ＳｎがＬ’のトポロジカルソートである、すなわちｊ＜ｉの場合Ｌ’には辺（Ｓｉ，Ｓｊ）が存在しないことを意味する。
【００３７】
アルゴリズム１．１は、若干の変更を加えることにより、選択された親候補に関してＧに対する最適ネットワークを計算することができる。アルゴリズム１．１を強連結成分Ｓｉに適用すべく、アルゴリズムを以下のように変更する。
１．ステップ１およびステップ２におけるＦの計算で、全てのｇ∈Ｓｉおよび全てのＡ⊆Ｃｇに対してＦ（ｇ，Ａ）だけを計算する。
２．ステップ４ａにおける項Ｆ（ｇ，Ａ−｛ｇ｝）（アルゴリズム１の定義参照）をＦ（ｇ，（Ｃｇ＝Ｓｉ）∪（Ｃｇ∩Ａ））で置き換える。
【００３８】
これら二つの変更により、アルゴリズムに親候補の制約を導入し、ｇ∈ＳｉがＳｉの外に親を持つことを許す。アルゴリズム１と同様に定理１．１（１５）の証明を変更されたアルゴリズム１．１にも行なうことができるため、これはアルゴリズム１．１の正しさに影響を及ぼさない。このように変更したアルゴリズム１．１を各々の強連結成分Ｓｉに任意の順序で適用して、各ｉ∈１，・・・，ｎに対して部分ネットワークＮｉ＝（Ｇ，Ｅｉ）を生成する。次いで、ネットワーク
【数５】
2007519988.tif 000007
を返す。定理１．１により、各部分ネットワークＮｉは最適である。Ｌ’の非巡回性から、Ｎも非巡回である。従って、定義ｓｃｏｒｅ（Ｎ）＝Σｇ∈Ｇｓ（ｇ，ＰＮ（ｇ））を用いて、Ｎは最適ネットワークである。
【００３９】
｜Ｓｉ｜≦ｃ且つ｜Ｃｇ｜≦ｍであるため、アルゴリズム１を一回呼び出す計算時間はある定数で制約される。従って、全ての強連結成分へのアルゴリズム１．１の適用はＯ（｜Ｇ｜）で実行可能であり、証明を完了する。
【００４０】
アルゴリズム１．１および証明で定義されたサブルーチンとしてのアルゴリズム１．１を、全ての最適ネットワークおよび準最適ネットワークを各自のランクの順序で列挙可能なように、適用できる点が注目される。計算時間は、計算対象であるネットワークの数に応じて増大するものの、依然として実行可能である。これは、スコアが若干変化しても推定されたネットワークの安定性を評価する上で貴重である。
【００４１】
１．２．２事前知識を用いない大規模遺伝子ネットワーク予測
定理１．２により、効果的なネットワーク予測を可能にすべく、あらゆる候補相互作用を含むグラフにおける強連結成分を定数で制約できるように、各々の遺伝子の親候補を選択することができる。遺伝子ネットワークは、相互には疎結合された（１０、１８）密結合ブロックで構成されていると仮定されるため、多くの研究設定においてこの制約は満たされると思える。従って、定理１．２は遺伝子ネットワーク推定の一般的手法の基礎を与える。
【００４２】
本研究において、定理１．２を利用すべく設計されたアルゴリズムの二つの例を与える。第一の例（アルゴリズム１．２）は、事前知識を前提とせず、クラスタリングにより密結合ブロックを検知する。第二のアルゴリズム（アルゴリズム１．３）は生物学的知識を応用して定理１．２の制約を実現する。
【００４３】
アルゴリズム１．２
ステップ１：どのクラスタもｃ遺伝子より大きくならいようにＧ内の遺伝子をクラスタリングする。
ステップ２：大きさの降順にクラスタをソートしてＣ１，・・・，Ｃｎを得る。
ステップ３：各ｉ∈｛１，・・・，ｎ｝および各ｇ∈Ｃｉに対して、Ｃ１∪．．．∪Ｃｎから最大ｍ個の親候補を選択する。
ステップ４：定理１．２を用いて最適遺伝子ネットワーク・モデルを計算する。
【００４４】
全ての遺伝子ｇに対する親候補は、ｇが含まれるクラスタまたは前のクラスタから選択されるため、Ｃｇ、ｇ∈Ｇにより生成されるグラフ内のどのサイクルも２個のクラスタにまたがることはできない。従って、全てのクラスタＣｉに対して｜Ｃｉ｜≦ｃである故、強連結成分はｃより大きくない。これはアルゴリズム１．２の正しさを示す。
【００４５】
同じクラスタに属する遺伝子は相関しており、従って、遺伝子ネットワークの密結合された部分を形成している筈であり、一方、異なるクラスタに属する遺伝子の相関はより低く、従って、遺伝子ネットワーク内で結合されている可能性は小さい。ピアソン相関を距離測度とするｋ-−平均クラスタリングを用いる。
【００４６】
大きさの降順にクラスタをソートする根拠は以下の通りである。ソーティングの開始時点におけるクラスタは、その前にクラスタが殆ど無いため、親候補の選択に最も強い制約を有する。親候補を選択できる遺伝子の数を最大にすべく、開始時点で大規模なクラスタを用意すべきである。
【００４７】
ステップ１およびステップ２により課される部分空間の制約は、遺伝子のいくつかのペアに対して辺の方向を制約するだけであり、依然として任意の２個の遺伝子は結合されてもよい点に注目する。
【００４８】
ステップ３において遺伝子ｇの親候補ｈの選択する多くの合理的な基準がある。そのような例として遺伝子の相関またはｓ（ｇ，｛ｈ｝）がある。本実施例において、後者の基準を利用した。
【００４９】
１．２．３  事前知識を用いた大規模遺伝子ネットワーク予測
遺伝子ネットワークに関して、密結合された要素および異なる要素における遺伝子の相互作用の方向に関する事前知識がある場合、以下のアルゴリズムを適用することができる。
【００５０】
アルゴリズム１．３
ステップ１：｜Ｃｉ｜≦ｃであるグループＣｉ内でＧの遺伝子をグループ分けして、生物学的知識に従いそれらをソートしてＣ１，・・・，Ｃｎを得る。
ステップ２：各ｉ∈｛１，・・・，ｎ｝および各ｇ∈Ｃｉに対して、Ｃ１∪．．．∪Ｃｉから最大ｍ個の親候補を選択する。
ステップ３：定理２を用いて最適遺伝子ネットワーク・モデルを計算する。
【００５１】
アルゴリズム１．２に関して、アルゴリズム１．３の正しさは、ステップ２において親候補を選択する仕方に直接追随する。密結合要素だけでなく、その順序が既知である例として、細胞周期の特定のフェーズにおいて活性化される遺伝子がある。密結合要素および／またはその順序に関する知識が部分的な知識である状況において、アルゴリズム１．２およびアルゴリズム１．３を組み合わせて用いることができる。
【００５２】
１．３  結果
アルゴリズム１．１（１５）の既存の実装および公的に利用可能なクラスタリング・ソフトウェア（８）を用いてアルゴリズム１．２およびアルゴリズム１．３を実装した。非線形遺伝子相互作用をモデル化して遺伝子発現データを離散化することなく扱うことができるため、スコア関数としてＢＮＲＣスコア（１１）を選択した。
【００５３】
１．３．１  枯草菌データへのアルゴリズム１．２の適用
アルゴリズム１．２を枯草菌データに適用した（９）。６個のシグマ因子に関するデータ、および選択されたシグマ因子（２１）により調節されることが知られている７９個のオペロンを選択した。オペロンの発現比率を、各自の遺伝子の発現比率の平均として定義した。従って、アルゴリズム１．２に予測能力があるならば、各オペロンに８５個の遺伝子を含むこのネットワークのうち７９個の既知の調節相互作用があるところ、これらの有意な部分を識別するであろう。選択された親候補の最大数であるパラメータｍを１５に、またクラスタの最大数であるパラメータｃを２５に設定した。ステップ１で計算されたクラスタの実際の大きさは８から２４まで変動した。この計算のために、各々９００ＭＨｚである９６個のＣＰＵを有するＳｕｎ  Ｆｉｒｅ  １５Ｋスーパーコンピュータを約１時間使用した。
【００５４】
推定されたネットワークの生物学的有意性を推定すべく、アルゴリズム１．２が予測したネットワークにおける各シグマ因子への各オペロンの距離を計算した。その際に、距離を、オペロンとシグマ因子を結合する無向経路上の辺の最小数として定義した。正しいシグマ因子が他の全てのシグマ因子よりもオペロンに近い場合、調整関係が検知されたものと推定した。これにより、無意味な予測手段の場合、調整相互作用を検知する確率は高々１／６であるため、この推定結果に対してｐ値を計算することができる。タイブレークを検知されていないとカウントするため、実際の確率は＜＜式：ＧＩＷ２０１４６＞＞より低い。
【００５５】
【表１】
2007519988.tif 000008
【００５６】
表１．１はこの評価スキームの結果を示す。推定された遺伝子ネットワークが２個のシグマ因子、２０１４８および２０１４９について生物学的知識と極めてよく合致することがわかった。このことは、推定結果が有益な生物学的情報を含むことを示す。検知されなかった調節関係の大多数はタイブレークであり、従って、オペロンから自身の正しいシグマ因子までの距離は依然として最小であったが、１個以上の他のシグマ因子も同程度近い点に注目した。２０１５０の有意性が特に強い点は、微分方程式ネットワークの場合[ＣＩＴＥ：ｄｅｈｏｏｎ０３ｂ]の結果と整合していることがわかった。従って、当該データ集合が、他の領域よりも２０１５１周辺にある遺伝子ネットワーク領域に関してより多くの情報をもたらす可能性がある。
【００５７】
発明者は、遺伝子ネットワーク推定の質を原理的に評価する問題に関してほとんど研究がされていない点に注目した。発明者の知る限り、推定の有意性を証明すべく明確なｐ値を与える公開文献が他に[ＣＩＴＥ：ｄｅｈｏｏｎ０３ｂ]一つしか存在しない。遺伝子ネットワーク推定の原理的な評価を行なう際に、いくつかの問題に遭遇する。そのような例として、遺伝子ネットワークに関する知識の不足、および遺伝子ネットワークのどの部分がどの実験条件の下で機能しているかが不確実であることが挙げられる。更に、真の遺伝子ネットワークと、推定されたネットワークの構造上の差異が発見された場合、かなり間違った差異と遷移辺のような受容できる差異とを区別しなければならない。推定されたネットワークにおける最短経路を解析する上述の方法はこれらの問題のいくつかに取り組む一つの道ではあるが、これらの問題の研究は更に追求しなければならない。
【００５８】
１．３．２  細胞周期データへのアルゴリズム１．３の適用
アルゴリズム１．３は、生物学的に有意味な部分空間を発見すべく事前知識の適用に基づいている。ここで、アルゴリズム１．３を、同期細胞培養（２２）を用いて経時的実験から得られた出芽酵母のＲＮＡマイクロアレイデータに適用する。このデータ集合は、イースト菌の細胞周期の間に作用する遺伝子ネットワークに関する何らかの情報を含んでいるものと期待できる。
【００５９】
３個の遺伝子ペアｃｌｎ１／ｃｌｎ２、ｃｌｂ５／ｃｌｂ６、およびｃｌｂ１／ｃｌｂ２は各々特に細胞周期Ｇ１／Ｓフェーズ、Ｓフェーズ、Ｍフェーズの間に活性化されることが知られている（２）。この事実に基づくクラスタリング解析により、いくつかの遺伝子がこれらのフェーズの一つの間に優先的に活性化されるものと予測された（１３）。４３個の遺伝子の集合を選択して、各々が細胞周期の一部に対応する３個のグループに分けた。これらのグループを、細胞周期フェーズの時間順にソートして（アルゴリズム１．３のステップ１）、親候補の最大数（パラメータ）を２０に設定し、遺伝子ネットワークを推定すべくアルゴリズム１．３を適用した。計算は、１．９ＧＨｚである単一のＣＰＵを１時間足らず用いて行なった。図１に推定結果を示す。
【００６０】
総計８７個の予測された遺伝子相互作用が観察された。表１．２に、遺伝子が属するグループ（行）および（列）により分割された、第一の遺伝子から第二の遺伝子への有向相互作用の個数を示す。８７個の相互作用のうち５２個がグループ内であって、細胞周期の同一フェーズの間に、優先的に活性化された遺伝子間で大部分の相互作用が起こるという予想をよく反映している。３５個の相互作用が、異なるグループから遺伝子ペアについて推定される。２個のグループの間の相互作用の個数は、時間的に隣接するＧ１／ＳとＳ／Ｇ２、およびＳ／Ｇ２とＭの各々についてそれぞれ１３であるが、Ｇ１／ＳフェーズからＭフェーズへの相互作用については９個しかない。従って、遺伝子ネットワーク推定は３個のグループ内の遺伝子の集合の分解に基づいているが、グループ間の相互作用をうまく説明できる。
【００６１】
【表２】
2007519988.tif 000009
【００６２】
各々の遺伝子について相互作用の個数を数えた。これは、２０１５８の親の個数に当該遺伝子の子の個数を加えたものである。推定された遺伝子ネットワーク（図１）の予測された相互作用が最大個数である遺伝子を表１．３に掲示するが、これはまた、サッカロミセス属ゲノムデータベースと注釈を付けられた通りに当該遺伝子の分子機能を与える。少なくとも７個の予測された相互作用を有する全ての遺伝子は、転写レベルで調整活動を有するか、または未知の機能を有することが観察された。６個の予測された相互作用を有する３個の遺伝子のうち、１個の遺伝子（ｈｃｍ１）は転写レベル上の調整機能を有するが、他の２個の遺伝子は細胞周期に依存するシグナリング・タンパク質である。予測された相互作用の個数がより多い全ての遺伝子は既知の活性化された調整機能を有するか、または未知の機能であると結論する。これにより、推定された遺伝子ネットワークは少なくともある程度の生物学的情報を含むことが示された。
【００６３】
重要な遺伝子ｃｌｎ２およびｃｌｂ５が表１．３の遺伝子のグループ内により低いランクを有するとの観察結果は、ｍＲＮＡ測定からいずれが観測可能であるかを間接的な仕方でのみシグナリングするそれらの役割により説明することができる。興味深いことに、ｓｗｉ５は、既知の調整関係[ＣＩＴＥ：ｚｈｕ００]であるｆｋｈ２により調整されることが予測される。従って、他の遺伝子の転写レベルでの働きをより容易に観測できると期待される。
【００６４】
【表３】
2007519988.tif 000010
【００６５】
これに対して表１．４は、予測された相互作用の最小数を有する遺伝子を示す。代謝でまたはシグナリング機能を有する遺伝子を数個、機能が未知の遺伝子を１個、ｍＲＮＡ結合遺伝子（ｓｔｏ１）を１個、および転写因子（ｓｗｉ４）を１個発見した。ｓｔｏ１が核スプライシングに関与することが知られているため、調整役割の可能性は低い。従って、予測された相互作用の最小数を有する１１個の遺伝子のただ１個のみが転写レベルでの既知の調整機能を有し、ｍＲＮＡ測定から多くのシグナリング・イベントが観測可能でないないため、これは生物学的にあり得る。遺伝子ネットワーク推定が遺伝子の機能に関して単なるクラスタリングの結果よりも多くの示唆を与えることができると結論付ける。
【００６６】
【表４】
2007519988.tif 000011
【００６７】
本明細書において、大規模遺伝子ネットワークを効率的に推定する一般的な方法に対する理論的基礎を与える。この方法に従うアルゴリズムは二つの部分から成る。第一の部分は、探索空間の生物学的に有意味な部分空間を識別するヒューリスティックまたは経験的方法である。第二の部分は、選択された部分空間における遺伝子ネットワークの最適推定である。
【００６８】
この方法により、有意味な遺伝子ネットワークを推定して、生物学的知識を適切且つ効果的に適用できる、一方、計算時間により遺伝子の個数が実際に制限されないことを示した。
【００６９】
本方法は遺伝子発現データ特定のスコアまたは特定の種類に依存しないため、一般的に適用可能である。シーケンス情報（２３）または構造情報（１３）等の事前知識を組み込んだ新たなスコアの開発は、予測能力を更に向上させるのに有用な方法である。部分空間選択をヒューリスティックに行なう別の方法には、ヒューリスティックにより推定された遺伝子ネットワーク周辺に部分空間を構築する、あるいはクラスタの遺伝子発現値の平均を求めてアルゴリズム１を適用して当該クラスタの最適遺伝子ネットワークを発見し、次いでこれを用いて（アルゴリズム１．２のステップ２のように単に大きさでソートするのではなく）クラスタの順序を発見することが含まれる。
【００７０】
例１の参考文献
1. Chickering, D. M., Learning Bayesian networks is NP-complete, in D. Fisher and H.-J. Lenz, editors,Learningfrom Data : Artcial Intelligence and Statistics V, Springer-Verlag, 1996. 
2. Cho, R. J. , Campbell, M. J.,Winzeler, E. A. , Steinmetz, L. , Conway, A. , Wodicka, L., Wolfsberg, T. G. , Gabrielian, A. E. , Landsman, D. , Lockhart, D. J. , Davis, R. W. , A genome-wide transcriptional analysis of themitotic cell cycle, Molecular Cell, 2: 65---73, 1998. 
3. Cooper, G. F. , Herskovits, E. , A Bayesian method for the induction of probabilistic networks from data. Machine Learning, 9: 309--347,1992. 
4. Friedman, N.,Goldszmidt, M. , Learning Bayesian networks with local structure, Jordan, M. I. (ed. ), Kluwer Academic Publishers, 421--459, 1998. 
5. Friedman, N. , Linial, M. , Nachman, I.,Pe'er, D. , Using Bayesian networks to analyze expression data, Journal of Computational Biology, 7: 601--620,2000. 
6. Hartemink, A. J., Gifford, D. K. , Jaakkola, T. S. , Young, R. A. , Using graphical models and genomic expression data to statistically validate models of genetic regulatory networks,Pacific Symposium on Biocomputing, World Scientific,422-433, 2001. 
7. Hartemink, A. J., Gifford, D. K. , Jaakkola, T. S. , Young, R. A. , Combining location and expression data for principled discovery of genetic regulatory network models, Pacific Symposium on Biocomputing, World Scientific, 7:437--449, 2002. 
8. de Hoon, M. J. L. , Imoto, S. , Nolan, J. , and Miyano, S. , Open source clustering software,Bioinformatics, 2003, in press. http://bonsai. ims. u-tokyo. ac.jp/ mdehoon/software/cluster/ 
9. de Hoon, M. J. L. , Ott, S. , Imoto, S. , Miyano, S., Validation of noisy dynamical system models of gene regulation inferred from time-course gene expression data at arbitrary time intervals, European Conference on Computational Biology, 2003. 
10. Hartwell, L. H., Hopfield, J. J. , Leibler, S. , Murray, A. W. , From molecular to modular cell biology, Nature, 402: C47--C52, 1999. 
11. Imoto, S. , Goto, T. , Miyano, S. , Estimation of genetic networks and functional structures between genes by using Bayesian networks and nonparametric regression, Pacific Symposium on Biocomputing, World Scientific, 175--186, 2002. 
12. Imoto, S., Kim, S. , Goto, T. , Aburatani, S. , Tashiro, K. , Kuhara, S. , Miyano, S. , Bayesian network and nonparametric heteroscedastic regression for nonlinear modeling of genetic network,Journal of Bioinformatics and Computational Biology, 1: 231--252,2003. 
13. Nariai, N., Kim, S. , Imoto, S. , Miyano, S. , Using protein-protein interactions for refining gene networks estimated from microarray data by Bayesian networks, Pacific Symposium on Biocomputing, World Scientific, in press, 2004. 
14. Ong, I. M., Glasner, J. D. , Page, D. , Modelling regulatory pathways inE. coli from time series expression profiles,Bioinformatics, 18: 241--248, 2002. 
15. Example 2. 
16. Pe'er, D. , Regev, A. , Elidan, G. , Friedman, N. , Inferring subnetworks from perturbedexpression profiles, Bioinformatics, 17: 215--224,2001. 
17. Robinson, R. W. , Counting labeled acyclic digraphs, New Directions in the Theory of Graphs, 239--273, 1973. 
18. http://www. yeastgenome. org/ 
19. Shen-Orr, S. S. , Milo, R. , Mangan, S. , Alon, U. , Network motifs in the transcriptional regulation networkof Escherichia coli, Nature Genetics, 31:64--68, 2002. 
20. Smith, V. A. , Jarvis, E. D., Hartemink, A. J. , Evaluating functional network inference usingsimulationsofcomplexbiological systems, Bioinformatics, 18:216--224, 2002. 
21. van Someren, E. P. , Wessels, L. F. A. , Backer, E. , Reinders, M. J. T., Genetic network modeling,Pharmacogenomics, 3 (4):507--525, 2002. 
22. Sonenshein, A. L. , Hoch, J. A. , Losick, R., Bacillus subtilis and its closest relatives : from genes to cells, ASM Press, Washington, D. C. , 2001. 
23.Spellman, P. , Sherlock, G. , Zhang, M.,lyer, V. , Anders, K. , Eisen, M. , Brown, P., Botstein, D. , Futcher, B. , Comprehensive identification of cell cycle-regulated genes of the yeast Saccharomyces cerevisiae by microarray hybridization, MolecularBiology ofthe Cell, 9: 3273--3297, 1998. 
24. Tamada, Y. , Kim, S. , Bannai, H. , Imoto, S. , Tashiro, K. , Kuhara, S. , Miyano, S., Estimating gene networks from gene expression data by combining Bayesian network model with promoter element detection, Bioinformatics, 19:227--236, 2003. 
25. Zhu, G.,Spellman, P. T., Volpe, T. , Brown, P. O. , Botstein, D. , Davis, T. N. , Futcher, B., Two yeast forkhead genes regulate the cell cycle and pseudohyphal growth, Nature, 406: 90--94, 2000.
【００７１】
例２：小規模遺伝子ネットワークの最適モデルを発見する
序論
遺伝子発現の測定から遺伝子ネットワークを推定することは、システムズ・バイオロジーにおける大きな課題である。遺伝子ネットワークが正しく推定できたならば、細胞プロセスのより良い理解につながり、従って、新薬発見、疾病研究その他の領域に応用できる。
【００７２】
ベイジアンネットワークは、遺伝子ネットワークのモデル化に広く利用されている方法である（文献３、４、７、９、１１、１３、１４、１７参照）。ベイジアンネットワークにおいて、遺伝子ネットワークの挙動は、全ての遺伝子に関する同時確率分布としてモデル化される。これにより、遺伝子相互作用の極めて一般的なモデリングが可能になる。同時確率分布は、いくつかの遺伝子ｇ１，・・・，ｇｎにより遺伝子ｇの調整を表わす条件付き確率Ｐ（Ｘｇ｜Ｘ１，・・・，Ｘｎ）の積に分解することができる。この分解は、有向非巡回グラフとして表わすことができる。ベイジアンネットワーク・モデルは、生物学的に可能な遺伝子ネットワークの発見が可能になることを示した（文献４、９参照）。
【００７３】
しかし、ベイジアンネットワークを調べることの難しさはその巨大な探索空間にある。ｎ個の遺伝子の遺伝子ネットワークの検索空間とは、ｎ個の頂点を有する有向非巡回グラフの空間である。ｎ個の頂点（ｃｎ）を有する有向非巡回グラフの個数の漸近式並びに再帰式は、ロビンソン（Ｒｏｂｉｎｓｏｎ）により導かれた（文献１５参照）。以下に漸近式を示す。
【数６】
2007519988.tif 000012
【００７４】
例えば、２０個の遺伝子には大体２．３４×１０７２個の可能なネットワーク化があり、３０個の遺伝子を含む遺伝子ネットワークの場合には約２．７１×１０１５８個の可能な解が存在する。遺伝子が９個（検索空間の大きさは大体１．２１×１０１５）しかない遺伝子ネットワークでさえ、力任せの方法ではスーパーコンピュータを使っても計算時間が何年もかかる。更に、最適ネットワークを発見する問題は、離散的スコアＢＤｅ（文献２、３参照）およびＭＤＬ（文献３参照）についてさえＮＰ困難（文献１参照）であることが知られている。従って、これまで研究者はベイジアンネットワーク（文献１８参照）を推定するためにシミュレーテッド・アニーリング（文献８参照）またはグリーディ・アルゴリズム（文献９参照）等のヒューリスティックな方法を使用した。
【００７５】
しかし、ヒューリスティックの精度が不確実なため、ヒューリスティック的に推定されたネットワークに基づいて結論を出すことは困難である。この問題を解決すべく、超指数探索空間の構造を解析して、指数時間の超指数探索空間内での最適解を発見するアルゴリズムを開発した。この方法は、使用する具体的な確率分布に応じて、２０個以上の遺伝子を含む遺伝子ネットワークについて実行可能である。更に、生物学的に正しいとされた仮定を加えることにより、最高４０個の遺伝子を含む遺伝子ネットワークについて最適ネットワークを推定することができる。
【００７６】
ヒューリスティックの不確実性を解決することで、生物学的に正確な遺伝子ネットワークを推定する能力に関して統計モデルを比較する可能性が開ける。また、この方法は、遺伝子の既知の官能基の遺伝子ネットワークを改良するための有益なツールである。
【００７７】
諸方法を第２．２節に示す。第２．３節において、本方法を適用した結果を提示し、これが生物学的に正確な遺伝子ネットワークを推定できることを示す。
【００７８】
２．２  諸方法
２．２．１
本節を通じて、遺伝子の集合をＧ、およびいくつかのグループで用いられるネットワークスコア関数、すなわち遺伝子
【数７】
2007519988.tif 000013
および親遺伝子の集合Ａ⊆Ｇにスコアを割り当てる関数ｓ：Ｇ×２Ｇ→ＩＲが与えられていると仮定する。所与のネットワークＮに対し、Ｎのスコアをｓｃｏｒｅ（Ｎ）＝ｄｅｆΣｇ∈Ｇｓ（ｇ，ＰＮ（ｇ））として定義し、ここに、ＰＮ（ｇ）はＮにおけるｇの親の集合を表わす。
【００７９】
例：
１．ＢＤｅスコア（文献２、３参照）。
所与のデータに対し、スコアはネットワークの事後確率に比例する。ＢＤｅスコアを用いた場合、マイクロアレイデータを離散化する必要がある。
２．ＭＤＬスコア（文献３参照）。ＭＤＬスコアは、最小記述長原理を利用し、離散化されたデータも使用する。
３．ＢＮＲＣスコア（文献９参照）。ＢＮＲＣスコアは、非線形遺伝子相互作用を取得すべくノンパラメトリック回帰を用いる。データを離散化する必要がないため、情報が失われない。ネットワークを推定する作業は各々の遺伝子に対して親遺伝子の集合を発見して、結果的に得られたネットワークが非巡回且つネットワークのスコアが最小であるようにすることである。
【００８０】
以下に、本アルゴリズムを記述するために必要ないくつかの表記を定義する。
【００８１】
定義２．１：Ｆ
Ｆ：Ｇ×２Ｇ→Ｒを、全ての
【数８】
2007519988.tif 000014
およびＡ⊆Ｇに対してＦ（ｇ，Ａ）＝ｄｅｆｍｉｎＢ⊆Ａｓ（ｇ，Ｂ）と定義する。
【００８２】
定義により、Ｆ（ｇ，Ａ）の意味は、親を部分集合Ａから選択する必要がある場合における遺伝子ｇの親の最適選択である。全ての非巡回グラフについて、全ての辺が順序付けの方向に向けられるように頂点の順序付けがある。逆に、Ｇの固定順序を所与の場合、次の定義でするように、所与の順序に従う全てのグラフの集合を考えることができる。
【００８３】
集合Ａ⊆Ｇの順序付けを、順列π：｛１，・・・，｜Ａ｜｝→Ａと記述することができる。ΠＡを用いてＡの全ての順列の集合を表わす。
【００８４】
定義２．２：  π線形性
Ａ⊆Ｇ且つ
【数９】
2007519988.tif 000015
とする。Ｎ⊆Ａ×Ａをネットワークとする。全ての（ｇ，ｈ）∈Ｎに対して、π−１（ｇ）＜π−１（ｈ）が成立すればＮはπ線形であると言う。
【００８５】
ここで、上の定義を用いて、関数ＱＡを定義することにより、以下に示すように、所与のπに対して、最適π線形ネットワークのスコアを計算することができる。
【００８６】
定義２．３：  ＱＡ
Ａ⊆Ｇとする。全ての
【数１０】
2007519988.tif 000016
に対して、ＱＡ：ΠＡ→Ｒを次式で定義する。
【数１１】
2007519988.tif 000017
【００８７】
関数ＦおよびＱを用いて、所与の順列πに対する最適π線形ネットワークを計算することができるならば、最適ネットワークを発見するためにすべきことは、グローバル最小値を与える最適順列πを発見することである。形式的に、このステップのための関数Ｍを定義する。
【００８８】
定義２．４：  Ｍ
全てのＡ⊆Ｇに対して、Ｍ：２Ｇ→Ａ⊆ＧΠＡを次式で定義する。
【数１２】
2007519988.tif 000018
【００８９】
アルゴリズム２．１
上の表記を用いて、アルゴリズム２．１を以下のように定義することができる。
ステップ１：全てのｇ∈Ｇに対して、Ｆ（ｇ，φ）＝ｓ（ｇ，φ）を計算する。
ステップ２：全てのＡ⊆Ｇ、Ａ≠φおよび全てのｇ∈Ｇに対して、Ｆ（ｇ，Ａ）を
【数１３】
2007519988.tif 000019
として計算する。
ステップ３：Ｍ（φ）＝φに設定する。
ステップ４：全てのＡ⊆Ｇ、Ａ＝φに対して、以下の二つのステップを実行する。
ステップ４ａ：ｇ＊＝ａｒｇ  ｍｉｎｇ⊆Ａ（Ｆ（ｇ，Ａ−｛ｇ｝）−ＱＡ−｛ｇ｝（Ｍ（Ａ−｛ｇ｝）））を計算する。
ステップ４ｂ：全てのｌ≦ｉ＜｜Ａ｜に対して、Ｍ（Ａ）（ｉ）＝Ｍ（Ａ−｛ｇ＊｝（ｉ）およびＭ（Ａ）（｜Ａ｜）＝ｇ＊に設定する。
ステップ５：ＱＧ（Ｍ（Ｇ））を返す。
【００９０】
ステップ２およびステップ４で与えられた再帰式において、濃度ｍ＝｜Ａ｜の部分集合Ａ⊆Ｇに対して関数ＦおよびＭを各々計算することが望まれており、濃度ｍ−１の部分集合について関数Ｆ  およびＭの各々の関数値を必要とする。従って、ステップ４だけでなくステップ２でも動的計画法を適用して、濃度が増大する部分集合Ａについて関数ＦおよびＭを各々計算することができる。ステップ４の再帰式において、最初に順列Ｍ（Ａ）の最後の要素ｇをステップ４ａで計算し、次いでＭ（Ａ）をステップ４ｂで設定する。
【００９１】
正しさと時間複雑度
アルゴリズム２．１のステップ２における再帰式の正しさは直接Ｆの定義に従う。従って、ステップ１およびステップ２を実行した後で、全ての遺伝子ｇおよび全ての部分集合Ａ□Ｇに対する関数Ｆの値がメモリに格納される。ステップ３およびステップ４へ進む前に、関数ＱＡの意味に関する補題を記述する。
【００９２】
補題２．１
Ａ⊆Ｇ且つπΠＡとする。Ｎ＊⊆Ａ×Ａを最小スコアを有するπ線形ネットワークとする。従って、ＱＡ（π）＝ｓｃｏｒｅ（Ｎ＊）が成り立つ。
【００９３】
証明：π線形グラフにおいて、遺伝子ｇは、πにより符号化された順序の上流にある親ｈだけ持つことができる。すなわちπ−１（ｈ）＜π−１（ｇ）である。従って、ｇの親を選択する場合、Ｂ＝｛ｈ∈Ａ｜π−１（ｈ）＜π−１（ｇ）｝に制約され、Ｆ（ｇ，Ｂ）がこの場合の最適な選択である。π線形グラフの全ての辺が、πにより符号化された順序に従うため、このように全ての遺伝子に対して独立に親を選択することができ、それにより補題２．１が証明される。
【００９４】
補題２．１を用いて、関数Ｍがステップ４で与えた式により計算できることを証明する。
【００９５】
補題２．２
Ａ⊆Ｇとする。ｇ＊＝ａｒｇ  ｍｉｎｇ∈Ａ（Ｆ（ｇ，Ａ＝｛ｇ｝）＋ＱＡ−｛ｇ｝（Ｍ（Ａ−｛ｇ｝）））とする。π（ｉ）＝Ｍ（Ａ−｛ｇ＊｝）（ｉ）、およびπ（｜Ａ｜）＝ｇによりπ∈ΠＡを定義する。
従って、π＝Ｍ（Ａ）である。
【００９６】
証明：π’∈ΠＡとする。Ｍの定義により、ＱＡ（π）ＱＡ（π’）を示す必要がある。Ｎ＊を最適π’線形ネットワーク、Ｍ＊を最適π’線形ネットワークとする。従って、補題２．１により、ＱＡ（π）ＱＡ（π'）はｓｃｏｒｅＮ＊ｓｃｏｒｅＭ＊に同値である。π’の最後の要素をｈ＝π’（｜Ａ｜）と表記する。任意のＢ⊆Ｇに対して、ＱＢ（Ｍ（Ｂ））が上記の定義によるＢに対するグローバル最適ネットワークのスコアである点に注意する。従って、次式を得る。
ｓｃｏｒｅ（Ｍ＊）＝ｓ（ｈ，ＰＭ＊（ｈ））＋Σｇ∈Ａ−｛ｈ｝ｓ（ｇ，ＰＭ＊（ｇ））
≧ｓ（ｈ，ＰＭ＊（ｈ））＋ＱＡ−｛ｈ｝（Ｍ（Ａ−｛ｈ｝））
≧Ｆ（ｈ，Ａ−｛ｈ｝）＋ＱＡ−｛ｈ｝（Ｍ（Ａ−｛ｈ｝））
≧ｍｉｎｈ∈ＡＦ（ｈ，Ａ−｛ｈ｝）＋ＱＡ−｛ｈ｝（Ｍ（Ａ−｛ｈ｝））
＝Ｆ（ｇ＊，Ａ−｛ｇ＊｝）＋ＱＡ−｛ｇ＊｝（Ｍ（Ａ−｛ｇ＊｝））ｓｃｏｒｅ（Ｎ＊）、
これは、補題２．２を示す。
【００９７】
Ｆを用いてＱを直接計算できるため、本アルゴリズムはステップ５でＱＧ（Ｍ（Ｇ））を計算することができる。最後に、ＱＧ（Ｍ（Ｇ））は定義により最適ベイジアンネットワークのスコアであり、その正しさを示す。
【００９８】
全ての遺伝子ｇおよび全ての部分集合Ａ⊆Ｇに対して最適な親の情報がＦ（ｇ，Ａ）と共に格納されていれば、最適ネットワークはＱＧ（Ｍ（Ｇ））を計算する間に構築することができる。
【００９９】
定理２．１  制約の無い最適ネットワーク
アルゴリズムは、Ｏ（ｎ・２ｎ）の動的計画法ステップを用いて最適ネットワークを発見する。
【０１００】
証明：ステップ１およびステップ２における動的計画法はＯ（ｎ・２ｎ）（ｎ＝｜Ｇ｜）ステップを要し、各ステップで１個のスコアが計算される。ステップ３およびステップ４の動的計画法において、Ｏ（２ｎ）ステップを要するが、各ステップは既に格納された若干のスコアを参照することを含んでいる。関数ＱＡをステップ４ａで実際に計算する必要がない点に注意されたい。その理由は、以前のステップでＱＡ−｛ｇ｝をＭ（Ａ−｛ｇ｝）と共に格納できるからである。従って、全体的な時間複雑度はＯ（ｎ・２ｎ）である。
【０１０１】
生物学的現実では、調整遺伝子の子の数が極めて多いが、親の数は制約されると仮定することができる。親の数を制約する場合、スコア計算回数は実質的に減少するため、より大きいネットワークの計算が可能になる。
【０１０２】
次に示す系は実際に極めて有用である（第２．３節参照）。
【０１０３】
系２．１  ｍ∈ＩＮを定数とする。ｍ個を超える親を持つ遺伝子が存在しない最適ネットワークをＯ（ｎ・２ｎ）の動的計画法ステップで発見することができる。［引用」親の数を定数で制約したくなく、その代わり各遺伝子に対して一定数の親候補を選択できる場合、複雑度は以下のように変化する。
【０１０４】
系２．２  ｍ∈ＩＮを定数とする。各々のｇ∈Ｇに対して、Ｃｇ∈ＧをＣｇ≦ｍである集合とする。各遺伝子ｇがＣｇ内だけに親を持つ最適ネットワークをＯ（２ｎ）の動的計画法ステップで発見することができる。
【０１０５】
証明：各遺伝子の親が一定サイズの集合から選択されるため、ステップ１およびステップ２における動的計画法の複雑度は一定になる。従って、全体的な複雑度はＯ（２ｎ）になる。
【０１０６】
本アルゴリズムにおける動的計画法の二通りの適用は、動的計画法の単一の適用として行なうことができる点に注目する。これは、サイズｍの集合に対して関数Ｍを計算する際に、サイズｍ−１の集合に対する関数Ｆの関数値だけを必要とするからである。従って、サイズｍ−１およびｍの集合に対する関数ＦおよびＭの関数値だけを同時にメモリに格納すればよい。これは、必要なメモリ容量を減らすために現実的に有意味である。
【０１０７】
本アルゴリズムが、準最適解の計算用にも修正可能である点にも注目する。二番目に最適または三番目に最適なネットワークを計算することは、限界付近でスコアが変化する場合に推定されたネットワークの安定性を評価するために重要であり得る。
【０１０８】
結果
上述のアルゴリズム２．１は、Ｃ＋＋プログラムとして実装された。スコアリング関数として、ＢＮＲＣスコア、ＢＤｅスコア、およびＭＤＬスコアの既存の実装例を用いる。３方法（定理２．１、系２．１および２．２）の全てを実装した。
【０１０９】
本プログラムを１７３個のマイクロアレイのデータ集合に適用して、各種のストレス条件に対する出芽酵母の反応を測定した。（文献５参照）。
【０１１０】
熱ショック・データへの適用
データ集合から、２５℃〜３７℃にわたる熱ショック実験から１５個のマイクロアレイを、また様々な温度から３７℃にわたる熱ショック実験から５個のマイクロアレイを選択した。次いで、熱ショック反応に関与または関与した推測される９個の遺伝子の集合を選択した。図２に、ＢＮＲＣスコアに関して最適ネットワークを示す。
【０１１１】
転写因子ＭＣＭ１が３個の他の遺伝子を調整するものと推定されるが、一方、当該因子はこの集合の遺伝子の１個により調整されないことが観察されたが、これは可能である。当該遺伝子集合の第二の転写因子ＨＳＦ１は、３個の他の熱ショック遺伝子を調整するものと推定される。これはまた、既に報告された（文献１６参照）ように、ＨＳＰ７０タンパク（ＳＳＡ１）により調整されると推定される。これらの遺伝子中の他のシャペロンＳＳＡ３もまた、熱ショック反応において能動的な役割を果たしてＳＳＡ１およびＨＳＰ１０４と相互作用するように見えるが、これはグローバー（Ｇｌｏｖｅｒ）とリンドクイスト（Ｌｉｎｄｑｕｉｓｔ）による報告（文献９８参照）と一致する。
【０１１２】
表２．１に上述の解析結果を示す。総括するに、この結果は生物学的に可能であり、熱ショック反応の間におけるシャプロンＳＳＡ１およびＳＳＡ３の能動的な役割を示唆する。最適に推定された遺伝子ネットワークが有意味且つ遺伝子調整の解明に役立つものと結論する。
【０１１３】
【表５】
2007519988.tif 000020
【０１１４】
計算可能性および制限
第２．３．１節において推定されたネットワークのような小規模ネットワークでさえ力任せの方法（式１）では推定できないが、発明者のプログラムにより１．９ＧＨｚで動作する単一のＰｅｎｔｉｕｍ（登録商標）ＣＰＵを用いて約１０分間で最適に推定することができる。本方法の現実的可能性を推定すべく、データ集合から遺伝子調整における既知の能動的な役割（文献１２参照）を有する２０個の遺伝子を選択して、１７３個のマイクロアレイの全てを用いて最適ＢＮＲＣスコアを有するネットワークを推定した。計算は、各々９００ＭＨｚで動作する９６個のＣＰＵを備えたＳｕｎＦｉｒｅ１５Ｋスーパーコンピュータを用いて約５０時間以内に完了した。この計算実験の結果として、本方法が２０個の遺伝子を含む遺伝子ネットワークに対し、たとえ一切の制約を加えず、ＢＮＲＣスコアのように複雑なスコアリング方式を用いる場合でも、実行可能であると結論する。離散的なスコアＢＤｅおよびＭＤＬの場合、はるかに高速に計算できるため、２０個を超える遺伝子のネットワークでさえ制約無しに最適に推定することができる。
【０１１５】
親の数が約６に制限されている（系２．１）か、あるいは、約２０個の親候補の集合が事前に選択されている（系２．２）場合、たとえＢＮＲＣスコアを用いるにせよ３０個を超える遺伝子の遺伝子ネットワークを最適に推定することができる。しかし、本方法は現状では、約４０個を超える遺伝子のネットワークを推定することはできない。
【０１１６】
系２．２で与えられる方法の理論上の時間複雑度は、系２．１で与えられる方法の時間複雑度を下回るものの、現実的には後者がより重要であると主張する。第一に、親の数を定数で制限することは容易に行なえ、且つ生物学的に正当化される一方、各遺伝子に対して親候補の集合を選択することは遺伝子選択の方法を必要とするため、潜在的に計算結果を偏らせる恐れがある。第二に、関数Ｆの計算における各々の動的計画法のステップは１個スコアの計算を要するのに対し、関数Ｍに対して１個の動的計画法のステップはいくつかの以前の結果を参照するだけで済むことを考慮しなければならない。親の数が系２．１に記載のように制限されている場合、必要なスコア計算の回数は多項式となり、系２．２の方法が理論的に優れているにもかかわらず、現実的な適用では当該方法の方がより早く計算できる。
【０１１７】
使用する確率分布および追加的な仮定をするか否かに依存して、２０〜４０個の遺伝子の遺伝子ネットワークを最適に推定することが可能な方法を提示した。これにより、最適解が得られるため、異なるスコアリング方式を比較し、所与のスコアリング方式に対して最適なパラメータを評価して、所与のマイクロアレイデータの有用性を評価することが可能になる。また、本方法は特に、研究者が遺伝子の特定のグループに注目して、これらの遺伝子に関する遺伝子発現測定を最大限に利用したい場合の設定に役立つ。
【０１１８】
ヒューリスティックな方法とは対照的に、結果が満足できないか、生物学的知識に矛盾する場合、統計モデルが不正確であるか、またはデータが不十分であると結論付けることができる。２０個の遺伝子のネットワークの場合でさえ、所与の巨大な探索空間から最適なネットワークを知るには膨大な情報量を要する。
【０１１９】
本方法が特定のスコアリング方式、または特定の種類の遺伝子発現測定に依存しない点に注目する。これは、第２．２節で定義したスコアが与えられる任意の設定に適用可能である。例えば、シーケンス情報（文献１９参照）、タンパク質相互作用データ（文献１０参照）、または他の知識がスコア関数に取り込まれている場合、本方法もまた適用可能である。
【０１２０】
４０個を超える遺伝子を含む遺伝子ネットワークを発見すべく、二つの方向を利用することができる。第一に、アルゴリズムが実際の検索を実行する対象である部分集合の集合の一部が剪定できる場合、実行現可能性の制限が増大することがある。第二に、遺伝子ネットワークの区分化（文献１８参照）を用いて、より巨大ネットワークをより小さい部分に分解して、各々の部分的ネットワークを最適に推定することができる。
【０１２１】
参考文献
1. D. M. Chickering. Learning Bayesian networks is NP-complete. In D. Fisher and H.-J. Lenz, editors,Learningfrom Data :Artificial Intelligence and Statistics V, Springer-Verlag, 1996. 
2. G. F. Cooper, E. Herskovits. A Bayesian method for the induction of probabilistic networks from data. Machine Learning, 9: 309--347,1992. 
3. N. Friedman, M.Goldszmidt. Learning Bayesian networks with local structure. Jordan, M. I. (ed. ), Kluwer Academic Publishers, pp. 421--459, 1998. 
4. N. Friedman, M. Linial,1. Nachman, D. Pe'er. Using Bayesian networks to analyze expression data. Journal of Computational Biology, 7:601--620, 2000. 
5. A. P. Gasch, et. al. Genomic expression programs in the response of yeast cells to environmental changes. Molecular Biology of the Cell, 11: 4241--4257, 2000. 
6. J. R. Glover, S. Lindquist.Hspl 04, Hsp70, and Hsp40 : a novel chaperone system that rescues previously aggregate proteins. Cell, 94: 73--82, 1998. 
7. A. J. Hartemink, D. K.Gifford, T. S. Jaakkola, R. A. Young. Using graphical models and genomic expression data to statistically validate models of genetic regulatory networks. Pacific Symposium on Biocomputing, 6: 422-433, 2001. 
8. A. J. Hartemink, D. K. Gifford, T. S. Jaakkola, R. A. Young. Combining location and expression data for principled discovery of genetic regulatory network models. Pacific Symposium on Biocomputing, 7:437--449, 2002. 
9. S. Imoto, T. Goto, S. Miyano. Estimation of genetic networks and functional structures between genes by using Bayesian networks and nonparametric regression. Pacific Symposium on Biocomputing, 7: 175--186, 2002. 
10. T. Ito, T. Chiba, R. Ozawa, M. Yoshida, M. Hattori, Y. Sakaki. A comprehensive two-hybrid analysis to explore the yeast protein interactome. Proceedings of the National Academy of   Sciences USA, 97:4569-4574, 2001. 
11. S. Imoto, S. Kim, T. Goto, S. Aburatani, K. Tashiro, S. Kuhara, S. Miyano. Bayesian network and nonparametric heteroscedastic regression for nonlinear modeling of geneticnetwork. Journal ofBioinformatics and Computational Biology, in press, 2003; U. S. Patent Application No : 10/259,723. 
12. T. I. Lee, N. J. Rinaldi, F. Robert, et. al. Transcriptional regulatory networks in Saccharomyces cerevisiae. Science, 298: 799--804, 2002. 
13.I. M. Ong, J. D. Glasner, D. Page.Modellingregulatorypathwaysin E. coli from time series expression profiles. Bioinformatics, 18: 241--248,2002. 
14. D. Pe'er, A. Regev, G. Elidan, N. Friedman. Inferring subnetworks from perturbed expression profiles. Bioinformatics, 17: 215--224,2001. 
15. R. W. Robinson. Counting labeled acyclic digraphs. New Directions in the Theory of Graphs, pp. 239--273, 1973. 
16. Y. Shi, D. D. Mosser, R. I. Morimoto. Molecular chaperones asHSF1-specific transcriptional repressors. Genes & Development, 12 : 654--666, 1998. 
17. V. A. Smith, E. D. Jarvis, A. J. Hartemink. Evaluating functional network inference using simulations of complex biological systems. Bioinformatics, 18:216--224, 2002. 
18. E. P. van Someren, L. F. A. Wessels, E. Backer, M. J. T. Reinders. Genetic network modeling. Pharmacogenomics, 3 (4):507--525, 2002. 
19. Y. Tamada, S. Kim, H. Bannai, S. Imoto, K. Tashiro, S. Kuhara, S. Miyano.
【０１２２】
例３：最適なものより良い遺伝子ネットワーク
３．１  序論
発現レベル測定から遺伝子ネットワークの推定は、昨今の生物情報科学研究の一つの焦点である（１、５、８、１３、３１）。発現データから遺伝子ネットワークの構造を知ることは、細胞の自身の環境に対する反応の理解および遺伝子ネットワークの構造原理に関する知識（１０、２９）を深めることができ、いくつかの分野において応用が可能になろう。
【０１２３】
遺伝子ネットワークのモデル化に広く利用されている方法として、ベイジアンネットワークがあり（３、４、９、１３、１５、２６、３０）、遺伝子の発現レベルを確率変数として、また遺伝子ネットワークを同時確率分布としてモデル化する。これらの分布は有向非巡回グラフを用いて分解され、これをネットワークと呼ぶ。ネットワークは、データを与えられたならば、ネットワークの尤度に基づくスコア関数を用いて記録される。
【０１２４】
最近の結果は、約３０個以下の遺伝子の遺伝子ネットワークに対していくつかのデータに関する最適遺伝子ネットワーク・モデルの知見が得られることを示す（２５）。この結果は、遺伝子ｇおよびｇの親の集合にスコアを割り当てる任意のスコア関数ｓ：Ｇ×２Ｇについて成立する。しかし、最適遺伝子ネットワーク・モデルは最も可能性のあるモデルであるが、いくつかのデータに関してほぼ同じ可能性を有する非常に異なるモデルが依然として有り得る。これは特に、遺伝子発現測定が一般に遺伝子ネットワークに関して部分的な情報しか提供しないためである。更に、遺伝子を１０個しか含んでいない遺伝子ネットワークでさえ、４．１７＊１０１８通りの可能なネットワーク・モデルが存在する（１７）。従って、最適な遺伝子ネットワーク・モデルでも一般に対象遺伝子ネットワークに合致しない。
【０１２５】
この問題に取り組む発明者の努力は三通りある。第一に、最適および準最適ネットワークをその尤度順に列挙するための理論的基礎およびアルゴリズムを提供する。第二に、最適遺伝子ネットワーク・モデルを既存の遺伝子ネットワークに関して利用可能な知識と厳密に比較する。第三に、最適遺伝子ネットワーク・モデルの信頼性が高い部分を抽出する方法を提示し、関連する種の遺伝子ネットワーク・モデルを比較すべく、この方法を枯草菌および大腸菌の利用可能なデータに適用する。
【０１２６】
この結果は、本方法により識別される部分的ネットワークが、最適ネットワーク・モデル自身よりも顕著に当該知識とより良く合致することを示す。本方法の基礎を最適なｎ遺伝子ネットワーク・モデルに置いているため、これらの推定から、ヒューリスティックに依存する方法よりも高い信頼性を以って結論を導くことができる。
【０１２７】
この理論的結果を第３．２節に示した後で、推定と知識を比較した結果を第３．３節に記述して議論する。次いで、本発明の遺伝子ネットワーク推定方法を第３．４節に示し、第３．５節においてこれを評価して適用する。
【０１２８】
３．２  遺伝子ネットワークの推定および列挙
本例を通じて、遺伝子Ｇの集合と、遺伝子ｇにスコアを割り当てるスコア関数ｓ：Ｇ×２Ｇ→ＩＲ、およびｇの親の集合が与えられていると仮定する。ネットワークＮを与えられたならば、Ｎのスコアをｓｃｏｒｅ（Ｎ）＝Σｇ０Ｇｓ（ｇ，ＰＮ（ｇ））として定義する。ここに、ＰＮ（ｇ）はＮにおけるｇの親の集合を表わす。与えられたデータを最も良く説明する同時確率分布を発見するには、ｓｃｏｒｅ（Ｎ）を最小化する有向非巡回グラフＮを発見する必要がある。
【０１２９】
この最適化問題はＮＰ困難であり（２）、探索空間が超指数サイズを有する（２７）ため、研究者は往々にしてグリーディ・アルゴリズム、シミュレーテッド・アニーリング、あるいは遺伝的アルゴリズム等のヒューリスティックな方法を適用してきた（５、１３、３１）。しかし、最近以下の結果が導かれて、実際の遺伝子ネットワークの遺伝子は親の数に制限がある事実を利用すれば、制約のない約２０個の遺伝子の遺伝子ネットワーク、および３０個または若干これより多い遺伝子を含むネットワークを推定することが可能なアルゴリズムが得られる。
【０１３０】
定理３．１（２５）
最適ネットワークは、Ｏ（ｎ・２ｎ）の動的計画法のステップを用いて発見することができる。
【０１３１】
より巨大遺伝子ネットワークの場合、経験的またはヒューリスティックな方法を用いて探索空間の部分空間を選択することができる。これを（２４）に記述されているように実行したならば、（２５）のアルゴリズムを用いて選択された部分空間内で最適解を発見することができ、決定論的且つヒューリスティックな技術を組み合わせた方法が得られる。
【０１３２】
しかし、異なるネットワークが同一発現データを同等またはほぼ同等の確率で与える場合がある。従って、所与のデータが遺伝子ネットワークに関する完全な情報を含むと仮定することはできない。しかし、所与のデータ集合から導くことができる情報が極めて部分的な情報であり得ることに留意しなければならない。このように、たとえ最適ネットワークが発見されても、それは一般に間違った辺を含んでいよう。この問題を解決すべく、発明者は（２５）のアルゴリズムを拡張して以下の結果を証明した。
【０１３３】
定理３．２
ｍ０ＩＮとする。最適なｍネットワークをＯ（ｎ・２ｎ）の動的計画法ステップを用いて発見することができる。
【０１３４】
３．２．１  証明−最適遺伝子ネットワークの列挙
３．２．１．１  アルゴリズム３．１
以下に、最適および準最適ネットワークの列挙が可能にすべく、（２５）からアルゴリズムを拡張する方法を示す。最初にいくつかの関数を定義し、次いで、これらの関数が相当な規模の遺伝子ネットワークについてどのように計算できるかを示す。第３．２節と同様に、遺伝子Ｇの集合およびスコア関数ｓ：Ｇ×２Ｇ！ＩＲが与えられると仮定する。ネットワークＮが与えられたならば、Ｎのスコアをｓｃｏｒｅ（Ｎ）＝３ｇ０Ｇｓ（ｇ，ＰＮ（ｇ））として定義される。ここに、ＰＮ（ｇ）はＮにおけるｇの親の集合を表わす。本節を通じて、同一スコアを有するネットワークは同様にソートされるとみなし、従って「第ｎ番目に最適なネットワーク」との表記を許す。
【０１３５】
定義３．１：Ｆｍ
ｍ∈ＩＮとする。誘導的にＦｍ：Ｇ×２Ｇ×Ｉ→→２Ｇを定義する。最初に、全てのｇ∈ＧおよびＡ⊆Ｇに対して、以下を定義する。
Ｆｍ（ｇ，Ａ，ｎ）＝ｄｅｆａｒｇ  ｍｉｎ  ｓ（ｇ，Ｂ）
Ｂ⊆Ａ    （２）
次いで、先の全ての解の集合｛Ｆｍ（ｇ，Ａ，ｐ）｜ｐ＜ｎ｝をＪ（ｎ）と表記すれば、
全ての１＜ｎ≦ｍに対して以下を得る。
Ｆｍ（ｇ，Ａ，ｎ）＝ｄｅｆａｒｇ  ｍｉｎ  ｓ（ｇ，Ｂ）
Ｂ⊆Ａ    （３）
Ｂ∈Ｊ（ｎ）
【０１３６】
定義３．２：Ｓｍ
ｍ∈ＩＮとする。全てのｇ∈Ｇ、Ａ⊆Ｇ、およびｎ∈ＩＮ≦ｍに対してＳｍ：Ｇ×２Ｇ×ＩＮ≦ｍ→ＩＲを以下のように定義する。
Ｓｍ（ｇ，Ａ，ｎ）＝ｄｅｆｓ（Ｆｍ（ｇ，Ａ，ｎ））  （４）
【０１３７】
定義により、親をＡから選択しなければならず、Ｓｍ（ｇ，Ａ，ｎ）がこの選択のスコアである場合、Ｆｍ（ｇ，Ａ，ｎ）が遺伝子ｇの親のｎ番目に最適な選択である。文脈からｍが明らかな場合、ＦｍおよびＳｍの代わりに各々ＦおよびＳを用いる。ｍがＡの部分集合の個数より大きい場合もあるため、ＦｍおよびＳｍが部分的に定義された関数である点に注意されたい。 
【０１３８】
集合ＡｆＧの順序付けを順列Ｂ：｛１，・・・，｜Ａ｜｝→Ａと記述することができる。ΠＡを用いてＡの全ての順列の集合を表わす。全ての辺が順列πにより与えられる方向に従うネットワークを示す以下の表記を必要とする。
【０１３９】
定義３．３：  π線形性
Ａ⊆Ｇ且つπ∈ΠＡとする。Ｎ⊆Ａ×Ａをネットワークとする。全ての（ｇ，ｈ）∈Ｎに対してπ−１（ｇ）＜π−１（ｈ）が成立するならばＮはπ線形であると言う。
【０１４０】
本アルゴリズムの基本戦略は、全てのπ∈ΠＡに対して、π線形ネットワークの部分集合内の集合Ａ⊆Ｇ上の全ての有向非巡回グラフの空間を分割して、最適および準最適ネットワークを発見する問題を以下の二つの問題に分解する。
１．探索される（準）最適ネットワークを含む探索空間の部分集合を発見する。（補題２）
２．選択された部分空間内で（準）最適ネットワークを発見する。
【０１４１】
（補題１および補題２）
探索空間の部分空間をπを用いて表記する。所与の順列πに対する最適および準最適ネットワークを発見すべく、以下の関数ＱＡを必要とする。所与の遺伝子ｇに対して、π内でｇに先行する全ての遺伝子の集合をＶ（π，ｇ）＝ｄｅｆ｛ｈ｜π−１（ｈ）＜π−１（ｇ）｝と表記する。
【０１４２】
定義３．４：ＱＡ
Ａ⊆Ｇとする。全てのπ∈ΠＡに対して、ＱＡ：ΠＡ×ＩＮ｜Ａ｜→２Ａ×Ａを
ＱＡ（π，ｖ）＝ｄｅｆ｛（ｈ，ｇ）｜ｈ∈Ｆ（ｇ，Ｖ（π、ｇ）、ｖπ−１（ｇ））｝（５）
と定義する。
【０１４３】
定義４において、特定の遺伝子の親の選択のランクを決定するために、ベクトルｖ∈ＩＮ｜Ａ｜を用いた。補題１および補題２において、ＱＡ（π，ｖ）が最適または準最適なπ線形ネットワークであり、そのランクはｖに依存することが示される。次に、部分空間を指定する２個の関数ＭｍおよびＤｍを定義し、これらは各々、内部で（次善）最適ネットワークを発見することが可能な部分空間、および当該部分空間からのネットワークの選択を指定する。
【０１４４】
定義３．５：Ｍｍ，Ｄｍ
ｍ∈ＩＮとする。Ｍｍ：２Ｇ×ＩＮ≦ｍ→∪Ａ⊆ＧΠＡおよび
【数１４】
2007519988.tif 000021
を、それらの第二パラメータについて誘導的に定義する。Ａ⊆Ｇとする。最初に、次式を定義する。
Ｄｍ（Ａ，１）＝ｄｅｆ（１，・・・，１）∈ＩＮ｜Ａ｜    （６）
および
Ｍｍ（Ａ，１）＝ｄｅｆａｒｇ  ｍｉｎ  ｓｃｏｒｅ（ＱＡ（π，Ｄｍ（Ａ，１）））  （７）
π∈ΠＡ
ｎ＞１のときπ∈ＩＮ≦ｍとし、ＮをＡに関するネットワークであってネットワークの最適スコアが｛ＱＡ（Ｍｍ（Ａ，ｐ），ＤＭ（Ａ，ｐ））｜ｐ＜ｎ）｝に含まれないものとする。π＊∈ΠＡを、Ｎがπ＊線形であるような順列とする。次式を定義する。
Ｍｍ（Ａ，ｎ）＝ｄｅｆπ＊      （８）
ｖ＊∈ＩＮ｜Ａ｜を、全てのｇ∈Ａに対してｇの親の集合ＰＮ（ｇ）が以下に等しいようにする。
Ｆｍ（ｇ，Ｖ，（π＊，ｇ），ｖ＊π＊−１（ｇ））５    （９）
次式を定義する。
Ｄｍ（Ａ，ｎ）＝ｄｅｆｖ＊    （１０）
ＦｍとＳｍ、およびＭｍとＤｍは部分関数である。表３．４において、上の表記を要約する。これらの表記を用いて、本アルゴリズムは所与のｍ０ＩＮに対してフローとして定義することができる。
【０１４５】
【表６】
2007519988.tif 000022
【０１４６】
これらの表記法用いて、アルゴリズム３．１は以下のように定義することができる。
【表７】
2007519988.tif 000023
【０１４７】
アルゴリズム３．１は、ステップ１およびステップ２において全てのｇ∈Ｇ、Ａ⊆Ｇ、およびｎ≦ｍに対して関数ＦｍおよびＳｍを計算する。ステップ２ａでＢ？を選択するために、より低濃度の集合Ａまたはより小さいｎに対するＦｍの関数値だけを必要とする。従って、ＦｍおよびＳｍは動的計画法を適用して計算することができる。
【０１４８】
ステップ３およびステップ４において、関数ＭｍおよびＤｍは動的計画法を用いて同様に計算することができる。その理由は、ステップ４ａにおける三つ組を選択する場合、Ｍｍ（Ａ，ｎ）およびＤｍ（Ａ，ｎ）を計算するために低い濃度の集合Ａまたはより低い濃度のｎに対するＭｍおよびＤｍの関数値だけを必要とするからである。三つ組（ｇ，ｐ，ｑ）の意味は次の通りである。ｇは探索される順列の最後の要素候補である。ｇが最後の要素になった場合、残りの順列はＡ−｛ｇ｝の既に計算された最大ｍ個の順列から選択することができる。選択された残りの順列はｐにより指定される。次いで、結果的に生じる順列により定義される部分空間内にネットワークを形成するために、ｇの親のｑ番目に最適な選択を用いる一方、他の遺伝子に対してはＤｍ（Ａ−｛ｇ｝，ｐ）で示すように親が選択される。
【０１４９】
４個の関数Ｆｍ、Ｓｍ、ＭｍおよびＤｍは全て部分的に定義されているため、低い濃度の集合Ａに対してｎ個の関数値の全てが計算できる訳ではない。簡便のため、本アルゴリズムの定式化に際してこのことを明示しなかった。
【０１５０】
アルゴリズム３．１は上述のように、集合Ａサイズに関係なく、一定数（ｍ）の解を計算する。現実的な応用において、計算時間は、ｌ個の要素を含むＧの多数の部分集合を用いて層ｌに対してより少ない数の（準）最適解を計算することにより最適化することができる。次いで、Ａの濃度が高くなると部分集合の個数が下がる場合、ｍを増加させることができる。低位層についてｍ個の最適解しか知られていない場合はｍ個を超える（準）最適解が導かれる保証はないが、試験的な計算ではうまく機能した。
【０１５１】
計算を行なう間、格納しておく必要があるネットワークの個数が極めて大きくなる可能性があるため、ネットワークを効率的にメモリに格納する一方、格納されたネットワークを高速に復号できることが重要である。本発明の実装では、上述の定式化に基づいて、順列πを格納して、π線形ネットワークにおける辺の制約を利用することにより、ネットワークのメモリ効率および時間効率に優れた符号化が得られる。
【０１５２】
選択されたネットワークが既に選択されたネットワークとは異なるようにするとのステップ４ａの要件を満たすために、動的計画法のステップを実行する間にネットワークを比較しなければならない。しかし、異なるスコアを有するネットワークは異なるため、スコアが等しいネットワークだけを比較すればよい。従って、実用的な応用においてネットワーク比較の回数を最小化することができる。
【０１５３】
Ｆｍ、Ｍｍ、またはＤｍの１個の層を計算する際に、層内の集合について任意の順序でこれを実行してよく、計算は下位層の値だけに依存する。従って、上記のアルゴリズムは、相当に並列化可能である。
【０１５４】
３．２．３  正しさおよび複雑度
集合Ａ⊆Ｇにおけるｎ番目に最適なネットワークをＮ＊Ａ，ｎで表わす。（２５）から再定式化された二つの補題を最初に記述する。
【０１５５】
補題３．１  Ａ⊆Ｇ且つπ∈ΠＡとする。Ｎ＊⊆Ａ×Ａを最小スコアのπ線形ネットワークとする。従って、ＱＡ（π，（１，・・・，１））＝ｓｃｏｒｅ（Ｎ＊）が成り立つ。
【０１５６】
補題３．２  Ａ⊆Ｇ且つｍ∈ＩＮとする。ｇ＊＝ａｒｇ  ｍｉｎｇ∈Ａ（Ｓｍ（ｇ，Ａ＝｛ｇ｝，１）＋Ｎ＊Ａ−｛ｇ｝，ｎ）とする。π∈ΠＡをπ（ｉ）＝Ｍ（Ａ−｛ｇ＊｝，１）（ｉ）およびπ｜Ａ｜＝ｇ＊により定義する。従って、π＝Ｍｍ（Ａ，１）である。
【０１５７】
これらを用いて、上で定義されたアルゴリズムの正しさを証明することができる。
【０１５８】
定理３．２  ｍ∈ＩＮとする。
最適なｍ個のネットワークをＯ（ｎ・２ｎ）の動的計画法ステップを用いて発見することができる。
【０１５９】
証明：アルゴリズムの出力、すなわちＱＧ（Ｍｍ（Ｇ，ｉ），Ｄｍ（Ｇ，ｉ））、ｉ≦ｍが定義により、Ｇに最適なｍ個のネットワークである。本アルゴリズムで与えられた再帰式が正しいことだけを証明すればよい。ステップ１で与えられる式は、ＦｍおよびＳｍの定義から正しい。ステップ２において集合Ａ？Ｇの部分集合を選択する場合、基本的に全体集合Ａまたは真の部分集合のどちらか２択が可能である。前者の場合、選択のスコアを直接計算することができ、後者の場合、既に計算されたＦｍおよびＳｍの値を使用することができ、これはステップ２の正しさを与える。
【０１６０】
ステップ２の実行後、ＦｍおよびＳｍの全ての値が計算されている。これらの値を用いて、関数Ｑを直接計算することができる。従って、ステップ５の出力を生成するために関数ＭｍおよびＤｍを計算するだけでよい。ステップ３の式は、ここでも定義により正しい。補題３．１を定義３．５と組み合わせれば、以下の式が成り立つことがわかる。
Ｎ＊ａ，ｎ＝ＱＡ（Ｍｍ（Ａ，ｎ），Ｄｍ（Ａ，ｎ））
【０１６１】
上式および補題３．２から、ステップ４の再帰がｎ＝１に対して正しいことがわかる。ｎ＞１の場合、既に選択されていないネットワークに限って、準最適順列Ｍｍ（Ａ，ｎ）および親Ｄｍ（Ａ，ｎ）の準最適選択を同様に計算する。
【０１６２】
動的計画法の１ステップに要する時間はｍに依存するが、本研究の計算で選択されたｍに対して一定であると見なすことができる。更に、付録に記載のプログラム技術を用いれば、現実の計算では大多数の動的計画法ステップでｍ個を大幅に下回るネットワークを計算すれば十分である。
【０１６３】
本アルゴリズムを用いて、最適および準最適ネットワークを尤度順に列挙することができる。そのようなネットワークの共通要素を見出せたならば、共通要素の尤度はそれを含んでいるネットワークの尤度の合計である。従って、遺伝子ネットワークモチーフと呼ばれる列挙されたネットワークの共通要素、すなわち、遺伝子ネットワークモチーフの定義は、最初は、例えば（２１）で用いた表記とは異なるが、緊密に関係していることがわかるであろう。従って、ネットワークモチーフ解析は、１回のネットワーク推定よりも信頼性が高いと期待される。
【０１６４】
３．３．  遺伝子ネットワーク推定の信頼性評価
列挙された最適および準最適遺伝子ネットワーク・モデルが、実際の遺伝子ネットワークに関する有益な情報を含んでいることを検証すべく、アルゴリズムを実装してＲＮＡ  マイクロアレイデータに適用した。細菌内での転写と翻訳の緊密な関係を考慮して、ＲＮＡデータに基づく遺伝子ネットワーク推定だけが、真核生物よりも細菌により適していて、転写および翻訳が異なる場所および異なる時間に生起するものと予想する。従って、遺伝子ネットワークに関する知識だけでなくマイクロアレイデータも利用可能な枯草菌および大腸菌を対象に選択した。
【０１６５】
３．３．１  データおよびソフトウェア
大腸菌について、遺伝子発現オムニバス（ＧＥＯ）から、データ集合ＧＤＳ９５〜ＧＤＳ１００を選択した（３５）。遺伝子発現レベルにおける変化は、トリプトファン代謝、ＵＶ露出、およびノボビオシン処理の摂動により導出された（１７、１８、１９）。次いで、推定結果と比較すべく調整ＤＢ（２８）から、大腸菌における既知の転写調整に関するデータを受け取った。
【０１６６】
枯草菌について、各種の処理の下で経時的実験から７０個のマイクロアレイ、および遺伝子破壊株実験から９９個のマイクロアレイを含むいくつかのデータ集合を用いた。データは未だ公開されていないが、これらのデータを用いて生物学的に有意味な推定を行なうことが可能なことが確認された（１２）。次いで、ＤＢＴＢＳから既知の転写調整のデータ集合を受け取った（１６、２０）。
【０１６７】
本発明の実装で用いた、ＭＤＬスコア（４）、ＢＤｅスコア（３、４）、およびＢＮＲＣスコア（１３）を含むスコア関数ｓ：Ｇ×２Ｇ→ＩＲ）のうち、以下の理由により本研究での計算にＢＮＲＣスコアを選択した。第一に、ＢＮＲＣスコアは、追加的なパラメータおよび情報の消失を避けながら、離散化せずにデータを使用する。第二に、遺伝子相互作用はＢスプラインを用いてモデル化されるため、線形関係に制約されない一般的なモデリングが可能になる。第三に、枯草菌データ集合の３種のスコア関数の全てを用いた予備計算実験において、ＢＮＲＣスコアを用いた推定が他のものより明らかに教科書的知識と良く一致した（３２）。
【０１６８】
３．３．２  対象ネットワークの選択
大腸菌に関する既知の調整関係のデータ集合から、実験的証拠が提供される全ての関係を選択した。これにより８９９個の既知の関係の集合が得られた。枯草菌のデータ集合から、文献に証拠がある８４０個の調整関係を選択した。これらの巨大ネットワークの一部を選択すべく、ランダムな手順を適用した。選択された遺伝子同士に何らかの既知の調整関係が存在するように遺伝子を選択する必要があるため、最初の数個の遺伝子を無作為に選択し、次いで既に選択された遺伝子に結合している遺伝子を繰り返し選択することにより、各ステップにおいて選択された部分ネットワークを１個の遺伝子および少なくとも１個の辺だけ拡張する。各々の反復において、部分ネットワークの結合された要素を同じ確率で選択し、次いで、当該要素の少なくとも１個の遺伝子と既知の関係を有する遺伝子が存在する場合、そのような遺伝子を無作為に選択した。対象ネットワークの無意味な選択を避ける必要があるため、５個の結合された遺伝子が続けて選択されている場合、既に選択された遺伝子に結合していない遺伝子を選択した。
【０１６９】
選択手順により、有意味な構造を有する既知の遺伝子ネットワークＮが得られた。Ｎを行列で表現した。既知の関係を有する遺伝子の各ペアを行列の１要素で表わす。未知の遺伝子のペアについては、大部分のペアをゼロ入力で表わすが、以下の４条件の１個以上が成り立つペア（ｇ，ｈ）には０．５入力で表わす。
１．ｇはｈにより調整される
２．ｇおよびｈを調整する遺伝子ｉが対象ネットワークに存在する
３．ｇ（ｈ）により調整され、且つｈ（ｇ）を調整する遺伝子ｉが対象ネットワークに存在する
４．対象ネットワークに含まれない遺伝子ｉに対して条件２または条件３が成り立つ
【０１７０】
これらの条件を用いて、ほぼ正しい予測が間違った予測から区別された。辺（ｇ，ｈ）が予測されていて、且つ（ｈ，ｇ）が既知の調整関係である場合、少なくともこれら２個の遺伝子が相互作用する事実は正しく予測された。同様に、間接的な調整関係が、おそらくは対象ネットワーク自体には含まれない遺伝子を介して予測された場合、これも完全に間違っている訳はない。
【０１７１】
３．３．３  評価結果
第３．３．２節で述べたように、１０個の遺伝子の集合３０個を無作為に選択して、各対象ネットワークに本方法を適用した。次いで、上位５００個の最適ネットワーク構造において各々の関係が推定された回数を数えて、推定の頻度と生物学的知識との関係を調べた。本発明の方法の結果得られるネットワークは有向である点に注意されたい。しかし、上の節において述べたように、間違った方向を示す推定された辺も情報を含んでいる。従って、有向および無向の両方のケースを考察した。
【０１７２】
各々の推定された辺が正しいか否かを知ることは、未発見の関係があり得るため、困難である。更に、既知の遺伝子関係が当該データで表わされていないかもしれない。従って、推定された辺をデータベース情報に基づいて３個のグループ（１、０．５および０）に分類し、それらの違いを調べた。良い推定はこれらのグループをうまく分離できるものと考察した。
【０１７３】
図２（ａ）、２（ｂ）は各々、枯草菌の経時的および破壊株データの結果を示す。ｘ軸は、５００個のネットワークに各々の辺が出現する百分率を示す。ｙ軸は、当該区間に含まれる対応グループの辺の割合を示す。ｘ軸の両端で２個のピークが見られるため、これら二つの領域を調べた。
【０１７４】
【表８】
2007519988.tif 000024
【０１７５】
表３．２は、出現頻度の百分率を示す。表３．２から、３グループの傾向を観察した。０％付近では、グループ０は他より１０％よりも高い出現頻度を示した。これに対して、１００％付近では、出現頻度は最低値まで下がった。一方、グループ０．５はグループ１により近い値を示した。経時的データと破壊株データを比較して、破壊株データを使用する間、グループは相当うまく良く区別された。
【０１７６】
【表９】
2007519988.tif 000025
【０１７７】
最後に、大腸菌マイクロアレイデータを解析して、その結果を図３および表３．３に示した。興味深いことに、この結果においてグループ間の差異は若干小さい。提案した方法がスコアリング関数の長所だけでなく、データの正確さも評価できると考察する。
【０１７８】
３．４  遺伝子ネットワークモチーフの抽出
前節における評価は辺の数え上げのみに基づいているが、ここで、列挙されたネットワークの完全な情報を利用することを考える。
【０１７９】
３．４．１  モチーフ抽出問題
定理３．２を本例で述べている技術と組み合わせることにより、かなり大きい遺伝子ネットワークの最適および準最適ネットワークを列挙できる可能性を示す。列挙されたネットワークが与える情報を利用する直接的な方法は、ネットワークにおける各辺の出現回数を数え、ある閾値を超えるカウント数を有する辺だけを選択して、これらの辺から部分ネットワークを構成することである。この方法は、ブートストラップ法により列挙されたネットワークから部分ネットワークを抽出すべく（２６）で用いられた。辺の信頼レベルが相互に独立しているという仮定の下で、遺伝子ネットワーク推定において高い信頼レベルを有するかなり多くの辺の領域が発見できることを示す。しかし、この仮定は、同一遺伝子ｇに結合している全ての辺の信頼レベルが、ｇの発現寸法に依存するという事実を説明しない。
【０１８０】
しかし、たとえ各々の単一辺がある閾値を超えるカウント数を有していても、これらの辺で構成された部分ネットワークは必ずしも手元のネットワークに頻繁に出現する必要はない。従って、可能性の高い辺で構成された部分ネットワーク自体は可能性が低い場合がある。これは以下の問題、すなわち遺伝子ネットワークのモチーフ抽出問題と呼ぶ問題をもたらす。
所与のグラフＮ１＝（Ｇ，Ｅ１），・・・，Ｎｍ＝（Ｇ，Ｅｍ）、およびｋ∈ＩＲに対して、集合Ｍ⊆Ｇ×Ｇを発見する。ここに、｜Ｍ｜＝ｋはＭを含む、すなわちＭ⊆ＥｉであるグラフＮｉの個数を最大にする。
【０１８１】
モチーフ抽出問題は、データマイニングから最大出現頻度の項目の集合を発見するよく知られた問題を同値である。バランスの取れた完全二分部分グラフを発見する問題（６）は両方の問題に分解することができ、従ってＮＰ困難である。しかし、この問題は、現実的な事例については以下に述べる全件探索戦略を用いて解くことができる。
【０１８２】
３．４．２  モチーフ抽出アルゴリズム（アルゴリズム３．２）
所与の遺伝子発現測定に含まれる遺伝子ネットワークに関する部分情報を抽出すべく、３個のパラメータｎ、ｃ、ｋ、∈、ＩＮを用いる以下の戦略を提案する。
【０１８３】
【表１０】
2007519988.tif 000026
【０１８４】
本アルゴリズムは、少なくともｃ回出現するモチーフＭの全ての辺自身が少なくともｃ回出現する筈であるという事実を利用する。従って、ｃ回より多く出現するモチーフを構成する候補辺の数は、実際には実質的に減らすことができる。この全件探索戦略の実行時間は、閾値ｃを小さく設定し過ぎた場合実行不可能になるが、信頼性の高い遺伝子ネットワークのモチーフに関心があるため、今回の計算には実際的な制約を課すものではない。
【０１８５】
３．５．  モチーフ抽出結果
第３節と同様に、本発明の推定の有望な対象として細菌を選択した。最初にモチーフ抽出の信頼性を評価し、次いで、枯草菌並びに大腸菌のトリプトファン代謝に関与する遺伝子のグループ本方法を適用する。
【０１８６】
３．５．１  細菌遺伝子ネットワークにおけるモチーフの発見
モチーフ抽出アルゴリズムの予測能力を評価して、最尤遺伝子ネットワーク・モデルだけを用いる予測能力と比較すべく、枯草菌の破壊株データ集合に適用した。第３．３節で述べたように、無作為に８個の遺伝子を含む２５個の対象遺伝子ネットワークＮｉを選択し、最も可能性のある１００個のネットワークを列挙して、閾値（本アルゴリズムのパラメータｃ）として９０を用いて、２個以上の辺を有するモチーフを抽出した。この計算は、１．９ＧＨｚの単一ＣＰＵを用いて約３時間行なわれた。スーパーコンピュータを用いたならば、最高約３０個の遺伝子の遺伝子ネットワークの列挙を行なうことができる点に注目する。
【０１８７】
次いで最適ネットワーク・モデルの一辺、および最も多くの辺を有するモチーフのうちスコアが最も高いモチーフの一辺を無作為に選択した。次いで、ＤＢＴＢＳデータを用いて、両方の辺の正しさを調べて、既知の調整関係を有する辺の数をＮｉにおける可能な辺の総数で除算することによりネットワークから単一の辺を無作為且つ正しく推測する確率ｐｉを計算した。第３節の結果によれば、１エントリおよび０．５エントリを正しいものと判断した。ｎ個のネットワークの中で少なくともｋ個の単一辺を正しく推測する確率Ｐ（ｎ，ｋ）の上限を、以下の不等式を用いて計算した。
【数１５】
2007519988.tif 000027
ここに、ｐはｍａｘｉ＝１２０ｐｉを表わす。この計算の結果を表３．３に示す。モチーフ抽出方法の結果が、有向モチーフの場合には最適遺伝子ネットワークの結果よりはるかに有意であり、無向モチーフの場合には同程度に有意であることが観察された。
【０１８８】
これは、最適ネットワーク・モデルが空グラフであったケースは無いものの、いくつかの対象ネットワークについては閾値を超えるモチーフが無かったという事実による。この対象ネットワークの無作為選択は、データに表われないいくつかのネットワークを選択する可能性があるため、どの遺伝子の集合に対してもスコアの高いモチーフを期待することはできない。従って、調整関係を予測することはある種の対象ネットワークには適していないかもしれない。評価のために任意の辺を選択したため、今回の結果も全体として予測に対して成り立つ点に注目する。
【０１８９】
【表１１】
2007519988.tif 000028
【０１９０】
３．５．２  関連する種の比較遺伝子ネットワークの比較
上述の本方法は、例えば大腸菌（トリプトファン代謝の実験）および枯草菌の両方の経時的実験並びに研究が進んでいるトリプトファン・ネットワークに関与することが知られている遺伝子から得られたマイクロアレイデータに上述のアルゴリズムを適用することにより、グラム陽性桿菌様枯草菌およびグラム陰性桿菌様大腸菌等の関連する種の遺伝子調整ネットワークの相違および類似点の解明に役立つ。閾値として５０を使用して、有向モチーフに基づく１００個の列挙された最適ネットワーク・モデルを抽出した。出力ファイルからデータ集合内で発見された最高ヒット数を有する最大モチーフを示すグラフを導出し、遺伝子をノードで表わし、各辺にはその重みをラベル付けしている。
【０１９１】
大腸菌のデータ集合から得られた最大モチーフは５４回発見されて重みが７９〜１００の範囲にある１３個の辺を含み（図４）、枯草菌データの最大モチーフは６１回発見されて重みが７７〜１００の範囲にある１５個の辺を含む（図５）。
【０１９２】
ＴｒｐＡ、ｔｒｐＢ、ｔｒｐＣ、ｔｒｐＤ、ｔｒｐＥは両方のモチーフ内で重み付けの高い辺により結合されて、この調整ネットワークの中核を形成しており、これらが両方の種のｔｒｐオペロンにおいて互いに緊密に配置されている事実に対応する。ｔｒｐオペロン内における位置の順序でさえも当該グラフで部分的に認識できる。しかし、大腸菌では７個の遺伝子がｔｒｐ生合成経路内の酵素をコード化するのに対し、枯草菌の場合は５個である。２個の余分な遺伝子ｔｒｐＦ、ｐａｂＡもまた、導出されたモチーフに含まれている。ｔｒｐＦはｔｒｐπおよびｔｒｐＤに結合していて、ｔｒｐオペロ上の近似的な位置に対応している。ＰａｂＡは、枯草菌内ではｔｒｐＧとも呼ばれており、ｔｒｐオペロン内にはなく、葉酸オペロンｃｉｔｅｈｅｎｎｅｒ９３，ｇｏｌｌｎｉｃｋ０２内に位置するものの、ｔｒｐＡおよびｔｒｐＢに緊密に結合していることがわかった。この緊密な結合は、トリプトファン生合成経路における緊密な機能的関係を示唆している。大腸菌には、ｐａｂＡ遺伝子も存在するが、結果はトリプトファン・ネットワーク内のどの遺伝子とも結合を示しておらず、これが関与することが知られていない事実と符合する。
【０１９３】
一方、ｍｔｒＢは枯草菌のグラフ内でｔｒｐＥおよびｔｒｐＦに緊密に結合されている。このことは、その遺伝子産物であるトリプトファンＲＮＡ結合減衰タンパク質（ＴＲＡＰ）が、枯草菌におけるトリプトファン生合成の主要調整剤に含まれることが判明している事実により説明できる（２４）。興味深いことに、大腸菌におけるトリプトファン固有輸送タンパク質をコード化するｍｔｒはまた、グラフ内のｔｒｐＢおよびｔｒｐＥ等のトリプトファン生合成の中核遺伝子に関連付けられている。ネットワークのモチーフの構造は、大腸菌におけるｍｔｒの機能が枯草菌におけるｍｔｒＢの機能と類似している可能性を示唆する。
【０１９４】
昨今、ＴＲＡＰ阻害タンパク質（反ＴＲＡＰ、ＡＴ）が発見され、ｙｃｚＡの遺伝子産物として確認された（３４）。しかし、得られた結果から、ｙｃｚＡと調査された他の遺伝子と間には関連が見られなかった。しかし、これはｙｃｚＡ−ｍＲＮＡレベルが低過ぎるためかもしれない。あるいは、ＡＴは、ｙｃｚＡ−ｍＲＮＡからの新たな翻訳を一切必要とせずに、細胞全体に存在している可能性もある。これは、ＡＴが単独で作用せず、細胞内のトリプトファン・レベルに依存し、Ｔｒｐ活性化ＴＲＡＰだけと結合するため、あり得る（３４）。これまで大腸菌にはｙｃｚＡの相同物が発見されていないが、ＢＬＡＳＴＰの探索を通じてＡＴのアミノ酸配列がシャペロンｄｎａＪの多システイン領域のそれと極めて高い類似性を示すこと明らかになった（３４）。従って、ｄｎａＪをこの計算に含めたが、結果は大腸菌におけるトリプトファン・ネットワークとの関連に示唆し、さらに枯草菌においてより顕著であった。その理由は恐らく、シャペロンがあらゆる種類の経路と相互作用できるからであろう。しかし、これはｙｃｚＡ−ｍＲＮＡとの交差交配によっても生じるかもしれない。あるいは、ｄｎａＪがｔｒｐネットワークにおいて実際に役割を果たしている。タンパク質レベルでの類似性のため、ＴＲＡＰ調整剤としても機能する可能性がある。これは、枯草菌のグラフにおけるｄｎａＪとｍｔｒＢの強い結合を説明する。
【０１９５】
トリプトファニル−ｔＲＮＡシンセターゼをコード化するＴｒｐＳは、ｔＲＮＡＴｒｐが、両方の細菌におけるトリプトファン生合成の調整において重要だが異なる役割を果たすことがわかったため、上述この解析に含まれている。各々のネットワーク・モチーフの構造は、実際に異なる。大腸菌のグラフは、ｔｒｐＳがｍｔｒおよびｗｒｂＡに結合されている様子を示す。上述のように、ｍｔｒはトリプトファン固有輸送タンパク質をコード化し（２３）、ｔｒｐＬによりコード化されてグラフ内でｔｒｐＬからｍｔｒへの辺に対応するリーダーペプチドで滞留しているｔＲＮＡＴｒｐが豊富なことで確認できる転写終了に関与している可能性がある（３４）。ｗｒｂＡは、トリプトファン・リプレッサ結合タンパク質をコード化するため、当該辺はトリプトファン・リプレッサを活性化させるＴｒｐに対応する。しかし、遺伝子はトリプトファン・リプレッサをコード化しており、ｔｒｐＳとｔｒｐＲの間に直接的な関連は無い。このことは、アレイがｍＲＮＡ発現を測定し、トリプトファン・リプレッサ・タンパク質は転写のレベルではなく、Ｔｒｐが結合する際の配座変化を通じて作用するという事実により説明できる。しかし、グラフはｔｒｐＲのｔｒｐＤとの関連を示し、ｔｒｐＲがｔｒｐオペロン内には存在しないものの、トリプトファン生合成酵素と共にトリプトファン・リプレッサ・タンパク質の遺伝子産物を示唆する。枯草菌のグラフにおいて、ｔｒｐＳ発現は、ｐａｂＡおよびｄｎａＪに関連付けられている。ｄｎａＪがｙｃｚＡとして機能する場合、これはｔＲＮＡＴｒｐがＡＴ不活性化ＴＲＡＰの形成を伴うことが判明しているとの知見と整合している（３４）。
【０１９６】
しかし、これには制限がある。マイクロアレイ・データは、ノイズの影響を受けやすく、ｍＲＮＡレベルでしかモニタしない。大腸菌データ集合は僅か（１８個のアレイ）しかなかった。枯草菌については、トリプトファン・ネットワークの解析用に設計された実験からの実験データが発明者の手元に無かった。しかし、大腸菌実験はトリプトファン代謝を調べる目的で設計されており、枯草菌データ集合は各種の培地、すなわち各種の栄養条件下で成長している細菌から得られた。従って、トリプトファン遺伝子調整ネットワークに影響を及ぼすような培地でのトリプトファン・レベルの差異が予想される。このように、上述の計算方法を評価するには選択された遺伝子の集合が最適且つ適当でないかもしれないが、これらのデータを与えられてトリプトファン・ネットワークは最適な対象であった。
【０１９７】
３．６．結論
相当巨大な遺伝子ネットワークに対して最適および準最適な遺伝子ネットワーク・モデルを列挙する理論上の基礎を与え、最適モデルを既存知識と総合的に比較した結果を提示して、最適ネットワーク・モデルの信頼性が高い部分を抽出する方法を紹介し、その方法を適用した。
【０１９８】
最先端のスコアに基づく最適推定から信頼性が最も高いモチーフを抽出することは、それ自体が重要であり、また一般的な遺伝子ネットワーク・モチーフの解析（ネットワーク・モチーフの意味で）、および関連する種同士での遺伝子ネットワークの比較等への道を開いた。
【０１９９】
本研究で記述したアルゴリズム的方法論は、特定のスコアリング方式または特定種類の遺伝子発現測定に依存しない。従って、これらの技術は、対応関係ｓ：Ｇ×２Ｇ！ＩＲを有するスコアが与えられているあらゆる状況において一般に適用可能である。これはベイジアンネットワークの枠組みの範囲内のあらゆるスコアを始め、他の大多数のスコア関数に対しても当てはまる。これは、シーケンス情報（３３）またはタンパク質相互作用データ（１４、２２）等の事前知識を取り込んだ新規スコアの研究が進行中であるため、遺伝子ネットワーク推定技術における重要な特性である。
【０２００】
本発明のアルゴリズムの計算面での限界を超える大きさの遺伝子ネットワークに対しては、生物学的に有意味なネットワークを受容して、選択された部分空間内で最適ネットワーク・モデルを発見できそうな探索空間の部分に制限すべく、（２４）に記載されているような技術を適用することができる。この場合でも、遺伝子ネットワークのモチーフを発見する本発明の方法が小規模遺伝子ネットワークに限定されないため、本発明の列挙アルゴリズムを同様に適用されることができる点に注目している。
【０２０１】
第３．３節で行なったように、利用できる知識を用いて遺伝子ネットワーク推定の精度を厳密に評価することは、遺伝子ネットワーク推定方法の性能を比較する基準を策定する有望な方法である。
【０２０２】
本明細書で引用した各々の参考文献は、その全体を以下に掲載している。
【０２０３】
参考文献
1. T. Chen, H. L. He, G. M. Church. Modeling gene expression with differential equations. Pacific Symposium on Biocomputing, 4: 29--40, 1999. 
2. D. M. Chickering. Learning Bayesian networks is NP-complete. In D. Fisher and H.-J. Lenz, editors, Learningfrom Data :Artificial Intelligence and Statistics V, Springer-Verlag, 1996. 
3. G. F. Cooper, E. Herskovits. A Bayesian method for the induction of probabilistic networks from data. Machine Learning, 9: 309--347, 1992. 
4. N. Friedman, M. Goldszmidt. Learning Bayesian networks with local structure. Jordan, M. I. (ed. ), Kluwer Academic Publishers, pp. 421--459, 1998. 
5.N. Friedman, M. Linial,1. Nachman, D. Pe'er. Using Bayesian networks to analyze expression data. Journal of Computational Biology, 7:601--620, 2000. 
6. M. R. Garey, D. S. Johnson. Computers and intractability, W. H. Freeman and Company, 1979. 
7. P. Gollnick, P. Babitzke, E. Merino, C. Yanofsky. Bacillus subtilis and itsclosest relatives : from genes to cells. A. L. Sonenshein, J. A. Hoch, R. Losick. (Eds.), American Society for Microbiology, 2002. 
8. A. J. Hartemink, D. K. Gifford, T. S. Jaakkola, R. A. Young. Using graphical models and genomic expression data to statistically validate models of genetic regulatory networks. Pacific Symposium on Biocomputing, 6: 422--433, 2001. 
9. A. J. Hartemink, D. K. Gifford, T. S. Jaakkola, R. A. Young. Combining location and expression data for principled discovery of genetic regulatory network models. Pacific Symposium on Biocomputing, 7: 437-449, 2002. 
10. L. H. Hartwell, J. J. Hopfield, S. Leibler, A. W. Murray. From molecular to modular cell biology. Nature, 402: C47--C52, 1999. 
11. D. Henner, C. Yanofsky. Bacillus subtilis and other gram-positive bacteria : biochemistry, physiology and molecular genetics. A. L. Sonenshein, J. A. Hoch,R Losick. (Eds. ), American Society for Microbiology, pp.269-280, 1993. 
12. M. J. L. de Hoon, S. Ott, S. Imoto, S. Miyano. Validation of noisy dynamical system models of gene regulation inferred from time-course gene expression data at arbitrary time intervals. European Conference on Computational Biology, 2003. 
13. S. Imoto, T. Goto, S. Miyano. Estimation of genetic networks and functional structures between genes by using Bayesian networks and nonparametric regression. Pacific Symposium on Biocomputing, 7: 175--186, 2002. 
14. S. Imoto, T. Higuchi, T. Goto, K. Tashiro, S. Kuhara, S. Miyano. Combining microarrays and biological knowledge for estimating gene networks via Bayesian networks. Proc. 2ndComputational Systems Bioinformatics, 104--113, 2003. 
15. S. Imoto, S. Kim, T. Goto, S. Aburatani, K. Tashiro, S. Kuhara, S. Miyano. Bayesian network and nonparametric heteroscedastic regression for nonlinear modeling of geneticnetwork. Journal ofBioinformatics and Computational Biology, in press, 2003; U. S. Patent Application No : 10/259,723. 
16. T. Ishii, K. Yoshida, G. Terai, Y. Fujita, K. Nakai. DBTBS: A database of Bacillus subtilis promoters and transcription factors. Nucleic Acids Research, 29: 278--280, 2001. 
17. A. B. Khodursky, B. J. Peter, M. B.Schmid, J. DeRisi, D. Botstein, P. O. Brown. Analysisof topoisomerase function in bacterial replication fork movement: Use of DNA microarrays. Proceedings of the National Academyof Sciences, 97: 9419--9424,2000. 
18. A. B. Khodursky, B. J. Peter, N. R. Cozzarelli, D. Botstein, P. O. Brown, C. Yanofsky. DNA microarray analysis of gene expression in response to physiological and genetic changes that affect tryptophan metabolism inEscherichia coli. Proceedings of the National Academyof Sciences, 97 :12170--12175, 2000. 
19. J. Courcelle, A. Khodursky, B. Peter, P. O. Brown, P. C. Hanawalt. Comparative gene expression profiles following LJV exposure in wild-type and SOS-deficient Escherichia coli. Genetics, 158: 41--64, 2001. 
20. Y. Makita, M. Nakao, N. Ogasawara, K. Nakai. DBTBS: Database of transcriptional regulation in Bacillus subtilis and its contribution to comparative genomics.submittedfor publication. 
21. R. Milo, S. Shen-Orr, S. Itzkovitz, N. Kashtan, D. Chklovskii, U. Alon. Network motifs: Simple building blocks of complex networks. Science, 298: 824--827,2002. 
22. N. Nariai, S. Kim, S. Imoto, S. Miyano. Using protein-protein interactions for refining gene networks estimated from microarray data by Bayesian networks. Pacific Symposium on Biocomputing, in press, 2004. 
23. I. M. Ong, J. D. Glasner, D. Page. Modelling regulatory pathways in E. coli from time series expression profiles.Bioinformatics, 18: 241--248, 2002. 
24. S. Ott, S. Miyano, Finding optimal gene networks using biological constraints, submittedforpublication. 
25. S. Ott, S. Imoto, S. Miyano. Finding optimal models for small gene networks. Pacific Symposium on Biocomputing, in press, 2004. 
26. D.Pe'er, A. Regev, G. Elidan, N. Friedman. Inferring subnetworks from perturbedexpressionprofiles. Bioinformatics, 17: 215--224,2001. 
27. R. W. Robinson. Counting labeled acyclic digraphs. New Directions in the Theory of Graphs, pp. 239-273,1973. 
28. H. Salgado, A. Santos-Zavaleta, S. Gama-Castro, D.Millan-Zarate, E. D"EDaz-Peredo, F.Sanchez-Solano, E.Perez-Rueda, C. Bonavides-Mart"EDnez, J. Collado-Vides. RegulonDB (version 3.2) : transcriptional regulation and operon organization inEscherichia coli K-12. Nucleic Acids Research, 29: 72--74, 2001. 
29. S. S. Shen-Orr, R. Milo, S. Mangan, U. Alon. Network motifs in the transcriptional regulation network ofEscherichia coli. Nature Genetics, 31: 64--68,2002. 
30. V. A. Smith, E. D. Jarvis, A. J. Hartemink. Evaluating functional network inference using simulations of complex biological systems. Bioinformatics, 18 : 216--224,2002. 
31. E. P. van Someren, L. F. A. Wessels, E. Backer, M. J. T. Reinders. Genetic network modeling. Pharmacogenomics, 3 (4):507--525, 2002. 
32. A. L. Sonenshein, J. A. Hoch, R. Losick.Bacillus subtilis and its closest relatives : from genes to cells. ASM Press, Washington, D. C. , 2001. 
33. Y. Tamada, S. Kim, H. Bannai, S. Imoto, K. Tashiro, S. Kuhara, S. Miyano. Estimating gene networks from gene expression data by combining Bayesian network model with promoter element detection.Bioinformatics, in press, 2003. 
34. A. Valbuzzi, C. Yanofsky. Inhibition of the B. subtilis regulatory protein TRAP by the TRAP-inhibitory protein, AT. Science, 293: 2057--2059, 2001. 
35. http://www. ncbi. nlm. nih. gov/geo/
【図面の簡単な説明】
【０２０４】
図面の簡単な説明
【図１】例１のアルゴリズム３により予測される遺伝子ネットワークを示す。
【図２ａ】枯草菌に関して得られた経時的データを示す。
【図２ｂ】枯草菌に関して得られた破壊株データを示す。
【図３】大腸菌から得られたデータを示す。
【図４】大腸菌に関して抽出された遺伝子ネットワーク関係モチーフを示す。
【図５】枯草菌に関して抽出された遺伝子ネットワーク関係モチーフを示す。
【図１】
2007519988.tif 000031
【図２ａ】
2007519988.tif 000032
【図２ｂ】
2007519988.tif 000033
【図３】
2007519988.tif 000034
【図４】
2007519988.tif 000035
【図５】
2007519988.tif 000036
【国際調査報告】
2007519988.tif 000037
2007519988.tif 000038