JPWO2002038749A1 - 遺伝子ネットワークの推定方法、推定システム及び記録媒体 - Google Patents
遺伝子ネットワークの推定方法、推定システム及び記録媒体 Download PDFInfo
- Publication number
- JPWO2002038749A1 JPWO2002038749A1 JP2002542065A JP2002542065A JPWO2002038749A1 JP WO2002038749 A1 JPWO2002038749 A1 JP WO2002038749A1 JP 2002542065 A JP2002542065 A JP 2002542065A JP 2002542065 A JP2002542065 A JP 2002542065A JP WO2002038749 A1 JPWO2002038749 A1 JP WO2002038749A1
- Authority
- JP
- Japan
- Prior art keywords
- gene
- network
- expression level
- genes
- estimation
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
Images
Classifications
-
- C—CHEMISTRY; METALLURGY
- C12—BIOCHEMISTRY; BEER; SPIRITS; WINE; VINEGAR; MICROBIOLOGY; ENZYMOLOGY; MUTATION OR GENETIC ENGINEERING
- C12N—MICROORGANISMS OR ENZYMES; COMPOSITIONS THEREOF; PROPAGATING, PRESERVING, OR MAINTAINING MICROORGANISMS; MUTATION OR GENETIC ENGINEERING; CULTURE MEDIA
- C12N15/00—Mutation or genetic engineering; DNA or RNA concerning genetic engineering, vectors, e.g. plasmids, or their isolation, preparation or purification; Use of hosts therefor
-
- G—PHYSICS
- G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
- G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
- G16B25/00—ICT specially adapted for hybridisation; ICT specially adapted for gene or protein expression
- G16B25/10—Gene or protein expression profiling; Expression-ratio estimation or normalisation
-
- G—PHYSICS
- G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
- G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
- G16B5/00—ICT specially adapted for modelling or simulations in systems biology, e.g. gene-regulatory networks, protein interaction networks or metabolic networks
-
- G—PHYSICS
- G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
- G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
- G16B25/00—ICT specially adapted for hybridisation; ICT specially adapted for gene or protein expression
Landscapes
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Engineering & Computer Science (AREA)
- Genetics & Genomics (AREA)
- Physics & Mathematics (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Biotechnology (AREA)
- General Health & Medical Sciences (AREA)
- Biophysics (AREA)
- Molecular Biology (AREA)
- Medical Informatics (AREA)
- Zoology (AREA)
- Bioinformatics & Computational Biology (AREA)
- Spectroscopy & Molecular Physics (AREA)
- Theoretical Computer Science (AREA)
- Organic Chemistry (AREA)
- Biomedical Technology (AREA)
- General Engineering & Computer Science (AREA)
- Chemical & Material Sciences (AREA)
- Evolutionary Biology (AREA)
- Wood Science & Technology (AREA)
- Microbiology (AREA)
- Plant Pathology (AREA)
- Biochemistry (AREA)
- Physiology (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
- Information Retrieval, Db Structures And Fs Structures Therefor (AREA)
Abstract
複数の遺伝子のうち一の遺伝子の発現量を２つの条件下でそれぞれ発現させたときの他の遺伝子の発現量をそれぞれ検出し、得られる検出値の差を求め、当該差を指標として前記一の遺伝子と他の遺伝子との因果関係を導くことを特徴とする遺伝子ネットワークの推定方法。
Description
技術分野
本発明は、遺伝子ネットワークの推定方法、推定システム及び記録媒体に関する。
背景技術
最近の分子生物学分野の進歩により、得られる遺伝子情報量は膨大である。従って、次々と明らかになり急激に増加する配列データ数又は膨大な遺伝子発現データから情報を取り出すには、コンピュータを利用して努力する必要がある。これまでは、ホモロジー検索、タンパク質分類、遺伝子集積など各種のコンピュータツールを開発することに関心が注がれてきた。
こうした試みの中で、遺伝子発現データから遺伝子調節ネットワークを推論する方法として、数例ではあるがいくつかの研究が知られている。遺伝子発現データは、時系列データ（時系列データとは、対象遺伝子群の遺伝子発現量を経時的に測定したもの）又は定常状態データ（定常状態データとは、対象遺伝子群の遺伝子発現量を複数の異なる実験条件下（例えば、遺伝子変異や薬剤投与）で測定したもの）のどちらかの形で得られる。
時系列を分析する方法の場合は、種々の方法、例えば情報理論、遺伝アルゴリズム又はシミュレーテッドアニーリング法を用いてネットワークを推定することができる（Ｌｉａｎｇ，Ｓ．ｅｔ ａｌ．，Ｐｒｏｃ．Ｐａｃｉｆｉｃ Ｓｙｍｐ．Ｂｉｏｃｏｍｐｕｔｉｎｇ ’９８，Ｗｏｒｌｄ Ｓｃｉｅｎｔｉｆｉｃ，１８−２９，１９９８．；Ｍｏｒｏｈａｓｈｉ，Ｍ．ａｎｄ Ｋｉｔａｎｏ，Ｈ．，Ｐｒｏｃ．５ｔｈ Ｅｕｒｏ．Ｃｏｎｆ．Ａｒｔｉｆｉｃｉａｌ Ｌｉｆｅ，Ｓｐｒｉｎｇｅｒ，４７７−４８６，１９９９．；Ｍｊｏｌｓｎｅｓｓ，Ｅ．，ｅｔ ａｌ．，Ｔｅｃｈ．Ｒｅｐｔ．ＪＰＬ−ＩＣＴＲ−９９−４，Ｊｅｔ Ｐｒｏｐｕｌｓｉｏｎ Ｌａｂ．，ＮＡＳＡ，１９９９．）。しかし、時系列アプローチでは、試験結果を非常に短いインターバルで取得し、かつ試験的ノイズのない状態で取得することが必要とされる。これは現在の技術では非常に困難である。
一方、定常状態を分析する方法は、すでにいくつかの方法が提唱されている。この定常状態データは、特定の遺伝子活性を変異させて、例えば遺伝子を欠失させたりあるいは過剰発現させることにより得ることができる。欠失は、酵母ゲノム欠失コンソーシアム（Ｙｅａｓｔ Ｇｅｎｏｍｅ Ｄｅｌｅｔｉｏｎ Ｃｏｎｓｏｒｔｉｕｍ）等で現在大規模に行われており、これによって各種遺伝子欠失型の発現プロファイルを近い将来容易に入手することができる（Ｗｉｎｚｅｌｅｒ，Ｅ．Ａ．ｅｔ ａｌ．，Ｓｃｉｅｎｃｅ，２８５（５４２９）：９０１−９０６，１９９９．）。
Ａｋｕｔｓｕらは、遺伝子ネットワークがブーリアンネットワークモデルで表現された場合（各遺伝子の遺伝子発現が活性、不活性の２状態で表現され、それら遺伝子間の関係がブール代数によって表現された場合）に、必要となる試験数の上限と下限を計算した（Ａｋｕｔｓｕ，Ｔ．ｅｔ ａｌ．，Ｐｒｏｃ．９ｔｈ ＡＣＭ−ＳＩＡＭ Ｓｙｍｐ．Ｄｉｓｃｒｅｔｅ Ａｌｇｏｒｉｔｈｍｓ，６９５−７０２，１９９８．）。さらに最近、Ｉｄｅｋｅｒらは、Ｐｒｅｄｉｃｔｏｒ法と呼ばれる推定方法を提唱した（Ｉｄｅｋｅｒ，Ｔ．Ｅ．ｅｔ ａｌ．，Ｐｒｏｃ．Ｐａｃｉｆｉｃ Ｓｙｍｐ．Ｂｉｏｃｏｍｐｕｔｉｎｇ ２０００，Ｗｏｒｌｄ Ｓｃｉｅｎｔｉｆｉｃ，３０５−３１６，２０００．）。この方法は、組み合わせに関する最適化法を用いて候補ネットワーク（発現データと一致するブーリアンネットワークモデル）を提供するものである。
しかし、これらの方法は、遺伝子発現レベルを２値として遺伝子ネットワークを現すものである。従って、試験データは通常は連続値（多値）を有するにもかかわらず、この方法を適用するにはデータを２値に補正する必要がある。こうした補正によって、ネットワークを推定するのに必要な情報が失われる可能性がある。例えば２値を利用して推定をすると、３つの状態（例えば野生型、欠失、過剰発現）の遺伝子発現レベルを正確に表すことが困難であり、調節関係が正確に反映されない結果となる。
発明の開示
本発明は、遺伝子ネットワークの推定方法、推定システム及び記録媒体を提供することを目的とする。
本発明者は、上記課題を解決するため鋭意研究を行った結果、突然変異を誘発したときに得られる遺伝子発現プロファイル（検出値）を用いて、遺伝子発現量を単純化（二値化）することなく遺伝子ネットワークの推定を行うことに成功し、本発明を完成するに至った。
すなわち、本発明は、複数の遺伝子のうち一の遺伝子の発現量を２つの条件下でそれぞれ発現させたときの他の遺伝子の発現量をそれぞれ検出し、得られる検出値の差を求め、当該差を指標として前記一の遺伝子と他の遺伝子との因果関係を導くことを特徴とする遺伝子ネットワークの推定方法である。
さらに、本発明は、複数の遺伝子のうち一の遺伝子を変異させたときの他の遺伝子の発現レベルと、当該一の遺伝子を変異させないときの当該他の遺伝子の発現レベルとをそれぞれ検出し、得られる検出値の差を求め、当該差を指標として前記一の遺伝子と他の遺伝子との因果関係を導くことを特徴とする遺伝子ネットワークの推定方法である。上記遺伝子ネットワークとしては、単細胞生物内のもの又は多細胞生物内のものが挙げられる。遺伝子の変異としては遺伝子の破壊又は過剰発現が挙げられる。上記推定方法においては、得られた因果関係の矛盾の有無を検定することが好ましい。
さらに、本発明は、複数の遺伝子のうち一の遺伝子の発現量を２つの条件下で発現させたときの他の遺伝子の発現量をそれぞれ検出する手段と、得られる検出値の差を求める手段と、当該差を指標として前記一の遺伝子と他の遺伝子との因果関係を導く手段とを含んでなる、遺伝子ネットワークの推定システムである。
さらに、本発明は、複数の遺伝子のうち一の遺伝子を変異させたときの他の遺伝子の発現レベルと、当該一の遺伝子を変異させないときの当該他の遺伝子の発現レベルとをそれぞれ検出する手段と、得られる検出値の差を求める手段と、当該差を指標として前記一の遺伝子と他の遺伝子との因果関係を導く手段とを含んでなる、遺伝子ネットワーク推定システムである。本発明の推定システムにおいては、得られた因果関係の矛盾の有無を検定する手段を含めることができる。
さらに、本発明は、複数の遺伝子のうち一の遺伝子を変異させたときの他の遺伝子の発現レベルと、当該一の遺伝子を変異させないときの当該他の遺伝子の発現レベルとをそれぞれ検出する手段と、得られる検出値の差を求める手段と、当該差を指標として前記一の遺伝子と他の遺伝子との因果関係を導く手段とを含んでなる、コンピュータを遺伝子ネットワーク推定システムとして機能させるためのプログラムを記録したコンピュータ読み取り可能な記録媒体である。
以下、本発明を詳細に説明する。本明細書は、本願の優先権の基礎である日本国特許出願２０００−３４５９８２号の明細書及び／又は図面に記載される内容を包含する。
本発明は、複数の遺伝子のうちの一の遺伝子（遺伝子ａ）の発現量を、２つの条件下（条件１、条件２）でそれぞれ発現させたときに、他の遺伝子（遺伝子ｂ）の発現量が、条件１のときと条件２のときにおいてどのように変化するかに着目し、それぞれの条件における遺伝子ｂの発現量をそれぞれ調べ（検出し）、その調べた値（検出値）の差を求め、当該差を指標として一の遺伝子ａと他の遺伝子ｂとの因果関係を導くことを特徴とするものである。
例えば、遺伝子ａの野生型の場合（条件１）、及び当該遺伝子ａを変異させた場合（条件２）の２条件を設定し、これらの２条件下における当該他の遺伝子ｂの発現量の差を検出することにより、遺伝子ａと遺伝子ｂとの間の因果関係を求め、これらの遺伝子間に生ずる遺伝子ネットワークを推定する。
本発明においては、ＯＮとＯＦＦとの関係のように二値を扱うのではなく、変化する値の全部を連続値（多値）として扱うことを特徴とするため、二値よりも更に高精度にネットワークを推定することが可能である。
１．遺伝子ネットワーク推定方法及び推定システム
本発明の推定システムは、
（ｉ）ある一の遺伝子の発現量を２条件に設定した場合に他の遺伝子の発現レベル（発現量）をそれぞれ検出する手段（「検出エンジン」ともいう）、
（ｉｉ）得られる検出値を比較してその差を求める手段（「比較エンジン」ともいう）、
（ｉｉｉ）当該差を指標として前記一の遺伝子と他の遺伝子との因果関係を求める手段（「因果関係作成エンジン」ともいう）、
（ｉｖ）当該因果関係の中から冗長な因果関係を除去する手段（「冗長関係除去エンジン」ともいう）
を含む装置により構成される。
ここで、本発明の推定システムの構成例を示すブロック図を示す（図１）。
図１に示す推定システムは、ＣＰＵ１０１ ＲＯＭ１０２、ＲＡＭ１０３、入力部１０４、送信／受信部１０５、出力部１０６、ハードディスクドライブ（ＨＤＤ）１０７及びＣＤ−ＲＯＭドライブ１０８を備える。
ＣＰＵ１０１は、ＲＯＭ１０２、ＲＡＭ１０３又はＨＤＤ１０７に記憶されているプログラムに従って、遺伝子ネットワークシステム全体を制御し、後述するネットワーク推定処理を実行する。ＲＯＭ１０２は、遺伝子ネットワーク推定システムの動作に必要な処理を命令するプログラム等を格納する。ＲＡＭ１０３は、ネットワーク推定処理を実行する上で必要なデータを一時的に格納する。入力部１０４は、キーボードやマウス等であり、ネットワーク推定処理を実行する上で必要な条件を入力するとき等に操作される。送信／受信部１０５は、ＣＰＵ１０１の命令に基づいて、通信回線を介してパブリックデータベース１１０等との間でデータの送受信処理を実行する。出力部１０６は、遺伝子の発現量、入力部１０４から入力された各種条件、遺伝子の塩基配列、ネットワークの推定結果等を、ＣＰＵ１０１からの命令に基づいて表示処理を実行する。なお、出力部１０６としては、コンピュータのディスプレイ、又はプリンターなどが例示される。ＨＤＤ１０７は、遺伝子ネットワーク推定プログラム、遺伝子発現量、塩基配列等を格納し、ＣＰＵ１０１の命令に基づいて、格納しているプログラム又はデータ等を読み出し、例えばＲＡＭ１０３に格納する。ＣＤ−ＲＯＭドライブ１０８は、ＣＰＵ１０１の指示に基づいて、ＣＤ−ＲＯＭ１０９に格納されている遺伝子ネットワーク推定プログラム又は発現量等から、プログラム又はデータ等を読み出し、例えばＲＡＭ１０３に格納する。
ＣＰＵ１０１は、入力部などから受け取ったデータを出力部１０６に供給するとともに、データベースから受け取ったデータに基づいて遺伝子ネットワークの推定処理を実行する。
ここで、遺伝子ネットワークを推定するための処理プログラムは、下記のごとく、遺伝子発現量の検出エンジン、比較エンジン及び因果関係作成エンジンから構成される。
（１）遺伝子発現量の検出エンジン
本発明において、遺伝子発現の検出は、ＤＮＡマイクロアレイ、オリゴヌクレオチドチップ、ＲＴ−ＰＣＲ、遺伝子発現の連続分析、定常発現レベル、プロテオミクス等を使用又は適用することにより行うことができ、これらの結果を遺伝子発現データ（発現量）として使用することができる。また、蛋白のリン酸化の状態や染色体ＤＮＡのメチル化の状態等、遺伝子機能発現を調節するあらゆる状態変化も遺伝子発現として検出できる。上記データは、本発明では検出値という。検出値は、絶対値であっても相対値であってもよい。また、実験により得られた値のみならず、文献的に又は遺伝子データベース等から得られた値でもよい。例えば、ＤＮＡマイクロアレイから得られる発現量は、マイクロアレイ上に発せられる蛍光強度等を測定し、その蛍光強度を数値化することができる。
データベースの具体例としては、インターネットを通じて利用できるデータベースを例示でき、より具体的には、ＧｅｎＢａｎｋ、ＤＤＢＪ、ＥＭＢＬ、ＮＣＩ６０等を例示できる。
（２）比較エンジン
比較エンジンでは、ある遺伝子の遺伝子発現量を２条件に設定した場合のこれらの２条件下における対象となる遺伝子の発現レベルの差の情報、並びに上記遺伝子間の機能（例えば活性化又は阻害）に関する情報を蓄積する。上記遺伝子発現量の２条件設定には遺伝子の変異が含まれ、遺伝子変異は、塩基配列の欠失、置換又は付加によりもとの遺伝子の機能を失わせることを意味し、遺伝子の破壊も変異に含まれる。
遺伝子に欠失、置換、付加等の変異を導入するには、Ｋｕｎｋｅｌ法、Ｇａｐｐｅｄ ｄｕｐｌｅｘ法等の公知の手法又はこれに準ずる方法を採用することができる。例えば部位特異的突然変異誘発法を利用した変異導入用キット（例えばＭｕｔａｎ−Ｋ（ＴＡＫＡＲＡ社製）やＭｕｔａｎ−Ｇ（ＴＡＫＡＲＡ社製））、ＴＡＫＡＲＡ社のＬＡ ＰＣＲ ｉｎ ｖｉｔｒｏ Ｍｕｔａｇｅｎｅｓｉｓシリーズキットなどを用いて、あるいは、遺伝子の一部又は全部を切断することにより、変異の導入が行われる。あるいは、遺伝子を変異したときの発現量データをデータベースから入手することもできる。
（３）因果関係作成エンジン
因果関係作成エンジンは、比較エンジンにより得られた発現量の差に基づいて、遺伝子間の因果関係を導き、最終的にネットワークを推定する手段である。この手段は、例えば遺伝子ａと遺伝子ｂとの因果関係を調べる場合において、遺伝子ａを変異させたときの遺伝子ｂの発現量が、遺伝子ａを変異させないとき（野生型のとき）の遺伝子ｂの発現量と比較してどの程度変化したときに、その遺伝子ａは遺伝子ｂに対して活性化するのか、抑制するのか、あるいは無関係であるのか、という判断を実行する手段である。
本発明において、便宜上、遺伝子ネットワークを「有向グラフ」と呼ばれる図に置き換えて考えてみる（図２）。図２において、アルファベットは遺伝子を表し、矢印（→）は遺伝子ａが遺伝子ｂの発現を促進する制御を、Ｔ字型矢印（「┤」）は遺伝子ｂが遺伝子ｃの発現を抑える制御を意味するものとする。個々の矢印を、本発明では「エッジ」という。遺伝子ａと遺伝子ｂとの間におけるネットワーク（因果関係）を推定する場合は、まず、遺伝子ａの発現量を２条件に設定した場合のこれらの２条件における遺伝子ｂの発現量、例えば、遺伝子ａを何ら操作しない定常型又は野生型の状態で検出した遺伝子ｂの発現量と、遺伝子ａを変異（例えば破壊）したときの遺伝子ｂの発現量とを検出し、次いで、それぞれ得られたデータ（検出値）の差を求める。この差を指標として遺伝子間の因果関係を導く。例えば、差に変化がなければ遺伝子ａは遺伝子ｂを活性化していない、又は無関係であると判断することができ、差に一定の変化があれば遺伝子ａは遺伝子ｂを活性化又は抑制していると判断することができる。ここで、「一定の変化」とは、いわゆる閾値を超える変化を示すものであり、閾値は、ネットワークを推定する際に目的となる遺伝子に応じて適宜設定できる。但し、閾値は、誤差検定を行い、例えば有意差危険度０．０１以下となるような基準で設定し、採用することが好ましい。
遺伝子の発現量は、一般に数値として得られる場合が多いため、仮に遺伝子間の相互関係を図３のように表すことができる。図３において、ａｎ（ｎ＝０，１，２，３）は「ノード（ｎｏｄｅ）」と呼ばれるものであり、遺伝子ネットワークを推定したときに作成されるエッジとエッジとの間の点（交点）を意味する。ノードは、遺伝子名で表示することもできることから、ノードの数と遺伝子の数（ｎの数）とは一致する。ｘｎ（ｎ＝０，１，２，３）は対応するノード（遺伝子）の発現量を意味する。ｗｔは野生型を意味する。
ある遺伝子ａ０〜ａ３の発現量の結果が図３Ｂのように得られたものとする。ａｎ −（ｎ＝０，１，２，３）は各ノードにおける遺伝子が変異（破壊）されたことを意味する。ここで、遺伝子ａ０を破壊した場合（ａ０ −）を考えてみる。
まず、ａ０ −における遺伝子ａ１の発現量と、野生型におけるａ１遺伝子の発現量との差は不変であるため（３．７５０−３．７５０＝０）、ａ０はａ１とは無関係、すなわちａ０からａ１へのネットワークは無いものと推定することができ、ノード間のエッジは無い（図４Ａ）。次に、ａ０ −における遺伝子ａ２の発現量と、野生型における遺伝子ａ２の発現量との差は−０．１７であり（８．７６９−８．９３９＝−０．１７）、遺伝子ａ０を破壊するとａ２の発現量は減少する。従って、閾値を無視して考えれば遺伝子ａ０はａ２を活性化させると推定することができ、ノードａ０からａ２に向かうエッジを引く（つまり矢印を表示する）ことができる（図４Ｂ）。さらに、ａ０ −における遺伝子ａ３の発現量と、野生型における遺伝子ａ３の発現量との差は−０．０６７であり（０．０１１−０．０７８＝−０．０６７）、遺伝子ａ０を破壊するとａ３の発現量は減少する。従って、遺伝子ａ０はａ３を活性化させるものと推定することができ、ａ０からａ３に向かってエッジを引くことができる（図４Ｃ）。
遺伝子ａ１を破壊した場合、遺伝子ａ２を破壊した場合及び遺伝子ａ３を破壊した場合についても上記と同じ要領で処理することにより、４つの遺伝子間のネットワークを推定することができる（図４Ｄ）。特に、遺伝子ａ２を破壊した場合における遺伝子ａ３の発現量と野生型における遺伝子ａ３の発現量との差は５．３９８であり（５．４７６−０．０７８＝５．３９８）、遺伝子ａ２を破壊するとａ３の発現量が上昇する。従って、遺伝子ａ２は遺伝子ａ３を抑制していると推定することができ、ノード間のエッジは抑制を示す矢印（Ｔ字型矢印）となる（図４Ｄ）。
以上のネットワークの関係をまとめると、図３Ａのようになる。図３Ａにおいて、「＋」は活性化（→）を、「−」は抑制（┤）を意味する。
ある遺伝子ａを欠失又は過発現させたときの遺伝子ｂの発現レベルの増加及び減少との関係を表１に示す。この比較プロセスの計算量はＯ（ｎ２）である。
（４）冗長関係除去エンジン
図５に示すように、遺伝子ａは遺伝子ｂを活性化し、活性化された遺伝子ｂは遺伝子ｃを活性化する遺伝子ネットワークを考える。この場合において、遺伝子ａを不活性化（例えば破壊）すると、遺伝子ｂは活性化されず、その結果遺伝子ｃも活性化されなくなる。従って、実際には遺伝子ａが遺伝子ｃを活性化する能力はなくても、結果として遺伝子ａは遺伝子ｃを間接的に活性化する能力があると誤って認識される場合がある。このように、誤った認識により作成される因果関係を、本発明では「冗長な因果関係」という。本発明においては、そのような間接的なエッジ（図５、破線の矢印）を見つけることができる。間接的なエッジ及び直接的なエッジが存在する遺伝子調節ネットワークにおいて、本発明では、各遺伝子が他の遺伝子に間接的な効果を及ぼすかどうかを探り、そのネットワーク内の間接的エッジをチェックすることができる。この間接的効果は、エッジの道筋（ルート）に関与する負の調節数の偶奇性にのみ依存する（Ｔｈｉｅｆｆｒｙ，Ｄ．，ａｎｄ Ｔｈｏｍａｓ，Ｒ．，Ｐｒｏｃ．Ｐａｃｉｆｉｃ Ｓｙｍｐ．Ｂｉｏ−ｃｏｍｐｕｔｉｎｇ ’９８，Ｗｏｒｌｄ Ｓｃｉｅｎｔｉｆｉｃ，７７−８８，１９９８．）。
本発明においては、冗長な因果関係、すなわち間接的な遺伝子調節を除去してネットワークを推定する。こうした冗長な因果関係は、グラフ理論（Ｇｒｏｓｓ，Ｊ．，ａｎｄ Ｙｅｌｌｅｎ，Ｊ．，ＣＲＣ Ｐｒｅｓｓ，１９９９．）の分野のＷａｒｓｈａｌｌのアルゴリズムを修正することにより見つけることができる。このような冗長な因果関係を除外する処理においては、冗長な因果関係を見つけるためのプログラムを実行する。
修正したＷａｒｓｈａｌｌのアルゴリズムを実行するためのフローチャートを図６に示す。
図６におけるステップは以下の通りである。すなわち、冗長な因果関係を含む遺伝子ネットワーク（冗長な因果関係が除去されていないネットワーク）を初期値Ｇ０とし、全ての遺伝子に１からｎ（遺伝子の総数）の通し番号をつける。次に遺伝子ｘから遺伝子ｙへの関係をマトリクスの要素（ａｘ，ａｙ）と表現し、活性（正）か抑制（負）の関係が存在する場合に（ａｘ，ａｙ）を１にセットする（存在しない場合には０にセットする）。また、その関係が抑制（負の制御）の関係の場合には、マトリクスの各要素が持つカウンタ（ａｘ，ａｙ）ｎｅｇを１にセットし、その他の場合には０をセットする（Ｓ１）。まず、負の制御の総数をカウントするカウンタｔｎを初期化する（Ｓ２）。次に図６のフローチャートに従って処理を進めていく。Ｓ３は、インデックスｉを０に初期化するという処理を行うステップである。Ｓ４は、インデックスｉをｉ＋１にセットするという処理を行うステップである。Ｓ５は、インデックスｊを０に初期化するという処理を行うステップである。Ｓ６は、インデックスｊをｊ＋１にセットするという処理を行うステップである。図６で（ａｊ，ａｉ）は遺伝子ｊから遺伝子ｉへの関係を意味しており、（ａｊ，ａｉ）が存在するという意味は、（ａｊ，ａｉ）が１であることを示しており、遺伝子ｊから遺伝子ｉへの関係が存在するということである。（ａｊ，ａｉ）が存在する場合（ｙｅｓ）は、Ｓ８においてインデックスｋを０に初期化する処理を行う（Ｓ７）。（ａｊ，ａｉ）が存在しない場合（ｎｏ）は、Ｓ６の処理を繰り返す。
次に、Ｓ８の処理の後に、Ｓ９においてインデックスｋをｋ＋１にセットするという処理を行い、Ｓ１０に進む。Ｓ１０の「（ａｉ，ａｋ）が存在」とは、（ａｉ，ａｋ）が１であることを示しており、遺伝子ｉから遺伝子ｋへの関係の存在を意味する。Ｓ１０において、（ａｉ，ａｋ）が存在する場合（ｙｅｓ）は、ｔｎに（ａｊ，ａｉ）ｎｅｇと（ａｉ，ａｋ）ｎｅｇの値を足したものをセットするという処理を行う（Ｓ１１）。（ａｉ，ａｋ）が存在しない場合（ｎｏ）は、Ｓ９の処理を行う。（ａｊ，ａｉ）ｎｅｇは、遺伝子ｊから遺伝子ｉへの経路の中に含まれる負の制御の数を示している。
Ｓ１２は、（ａｊ，ａｋ）とｔｎの値を比較するステップであり（Ｓ１２）、（ａｊ，ａｋ）ｎｅｇが偶数であり、かつ、ｔｎが偶数の場合（ｙｅｓ）は（ａｊ，ａｋ）を０にして遺伝子ｊから遺伝子ｋへの関係を削除し（ａｊ，ａｋ）ｎｅｇをｔｎにセットする処理を行う（Ｓ１３）。（ａｊ，ａｋ）ｎｅｇ及びｔｎの少なくとも一方が奇数の場合（両方が偶数でない場合）（ｎｏ）は、（ａｊ，ａｋ）ｎｅｇが奇数であり、かつ、ｔｎが奇数であるかどうかを比較する処理を行う（Ｓ１４）。Ｓ１４において、（ａｊ，ａｋ）ｎｅｇが奇数であり、かつ、ｔｎが奇数の場合（ｙｅｓ）は、（ａｊ，ａｋ）を０にして遺伝子ｊから遺伝子ｋへの関係を削除し（ａｊ，ａｋ）ｎｅｇをｔｎにセットする処理を行う（Ｓ１５）。Ｓ１４において、（ａｊ，ａｋ）ｎｅｇ及びｔｎの少なくとも一方が偶数の場合（ｎｏ）は、インデックスｋがｎ（遺伝子の総数を示している）かどうか比較する処理を行う（Ｓ１５）。
次に、Ｓ１６においてインデックスｋがｎかどうか比較する処理を行う。ｋ＝ｎの場合（ｙｅｓ）はＳ１７においてインデックスｊがｎかどうか比較する処理を行う。ｋ＜ｎの場合（ｎｏ）はＳ９の処理を行う。
Ｓ１７においてｊ＝ｎの場合（ｙｅｓ）は、Ｓ１８においてインデックスｉがｎかどうか比較する処理を行う。Ｓ１７の処理でｊ＜ｎの場合（ｎｏ）は、Ｓ４の処理を行う。
このようにして全ての遺伝子に対して冗長であると考えられる遺伝子間の制御を削除することによって、遺伝子ネットワークＧｎが得られる（Ｓ１９）。
上記フローチャートによれば、冗長な可能性があるネットワークの全てをチェックして排除することができる。このアルゴリズムの計算量はＯ（ｎ３）である。
本発明の推定方法は、単細胞生物のみならず、細胞生物における遺伝子ネットワークについても適用することができる。多細胞生物における遺伝子ネットワークは、上記推定方法に、さらに以下の要素を加えて推定することができる。すなわち、単細胞生物の場合と同様に（１）遺伝子発現量の検出エンジン、（２）比較エンジン、（３）因果関係作成エンジン、を適用した後、各細胞あるいは組織ごとに（４）冗長関係除去エンジンを適用する。
２．プログラムを記録した記録媒体
本発明の遺伝子ネットワーク推定方法においては、遺伝子の発現量（検出値）を読み取るプログラム、検出値から差を求めるプログラム及び因果関係作成プログラムを記録したコンピュータ読み取り可能な記録媒体を使用するのが好ましい。このプログラムは、別の記録媒体に記録されていてもよい。また、記録媒体には、ＣＤ−ＲＯＭ、ハードディスク、ＲＯＭ、ＲＡＭ等が含まれる。
発明を実施するための最良の形態
以下、実施例により本発明をさらに具体的に説明する。但し、本発明はこれら実施例にその技術的範囲が限定されるものではない。
〔実施例１〕遺伝子ネットワークの推定
本実施例は、本発明の方法の効果を確認するためのモデル実験であり、仮想ネットワークを作成し、その仮想ネットワークに対して本発明手法がどの程度の割合で遺伝子ネットワークを推定できるかを確認するものである。
コンピュータ上で、遺伝子数ＮをＮ＝１０、２０、５０又は１００に設定し、最大入次数ｋを２、４又は８に設定した。最大入次数とは、ある遺伝子が最大いくつの遺伝子から制御を受けているかを意味する。この設定を基に乱数表を使用して、ランダムに、それぞれの数に設定した遺伝子群について１００通りのパターンの有向グラフで構成されるネットワークＴをシミュレートした。なお、各ネットワークは、環状調節ネットワークを含むが、自己調節ネットワークは含まないようプログラムした。
モデル式のパラメータ及び調節タイプ（式（Ｉ）において、Ｒａ、Ｗａｂ及びλａの値）、すなわち各遺伝子がそれぞれの遺伝子とは無関係に調節されるのか、あるいは協同的に調節されるのかについてのパラメータは、ランダムに決定した。
（式Ｉ中、Ｒａは、遺伝子ａから微小時間ｄｔの間に生成される最大量を決定するレートを表わし、ｇはシグモイド関数を表わし、Ｗａｂはｖｂから受ける影響の重みを表わし、ｖｂ遺伝子ｂの遺伝子発現レベルを表わし、ｈａは遺伝子ａの一般的な転写因子からの影響を表わし、λａは遺伝子ａの遺伝子産物の減衰率を表わし、ｖａは遺伝子ａの遺伝子発現レベルを表わす。）
野生型遺伝子ネットワークの各遺伝子の発現量は、モデル式で１０００シミュレーションステップ計算し、安定化させることによりシミュレートし、変異型遺伝子ネットワーク（ネットワーク中の特定の遺伝子の機能を破壊したもの）の発現量は対応する遺伝子のパラメーターＲａをゼロにした後、野生型遺伝子ネットワークと同様にしてシミュレートした。各遺伝子を単独で欠失させたネットワーク（単一欠失ネットワーク）を、ネットワーク上に存在する全ての遺伝子それぞれについて作成し、シミュレートした。
上記のようにしてシミュレートされたネットワーク（標的ネットワーク）に対し、本発明の方法により推定したネットワークが、標的ネットワークとどの程度相違するか、その比率を求めた。
推定したネットワークと標的ネットワークとの間の相違点は、２つの基準（感度及び特異性）によって評価した。「感度」は、標的ネットワーク中のエッジ数に対する、正しく推定されたエッジ数の割合（％）であると定義する。この感度は、シミュレートした標的ネットワーク中のエッジの総数に対して、推定したネットワーク中に含まれる標的ネットワークに含まれているエッジの数の割合を示すものである。例えば、標的ネットワーク中のエッジが１０本であり、推定ネットワーク中に存在し、かつ標的ネットワーク中にも存在するエッジが８本であれば、感度は８０％となる。
また、特異性は、推定ネットワーク中に存在するエッジの総数に対する、正しく推定されたエッジ数の割合（％）を意味する。
この特異性は、推定したネットワーク中のエッジのうち何本が、標的ネットワーク中のエッジと一致するかを示すものである。上述の例によれば、推定ネットワーク中に存在するエッジ８本のうち、７本が標的ネットワーク中に存在するエッジと一致したとすれば、特異性は８７．５％となる。
試験結果を表２に示す。１０個の遺伝子を想定してシミュレートしたときの標的ネットワークは１００通り作成したので、推定操作は１００回行うことになる（２０個、５０個及び１００個の遺伝子の場合も同様）。従って、感度、特異性等のデータ（表２）は、１００回の推定操作の平均値を示す。
各測定値は、シミュレートした１００個の標的ネットワークの平均値であり、標準誤差は（）内に記す。エッジの総数とは、ネットワーク中の２つのノード間を結ぶ矢印（活性化又は抑制を示す矢印）の本数の平均値±標準誤差を意味する。
表２から明らかなように、特異性の平均値は常に感度の平均値よりも高く、また感度はネットワークサイズＮに比例して増加した。
特異性の平均値は、Ｎとは無関係に、入次数ｋ＝２でおよそ９０％である。すなわち、遺伝子数が２０個のときの標的ネットワークを考えると、標的ネットワーク中に引かれるエッジ数は平均３０．５本であり、推定ネットワーク中に引かれるエッジ数は２０．６本である。この２０．６本のうち、標的ネットワーク中のエッジと一致するエッジ数は１８．６本である。
感度が低くても特異性が高い場合は、標的ネットワーク中のエッジ数に見合う数のエッジ数が推定ネットワーク中に存在しなくても、エッジを引いた分については標的ネットワーク中のエッジと一致する確率が高いことを意味する。従って、本発明の方法は、極めて高い割合でネットワークを推定することができる。
なお、感度と特異性の平均値は、ｋの増加に比例して減少する。
〔実施例２〕連続値データと２値データとの比較
本発明の方法の主な特徴の一つは、発現データを連続値として適用できる点である。連続値発現データを使うことの利点を明らかにすべく、本発明者は、連続値データを使う場合及び２値データを使う場合について、感度及び特異性の比較を行った。２値は、閾値に従って連続値発現データを２値データに補正した。この閾値は、最小発現レベルｘｍｉｎと最大発現レベルｘｍａｘとの間の中間値とした。但し、最小発現レベルｘｍｉｎは、０（ゼロ）と定義した（ｘｍａｘの１／２を閾値とした）。遺伝子数、発現量の検出及び算定並びに推定処理は実施例１と同様にして行った。
結果を表３に示す。各測定値は、シミュレートした１００個の標的ネットワークの平均値である。連続値を用いたときの推定結果は、感度及び特異性のいずれも、２値データを用いたときの推定結果よりもかなり上回った。
〔実施例３〕従来法との比較
従来知られている遺伝子ネットワーク推定法は、２値データを分析するように設計してあるため、連続値データには適用できない。
本実施例では、本発明の方法を用いた推定方法と、従来法による推定方法とを比較した。従来法では、実施例（１）と同様にして得られたオリジナルデータを実施例２と同様にして２値データに補正したしたものを使用し、
本発明の方法は、上記オリジナルデータをそのまま使用した。遺伝子数、発現量の検出及び算定並びに推定処理は実施例１と同様にして行った。
その結果、本発明の方法は、従来法と比較して優れていることが示された（表４）。
各測定値は、シミュレートした１００種類の標的ネットワークから得られた値の平均値である。シンボル「＊」は、計算実行不可能を意味する。
遺伝子数が２０（Ｎ＝２０）であり、ｋ＝８の場合、従来法の感度は本発明の方法をわずかに上回るが、この２つの方法間に統計的に有意な相違はない（ｐ＜０．０５）。この結果は、従来法が、ブーリアンネットワークモデルを使って生成した人工遺伝子発現データからエッジを推定できるとしても、実際の遺伝子発現データには適さないことを意味する。
本発明の方法における計算量はＯ（ｎ３）であるので、Ｎ＝１００、ｋ＝２の場合であっても、遺伝子調節ネットワークを推定することができた。
〔実施例４〕ＤＮＡマイクロアレイから得られた発現量を用いた遺伝子ネットワークの推定
本実施例では、実際の生体材料として酵母（Ｓａｃｃｈａｒｏｍｙｃｅｓ ｃｅｒｅｖｉｓｉａｅ）の遺伝子（約２５０種類）をＤＮＡマイクロアレイを用いて、発現量を計測した遺伝子発現データを用いて遺伝子ネットワークを推定した。
用いた遺伝子発現データは、Ｈｕｇｈｅｓらの論文（Ｈｕｇｈｅｓ，Ｔ．Ｒ．ｅｔ ａｌ．，Ｃｅｌｌ，１０２：１０９−１２６，２０００．）にて公開されているものであり、野生型と突然変異型の遺伝子発現レベルの比をマイクロアレイにより計測したものである。このデータのことをＣｏｍｐｅｎｄｉｕｍと呼んでおり、遺伝子発現レベルの比の他にもｐ−ｖａｌｕｅという突然変異型の発現レベルが野生型と比べてどれくらい有意的に差が生じているかを示した値も遺伝子発現レベルの比に付随している。この値が低ければ低いほど、野生型と比較して有意的に遺伝子発現レベルに差が生じていることになる。本実施例では、ｐ−ｖａｌｕｅが０．０１以下のものに発現レベルに差があるとみなし、その他のものには野生型と比較して差を生じていないとした。このデータの中から単一遺伝子の欠失型の遺伝子発現データのみを取り出し、さらに欠失した遺伝子の発現データのみを選び出し、約２５０個の遺伝子欠失型に対する約２５０個の遺伝子の発現データを取り出して本発明の方法に適用した。その結果、図７に示した遺伝子ネットワークを推定することができた。
また、文献的に知られているネットワークについて、本発明の方法により推定を行った。結果を図８に示す。図８左パネルは文献的に知られているネットワーク（Ｒｏｂｅｒｔｓ，Ｃ．Ｊ．ｅｔ ａｌ．，Ｓｃｉｅｎｃｅ，２８７（５４５４）：８７３−８８０，２０００．；Ｏｅｈｌｅｎ，Ｌ．Ｊ．ｅｔ ａｌ．，Ｍｏｌ．Ｃｅｌｌ．Ｂｉｏｌ．，１６（６）：２８３０−２８３７，１９９６．；Ｄｉｅｔｚｅｌ，Ｃ．ａｎｄ Ｋｕｒｊａｎ，Ｊ．，Ｍｏｌ．Ｃｅｌｌ．Ｂｉｏｌ．，７（１２）：４１６９−４１７７，１９８７．；Ｅｒｒｅｄｅ，Ｂ．ａｎｄ Ａｍｍｅｒｅｒ，Ｇ．，Ｇｅｎｅｓ ａｎｄ Ｄｅｖ．，３（９）：１３４９−１３６１，１９８９．；Ｏｅｈｌｅｎ，Ｌ．ａｎｄ Ｃｒｏｓｓ，Ｆ．Ｒ．，ＦＥＢＳ Ｌｅｔｔ．，４２９（１）：８３−８８，１９９８．）を、右パネルは本発明の推定方法により推定されたネットワークである。転写因子であるｓｔｅ１２より上流の遺伝子に関しては、上記文献によれば、ｓｔｅ２からｓｔｅ１２へのリン酸化リレーの関係が特定されているが、本発明の方法によれば、マイクロアレイの発現データがｍＲＮＡの発現レベルを計測していることとｓｔｅ１２が自己制御されているため、ｓｔｅ１２の上流の遺伝子からｓｔｅ１２への関係が推定されている。ｓｔｅ１２が制御する遺伝子に関しては、上記文献によれば、ｆｕｓ３，ｔｅｃ１，ｆａｒ１，ｓｔｅ２，ｓｓｔ２が存在するが、本発明の方法によれば、マイクロアレイの発現データがｍＲＮＡの発現レベルを計測していることから、文献のものと一致したｆｕｓ３，ｔｅｃ１，ｆａｒ１，ｓｔｅ２，ｓｓｔ２への関係が推定された。
従って、本発明の方法は、マイクロアレイのデータから推定されるネットワークとして既知の遺伝子ネットワークと比較して妥当な遺伝子ネットワークが推定できた。
本明細書で引用した全ての刊行物、特許及び特許出願は、そのまま参考として本明細書に取り入れるものとする。
産業上の利用可能性
本発明により、遺伝子ネットワークの推定方法が提供される。本発明の方法により得られたネットワークは、例えば医薬品を開発（ドラッグデザイン等）するために有用である。
【図面の簡単な説明】
図１は、本発明の推定装置のブロック図である。
図２は、遺伝子ネットワークを示す図である。
図３は、遺伝子の発現量と遺伝子調節との関係を示す図である。
図４は、本発明の推定方法を示す有向グラフの図である。
図５は、本発明の推定方法を示す有向グラフの図である。
図６は、冗長な因果関係を排除するプログラムを実行するためのフローチャートの図である。
図７は、本発明の方法により推定された遺伝子ネットワークを示す図である。
図８は、本発明の方法により推定された遺伝子ネットワークを示す図である。
符号の説明
１０１…ＣＰＵ、１０２…ＲＯＭ、１０３…ＲＡＭ、１０４…入力部、１０５…送信／受信部、１０６…表示部、１０７…ハードディスクドライブ、１０８…ＣＤ−ＲＯＭドライブ、１０９…ＣＤ−ＲＯＭ、１１０…パブリックデータベース
本発明は、遺伝子ネットワークの推定方法、推定システム及び記録媒体に関する。
背景技術
最近の分子生物学分野の進歩により、得られる遺伝子情報量は膨大である。従って、次々と明らかになり急激に増加する配列データ数又は膨大な遺伝子発現データから情報を取り出すには、コンピュータを利用して努力する必要がある。これまでは、ホモロジー検索、タンパク質分類、遺伝子集積など各種のコンピュータツールを開発することに関心が注がれてきた。
こうした試みの中で、遺伝子発現データから遺伝子調節ネットワークを推論する方法として、数例ではあるがいくつかの研究が知られている。遺伝子発現データは、時系列データ（時系列データとは、対象遺伝子群の遺伝子発現量を経時的に測定したもの）又は定常状態データ（定常状態データとは、対象遺伝子群の遺伝子発現量を複数の異なる実験条件下（例えば、遺伝子変異や薬剤投与）で測定したもの）のどちらかの形で得られる。
時系列を分析する方法の場合は、種々の方法、例えば情報理論、遺伝アルゴリズム又はシミュレーテッドアニーリング法を用いてネットワークを推定することができる（Ｌｉａｎｇ，Ｓ．ｅｔ ａｌ．，Ｐｒｏｃ．Ｐａｃｉｆｉｃ Ｓｙｍｐ．Ｂｉｏｃｏｍｐｕｔｉｎｇ ’９８，Ｗｏｒｌｄ Ｓｃｉｅｎｔｉｆｉｃ，１８−２９，１９９８．；Ｍｏｒｏｈａｓｈｉ，Ｍ．ａｎｄ Ｋｉｔａｎｏ，Ｈ．，Ｐｒｏｃ．５ｔｈ Ｅｕｒｏ．Ｃｏｎｆ．Ａｒｔｉｆｉｃｉａｌ Ｌｉｆｅ，Ｓｐｒｉｎｇｅｒ，４７７−４８６，１９９９．；Ｍｊｏｌｓｎｅｓｓ，Ｅ．，ｅｔ ａｌ．，Ｔｅｃｈ．Ｒｅｐｔ．ＪＰＬ−ＩＣＴＲ−９９−４，Ｊｅｔ Ｐｒｏｐｕｌｓｉｏｎ Ｌａｂ．，ＮＡＳＡ，１９９９．）。しかし、時系列アプローチでは、試験結果を非常に短いインターバルで取得し、かつ試験的ノイズのない状態で取得することが必要とされる。これは現在の技術では非常に困難である。
一方、定常状態を分析する方法は、すでにいくつかの方法が提唱されている。この定常状態データは、特定の遺伝子活性を変異させて、例えば遺伝子を欠失させたりあるいは過剰発現させることにより得ることができる。欠失は、酵母ゲノム欠失コンソーシアム（Ｙｅａｓｔ Ｇｅｎｏｍｅ Ｄｅｌｅｔｉｏｎ Ｃｏｎｓｏｒｔｉｕｍ）等で現在大規模に行われており、これによって各種遺伝子欠失型の発現プロファイルを近い将来容易に入手することができる（Ｗｉｎｚｅｌｅｒ，Ｅ．Ａ．ｅｔ ａｌ．，Ｓｃｉｅｎｃｅ，２８５（５４２９）：９０１−９０６，１９９９．）。
Ａｋｕｔｓｕらは、遺伝子ネットワークがブーリアンネットワークモデルで表現された場合（各遺伝子の遺伝子発現が活性、不活性の２状態で表現され、それら遺伝子間の関係がブール代数によって表現された場合）に、必要となる試験数の上限と下限を計算した（Ａｋｕｔｓｕ，Ｔ．ｅｔ ａｌ．，Ｐｒｏｃ．９ｔｈ ＡＣＭ−ＳＩＡＭ Ｓｙｍｐ．Ｄｉｓｃｒｅｔｅ Ａｌｇｏｒｉｔｈｍｓ，６９５−７０２，１９９８．）。さらに最近、Ｉｄｅｋｅｒらは、Ｐｒｅｄｉｃｔｏｒ法と呼ばれる推定方法を提唱した（Ｉｄｅｋｅｒ，Ｔ．Ｅ．ｅｔ ａｌ．，Ｐｒｏｃ．Ｐａｃｉｆｉｃ Ｓｙｍｐ．Ｂｉｏｃｏｍｐｕｔｉｎｇ ２０００，Ｗｏｒｌｄ Ｓｃｉｅｎｔｉｆｉｃ，３０５−３１６，２０００．）。この方法は、組み合わせに関する最適化法を用いて候補ネットワーク（発現データと一致するブーリアンネットワークモデル）を提供するものである。
しかし、これらの方法は、遺伝子発現レベルを２値として遺伝子ネットワークを現すものである。従って、試験データは通常は連続値（多値）を有するにもかかわらず、この方法を適用するにはデータを２値に補正する必要がある。こうした補正によって、ネットワークを推定するのに必要な情報が失われる可能性がある。例えば２値を利用して推定をすると、３つの状態（例えば野生型、欠失、過剰発現）の遺伝子発現レベルを正確に表すことが困難であり、調節関係が正確に反映されない結果となる。
発明の開示
本発明は、遺伝子ネットワークの推定方法、推定システム及び記録媒体を提供することを目的とする。
本発明者は、上記課題を解決するため鋭意研究を行った結果、突然変異を誘発したときに得られる遺伝子発現プロファイル（検出値）を用いて、遺伝子発現量を単純化（二値化）することなく遺伝子ネットワークの推定を行うことに成功し、本発明を完成するに至った。
すなわち、本発明は、複数の遺伝子のうち一の遺伝子の発現量を２つの条件下でそれぞれ発現させたときの他の遺伝子の発現量をそれぞれ検出し、得られる検出値の差を求め、当該差を指標として前記一の遺伝子と他の遺伝子との因果関係を導くことを特徴とする遺伝子ネットワークの推定方法である。
さらに、本発明は、複数の遺伝子のうち一の遺伝子を変異させたときの他の遺伝子の発現レベルと、当該一の遺伝子を変異させないときの当該他の遺伝子の発現レベルとをそれぞれ検出し、得られる検出値の差を求め、当該差を指標として前記一の遺伝子と他の遺伝子との因果関係を導くことを特徴とする遺伝子ネットワークの推定方法である。上記遺伝子ネットワークとしては、単細胞生物内のもの又は多細胞生物内のものが挙げられる。遺伝子の変異としては遺伝子の破壊又は過剰発現が挙げられる。上記推定方法においては、得られた因果関係の矛盾の有無を検定することが好ましい。
さらに、本発明は、複数の遺伝子のうち一の遺伝子の発現量を２つの条件下で発現させたときの他の遺伝子の発現量をそれぞれ検出する手段と、得られる検出値の差を求める手段と、当該差を指標として前記一の遺伝子と他の遺伝子との因果関係を導く手段とを含んでなる、遺伝子ネットワークの推定システムである。
さらに、本発明は、複数の遺伝子のうち一の遺伝子を変異させたときの他の遺伝子の発現レベルと、当該一の遺伝子を変異させないときの当該他の遺伝子の発現レベルとをそれぞれ検出する手段と、得られる検出値の差を求める手段と、当該差を指標として前記一の遺伝子と他の遺伝子との因果関係を導く手段とを含んでなる、遺伝子ネットワーク推定システムである。本発明の推定システムにおいては、得られた因果関係の矛盾の有無を検定する手段を含めることができる。
さらに、本発明は、複数の遺伝子のうち一の遺伝子を変異させたときの他の遺伝子の発現レベルと、当該一の遺伝子を変異させないときの当該他の遺伝子の発現レベルとをそれぞれ検出する手段と、得られる検出値の差を求める手段と、当該差を指標として前記一の遺伝子と他の遺伝子との因果関係を導く手段とを含んでなる、コンピュータを遺伝子ネットワーク推定システムとして機能させるためのプログラムを記録したコンピュータ読み取り可能な記録媒体である。
以下、本発明を詳細に説明する。本明細書は、本願の優先権の基礎である日本国特許出願２０００−３４５９８２号の明細書及び／又は図面に記載される内容を包含する。
本発明は、複数の遺伝子のうちの一の遺伝子（遺伝子ａ）の発現量を、２つの条件下（条件１、条件２）でそれぞれ発現させたときに、他の遺伝子（遺伝子ｂ）の発現量が、条件１のときと条件２のときにおいてどのように変化するかに着目し、それぞれの条件における遺伝子ｂの発現量をそれぞれ調べ（検出し）、その調べた値（検出値）の差を求め、当該差を指標として一の遺伝子ａと他の遺伝子ｂとの因果関係を導くことを特徴とするものである。
例えば、遺伝子ａの野生型の場合（条件１）、及び当該遺伝子ａを変異させた場合（条件２）の２条件を設定し、これらの２条件下における当該他の遺伝子ｂの発現量の差を検出することにより、遺伝子ａと遺伝子ｂとの間の因果関係を求め、これらの遺伝子間に生ずる遺伝子ネットワークを推定する。
本発明においては、ＯＮとＯＦＦとの関係のように二値を扱うのではなく、変化する値の全部を連続値（多値）として扱うことを特徴とするため、二値よりも更に高精度にネットワークを推定することが可能である。
１．遺伝子ネットワーク推定方法及び推定システム
本発明の推定システムは、
（ｉ）ある一の遺伝子の発現量を２条件に設定した場合に他の遺伝子の発現レベル（発現量）をそれぞれ検出する手段（「検出エンジン」ともいう）、
（ｉｉ）得られる検出値を比較してその差を求める手段（「比較エンジン」ともいう）、
（ｉｉｉ）当該差を指標として前記一の遺伝子と他の遺伝子との因果関係を求める手段（「因果関係作成エンジン」ともいう）、
（ｉｖ）当該因果関係の中から冗長な因果関係を除去する手段（「冗長関係除去エンジン」ともいう）
を含む装置により構成される。
ここで、本発明の推定システムの構成例を示すブロック図を示す（図１）。
図１に示す推定システムは、ＣＰＵ１０１ ＲＯＭ１０２、ＲＡＭ１０３、入力部１０４、送信／受信部１０５、出力部１０６、ハードディスクドライブ（ＨＤＤ）１０７及びＣＤ−ＲＯＭドライブ１０８を備える。
ＣＰＵ１０１は、ＲＯＭ１０２、ＲＡＭ１０３又はＨＤＤ１０７に記憶されているプログラムに従って、遺伝子ネットワークシステム全体を制御し、後述するネットワーク推定処理を実行する。ＲＯＭ１０２は、遺伝子ネットワーク推定システムの動作に必要な処理を命令するプログラム等を格納する。ＲＡＭ１０３は、ネットワーク推定処理を実行する上で必要なデータを一時的に格納する。入力部１０４は、キーボードやマウス等であり、ネットワーク推定処理を実行する上で必要な条件を入力するとき等に操作される。送信／受信部１０５は、ＣＰＵ１０１の命令に基づいて、通信回線を介してパブリックデータベース１１０等との間でデータの送受信処理を実行する。出力部１０６は、遺伝子の発現量、入力部１０４から入力された各種条件、遺伝子の塩基配列、ネットワークの推定結果等を、ＣＰＵ１０１からの命令に基づいて表示処理を実行する。なお、出力部１０６としては、コンピュータのディスプレイ、又はプリンターなどが例示される。ＨＤＤ１０７は、遺伝子ネットワーク推定プログラム、遺伝子発現量、塩基配列等を格納し、ＣＰＵ１０１の命令に基づいて、格納しているプログラム又はデータ等を読み出し、例えばＲＡＭ１０３に格納する。ＣＤ−ＲＯＭドライブ１０８は、ＣＰＵ１０１の指示に基づいて、ＣＤ−ＲＯＭ１０９に格納されている遺伝子ネットワーク推定プログラム又は発現量等から、プログラム又はデータ等を読み出し、例えばＲＡＭ１０３に格納する。
ＣＰＵ１０１は、入力部などから受け取ったデータを出力部１０６に供給するとともに、データベースから受け取ったデータに基づいて遺伝子ネットワークの推定処理を実行する。
ここで、遺伝子ネットワークを推定するための処理プログラムは、下記のごとく、遺伝子発現量の検出エンジン、比較エンジン及び因果関係作成エンジンから構成される。
（１）遺伝子発現量の検出エンジン
本発明において、遺伝子発現の検出は、ＤＮＡマイクロアレイ、オリゴヌクレオチドチップ、ＲＴ−ＰＣＲ、遺伝子発現の連続分析、定常発現レベル、プロテオミクス等を使用又は適用することにより行うことができ、これらの結果を遺伝子発現データ（発現量）として使用することができる。また、蛋白のリン酸化の状態や染色体ＤＮＡのメチル化の状態等、遺伝子機能発現を調節するあらゆる状態変化も遺伝子発現として検出できる。上記データは、本発明では検出値という。検出値は、絶対値であっても相対値であってもよい。また、実験により得られた値のみならず、文献的に又は遺伝子データベース等から得られた値でもよい。例えば、ＤＮＡマイクロアレイから得られる発現量は、マイクロアレイ上に発せられる蛍光強度等を測定し、その蛍光強度を数値化することができる。
データベースの具体例としては、インターネットを通じて利用できるデータベースを例示でき、より具体的には、ＧｅｎＢａｎｋ、ＤＤＢＪ、ＥＭＢＬ、ＮＣＩ６０等を例示できる。
（２）比較エンジン
比較エンジンでは、ある遺伝子の遺伝子発現量を２条件に設定した場合のこれらの２条件下における対象となる遺伝子の発現レベルの差の情報、並びに上記遺伝子間の機能（例えば活性化又は阻害）に関する情報を蓄積する。上記遺伝子発現量の２条件設定には遺伝子の変異が含まれ、遺伝子変異は、塩基配列の欠失、置換又は付加によりもとの遺伝子の機能を失わせることを意味し、遺伝子の破壊も変異に含まれる。
遺伝子に欠失、置換、付加等の変異を導入するには、Ｋｕｎｋｅｌ法、Ｇａｐｐｅｄ ｄｕｐｌｅｘ法等の公知の手法又はこれに準ずる方法を採用することができる。例えば部位特異的突然変異誘発法を利用した変異導入用キット（例えばＭｕｔａｎ−Ｋ（ＴＡＫＡＲＡ社製）やＭｕｔａｎ−Ｇ（ＴＡＫＡＲＡ社製））、ＴＡＫＡＲＡ社のＬＡ ＰＣＲ ｉｎ ｖｉｔｒｏ Ｍｕｔａｇｅｎｅｓｉｓシリーズキットなどを用いて、あるいは、遺伝子の一部又は全部を切断することにより、変異の導入が行われる。あるいは、遺伝子を変異したときの発現量データをデータベースから入手することもできる。
（３）因果関係作成エンジン
因果関係作成エンジンは、比較エンジンにより得られた発現量の差に基づいて、遺伝子間の因果関係を導き、最終的にネットワークを推定する手段である。この手段は、例えば遺伝子ａと遺伝子ｂとの因果関係を調べる場合において、遺伝子ａを変異させたときの遺伝子ｂの発現量が、遺伝子ａを変異させないとき（野生型のとき）の遺伝子ｂの発現量と比較してどの程度変化したときに、その遺伝子ａは遺伝子ｂに対して活性化するのか、抑制するのか、あるいは無関係であるのか、という判断を実行する手段である。
本発明において、便宜上、遺伝子ネットワークを「有向グラフ」と呼ばれる図に置き換えて考えてみる（図２）。図２において、アルファベットは遺伝子を表し、矢印（→）は遺伝子ａが遺伝子ｂの発現を促進する制御を、Ｔ字型矢印（「┤」）は遺伝子ｂが遺伝子ｃの発現を抑える制御を意味するものとする。個々の矢印を、本発明では「エッジ」という。遺伝子ａと遺伝子ｂとの間におけるネットワーク（因果関係）を推定する場合は、まず、遺伝子ａの発現量を２条件に設定した場合のこれらの２条件における遺伝子ｂの発現量、例えば、遺伝子ａを何ら操作しない定常型又は野生型の状態で検出した遺伝子ｂの発現量と、遺伝子ａを変異（例えば破壊）したときの遺伝子ｂの発現量とを検出し、次いで、それぞれ得られたデータ（検出値）の差を求める。この差を指標として遺伝子間の因果関係を導く。例えば、差に変化がなければ遺伝子ａは遺伝子ｂを活性化していない、又は無関係であると判断することができ、差に一定の変化があれば遺伝子ａは遺伝子ｂを活性化又は抑制していると判断することができる。ここで、「一定の変化」とは、いわゆる閾値を超える変化を示すものであり、閾値は、ネットワークを推定する際に目的となる遺伝子に応じて適宜設定できる。但し、閾値は、誤差検定を行い、例えば有意差危険度０．０１以下となるような基準で設定し、採用することが好ましい。
遺伝子の発現量は、一般に数値として得られる場合が多いため、仮に遺伝子間の相互関係を図３のように表すことができる。図３において、ａｎ（ｎ＝０，１，２，３）は「ノード（ｎｏｄｅ）」と呼ばれるものであり、遺伝子ネットワークを推定したときに作成されるエッジとエッジとの間の点（交点）を意味する。ノードは、遺伝子名で表示することもできることから、ノードの数と遺伝子の数（ｎの数）とは一致する。ｘｎ（ｎ＝０，１，２，３）は対応するノード（遺伝子）の発現量を意味する。ｗｔは野生型を意味する。
ある遺伝子ａ０〜ａ３の発現量の結果が図３Ｂのように得られたものとする。ａｎ −（ｎ＝０，１，２，３）は各ノードにおける遺伝子が変異（破壊）されたことを意味する。ここで、遺伝子ａ０を破壊した場合（ａ０ −）を考えてみる。
まず、ａ０ −における遺伝子ａ１の発現量と、野生型におけるａ１遺伝子の発現量との差は不変であるため（３．７５０−３．７５０＝０）、ａ０はａ１とは無関係、すなわちａ０からａ１へのネットワークは無いものと推定することができ、ノード間のエッジは無い（図４Ａ）。次に、ａ０ −における遺伝子ａ２の発現量と、野生型における遺伝子ａ２の発現量との差は−０．１７であり（８．７６９−８．９３９＝−０．１７）、遺伝子ａ０を破壊するとａ２の発現量は減少する。従って、閾値を無視して考えれば遺伝子ａ０はａ２を活性化させると推定することができ、ノードａ０からａ２に向かうエッジを引く（つまり矢印を表示する）ことができる（図４Ｂ）。さらに、ａ０ −における遺伝子ａ３の発現量と、野生型における遺伝子ａ３の発現量との差は−０．０６７であり（０．０１１−０．０７８＝−０．０６７）、遺伝子ａ０を破壊するとａ３の発現量は減少する。従って、遺伝子ａ０はａ３を活性化させるものと推定することができ、ａ０からａ３に向かってエッジを引くことができる（図４Ｃ）。
遺伝子ａ１を破壊した場合、遺伝子ａ２を破壊した場合及び遺伝子ａ３を破壊した場合についても上記と同じ要領で処理することにより、４つの遺伝子間のネットワークを推定することができる（図４Ｄ）。特に、遺伝子ａ２を破壊した場合における遺伝子ａ３の発現量と野生型における遺伝子ａ３の発現量との差は５．３９８であり（５．４７６−０．０７８＝５．３９８）、遺伝子ａ２を破壊するとａ３の発現量が上昇する。従って、遺伝子ａ２は遺伝子ａ３を抑制していると推定することができ、ノード間のエッジは抑制を示す矢印（Ｔ字型矢印）となる（図４Ｄ）。
以上のネットワークの関係をまとめると、図３Ａのようになる。図３Ａにおいて、「＋」は活性化（→）を、「−」は抑制（┤）を意味する。
ある遺伝子ａを欠失又は過発現させたときの遺伝子ｂの発現レベルの増加及び減少との関係を表１に示す。この比較プロセスの計算量はＯ（ｎ２）である。
図５に示すように、遺伝子ａは遺伝子ｂを活性化し、活性化された遺伝子ｂは遺伝子ｃを活性化する遺伝子ネットワークを考える。この場合において、遺伝子ａを不活性化（例えば破壊）すると、遺伝子ｂは活性化されず、その結果遺伝子ｃも活性化されなくなる。従って、実際には遺伝子ａが遺伝子ｃを活性化する能力はなくても、結果として遺伝子ａは遺伝子ｃを間接的に活性化する能力があると誤って認識される場合がある。このように、誤った認識により作成される因果関係を、本発明では「冗長な因果関係」という。本発明においては、そのような間接的なエッジ（図５、破線の矢印）を見つけることができる。間接的なエッジ及び直接的なエッジが存在する遺伝子調節ネットワークにおいて、本発明では、各遺伝子が他の遺伝子に間接的な効果を及ぼすかどうかを探り、そのネットワーク内の間接的エッジをチェックすることができる。この間接的効果は、エッジの道筋（ルート）に関与する負の調節数の偶奇性にのみ依存する（Ｔｈｉｅｆｆｒｙ，Ｄ．，ａｎｄ Ｔｈｏｍａｓ，Ｒ．，Ｐｒｏｃ．Ｐａｃｉｆｉｃ Ｓｙｍｐ．Ｂｉｏ−ｃｏｍｐｕｔｉｎｇ ’９８，Ｗｏｒｌｄ Ｓｃｉｅｎｔｉｆｉｃ，７７−８８，１９９８．）。
本発明においては、冗長な因果関係、すなわち間接的な遺伝子調節を除去してネットワークを推定する。こうした冗長な因果関係は、グラフ理論（Ｇｒｏｓｓ，Ｊ．，ａｎｄ Ｙｅｌｌｅｎ，Ｊ．，ＣＲＣ Ｐｒｅｓｓ，１９９９．）の分野のＷａｒｓｈａｌｌのアルゴリズムを修正することにより見つけることができる。このような冗長な因果関係を除外する処理においては、冗長な因果関係を見つけるためのプログラムを実行する。
修正したＷａｒｓｈａｌｌのアルゴリズムを実行するためのフローチャートを図６に示す。
図６におけるステップは以下の通りである。すなわち、冗長な因果関係を含む遺伝子ネットワーク（冗長な因果関係が除去されていないネットワーク）を初期値Ｇ０とし、全ての遺伝子に１からｎ（遺伝子の総数）の通し番号をつける。次に遺伝子ｘから遺伝子ｙへの関係をマトリクスの要素（ａｘ，ａｙ）と表現し、活性（正）か抑制（負）の関係が存在する場合に（ａｘ，ａｙ）を１にセットする（存在しない場合には０にセットする）。また、その関係が抑制（負の制御）の関係の場合には、マトリクスの各要素が持つカウンタ（ａｘ，ａｙ）ｎｅｇを１にセットし、その他の場合には０をセットする（Ｓ１）。まず、負の制御の総数をカウントするカウンタｔｎを初期化する（Ｓ２）。次に図６のフローチャートに従って処理を進めていく。Ｓ３は、インデックスｉを０に初期化するという処理を行うステップである。Ｓ４は、インデックスｉをｉ＋１にセットするという処理を行うステップである。Ｓ５は、インデックスｊを０に初期化するという処理を行うステップである。Ｓ６は、インデックスｊをｊ＋１にセットするという処理を行うステップである。図６で（ａｊ，ａｉ）は遺伝子ｊから遺伝子ｉへの関係を意味しており、（ａｊ，ａｉ）が存在するという意味は、（ａｊ，ａｉ）が１であることを示しており、遺伝子ｊから遺伝子ｉへの関係が存在するということである。（ａｊ，ａｉ）が存在する場合（ｙｅｓ）は、Ｓ８においてインデックスｋを０に初期化する処理を行う（Ｓ７）。（ａｊ，ａｉ）が存在しない場合（ｎｏ）は、Ｓ６の処理を繰り返す。
次に、Ｓ８の処理の後に、Ｓ９においてインデックスｋをｋ＋１にセットするという処理を行い、Ｓ１０に進む。Ｓ１０の「（ａｉ，ａｋ）が存在」とは、（ａｉ，ａｋ）が１であることを示しており、遺伝子ｉから遺伝子ｋへの関係の存在を意味する。Ｓ１０において、（ａｉ，ａｋ）が存在する場合（ｙｅｓ）は、ｔｎに（ａｊ，ａｉ）ｎｅｇと（ａｉ，ａｋ）ｎｅｇの値を足したものをセットするという処理を行う（Ｓ１１）。（ａｉ，ａｋ）が存在しない場合（ｎｏ）は、Ｓ９の処理を行う。（ａｊ，ａｉ）ｎｅｇは、遺伝子ｊから遺伝子ｉへの経路の中に含まれる負の制御の数を示している。
Ｓ１２は、（ａｊ，ａｋ）とｔｎの値を比較するステップであり（Ｓ１２）、（ａｊ，ａｋ）ｎｅｇが偶数であり、かつ、ｔｎが偶数の場合（ｙｅｓ）は（ａｊ，ａｋ）を０にして遺伝子ｊから遺伝子ｋへの関係を削除し（ａｊ，ａｋ）ｎｅｇをｔｎにセットする処理を行う（Ｓ１３）。（ａｊ，ａｋ）ｎｅｇ及びｔｎの少なくとも一方が奇数の場合（両方が偶数でない場合）（ｎｏ）は、（ａｊ，ａｋ）ｎｅｇが奇数であり、かつ、ｔｎが奇数であるかどうかを比較する処理を行う（Ｓ１４）。Ｓ１４において、（ａｊ，ａｋ）ｎｅｇが奇数であり、かつ、ｔｎが奇数の場合（ｙｅｓ）は、（ａｊ，ａｋ）を０にして遺伝子ｊから遺伝子ｋへの関係を削除し（ａｊ，ａｋ）ｎｅｇをｔｎにセットする処理を行う（Ｓ１５）。Ｓ１４において、（ａｊ，ａｋ）ｎｅｇ及びｔｎの少なくとも一方が偶数の場合（ｎｏ）は、インデックスｋがｎ（遺伝子の総数を示している）かどうか比較する処理を行う（Ｓ１５）。
次に、Ｓ１６においてインデックスｋがｎかどうか比較する処理を行う。ｋ＝ｎの場合（ｙｅｓ）はＳ１７においてインデックスｊがｎかどうか比較する処理を行う。ｋ＜ｎの場合（ｎｏ）はＳ９の処理を行う。
Ｓ１７においてｊ＝ｎの場合（ｙｅｓ）は、Ｓ１８においてインデックスｉがｎかどうか比較する処理を行う。Ｓ１７の処理でｊ＜ｎの場合（ｎｏ）は、Ｓ４の処理を行う。
このようにして全ての遺伝子に対して冗長であると考えられる遺伝子間の制御を削除することによって、遺伝子ネットワークＧｎが得られる（Ｓ１９）。
上記フローチャートによれば、冗長な可能性があるネットワークの全てをチェックして排除することができる。このアルゴリズムの計算量はＯ（ｎ３）である。
本発明の推定方法は、単細胞生物のみならず、細胞生物における遺伝子ネットワークについても適用することができる。多細胞生物における遺伝子ネットワークは、上記推定方法に、さらに以下の要素を加えて推定することができる。すなわち、単細胞生物の場合と同様に（１）遺伝子発現量の検出エンジン、（２）比較エンジン、（３）因果関係作成エンジン、を適用した後、各細胞あるいは組織ごとに（４）冗長関係除去エンジンを適用する。
２．プログラムを記録した記録媒体
本発明の遺伝子ネットワーク推定方法においては、遺伝子の発現量（検出値）を読み取るプログラム、検出値から差を求めるプログラム及び因果関係作成プログラムを記録したコンピュータ読み取り可能な記録媒体を使用するのが好ましい。このプログラムは、別の記録媒体に記録されていてもよい。また、記録媒体には、ＣＤ−ＲＯＭ、ハードディスク、ＲＯＭ、ＲＡＭ等が含まれる。
発明を実施するための最良の形態
以下、実施例により本発明をさらに具体的に説明する。但し、本発明はこれら実施例にその技術的範囲が限定されるものではない。
〔実施例１〕遺伝子ネットワークの推定
本実施例は、本発明の方法の効果を確認するためのモデル実験であり、仮想ネットワークを作成し、その仮想ネットワークに対して本発明手法がどの程度の割合で遺伝子ネットワークを推定できるかを確認するものである。
コンピュータ上で、遺伝子数ＮをＮ＝１０、２０、５０又は１００に設定し、最大入次数ｋを２、４又は８に設定した。最大入次数とは、ある遺伝子が最大いくつの遺伝子から制御を受けているかを意味する。この設定を基に乱数表を使用して、ランダムに、それぞれの数に設定した遺伝子群について１００通りのパターンの有向グラフで構成されるネットワークＴをシミュレートした。なお、各ネットワークは、環状調節ネットワークを含むが、自己調節ネットワークは含まないようプログラムした。
モデル式のパラメータ及び調節タイプ（式（Ｉ）において、Ｒａ、Ｗａｂ及びλａの値）、すなわち各遺伝子がそれぞれの遺伝子とは無関係に調節されるのか、あるいは協同的に調節されるのかについてのパラメータは、ランダムに決定した。
野生型遺伝子ネットワークの各遺伝子の発現量は、モデル式で１０００シミュレーションステップ計算し、安定化させることによりシミュレートし、変異型遺伝子ネットワーク（ネットワーク中の特定の遺伝子の機能を破壊したもの）の発現量は対応する遺伝子のパラメーターＲａをゼロにした後、野生型遺伝子ネットワークと同様にしてシミュレートした。各遺伝子を単独で欠失させたネットワーク（単一欠失ネットワーク）を、ネットワーク上に存在する全ての遺伝子それぞれについて作成し、シミュレートした。
上記のようにしてシミュレートされたネットワーク（標的ネットワーク）に対し、本発明の方法により推定したネットワークが、標的ネットワークとどの程度相違するか、その比率を求めた。
推定したネットワークと標的ネットワークとの間の相違点は、２つの基準（感度及び特異性）によって評価した。「感度」は、標的ネットワーク中のエッジ数に対する、正しく推定されたエッジ数の割合（％）であると定義する。この感度は、シミュレートした標的ネットワーク中のエッジの総数に対して、推定したネットワーク中に含まれる標的ネットワークに含まれているエッジの数の割合を示すものである。例えば、標的ネットワーク中のエッジが１０本であり、推定ネットワーク中に存在し、かつ標的ネットワーク中にも存在するエッジが８本であれば、感度は８０％となる。
また、特異性は、推定ネットワーク中に存在するエッジの総数に対する、正しく推定されたエッジ数の割合（％）を意味する。
この特異性は、推定したネットワーク中のエッジのうち何本が、標的ネットワーク中のエッジと一致するかを示すものである。上述の例によれば、推定ネットワーク中に存在するエッジ８本のうち、７本が標的ネットワーク中に存在するエッジと一致したとすれば、特異性は８７．５％となる。
試験結果を表２に示す。１０個の遺伝子を想定してシミュレートしたときの標的ネットワークは１００通り作成したので、推定操作は１００回行うことになる（２０個、５０個及び１００個の遺伝子の場合も同様）。従って、感度、特異性等のデータ（表２）は、１００回の推定操作の平均値を示す。
表２から明らかなように、特異性の平均値は常に感度の平均値よりも高く、また感度はネットワークサイズＮに比例して増加した。
特異性の平均値は、Ｎとは無関係に、入次数ｋ＝２でおよそ９０％である。すなわち、遺伝子数が２０個のときの標的ネットワークを考えると、標的ネットワーク中に引かれるエッジ数は平均３０．５本であり、推定ネットワーク中に引かれるエッジ数は２０．６本である。この２０．６本のうち、標的ネットワーク中のエッジと一致するエッジ数は１８．６本である。
感度が低くても特異性が高い場合は、標的ネットワーク中のエッジ数に見合う数のエッジ数が推定ネットワーク中に存在しなくても、エッジを引いた分については標的ネットワーク中のエッジと一致する確率が高いことを意味する。従って、本発明の方法は、極めて高い割合でネットワークを推定することができる。
なお、感度と特異性の平均値は、ｋの増加に比例して減少する。
〔実施例２〕連続値データと２値データとの比較
本発明の方法の主な特徴の一つは、発現データを連続値として適用できる点である。連続値発現データを使うことの利点を明らかにすべく、本発明者は、連続値データを使う場合及び２値データを使う場合について、感度及び特異性の比較を行った。２値は、閾値に従って連続値発現データを２値データに補正した。この閾値は、最小発現レベルｘｍｉｎと最大発現レベルｘｍａｘとの間の中間値とした。但し、最小発現レベルｘｍｉｎは、０（ゼロ）と定義した（ｘｍａｘの１／２を閾値とした）。遺伝子数、発現量の検出及び算定並びに推定処理は実施例１と同様にして行った。
結果を表３に示す。各測定値は、シミュレートした１００個の標的ネットワークの平均値である。連続値を用いたときの推定結果は、感度及び特異性のいずれも、２値データを用いたときの推定結果よりもかなり上回った。
従来知られている遺伝子ネットワーク推定法は、２値データを分析するように設計してあるため、連続値データには適用できない。
本実施例では、本発明の方法を用いた推定方法と、従来法による推定方法とを比較した。従来法では、実施例（１）と同様にして得られたオリジナルデータを実施例２と同様にして２値データに補正したしたものを使用し、
本発明の方法は、上記オリジナルデータをそのまま使用した。遺伝子数、発現量の検出及び算定並びに推定処理は実施例１と同様にして行った。
その結果、本発明の方法は、従来法と比較して優れていることが示された（表４）。
遺伝子数が２０（Ｎ＝２０）であり、ｋ＝８の場合、従来法の感度は本発明の方法をわずかに上回るが、この２つの方法間に統計的に有意な相違はない（ｐ＜０．０５）。この結果は、従来法が、ブーリアンネットワークモデルを使って生成した人工遺伝子発現データからエッジを推定できるとしても、実際の遺伝子発現データには適さないことを意味する。
本発明の方法における計算量はＯ（ｎ３）であるので、Ｎ＝１００、ｋ＝２の場合であっても、遺伝子調節ネットワークを推定することができた。
〔実施例４〕ＤＮＡマイクロアレイから得られた発現量を用いた遺伝子ネットワークの推定
本実施例では、実際の生体材料として酵母（Ｓａｃｃｈａｒｏｍｙｃｅｓ ｃｅｒｅｖｉｓｉａｅ）の遺伝子（約２５０種類）をＤＮＡマイクロアレイを用いて、発現量を計測した遺伝子発現データを用いて遺伝子ネットワークを推定した。
用いた遺伝子発現データは、Ｈｕｇｈｅｓらの論文（Ｈｕｇｈｅｓ，Ｔ．Ｒ．ｅｔ ａｌ．，Ｃｅｌｌ，１０２：１０９−１２６，２０００．）にて公開されているものであり、野生型と突然変異型の遺伝子発現レベルの比をマイクロアレイにより計測したものである。このデータのことをＣｏｍｐｅｎｄｉｕｍと呼んでおり、遺伝子発現レベルの比の他にもｐ−ｖａｌｕｅという突然変異型の発現レベルが野生型と比べてどれくらい有意的に差が生じているかを示した値も遺伝子発現レベルの比に付随している。この値が低ければ低いほど、野生型と比較して有意的に遺伝子発現レベルに差が生じていることになる。本実施例では、ｐ−ｖａｌｕｅが０．０１以下のものに発現レベルに差があるとみなし、その他のものには野生型と比較して差を生じていないとした。このデータの中から単一遺伝子の欠失型の遺伝子発現データのみを取り出し、さらに欠失した遺伝子の発現データのみを選び出し、約２５０個の遺伝子欠失型に対する約２５０個の遺伝子の発現データを取り出して本発明の方法に適用した。その結果、図７に示した遺伝子ネットワークを推定することができた。
また、文献的に知られているネットワークについて、本発明の方法により推定を行った。結果を図８に示す。図８左パネルは文献的に知られているネットワーク（Ｒｏｂｅｒｔｓ，Ｃ．Ｊ．ｅｔ ａｌ．，Ｓｃｉｅｎｃｅ，２８７（５４５４）：８７３−８８０，２０００．；Ｏｅｈｌｅｎ，Ｌ．Ｊ．ｅｔ ａｌ．，Ｍｏｌ．Ｃｅｌｌ．Ｂｉｏｌ．，１６（６）：２８３０−２８３７，１９９６．；Ｄｉｅｔｚｅｌ，Ｃ．ａｎｄ Ｋｕｒｊａｎ，Ｊ．，Ｍｏｌ．Ｃｅｌｌ．Ｂｉｏｌ．，７（１２）：４１６９−４１７７，１９８７．；Ｅｒｒｅｄｅ，Ｂ．ａｎｄ Ａｍｍｅｒｅｒ，Ｇ．，Ｇｅｎｅｓ ａｎｄ Ｄｅｖ．，３（９）：１３４９−１３６１，１９８９．；Ｏｅｈｌｅｎ，Ｌ．ａｎｄ Ｃｒｏｓｓ，Ｆ．Ｒ．，ＦＥＢＳ Ｌｅｔｔ．，４２９（１）：８３−８８，１９９８．）を、右パネルは本発明の推定方法により推定されたネットワークである。転写因子であるｓｔｅ１２より上流の遺伝子に関しては、上記文献によれば、ｓｔｅ２からｓｔｅ１２へのリン酸化リレーの関係が特定されているが、本発明の方法によれば、マイクロアレイの発現データがｍＲＮＡの発現レベルを計測していることとｓｔｅ１２が自己制御されているため、ｓｔｅ１２の上流の遺伝子からｓｔｅ１２への関係が推定されている。ｓｔｅ１２が制御する遺伝子に関しては、上記文献によれば、ｆｕｓ３，ｔｅｃ１，ｆａｒ１，ｓｔｅ２，ｓｓｔ２が存在するが、本発明の方法によれば、マイクロアレイの発現データがｍＲＮＡの発現レベルを計測していることから、文献のものと一致したｆｕｓ３，ｔｅｃ１，ｆａｒ１，ｓｔｅ２，ｓｓｔ２への関係が推定された。
従って、本発明の方法は、マイクロアレイのデータから推定されるネットワークとして既知の遺伝子ネットワークと比較して妥当な遺伝子ネットワークが推定できた。
本明細書で引用した全ての刊行物、特許及び特許出願は、そのまま参考として本明細書に取り入れるものとする。
産業上の利用可能性
本発明により、遺伝子ネットワークの推定方法が提供される。本発明の方法により得られたネットワークは、例えば医薬品を開発（ドラッグデザイン等）するために有用である。
【図面の簡単な説明】
図１は、本発明の推定装置のブロック図である。
図２は、遺伝子ネットワークを示す図である。
図３は、遺伝子の発現量と遺伝子調節との関係を示す図である。
図４は、本発明の推定方法を示す有向グラフの図である。
図５は、本発明の推定方法を示す有向グラフの図である。
図６は、冗長な因果関係を排除するプログラムを実行するためのフローチャートの図である。
図７は、本発明の方法により推定された遺伝子ネットワークを示す図である。
図８は、本発明の方法により推定された遺伝子ネットワークを示す図である。
符号の説明
１０１…ＣＰＵ、１０２…ＲＯＭ、１０３…ＲＡＭ、１０４…入力部、１０５…送信／受信部、１０６…表示部、１０７…ハードディスクドライブ、１０８…ＣＤ−ＲＯＭドライブ、１０９…ＣＤ−ＲＯＭ、１１０…パブリックデータベース
Claims (10)
- 複数の遺伝子のうち一の遺伝子の発現量を２つの条件下でそれぞれ発現させたときの他の遺伝子の発現量をそれぞれ検出し、得られる検出値の差を求め、当該差を指標として前記一の遺伝子と他の遺伝子との因果関係を導くことを特徴とする遺伝子ネットワークの推定方法。
- 複数の遺伝子のうち一の遺伝子を変異させたときの他の遺伝子の発現レベルと、当該一の遺伝子を変異させないときの当該他の遺伝子の発現レベルとをそれぞれ検出し、得られる検出値の差を求め、当該差を指標として前記一の遺伝子と他の遺伝子との因果関係を導くことを特徴とする遺伝子ネットワークの推定方法。
- さらに、得られた因果関係の矛盾の有無を検定することを含む、請求項１又は２記載の推定方法。
- 遺伝子の発現レベルが遺伝子のリン酸化の程度である請求項１又は２記載の推定方法。
- 遺伝子の変異が遺伝子の破壊である請求項２記載の推定方法。
- 遺伝子ネットワークが単細胞生物内のものである請求項１又は２記載の推定方法。
- 遺伝子ネットワークが多細胞生物内のものである請求項１又は２記載の推定方法。
- 複数の遺伝子のうち一の遺伝子の発現量を２つの条件下でそれぞれ発現させたときの他の遺伝子の発現量をそれぞれ検出する手段と、得られる検出値の差を求める手段と、当該差を指標として前記一の遺伝子と他の遺伝子との因果関係を導く手段とを含んでなる、遺伝子ネットワーク推定システム。
- さらに、得られた因果関係の矛盾の有無を検定することを含む、請求項８記載の推定システム。
- 複数の遺伝子のうち一の遺伝子を変異させたときの他の遺伝子の発現レベルと、当該一の遺伝子を変異させないときの当該他の遺伝子の発現レベルとをそれぞれ検出する手段と、得られる検出値の差を求める手段と、当該差を指標として前記一の遺伝子と他の遺伝子との因果関係を導く手段と、得られる因果関係の矛盾の有無を検定する手段とを含んでなる、コンピュータを遺伝子ネットワーク推定システムとして機能させるためのプログラムを記録したコンピュータ読み取り可能な記録媒体。
Applications Claiming Priority (3)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2000345982 | 2000-11-13 | ||
| JP2000345982 | 2000-11-13 | ||
| PCT/JP2001/009909 WO2002038749A1 (fr) | 2000-11-13 | 2001-11-13 | Procede d'anticipation de reseau de genes, systeme d'anticipation et support d'enregistrement |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| JPWO2002038749A1 true JPWO2002038749A1 (ja) | 2004-03-18 |
Family
ID=18819928
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP2002542065A Pending JPWO2002038749A1 (ja) | 2000-11-13 | 2001-11-13 | 遺伝子ネットワークの推定方法、推定システム及び記録媒体 |
Country Status (4)
| Country | Link |
|---|---|
| US (1) | US20040029154A1 (ja) |
| EP (1) | EP1342782A4 (ja) |
| JP (1) | JPWO2002038749A1 (ja) |
| WO (1) | WO2002038749A1 (ja) |
Families Citing this family (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP4617448B2 (ja) * | 2004-05-24 | 2011-01-26 | 独立行政法人科学技術振興機構 | 遺伝子ネットワークの推定方法及び推定プログラム |
| US7542854B2 (en) | 2004-07-22 | 2009-06-02 | International Business Machines Corporation | Method for discovering gene regulatory models and genetic networks using relational fuzzy models |
| JP5807336B2 (ja) * | 2011-02-08 | 2015-11-10 | 富士ゼロックス株式会社 | 情報処理装置および情報処理システム |
Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2001508303A (ja) * | 1997-01-13 | 2001-06-26 | デイヴィッド・エイチ・マック | 遺伝子機能同定のための発現モニタリング |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2000023933A1 (en) * | 1998-10-21 | 2000-04-27 | Bios Group Lp | Systems and methods for analysis of genetic networks |
-
2001
- 2001-11-13 EP EP01981092A patent/EP1342782A4/en not_active Withdrawn
- 2001-11-13 JP JP2002542065A patent/JPWO2002038749A1/ja active Pending
- 2001-11-13 US US10/416,406 patent/US20040029154A1/en not_active Abandoned
- 2001-11-13 WO PCT/JP2001/009909 patent/WO2002038749A1/ja active Application Filing
Patent Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2001508303A (ja) * | 1997-01-13 | 2001-06-26 | デイヴィッド・エイチ・マック | 遺伝子機能同定のための発現モニタリング |
Non-Patent Citations (2)
| Title |
|---|
| BIOINFOMATICS, vol. 16, no. 8, JPN4006004077, August 2000 (2000-08-01), pages 727 - 734, ISSN: 0000719857 * |
| BIOINFOMATICS, vol. 16, no. 8, JPNX007014471, August 2000 (2000-08-01), pages 727 - 734, ISSN: 0000833328 * |
Also Published As
| Publication number | Publication date |
|---|---|
| EP1342782A1 (en) | 2003-09-10 |
| WO2002038749A1 (fr) | 2002-05-16 |
| US20040029154A1 (en) | 2004-02-12 |
| EP1342782A4 (en) | 2006-04-12 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| de Koning et al. | Repetitive elements may comprise over two-thirds of the human genome | |
| Benidt et al. | SimSeq: a nonparametric approach to simulation of RNA-sequence datasets | |
| Markowetz et al. | Nested effects models for high-dimensional phenotyping screens | |
| US8332347B2 (en) | System and method for inferring a network of associations | |
| Forslund et al. | Evolution of protein domain architectures | |
| Carlborg et al. | Methodological aspects of the genetic dissection of gene expression | |
| Steenwyk et al. | Incongruence in the phylogenomics era | |
| US11322225B2 (en) | Systems and methods for determining effects of therapies and genetic variation on polyadenylation site selection | |
| JP2008546046A (ja) | マハラノビスの距離の遺伝的アルゴリズムの方法及びシステム | |
| Wang et al. | High-dimensional Bayesian network inference from systems genetics data using genetic node ordering | |
| CN115997255A (zh) | 从基因组预测细菌表型性状的分子技术 | |
| CN113488104A (zh) | 基于局部和全局的网络中心性分析的癌症驱动基因预测方法及系统 | |
| Pornputtapong et al. | KITSUNE: A tool for identifying empirically optimal K-mer length for alignment-free phylogenomic analysis | |
| US20200082910A1 (en) | Systems and Methods for Determining Effects of Genetic Variation of Splice Site Selection | |
| CN116959585B (zh) | 基于深度学习的全基因组预测方法 | |
| Lei et al. | An approach of gene regulatory network construction using mixed entropy optimizing context-related likelihood mutual information | |
| Lorena et al. | Evaluation of noise reduction techniques in the splice junction recognition problem | |
| JPWO2002038749A1 (ja) | 遺伝子ネットワークの推定方法、推定システム及び記録媒体 | |
| Anvar et al. | Interspecies translation of disease networks increases robustness and predictive accuracy | |
| CN116153396A (zh) | 一种基于迁移学习的非编码变异预测方法 | |
| Wang et al. | Network-adjusted Kendall’s Tau measure for feature screening with application to high-dimensional survival genomic data | |
| Sambhe et al. | Multi-objective classification and feature selection of Covid-19 proteins sequences using NSGA-II and MAP-elites | |
| JP4617448B2 (ja) | 遺伝子ネットワークの推定方法及び推定プログラム | |
| US11435357B2 (en) | System and method for discovery of gene-environment interactions | |
| Onami et al. | The DBRF method for inferring a gene network from large-scale steady-state gene expression data |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A711 | Notification of change in applicant |
Free format text: JAPANESE INTERMEDIATE CODE: A712 Effective date: 20031031 |
|
| RD03 | Notification of appointment of power of attorney |
Free format text: JAPANESE INTERMEDIATE CODE: A7423 Effective date: 20031216 |
|
| A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20040416 |
|
| A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20060307 |
|
| A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20060508 |
|
| A02 | Decision of refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A02 Effective date: 20070320 |