(57)【要約】【課題】マイクロアレイデータを用いた解析の精度および効率を高めるために、遺伝子の発現状況に関する真の信号と測定誤差の分離を行い、チャンネル毎に遺伝子が発現している確率を推定する。【解決手段】本発明の遺伝子発現状況推定システムは、マイクロアレイデータを入力する入力装置１と、プログラム制御により動作するデータ解析装置２と、出力装置３とを含む。データ解析装置２は、入力装置１から与えられた遺伝子発現強度データを用いて、混合正規分布のコンポーネント毎の分布パラメータおよび混合比パラメータを推定するパラメータ推定手段２１，２２と、推定された各パラメータを用いて、チャンネル毎に遺伝子の発現状況に関する事後確率を計算する事後確率計算手段２３を有し、計算された事後確率は出力装置３に出力する。【選択図】図２
2005038279.tif 000002
【特許請求の範囲】
【請求項１】
遺伝子発現強度データを送出する入力装置と、プログラム制御により動作するデータ解析装置と、出力装置とを含み、チャンネル毎に遺伝子が発現している確率を推定する遺伝子発現状況推定システムにおいて、
前記データ解析装置は、
前記入力装置から送出された遺伝子発現強度データを用いて、以下の数１に示される混合正規分布（尚、外１は平均０,分散σ２の１次元正規分布の密度関数を表し、外２−１及び外２−２はそれぞれ第１および第２のコンポーネントの平均と分散パラメータを表し,ξは混合率を表し、外３が満たされているものとする。）の分布パラメータを推定し、この推定された分布パラメータを送出する分布パラメータ推定手段と、
【数１】
2005038279.tif 000049
［外１］
2005038279.tif 000050
［外２−１］
2005038279.tif 000051
［外２−２］
2005038279.tif 000052
［外３］
2005038279.tif 000053
前記入力装置から送出された遺伝子発現強度データと、前記分布パラメータ推定手段から送出された分布パラメータとを用いて、混合正規分布の混合比パラメータを推定し、この推定された混合比パラメータを送出する混合比パラメータ推定手段と、
前記遺伝子発現強度データと、前記推定された分布パラメータ及び混合比パラメータを用いて、各チャンネル毎に各遺伝子の発現状況に関する事後確率を計算し、計算された事後確率を送出する事後確率計算手段を有し、
計算された前記事後確率が前記出力装置に出力されることを特徴とする遺伝子発現状況推定システム。
【請求項２】
前記分布パラメータ推定手段は、２つのチャンネルの遺伝子発現強度Ｘ，Ｙの差を示すＶが０近傍である遺伝子に関する発現量の和のデータに対して、２つのコンポーネントからなる前記混合正規分布を当てはめ、混合率（ξ）、平均（μ0，μ１）、及び分散外３−１を推定する
［外３−１］
2005038279.tif 000054
ことを特徴とする請求項１記載の遺伝子発現状況推定システム。
【請求項３】
前記遺伝子に関する発現量の和のデータは、遺伝子発現強度Ｘ，Ｙの差の絶対値外４の中央値をｃＭとするとき、外５で示される
［外４］
2005038279.tif 000055
［外５］
2005038279.tif 000056
ことを特徴とする請求項２記載の遺伝子発現状況推定システム。
【請求項４】
前記分布パラメータ推定手段は、推定された外６を用いて、外７を以下の数２、数３（尚、Ｎ０は外８を満たすデータのインデクッス集合とし、外９はその要素の個数とする。）、数４、及び数５に従って推定する
【数２】
2005038279.tif 000057
【数３】
2005038279.tif 000058
【数４】
2005038279.tif 000059
【数５】
2005038279.tif 000060
［外６］
2005038279.tif 000061
［外７］
2005038279.tif 000062
［外８］
2005038279.tif 000063
［外９］
2005038279.tif 000064
ことを特徴とする請求項３記載の遺伝子発現状況推定システム。
【請求項５】
前記混合比パラメータ推定手段は、前記入力装置から送出された遺伝子発現強度データ｛（ｕi，ｖi）｜ｉ＝１，．．．，ｎ｝に対して、前記分布パラメータ推定手段から与えられた外１０を用いて、以下の数６（ただし、以下の数７及び数８に示される関係を満たすとする。）に示される４つのコンポーネントからなる２変量混合正規分布を当てはめ、混合比パラメータｐ＝（ｐ00，ｐ10，ｐ01，ｐ11）（ｐ00は細胞１，２において、共に遺伝子が発現していない場合における混合比パラメータを表し、ｐ11は細胞１，２において、共に遺伝子が発現している場合における混合比パラメータを表し、ｐ10は細胞１において遺伝子が発現しており、細胞２において遺伝子が発現していない場合における混合比パラメータを表し、ｐ01は細胞１において遺伝子が発現していなく、細胞２において遺伝子が発現している場合における混合比パラメータを表す）を推定する
［外１０］
2005038279.tif 000065
【数６】
2005038279.tif 000066
【数７】
2005038279.tif 000067
【数８】
2005038279.tif 000068
ことを特徴とする請求項４記載の遺伝子発現状況推定システム。
【請求項６】
前記事後確率計算手段は、前記入力装置から送出された遺伝子発現強度データの対（ｕ，ｖ）に対して、細胞１および細胞２における任意の遺伝子が発現している事後確率を、以下の数９（尚、外１０−１は、以下の数１０によって与えられるものとする。）及び数１１（尚、τ１，τ２はそれぞれの細胞における遺伝子の真の発現の有無を１および０で表したものである。）にしたがって算出する
【数９】
2005038279.tif 000069
【数１０】
2005038279.tif 000070
【数１１】
2005038279.tif 000071
［外１０−１］
2005038279.tif 000072
ことを特徴とする請求項５記載の遺伝子発現状況推定システム。
【請求項７】
前記事後確率計算手段は、細胞１および細胞２で遺伝子の発現状態が異なる状態である事後確率を、以下の数１２（尚、τ１，τ２はそれぞれの細胞における遺伝子の真の発現の有無を１および０で表したものである。）にしたがって算出する
【数１２】
2005038279.tif 000073
ことを特徴とする請求項５記載の遺伝子発現状況推定システム。
【請求項８】
前記事後確率計算手段は、全ての遺伝子発現強度データの対（ｕ，ｖ）に基づいて、遺伝子が発現している事後確率を計算したかどうかを判定し、計算していれば終了し、計算していなければ、次の遺伝子に関する事後確率を計算し、
計算された各チャンネル毎の遺伝子が発現している事後確率は、前記出力装置へ送出され、
該出力装置は、各チャンネル毎の遺伝子が発現している事後確率を表示する
ことを特徴とする請求項６又は７記載の遺伝子発現状況推定システム。
【請求項９】
遺伝子発現強度データに基づいてチャンネル毎に遺伝子が発現している確率を推定する遺伝子発現状況推定方法において、
前記遺伝子発現強度データを用いて、以下の数１３に示される混合正規分布（尚、外１１は平均０,分散σ２の１次元正規分布の密度関数を表し、外１２−１及び外１２−２はそれぞれ第１および第２のコンポーネントの平均と分散パラメータを表し,ξは混合率を表し、外１３が満たされているものとする。）の分布パラメータを推定し、この推定された分布パラメータを送出するステップと、
【数１３】
2005038279.tif 000074
［外１１］
2005038279.tif 000075
［外１２−１］
2005038279.tif 000076
［外１２−２］
2005038279.tif 000077
［外１３］
2005038279.tif 000078
前記遺伝子発現強度データと、前記推定された分布パラメータとを用いて、混合正規分布の混合比パラメータを推定し、この推定された混合比パラメータを送出するステップと、
前記遺伝子発現強度データと、前記推定された分布パラメータ及び前記推定された混合比パラメータを用いて、各チャンネル毎に各遺伝子の発現状況に関する事後確率を計算し、計算された事後確率を送出するステップと、
計算された前記事後確率を出力するステップ
を有することを特徴とする遺伝子発現状況推定方法。
【請求項１０】
前記分布パラメータを推定するステップはさらに、２つのチャンネルの遺伝子発現強度Ｘ，Ｙの差を示すＶが０近傍である遺伝子発現強度データに対して、２つのコンポーネントからなる前記混合正規分布を当てはめ、混合率（ξ）、平均（μ0，μ１）、及び分散外１３−１を推定するステップを有する
［外１３−１］
2005038279.tif 000079
ことを特徴とする請求項９記載の遺伝子発現状況推定方法。
【請求項１１】
前記遺伝子に関する発現量の和のデータは、遺伝子発現強度Ｘ，Ｙの差の絶対値外１４−１の中央値をｃＭとするとき、外１４−２で示される
［外１４−１］
2005038279.tif 000080
［外１４−２］
2005038279.tif 000081
ことを特徴とする請求項１０記載の遺伝子発現状況推定方法。
【請求項１２】
前記分布パラメータを推定するステップにおいて、推定された外１５を用いて、外１６を以下の数１４、数１５（尚、Ｎ０は外１７を満たすデータのインデクッス集合とし、外１８はその要素の個数とする。）、数１６、及び数１７に従って推定する
【数１４】
2005038279.tif 000082
【数１５】
2005038279.tif 000083
【数１６】
2005038279.tif 000084
【数１７】
2005038279.tif 000085
［外１５］
2005038279.tif 000086
［外１６］
2005038279.tif 000087
［外１７］
2005038279.tif 000088
［外１８］
2005038279.tif 000089
ことを特徴とする請求項１１記載の遺伝子発現状況推定方法。
【請求項１３】
前記混合比パラメータを推定するステップにおいて、前記送出された遺伝子発現強度データ｛（ｕi，ｖi）｜ｉ＝１，．．．，ｎ｝に対して、前記送出された外１９を用いて、以下の数１８（ただし、以下の数１９及び数２０に示される関係を満たすとする。）に示される４つのコンポーネントからなる２変量混合正規分布を当てはめ、混合比パラメータｐ＝（ｐ00，ｐ10，ｐ01，ｐ11）（ｐ00は細胞１，２において、共に遺伝子が発現していない場合における混合比パラメータを表し、ｐ11は細胞１，２において、共に遺伝子が発現している場合における混合比パラメータを表し、ｐ10は細胞１において遺伝子が発現しており、細胞２において遺伝子が発現していない場合における混合比パラメータを表し、ｐ01は細胞１において遺伝子が発現していなく、細胞２において遺伝子が発現している場合における混合比パラメータを表す）を推定する
【数１８】
2005038279.tif 000090
【数１９】
2005038279.tif 000091
【数２０】
2005038279.tif 000092
［外１９］
2005038279.tif 000093
ことを特徴とする請求項１２記載の遺伝子発現状況推定方法。
【請求項１４】
前記事後確率を計算するステップでは、前記送出された遺伝子発現強度データの対（ｕ，ｖ）に対して、細胞１および細胞２における任意の遺伝子が発現している事後確率を、以下の数２１（尚、外１９−１は、以下の数２２によって与えられるものとする。）及び数２３（尚、τ１，τ２はそれぞれの細胞における遺伝子の真の発現の有無を１および０で表したものである。）にしたがって算出する
【数２１】
2005038279.tif 000094
【数２２】
2005038279.tif 000095
【数２３】
2005038279.tif 000096
［外１９−１］
2005038279.tif 000097
ことを特徴とする請求項１３記載の遺伝子発現状況推定方法。
【請求項１５】
前記事後確率を計算するステップでは、細胞１および細胞２で遺伝子の発現状態が異なる状態である事後確率を、以下の数２４（尚、τ１，τ２はそれぞれの細胞における遺伝子の真の発現の有無を１および０で表したものである。）にしたがって算出する
【数２４】
2005038279.tif 000098
ことを特徴とする請求項１３記載の遺伝子発現状況推定方法。
【請求項１６】
前記事後確率を計算するステップでは、全ての遺伝子発現強度データの対（ｕ，ｖ）に基づいて、遺伝子が発現している事後確率を計算したかどうかを判定し、計算していれば終了し、計算していなければ、次の遺伝子に関する事後確率を計算する
ことを特徴とする請求項１４又は１５記載の遺伝子発現状況推定方法。
【請求項１７】
遺伝子発現強度データに基づいてチャンネル毎に遺伝子が発現している確率を推定するためコンピュータに、
前記遺伝子発現強度データを用いて、以下の数２５に示される混合正規分布（尚、外２０は平均０,分散σ２の１次元正規分布の密度関数を表し、外２１−１及び外２１−２はそれぞれ第１および第２のコンポーネントの平均と分散パラメータを表し,ξは混合率を表し、外２２が満たされているものとする。）の分布パラメータを推定し、この推定された分布パラメータを送出するステップと、
【数２５】
2005038279.tif 000099
［外２０］
2005038279.tif 000100
［外２１−１］
2005038279.tif 000101
［外２１−２］
2005038279.tif 000102
［外２２］
2005038279.tif 000103
前記遺伝子発現強度データと、前記推定された分布パラメータとを用いて、混合正規分布の混合比パラメータを推定し、この推定された混合比パラメータを送出するステップと、
前記遺伝子発現強度データと、前記推定された分布パラメータ及び前記推定された混合比パラメータを用いて、各チャンネル毎に各遺伝子の発現状況に関する事後確率を計算し、計算された事後確率を送出するステップと、
計算された前記事後確率を出力するステップ
を実行させるための遺伝子発現状況推定プログラム。
【請求項１８】
遺伝子発現強度データに基づいてチャンネル毎に遺伝子が発現している確率を推定するためコンピュータに、
前記遺伝子発現強度データを用いて、以下の数２６に示される混合正規分布（尚、外２３は平均０,分散σ２の１次元正規分布の密度関数を表し、外２４−１及び外２４−２はそれぞれ第１および第２のコンポーネントの平均と分散パラメータを表し,ξは混合率を表し、外２５が満たされているものとする。）の分布パラメータを推定し、この推定された分布パラメータを送出するステップと、
【数２６】
2005038279.tif 000104
［外２３］
2005038279.tif 000105
［外２４−１］
2005038279.tif 000106
［外２４−２］
2005038279.tif 000107
［外２５］
2005038279.tif 000108
前記遺伝子発現強度データと、前記推定された分布パラメータとを用いて、混合正規分布の混合比パラメータを推定し、この推定された混合比パラメータを送出するステップと、
前記遺伝子発現強度データと、前記推定された分布パラメータ及び前記推定された混合比パラメータを用いて、各チャンネル毎に各遺伝子の発現状況に関する事後確率を計算し、計算された事後確率を送出するステップと、
計算された前記事後確率を出力するステップ
を実行させるための遺伝子発現状況推定プログラムを記録したコンピュータ読み取り可能な記録媒体。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、二色蛍光法によるｃＤＮＡマイクロアレイデータの統計的解析方法、解析システム、及び記録媒体に関し、特にチャンネル毎に遺伝子が発現している確率を推定するシステム、方法、プログラム、及び記録媒体に関するものである。
【背景技術】
【０００２】
現在、ゲノム研究は個々の遺伝子についての構造解析から体系的な遺伝子の機能解析へと展開しつつある。機能未知の遺伝子や総体としての遺伝子の機能解析のために、多数の遺伝子の発現強度を同時に定量化することのできるｃＤＮＡ（相補的なＤＮＡ）マイクロアレイを用いた実験はその有効性が大いに期待されている。
【０００３】
二色蛍光法によるｃＤＮＡマイクロアレイを用いた実験の目的は二種類の細胞の遺伝子発現強度の違いを検出することにある。ここで、二色蛍光法によるｃＤＮＡマイクロアレイの概要について述べる。まず、多数の遺伝子セットのｃＤＮＡを参照用のプローブとして、スライドグラス上にアレイ状に高密度に固定化する（マイクロアレイ）。
【０００４】
次に、条件の異なる２種類のサンプル、細胞１と細胞２（例えば正常細胞と癌細胞）から抽出したｍＲＮＡをそれぞれ波長の異なる蛍光色素でラベルし、ターゲットｃＤＮＡを合成する。そして、それらを等量混合したものをマイクロアレイに固定化された参照用のプローブｃＤＮＡに競合的にハイブリダイズさせる。ハイブリダイゼーション後、スキャナーでそれぞれの蛍光色素強度を測定する。細胞１にラベルされた蛍光色素をチャンネル1により、細胞２にラベルされた蛍光色素をチャンネル２により読み取り、それぞれを各細胞の遺伝子発現強度データ（マイクロアレイデータ）とする。
【０００５】
このように、マイクロアレイデータが得られるまでの過程は複雑であり、高度な実験技術が必要とされることから、実験の各段階において様々な実験誤差が生じると考えられる。このため、マイクロアレイデータから真に生物学的意味のあるデータを取り出すためには発現強度の分布と実験誤差の解析は解決すべき重要な課題である。
【０００６】
発現強度の分布に関しては、例えば、以下の非特許文献１を参照すると、Ｎｅｗｔｏｎ等は発現強度にガンマ分布関数を仮定し、発現強度比（チャンネル１とチャンネル２の発現強度の比）についての統計学的性質を考察している。
【０００７】
また、観測された発現強度データに対しては、例えば、以下の非特許文献２を参照すると、Ｌｅｅ等は真の発現強度を２個の水準値に分離できることおよび偶然誤差の存在を前提として、以下の数２７に示されるような混合正規分布を適用し、発現強度データについての統計学的考察を行っている。
【０００８】
【数２７】
2005038279.tif 000003
ここで、ｘはスキャナーなどによって得られる蛍光強度などの遺伝子発現強度（の対数値）を表し、右辺第１項の外２６は遺伝子が発現しているときの平均μ1, 分散外２７の正規分布、
［外２６］
2005038279.tif 000004
［外２７］
2005038279.tif 000005
また、同第２項の外２８は遺伝子が発現していないときの平均μ2、分散外２９の正規分布の密度関数を表し,ｐはその混合率を表す母数である。
【０００９】
［外２８］
2005038279.tif 000006
［外２９］
2005038279.tif 000007
【００１０】
実験誤差の解析については、系統誤差の除去、いわゆるノーマライゼーションの方法がいくつか提案されている。非特許文献３を参照すると、Ｃｈｅｎ等は二つの細胞の遺伝子発現強度の中央値は等しいとしてチャンネル１とチャンネル２で得られた測定値の補正を行っている。また，非特許文献４、５、６を参照すると、ＤｕｄoｉｔやＳｃｈｕcｈｈａｒｄｔやＹａｎｇは、系統誤差がスポットのスライドグラス上の位置や、二種類の蛍光色素の感度の違いによって生じたものと考え、それらを除去する方法を提案している。
【非特許文献１】Newton et. al、2001年、ジャーナル・オブ・コンピュテーショナル・バイオロジー、第8巻、37〜52頁（Journal of Computational Biology Vol. 8, pp. 37-52)
【非特許文献２】Lee et. al、2000年、プロシーディング・オブ・ザ・ナショナル・アカデミィー・オブ・サイエンシズ、第97巻、第18号、9834〜9839頁（Proceeding of the National Academy of Sciences Vol. 97, No 18, pp. 9834-9839）
【非特許文献３】Chen et. al、1997年、ジャーナル・オブ・バイオメディカル・オプティクス、第2号、364〜374頁（Journal of Biomedical Optics Vol. 2, pp. 364-374）
【非特許文献４】Dudoit et. al、2000. Statistical methods for identifying differentially expressed genes in replicated ｃＤＮＡmicroarray experiments. Technical~Report #578 2.
【非特許文献５】Schuchhardt et. al、2000年、ヌクレ・アシッド・リサーチ、第28巻、第10号（Nucleic Acids Research, 2000, Vol.28, No. 10）
【非特許文献６】Yang et. al、 2002年、ヌクレ・アシッド・リサーチ、第30巻、第4号（Nucleic Acids Research, 2002, Vol.30, No. 4）
【発明の開示】
【発明が解決しようとする課題】
【００１１】
上記した従来技術における問題点は、マイクロアレイデータの解析結果は再現性に乏しく解析の精度や効率が低いことにある。その原因として，従来の解析方法では，得られたマイクロアレイデータが真の信号と系統誤差および測定誤差に十分に分離できていないことが考えられる。したがって，系統誤差の除去および測定誤差の評価が重要な課題となる。
【００１２】
系統誤差の除去に関しては「ｃＤＮＡ マイクロアレイデータの補正法、システム、及び記録媒体」として，別途特許にて申請中である。本発明におけるマイクロアレイデータは、系統誤差はあらかじめ取り除かれているものとする。
【００１３】
従来のマイクロアレイデータを用いた解析では，各スポットのｃＤＮＡ量の均一性が保証されていないという理由で、２つのチャンネルの遺伝子発現強度の比（対数値の差）のみを扱っており、各チャンネル毎の遺伝子発現強度を扱っていない。従って、遺伝子の発現状況に関する真の信号と測定誤差の分離が十分に行われていない。
【００１４】
本発明の目的は、マイクロアレイデータを用いた解析の精度および効率を高めるために、遺伝子の発現に関する真の信号と測定誤差の分離を行い、さらに，各チャンネル毎に、遺伝子が発現している確率を推定する手法およびシステムを提供することにある。
【課題を解決するための手段】
【００１５】
本発明の遺伝子発現状況推定システムは、マイクロアレイデータを入力する入力装置と、プログラム制御により動作するデータ解析装置と、出力装置とを含む。前記データ解析装置は、前記入力装置から与えられた遺伝子発現強度データを用いて、混合正規分布のコンポーネント毎の分布パラメータおよび混合比パラメータを推定するパラメータ推定手段と、推定された各パラメータを用いて、各チャンネル毎に遺伝子の発現に関する事後確率を計算する事後確率計算手段を有し、計算された事後確率は前記出力装置に出力することを特徴とする。
【００１６】
このような構成を採用し、遺伝子発現状況を推定することにより、本発明の目的は達成することができる。
【発明の効果】
【００１７】
本発明による第１の効果は、マイクロアレイによって得られた遺伝子発現強度データに対して遺伝子の発現および非発現という概念を導入し，数理モデルを構築することにより遺伝子の発現に関する真の信号と実験誤差との分離が可能となったことにある。
【００１８】
本発明による第２の効果は、マイクロアレイによって得られた遺伝子発現強度データを、２つのチャンネルの遺伝子発現強度データの和と差のデータに変換することにより，２つのチャンネルの蛍光強度の感度情報が得られ易くなったことにある。この結果、感度の違いによる実験誤差の大きさを視覚的に示すことが可能となった。
【００１９】
本発明による第３の効果は、マイクロアレイによって得られた遺伝子発現強度データを、２つのチャンネルの遺伝子発現強度の和と差のデータに変換し，それらの２次元同時分布を記述することにより，２つのチャンネル毎に、各遺伝子の発現および非発現に関する事後確率を推定することが可能となったことにある。その結果、細胞１と細胞２の違いに関係する遺伝子を高い精度で検出することが可能となった。
【発明を実施するための最良の形態】
【００２０】
はじめに本発明におけるマイクロアレイによって得られた遺伝子発現強度データに対する数理モデルを説明する。Ｘをチャンネル１によって得られた細胞１の遺伝子発現強度とし、Ｙをチャンネル２によって得られた細胞２の遺伝子発現強度としたときに、それぞれの遺伝子発現データを数２８で示す。なお,ＸおよびＹは観察値に対して対数変換ないしは適当なべき変換および１次変換を含む適当な変換が施されている量である。
【００２１】
【数２８】
2005038279.tif 000008
【００２２】
ここで、τ１，τ２はそれぞれの細胞における遺伝子の真の発現の有無（ＯＮ／ＯＦＦ）を１および０で表す。また、αは遺伝子がＯＮの状態にある場合に産生されるｍＲＮＡの量およびスポットの状態に規定される遺伝子発現強度変量とし、βはチャンネル１およびチャンネル２に共通な測定誤差とし、εはチャンネル間で独立な測定誤差とする。なお、それぞれの確率変数の分布は数２９に従うものとする。
【００２３】
【数２９】
2005038279.tif 000009
【００２４】
ここで、Ｎ（μ，σ2）は平均μ，分散σ2の１次元正規分布とする。また、α，βおよびεはそれぞれ独立とする。この数理モデルにおいては、遺伝子が発現している場合（ＯＮ状態）の真の発現強度は非負値をとる確率変数であり、遺伝子が発現していない場合（ＯＦＦ状態）は、単なる実験誤差のみが観察されるとしている。さらに、非特許文献６を参照して、Ｙａｎｇ ＹＨらが導入したＭ−Ａ変換を変形した数３０で示すようなＳ−Ｄ変換を行う。
【００２５】
【数３０】
2005038279.tif 000010
【００２６】
すなわち、２つのチャンネルの遺伝子発現強度の和をＵとし、差をＶと変換する。このＳ−Ｄ変換を行ったときの本モデルの模式図を図１に示す。なお、以下ではこのプロットをＳ−Ｄプロットとよぶ。図１において、ｇ00はそれぞれの細胞において遺伝子が発現していない同時分布を表し、ｇ10は細胞１において遺伝子が発現しており、細胞２において遺伝子が発現していない同時分布を表し、ｇ01は細胞１において遺伝子が発現していなく、細胞２において遺伝子が発現している同時分布を表し、ｇ11はそれぞれの細胞において遺伝子が発現している同時分布を表している。ｇ00の密度関数を数３１に示す。
【００２７】
【数３１】
2005038279.tif 000011
ここで、φ（ｕ｜σ2）は平均０分散σ2の１次元正規分布の密度関数とする。ｇ10の密度関数を数３２に示す。
【００２８】
【数３２】
2005038279.tif 000012
ここで、φ２（ｕ，ｖ｜Σ）は平均ベクトル０、分散共分散行列Σの２次元正規分布の密度関数とし、Σ１０は２×２の分散共分散行列であり、数３３に示す。
【００２９】
【数３３】
2005038279.tif 000013
ｇ01の密度関数を数３４に示す。
【００３０】
【数３４】
2005038279.tif 000014
ここで、Σ０１は２×２の分散共分散行列であり、数３５に示す。
【００３１】
【数３５】
2005038279.tif 000015
ｇ11の密度関数を数３６に示す。
【００３２】
【数３６】
2005038279.tif 000016
以上の分布をもとに、細胞１および細胞２における遺伝子が発現している事後確率を数３７および数３８に示す。
【００３３】
【数３７】
2005038279.tif 000017
【数３８】
2005038279.tif 000018
ここで、ｆ（ｕ，ｖ｜p, θ）は数３９によって与えられるものとする。
【００３４】
【数３９】
2005038279.tif 000019
ただし、ただし、ｐ＝（ｐ00，ｐ10，ｐ01，ｐ11）は各分布の混合率を表すパラメータとする。
【００３５】
次に、本発明の実施の形態について図面を参照して詳細に説明する。図２を参照すると、本発明の第１の実施の形態は、遺伝子発現強度データに関する数理モデルの定式化およびそのデータ解析への適用による未知母数の推定する過程及び算出された母数の推定値を使用して、細胞１と細胞２の遺伝子の発現状態に関する事後推定を行うためのシステムである。同システムにはキーボード等の入力装置１と、プログラムの制御により動作するデータ解析装置２と、ディスプレイ装置や印刷装置等の出力装置３が含まれる。
【００３６】
データ解析装置２は、分布パラメータ推定手段２１と、混合比パラメータ推定手段２２と、事後確率計算手段２３とを備えている。分布パラメータ推定手段２１は、入力装置１から与えられた遺伝子発現強度データを用いて、混合正規分布の各コンポーネント毎の分布パラメータを推定する。推定された分布パラメータは、混合比パラメータ推定手段２２と、事後確率計算手段２３へ送られる。混合比パラメータ推定手段２２は、入力装置１から与えられた遺伝子発現強度データと、分布パラメータ推定手段２１から与えられた各コンポーネント毎の分布パラメータを用いて、混合正規分布の混合比パラメータを条件付き最尤法により推定する。推定された混合比パラメータは、事後確率計算手段２３へ送られる。事後確率計算手段２３は、入力装置１から与えられた遺伝子発現強度データと、分布パラメータ推定手段２１から与えられた各コンポーネント毎の分布パラメータと、混合比パラメータ推定手段２２から与えられた混合比パラメータを用いて、各チャンネル毎に遺伝子の発現状況に関する事後確率を計算する。計算された事後確率は出力装置３へ送られる。
【００３７】
次に、図２及び図３を参照して、遺伝子発現強度データに関する数理モデルの定式化およびそのデータ解析への適用による未知母数の推定する過程について詳細に説明する。入力装置１から与えられた遺伝子発現強度データ外３０は、分布パラメータ推定手段２１及び混合比パラメータ推定手段２２へ送られる。分布パラメータ推定手段２１は遺伝子発現強度の差の絶対値外３１の中央値をｃＭとするとき、Ｖ＝０近傍の遺伝子に関する発現量の和のデータ外３２に対して、以下の数４０に示す、２つのコンポーネントからなる混合正規分布を当てはめ、ξ，μ0，σ0，μ１，σ１を推定する（図３のステップＡ１）。
【００３８】
【数４０】
2005038279.tif 000020
［外３０］
2005038279.tif 000021
［外３１］
2005038279.tif 000022
［外３２］
2005038279.tif 000023
【００３９】
ここで、外３３は、平均０分散σ2の１次元正規分布の密度関数であり、外３４と外３５はそれぞれ第１および第２コンポーネントの平均と分散パラメータ、ξは混合率であり、外３６が満たされているものとする。
【００４０】
次に、分布パラメータ推定手段２１は、推定された外３７を用いて、外３８を以下の数４１、数４２、数４３、及び数４４に従って推定する（ステップＡ２）。
【００４１】
［外３３］
2005038279.tif 000024
［外３４］
2005038279.tif 000025
［外３５］
2005038279.tif 000026
［外３６］
2005038279.tif 000027
［外３７］
2005038279.tif 000028
［外３８］
2005038279.tif 000029
【数４１】
2005038279.tif 000030
【数４２】
2005038279.tif 000031
【数４３】
2005038279.tif 000032
【数４４】
2005038279.tif 000033
【００４２】
ここで、Ｎ０は外３９を満たすデータのインデックス集合とし、外４０はその要素の個数とする。
【００４３】
［外３９］
2005038279.tif 000034
［外４０］
2005038279.tif 000035
【００４４】
次に、混合比パラメータ推定手段２２は、入力装置１から与えられた遺伝子発現強度データ｛（ｕi，ｖi）｜ｉ＝１，．．．，ｎ｝に対して、分布パラメータ推定手段２１から与えられた各パラメータの推定値外４１を用いて、以下の数４５（尚、以下の数４６（ただし、外４２は数３３から導かれる２×２の分散共分散行列）及び数４７（ただし、外４３は数３５から導かれる２×２の分散共分散行列）に示される関係を満たすとする。）に示される４つのコンポーネントからなる２変量混合正規分布を当てはめ、混合比パラメータｐ＝（ｐ00，ｐ10，ｐ01，ｐ11）を条件付き最尤法により推定する（ステップＡ３）。
【００４５】
［外４１］
2005038279.tif 000036
【数４５】
2005038279.tif 000037
【数４６】
2005038279.tif 000038
【数４７】
2005038279.tif 000039
［外４２］
2005038279.tif 000040
［外４３］
2005038279.tif 000041
【００４６】
次に、算出された母数の推定値を使用して、細胞１と細胞２の遺伝子の発現状態に関する事後推定を行う過程を説明する。
【００４７】
各細胞の遺伝子毎の発現状況に関する事後確率計算手段２３は、入力装置１から与えられた遺伝子発現強度データの対（ｕ，ｖ）に対して、分布パラメータ推定手段２１及び混合比パラメータ推定手段２２から与えられた各パラメータの推定値外４４および外４５を用いて記述できる。
【００４８】
［外４４］
2005038279.tif 000042
［外４５］
2005038279.tif 000043
【００４９】
すなわち、細胞１および細胞２における任意の遺伝子の発現がＯＮになっている事後確率は、以下の数４８及び数４９より算出できる（ステップＡ４）。
【００５０】
【数４８】
2005038279.tif 000044
【数４９】
2005038279.tif 000045
【００５１】
全ての遺伝子発現強度データの対（ｕ，ｖ）に基づいて、遺伝子の発現がＯＮになっている事後確率を計算したかどうかを判定し（ステップＡ５）、計算していれば終了する。計算していなければ、次の遺伝子に関する事後確率を計算する。
【００５２】
計算された各チャンネル毎の遺伝子が発現している事後確率は、出力装置３へ送る。出力装置３は、各チャンネル毎の遺伝子が発現している事後確率をグラフで表示したり、印刷したりする。
【００５３】
次に、本実施の形態の効果について説明する。本実施の形態では、マイクロアレイで得られた遺伝子発現データに対して、遺伝子の発現および非発現を導入した数理モデルを構築し、真の信号と実験誤差との分離を行った。また、２つのチャンネルの遺伝子発現強度による和と差のデータを用いることにより、マイクロアレイデータでの各チャンネルの蛍光強度の感度情報が得られ易くなり、実験誤差の大きさをより的確に抽出することが可能となった。さらに、これらの和と差のデータに対して２次元同時分布を記述することにより、チャンネル毎に各遺伝子の発現に関して、高い精度の事後確率を推定することが可能となった。
【００５４】
また、細胞１と細胞２で遺伝子の発現状態が異なるという事象(ＯＮ−ＯＦＦが不一致である状態)（図３のステップＡ４）の事後確率は以下の数５０により算出される。
【００５５】
【数５０】
2005038279.tif 000046
【００５６】
本実施の形態の効果として、細胞１と細胞２において発現状態が異なる可能性の高い遺伝子の候補を検出することが可能となる。
【００５７】
次に、本発明の第２の実施の形態について図面を参照して詳細に説明する。
【００５８】
図４を参照すると、本発明の第２の実施形態は、本発明の第１の実施形態と同様に、入力装置、データ解析装置、出力装置を備え、更に、データ解析プログラムを記録した記録媒体４を備える。この記録媒体４は可搬形あるいは固定型のいずれであってもよく、磁気ディスク、半導体メモリ、ＣＤ−ＲＯＭその他の記録媒体であってもよい。また、本手法を実行できるコンピュータプログラムを、ネットワークに接続されたコンピュータの記録装置に格納しておき、ネットワークを介して他のコンピュータに転送することもできる。本アルゴリズムを実行するコンピュータプログラムを提供する提供媒体としては、様々な形式のコンピュータに読み出し可能な媒体として頒布可能であって、特定のタイプの媒体に限定されるものではない。
【００５９】
データ解析プログラムは記録媒体４からデータ解析装置５に読み込まれ、データ解析装置５の動作を制御し、入力装置１から入力されたデータファイルに対して第１の実施の形態におけるデータ処理装置２による処理と同一の処理を実行する。
【実施例１】
【００６０】
次に、本発明の実施例について説明する。例として用いたデータは，異なる２種類の癌細胞（細胞１、細胞２）の遺伝子発現状況の比較のために行われた実験から得られたものである。
【００６１】
一枚のチップ上に４８グリッド，１グリッドあたり４４１（２１×２１）スポット，計２１１６８の遺伝子の発現パターンについて調べたものである。
【００６２】
図５および図６に、細胞１および細胞２の遺伝子の発現状態が両者ともにＯＦＦ、あるいは両者ともにＯＮである場合(Ｖ＝０)、それらの発現強度の和の分布Ｕに混合正規分布と、その対比として単峰の正規分布を当てはめた場合の各分布の推定結果を示す。以下の表１にコンポーネント毎の分布パラメータの推定結果を示す（図３のステップＡ１の結果）。
【００６３】
【表１】
2005038279.tif 000047
【００６４】
図５に推定された累積分布関数を、図６に推定された密度関数を示す。細実線は混合正規分布を仮定した場合、二点鎖線はその第１コンポーネント（ＯＦＦ−ＯＦＦ）、太実線は第２コンポーネント（ＯＮ−ＯＮ）、そして点線は単峰の正規分布を仮定した場合の推定結果を示している。
【００６５】
図５における一点鎖線は観察されたデータに基づく経験累積分布関数を示しており、観察値が推定された混合正規分布（細実線）によく従っていることを示している。
【００６６】
図６において＊印（尚、左右両端は＊印が確認できるが、遺伝子発現強度が０付近から３０付近までは重なり合っているため＊印ではなく以下の（１）〜（５）に示すハッチングで表示することとした。）は観察データを示しており、それぞれのハッチング領域（（１）〜（５））は、第１コンポーネントに属している事後確率の大きさを示している。（１）の黒塗り領域は０〜０．２、（２）のハッチング領域は０．２〜０．４、（３）のハッチング領域は０．４〜０．６、（４）のハッチング領域は０．６〜０．８、（５）のハッチング領域は０．８〜１．０を表している。
【００６７】
図７に、遺伝子発現強度データのＳ−Ｄプロットを示す。横軸は細胞１と細胞２の遺伝子発現強度の対数値の和を、縦軸はそれらの差を示している。各マークの色は、細胞１と細胞２の遺伝子の発現状態が不一致（ＯＮ−ＯＦＦあるいはＯＦＦ−ＯＮ）である事後確率の大きさを示している。（１）の黒塗り領域は０〜０．２、（２）のハッチング領域は０．２〜０．４、（３）のハッチング領域は０．４〜０．６、（４）のハッチング領域は０．６〜０．８、（５）のハッチング領域は０．８〜１．０を表している。
【００６８】
以下の表２に、条件付き最尤法による混合比パラメータの推定結果を示す（図３のステップＡ２，Ａ３の結果）。図７における各プロットに対しては、それぞれ遺伝子が対応しているので、細胞１と細胞２の違いに関係する遺伝子の候補を容易に絞り込むことができる。
【００６９】
【表２】
2005038279.tif 000048
【図面の簡単な説明】
【００７０】
【図１】本発明における数理モデルのＳ−Ｄプロットによる模式図である。
【図２】本発明の第１の実施の形態の構成を示すブロック図である。
【図３】本発明の第１の実施の形態の動作を示すフローチャートである。
【図４】本発明の第２の実施形態の構成を示すブロック図である。
【図５】Ｖ＝０近傍の遺伝子発現強度データに対して推定した正規分布の累積分布グラフである。
【図６】Ｖ＝０近傍の遺伝子発現強度データに対して推定した正規分布の密度関数を示したグラフである。
【図７】遺伝子発現強度データのＳ−Ｄプロットを示した図である。
【符号の説明】
【００７１】
１  入力装置
２  データ解析装置
３  出力装置
４  記録媒体
５  データ解析装置
２１  分布パラメータ推定手段
２２  混合比パラメータ推定手段
２３  事後確率計算手段
【図１】
2005038279.tif 000109
【図２】
2005038279.tif 000110
【図３】
2005038279.tif 000111
【図４】
2005038279.tif 000112
【図５】
2005038279.tif 000113
【図６】
2005038279.tif 000114
【図７】
2005038279.tif 000115