(57)【要約】  【課題】  バイアスを含ませることなく打切りサンプルを処理して精度の高い生存時間解析モデルの構築を行えるようにする。  【解決手段】  ランダムさを持って実際の打切り時間に対応する仮イベント発生時間を仮設定し、生体状態の変化データを生存時間解析に使用する生命表へ適用する場合に、イベントサンプルのイベント発生時間と打切りサンプルの仮イベント発生時間の両方のイベント発生時間でのイベント非発生確率に基づいて目的変数を算出し、遺伝子発現のデータ、細胞内物質のデータの少なくとも一方のデータを説明変数として第１次生存時間解析モデルを最適化することにより重要な説明変数を選抜する。【選択図】  図２
2005222422.tif 000002
【特許請求の範囲】
【請求項１】
時間の経過とともに確率的に発生する生体状態の変化として所定のイベントを定義し、そのイベントが起こったイベントサンプルと打切りサンプルを含む生体状態の変化データを使用したコンピュータによるデータ解析方法であって、
打切りサンプルの実際には観測されなかったイベント発生時間を、統計的にランダムさを持って実際の打切り時間に対応する仮イベント発生時間に仮設定する打切りサンプル仮設定処理と、
生体状態の変化データを生存時間解析に使用する生命表へ適用する場合に、イベントサンプルのイベント発生時間と打切りサンプルの仮イベント発生時間の両方のイベント発生時間でのイベント非発生確率を算出し、その算出されたイベント非発生確率に基づいて個々のサンプルの相対的確率を算出し、その相対的確率を仮定した分布に基づいた変換又は仮定を前提としない変換をし、その変換結果を目的変数とし、
遺伝子発現のデータ、細胞内物質のデータの少なくとも一方のデータを説明変数とし、
前記目的変数と説明変数に基づいて、第１次生存時間解析モデルの第１次選抜基準を算出しつつ、所定選抜方法によって第１次生存時間解析モデルを最適化することにより重要な説明変数を選抜する第１次説明変数群選抜処理と、を含むことを特徴とする、データ解析方法。
【請求項２】
請求項１に記載のデータ解析方法において、前記打切りサンプル仮設定処理と第１次説明変数群選抜処理とを少なくとも含んで一組として、その一組を複数回繰り返す繰返し処理を行い、その結果得られた第１次説明変数群に基づいて重要な説明変数を選抜する、データ解析方法。
【請求項３】
請求項２に記載のデータ解析方法において、繰返し処理によって得られた第１次説明変数群から重複する説明変数を選択することにより、重要な説明変数を選抜する処理を含む、データ解析方法。
【請求項４】
請求項２ないし請求項３のいずれか一項に記載のデータ解析方法において、繰返し処理の際に、遺伝子発現のデータ、細胞内物質のデータの一部のデータのみ使用する、データ解析方法。
【請求項５】
請求項２ないし請求項４のいずれか一項に記載のデータ解析方法において、繰返し処理の際に、既に選抜された遺伝子発現のデータ、細胞内物質のデータを、その後の説明変数の候補として採用しないようにする、データ解析方法。
【請求項６】
請求項１ないし請求項５のいずれか一項に記載のデータ解析方法において、打切りサンプルの仮イベント発生時間をイベントサンプル及び打切りサンプルの最大観測時間を考慮した基準観測時間に基づいて乱数を用いて仮設定する、データ解析方法。
【請求項７】
請求項６に記載のデータ解析方法において、打切りサンプルの仮イベント発生時間を、乱数を用いて、実際の打切り時間を下限とし、基準観測時間を上限とする所定時間に設定する、データ解析方法。
【請求項８】
請求項１ないし請求項７のいずれか一項に記載のデータ解析方法において、前記第１次生存時間解析モデルが部分最小自乗モデルである、データ解析方法。
【請求項９】
請求項１ないし請求項８のいずれか一項に記載のデータ解析方法において、生体状態の変化データと、少なくとも選抜された第１次説明変数群に基づいて、第２次生存時間解析モデルの第２次選抜基準を算出しつつ、所定選抜方法によって第２次生存時間解析モデルを最適化する第２次生存時間解析モデル最適化処理を含む、データ解析方法。
【請求項１０】
請求項９に記載のデータ解析方法において、選抜された第１次説明変数群のみを第２次生存時間解析モデル構築の際の説明変数として採用する、データ解析方法。
【請求項１１】
請求項９に記載のデータ解析方法において、選抜された第１次説明変数群のみならず、他の説明変数を入れたものを第２次生存時間解析モデル構築の際の説明変数として採用する、データ解析方法。
【請求項１２】
請求項９ないし請求項１１のいずれか一項に記載のデータ解析方法において、最適化された第１次生存時間解析モデルの潜在変数を用いて第２次生存時間解析モデルを最適化する、データ解析方法。
【請求項１３】
請求項９ないし請求項１２のいずれか一項に記載のデータ解析方法において、前記第２次生存時間解析モデルが回帰モデルである、データ解析方法。
【請求項１４】
請求項１３に記載のデータ解析方法において、前記回帰モデルが比例ハザードモデル又は加速モデルである、データ解析方法。
【請求項１５】
請求項１ないし請求項１４のいずれか一項に記載のデータ解析方法において、選抜規準が残差平方和、相関係数、予測残差平方和、予測相関係数、赤池情報規準(AIC)、最小記述長(MDL)に基づいて算出される、データ解析方法。
【請求項１６】
予測対象サンプルの遺伝子発現のデータ、細胞内物質のデータの少なくとも一方のデータを入力して、請求項１ないし請求項８のいずれか一つに記載の発明によって決定された第１次生存時間解析モデルとその選抜された第１次説明変数群に基づいて生存時間解析に関する予測出力を得ることを特徴とする、データ解析方法。
【請求項１７】
予測対象サンプルの遺伝子発現のデータ、細胞内物質のデータの少なくとも一方のデータを入力して、請求項９ないし請求項１５のいずれか一つに記載の発明によって決定された第２次生存時間解析モデルとその選抜された第２次説明変数群に基づいて生存時間解析に関する予測出力を得ることを特徴とする、データ解析方法。
【請求項１８】
時間の経過とともに確率的に発生する生体状態の変化として所定のイベントを定義し、そのイベントが起こったイベントサンプルと打切りサンプルを含む生体状態の変化データを使用したコンピュータによるデータ解析システムであって、
打切りサンプルの実際には観測されなかったイベント発生時間を、統計的にランダムさを持って実際の打切り時間に対応する仮イベント発生時間に仮設定する打切りサンプル仮設定手段と、
生体状態の変化データを生存時間解析に使用する生命表へ適用する場合に、イベントサンプルのイベント発生時間と打切りサンプルの仮イベント発生時間の両方のイベント発生時間でのイベント非発生確率を算出し、その算出されたイベント非発生確率に基づいて個々のサンプルの相対的確率を算出し、その相対的確率を仮定した分布に基づいた変換又は仮定を前提としない変換をし、その変換結果を目的変数とし、
遺伝子発現のデータ、細胞内物質のデータの少なくとも一方のデータを説明変数とし、
前記目的変数と説明変数に基づいて、第１次生存時間解析モデルの第１次選抜基準を算出しつつ、所定選抜方法によって第１次生存時間解析モデルを最適化することにより重要な説明変数を選抜する第１次説明変数群選抜手段と、を含むことを特徴とする、データ解析システム。
【請求項１９】
請求項１８に記載のデータ解析システムにおいて、前記打切りサンプル仮設定手段が行う処理と第１次説明変数群選抜手段が行う処理とを少なくとも含んで一組として、その一組を複数回繰り返す繰返し処理手段とを備え、その結果得られた第１次説明変数群に基づいて重要な説明変数を選抜する、データ解析システム。
【請求項２０】
請求項１９に記載のデータ解析システムにおいて、繰返し処理によって得られた第１次説明変数群から重複する説明変数を選択することにより、重要な説明変数を選抜する重複選抜手段を含む、データ解析システム。
【請求項２１】
請求項１９ないし請求項２０のいずれか一項に記載のデータ解析システムにおいて、繰返し処理の際に、遺伝子発現のデータ、細胞内物質のデータの一部のデータのみ使用する繰り返し制御手段を含む、データ解析システム。
【請求項２２】
請求項１９ないし請求項２１のいずれか一項に記載のデータ解析システムにおいて、繰返し処理の際に、既に選抜された遺伝子発現のデータ、細胞内物質のデータを、その後の説明変数の候補として採用しないようにする繰り返し制御手段を含む、データ解析システム。
【請求項２３】
請求項１８ないし請求項２２のいずれか一項に記載のデータ解析システムにおいて、前記打切りサンプル仮設定手段が、打切りサンプルの仮イベント発生時間をイベントサンプル及び打切りサンプルの最大観測時間を考慮した基準観測時間に基づいて乱数を用いて仮設定する、データ解析システム。
【請求項２４】
請求項２３に記載のデータ解析システムにおいて、前記打切りサンプル仮設定手段が、打切りサンプルの仮イベント発生時間を、乱数を用いて、実際の打切り時間を下限とし、基準観測時間を上限とする所定時間に設定する、データ解析システム。
【請求項２５】
請求項１８ないし請求項２４のいずれか一項に記載のデータ解析システムにおいて、前記第１次生存時間解析モデルが部分最小自乗モデルである、データ解析システム。
【請求項２６】
請求項１８ないし請求項２５のいずれか一項に記載のデータ解析システムにおいて、生体状態の変化データと、少なくとも選抜された第１次説明変数群に基づいて、第２次生存時間解析モデルの第２次選抜基準を算出しつつ、所定選抜方法によって第２次生存時間解析モデルを最適化する第２次生存時間解析モデル最適化手段とを含む、データ解析システム。
【請求項２７】
請求項２６に記載のデータ解析システムにおいて、選抜された第１次説明変数群のみを第２次生存時間解析モデル構築の際の説明変数として採用する、データ解析システム。
【請求項２８】
請求項２６に記載のデータ解析システムにおいて、選抜された第１次説明変数群のみならず、他の説明変数を入れたものを第２次生存時間解析モデル構築の際の説明変数として採用する、データ解析システム。
【請求項２９】
請求項２６ないし請求項２８のいずれか一項に記載のデータ解析システムにおいて、第２次生存時間解析モデル最適化手段が、最適化された第１次生存時間解析モデルの潜在変数を用いて第２次生存時間解析モデルを最適化する、データ解析システム。
【請求項３０】
請求項２６ないし請求項２９のいずれか一項に記載のデータ解析システムにおいて、前記第２次生存時間解析モデルが回帰モデルである、データ解析システム。
【請求項３１】
請求項３０に記載のデータ解析システムにおいて、前記回帰モデルが比例ハザードモデル又は加速モデルである、データ解析システム。
【請求項３２】
請求項１８ないし請求項３１のいずれか一項に記載のデータ解析システムにおいて、選抜規準が残差平方和、相関係数、予測残差平方和、予測相関係数、赤池情報規準(AIC)、最小記述長(MDL)に基づいて算出される、データ解析システム。
【請求項３３】
予測対象サンプルの遺伝子発現のデータ、細胞内物質のデータの少なくとも一方のデータを入力する入力手段と、請求項１８ないし請求項２５のいずれか一つに記載の発明によって決定された第１次生存時間解析モデルとその選抜された第１次説明変数群に基づいて生存時間解析に関する予測出力を得る予測手段と、を有することを特徴とする、データ解析システム。
【請求項３４】
予測対象サンプルの遺伝子発現のデータ、細胞内物質のデータの少なくとも一方のデータを入力する入力手段と、請求項２６ないし請求項３２のいずれか一つに記載の発明によって決定された第２次生存時間解析モデルと、その選抜された第２次説明変数群に基づいて生存時間解析に関する予測出力を得る予測手段と、を有することを特徴とする、データ解析システム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明はデータ解析方法及びそのシステムに関し、特に生体状態の変化と、遺伝子発現データ、細胞内物質データの少なくとも一方を利用するデータ解析方法及びデータ解析システムに関する。
【背景技術】
【０００２】
２０００年６月のヒトゲノムの解読宣言以降、ゲノムに書かれた遺伝情報がどのように発現して機能しているかを解明するポストゲノム時代に突入したと言われている。ヒトゲノム計画の進展の中で、ゲノム発現状態を測定する方法論も進展してきた。
トランスクリプトーム（ｍＲＮＡの総称）測定手段としてオリゴヌクレオチドアレイやマイクロチップが知られている。またプロテオーム（蛋白質）測定手段として、以前からある２次元電気泳動に加えて、最近では質量分析の方法が進歩してきた。また抗体チップなどの先進の技術も注目されている。これらの測定技術は、生体の状態パラメータを短時間に一挙に測定できることがそれまでの技術と比較して画期的であるといえる。
【０００３】
遺伝子発現状態を効率的に測定する技術として次のものがあげられる。
トランスクリプトームを特定するものとして、基盤に複数種のＤＮＡを担持し、それに相補的なｍＲＮＡを検出するＤＮＡチップが知られている。代表的なＤＮＡチップには、遺伝子チップやＤＮＡマイクロアレイがある。また、プロテオーム（蛋白質の総体）を特定するものには、２次元電気泳動、抗体チップ、質量スペクトルを用いるものがある。またメタボローム（代謝中間体を含めた代謝産物の総体）を測定する手法も質量分析などによって試みられており、進展が見られる。
【０００４】
生体内の細胞の状態は遺伝子産物の発現によってよく記述されるため、従来の診断マーカーでは情報が不足している場面でも、精度のより高い診断が可能になると期待される。たとえば、非特許文献１において、P.0.Brownらは、ＤＮＡチップによってびまん性巨大Ｂ細胞リンパ種患者の細胞のトランスクリプトームを測定し、クラスター解析によってより悪性度の高いものとそうでないものとを別クラスターに分離した。しかし、この方法は分類分けしたにすぎず、どの遺伝子がどの程度、寄与することによって比較的良性のものと悪性の高いものに区分けされるのかという重要な情報は得られない。
【０００５】
また、特許文献１に示すように、D.M  RockeとD.V.  Nguyenによって遺伝子発現量のような膨大な数の説明変数であっても、部分最小自乗法が次元圧縮とモデルフィットを同時に行えることに着目して、ＤＮＡチップの解析に部分最小自乗法を用いる研究が公開されている。
【０００６】
一方、生存時間解析の説明変数に遺伝子産物の発現量を用いるというより直接的な研究も発表されるに至っている。生存時間解析は臨床の多くの場面で重要とされる情報処理であり、ポストゲノム科学の成果が応用されることの意義は大きい。ここで生存時間解析とは臨床医学などを中心に発展してきた統計学の一分野であり、時間を経過する中で確率的に発生するイベント(event)を取扱うものである。イベントは必ずしも死亡に限定されるものではなく、疾病の発症、再発、など多くの生命現象をカバーする用語である。生存時間解析の代表的な方法としてカプラン・マイヤー法（Kaplan-Meier）と比例ハザードモデルなどを挙げることができる。これらのモデルでは、いわゆる、打切り(censor)サンプルも統計的に整合性をもって処理される。
【０００７】
生存時間解析の説明変数に遺伝子産物の発現量を用いた従来技術として、本発明者らが先に提案した下記特許文献２が挙げられる。
この特許文献２は、カプラン・マイヤー法（Kaplan-Meier）で求めた生存率から得た生存率をlogit変換したものを目的変数とし、部分最小自乗法の手法を用いたモデル構築によって、最終的には比例ハザードモデルに用いる説明変数を効率的に絞り込むことができることを示した。これはＤＮＡチップなどで測定される遺伝子産物の発現量が膨大な数に上り、遺伝子産物の発現量を解析することによって病気の種類、進行度合などを判断できる可能性を高める上で、非常に意義深いことである。
【０００８】
【特許文献１】国際公開  WO  00/070340
【特許文献２】国際公開  WO  03/085548
【非特許文献１】Nature  403(3),  503-11  (2000)
【発明の開示】
【発明が解決しようとする課題】
【０００９】
しかしながら、上記特許文献２に開示された解析技術でも更なる改良を目指す点において、下記のような検討すべき課題がある。
打切りサンプルとは、死亡などのイベントに至ることなく、観測を打ち切ったサンプルのことであり、例えて言うと、イベント発生時間が値ではなく不等式によって制限されているようなサンプルのことである。あるいはより一般にはある一定の制約条件によって表現されるサンプルを一般に打切りサンプルと呼ぶ。
【００１０】
生存時間解析において打切りサンプルをどのように扱うかによって結果的にバイアスが入る危険性がある。打切りサンプルにおいて、その後、潜在的に観測されるイベント発生に関してなんらかの情報を持っている状態で打切りが行われた場合、例えば、極端な場合、患者が死亡しそうになったので打切りにしてしまうような打切りであると、バイアスが入ってしまうことになる。
また、打切りサンプルをイベントサンプルから除外して処理することは、それ自体がバイアスを入れることになってしまい、生存時間解析では好ましいことではないとされている。
【００１１】
つまり、遺伝子発現量のようにサンプル数に比べて極めて多量の説明変数を解析する場合、部分最小自乗法などの相関モデルを用いた従来技術であっても、打切りサンプルの扱いに着目し、かつ打切りサンプルをどのように扱って生存時間解析によるモデル構築を行えばよいかについてその検証を踏まえた具体的な改善方法の提案は報告されていないのである。
【００１２】
本発明は上記課題に鑑みてなされたものであり、本発明の目的は、上記課題を解決できる、データ解析方法及びデータ解析システムを提供することにある。
具体的な目的の一例を示すと、以下の通りである。
（ａ）遺伝子発現の量、細胞内物質の量のようにサンプル数に比べて極めて多量の説明変数を扱って生存時間解析する場合でも、バイアスを含ませる可能性を極力低くして打切りサンプルを処理して生存時間解析モデルの構築を行えるようにする。
（ｂ）生存時間解析において、遺伝子発現の量、細胞内物質の量のようにサンプル数に比べて極めて多量の説明変数を扱って生存時間解析する場合でも、重要な説明変数を漏らすことなく、多様性を持った説明変数をセットできる打切りサンプルの処理方法を提供する。
なお、上記に記載した以外の発明の課題及びその解決手段は、後述する明細書内の記載において詳しく説明する。
【課題を解決するための手段】
【００１３】
本発明は多面的に表現できるが、例えば、代表的なものを挙げると、次のように構成したものである。なお、下記第１発明〜第３４発明は代表的なものにすぎず、本明細書の他の記載によって他の発明の構成及び効果を示すこととする。
第１発明のデータ解析方法は、時間の経過とともに確率的に発生する生体状態の変化として所定のイベントを定義し、そのイベントが起こったイベントサンプルと打切りサンプルを含む生体状態の変化データを使用したコンピュータによるデータ解析方法であって、
打切りサンプルの実際には観測されなかったイベント発生時間を、統計的にランダムさを持って実際の打切り時間に対応する仮イベント発生時間に仮設定する打切りサンプル仮設定処理と、
生体状態の変化データを生存時間解析に使用する生命表へ適用する場合に、イベントサンプルのイベント発生時間と打切りサンプルの仮イベント発生時間の両方のイベント発生時間でのイベント非発生確率を算出し、その算出されたイベント非発生確率に基づいて個々のサンプルの相対的確率を算出し、その相対的確率を仮定した分布に基づいた変換又は仮定を前提としない変換をし、その変換結果を目的変数とし、
遺伝子発現のデータ、細胞内物質のデータの少なくとも一方のデータを説明変数とし、
前記目的変数と説明変数に基づいて、第１次生存時間解析モデルの第１次選抜基準を算出しつつ、所定選抜方法によって第１次生存時間解析モデルを最適化することにより重要な説明変数を選抜する第１次説明変数群選抜処理と、を含むことを特徴とする。
【００１４】
第２発明は、第１発明において、前記打切りサンプル仮設定処理と第１次説明変数群選抜処理とを少なくとも含んで一組（一括り）として、その一組を複数回繰り返す繰返し処理を行い、その結果得られた第１次説明変数群に基づいて重要な説明変数を選抜することを特徴とする。
第３発明は、第２発明において、繰返し処理によって得られた第１次説明変数群から重複する説明変数を選択することにより、重要な説明変数を選抜する処理を含むことを特徴とする。
【００１５】
第４発明は、第２発明ないし第３発明のいずれか一つに記載の発明において、繰返し処理の際に、遺伝子発現のデータ、細胞内物質のデータの一部のデータのみ使用することを特徴とする。
第５発明は、第２発明ないし第４発明のいずれか一つに記載の発明において、繰返し処理の際に、既に選抜された遺伝子発現のデータ、細胞内物質のデータを、その後の説明変数の候補として採用しないようにすることを特徴とする。
【００１６】
第６発明は、第１発明ないし第５発明のいずれか一つに記載の発明において、打切りサンプルの仮イベント発生時間をイベントサンプル及び打切りサンプルの最大観測時間を考慮した基準観測時間に基づいて乱数を用いて仮設定することを特徴とする。
第７発明は、第６発明において、打切りサンプルの仮イベント発生時間を、乱数を用いて、実際の打切り時間を下限とし、基準観測時間を上限とする所定時間に設定することを特徴とする。
【００１７】
第８発明は、第１発明ないし第７発明のいずれか一つに記載の発明において、前記第１次生存時間解析モデルが部分最小自乗モデルであることを特徴とする。
第９発明は、第１発明ないし第８発明のいずれか一つに記載の発明において、生体状態の変化データと、少なくとも選抜された第１次説明変数群に基づいて、第２次生存時間解析モデルの第２次選抜基準を算出しつつ、所定選抜方法によって第２次生存時間解析モデルを最適化する第２次生存時間解析モデル最適化処理を含むことを特徴とする。
【００１８】
第１０発明は、第９発明において、選抜された第１次説明変数群のみを第２次生存時間解析モデル構築の際の説明変数として採用することを特徴とする。
第１１発明は、第９発明において、選抜された第１次説明変数群のみならず、他の説明変数を入れたものを第２次生存時間解析モデル構築の際の説明変数として採用することを特徴とする。
【００１９】
第１２発明は、第９発明ないし第１１発明のいずれか一つに記載の発明において、最適化された第１次生存時間解析モデルの潜在変数を用いて第２次生存時間解析モデルを最適化することを特徴とする。
第１３発明は、第９発明ないし第１２発明のいずれか一つに記載の発明において、前記第２次生存時間解析モデルが回帰モデルであることを特徴とする。
【００２０】
第１４発明は、第１３発明において、前記回帰モデルが比例ハザードモデル又は加速モデルであることを特徴とする。
第１５発明は、第１発明ないし第１４発明のいずれか一つに記載の発明において、選抜規準が残差平方和、相関係数、予測残差平方和、予測相関係数、赤池情報規準(AIC)、最小記述長(MDL)に基いて算出されることを特徴とする。
【００２１】
第１６発明は、予測対象サンプルの遺伝子発現のデータ、細胞内物質のデータの少なくとも一方のデータを入力して、第１発明ないし第８発明のいずれか一つに記載の発明によって決定された第１次生存時間解析モデルとその選抜された第１次説明変数群に基づいて生存時間解析に関する予測出力を得ることを特徴とする。
第１７発明は、予測対象サンプルの遺伝子発現のデータ、細胞内物質のデータの少なくとも一方のデータを入力して、第９発明ないし第１５発明のいずれか一つに記載の発明によって決定された第２次生存時間解析モデルとその選抜された第２次説明変数群に基づいて生存時間解析に関する予測出力を得ることを特徴とする。
【００２２】
第１８発明のデータ解析システムは、時間の経過とともに確率的に発生する生体状態の変化として所定のイベントを定義し、そのイベントが起こったイベントサンプルと打切りサンプルを含む生体状態の変化データを使用したコンピュータによるデータ解析システムであって、
打切りサンプルの実際には観測されなかったイベント発生時間を、統計的にランダムさを持って実際の打切り時間に対応する仮イベント発生時間に仮設定する打切りサンプル仮設定手段と、
生体状態の変化データを生存時間解析に使用する生命表へ適用する場合に、イベントサンプルのイベント発生時間と打切りサンプルの仮イベント発生時間の両方のイベント発生時間でのイベント非発生確率を算出し、その算出されたイベント非発生確率に基づいて個々のサンプルの相対的確率を算出し、その相対的確率を仮定した分布に基づいた変換又は仮定を前提としない変換をし、その変換結果を目的変数とし、
遺伝子発現のデータ、細胞内物質のデータの少なくとも一方のデータを説明変数とし、
前記目的変数と説明変数に基づいて、第１次生存時間解析モデルの第１次選抜基準を算出しつつ、所定選抜方法によって第１次生存時間解析モデルを最適化することにより重要な説明変数を選抜する第１次説明変数群選抜手段と、を含むことを特徴とする。
【００２３】
第１９発明は、第１８発明において、前記打切りサンプル仮設定手段が行う処理と第１次説明変数群選抜手段が行う処理とを少なくとも含んで一組として、その一組を複数回繰り返す繰返し処理手段とを備え、その結果得られた第１次説明変数群に基づいて重要な説明変数を選抜することを特徴とする。
第２０発明は、第１９発明において、繰返し処理によって得られた第１次説明変数群から重複する説明変数を選択することにより、重要な説明変数を選抜する重複選抜手段を含むことを特徴とする。
【００２４】
第２１発明は、第１９発明ないし第２０発明のいずれか一つに記載の発明において、繰返し処理の際に、遺伝子発現のデータ、細胞内物質のデータの一部のデータのみ使用する繰り返し制御手段を含むことを特徴とする。
第２２発明は、第１９発明ないし第２１発明のいずれか一つに記載の発明において、繰返し処理の際に、既に選抜された遺伝子発現のデータ、細胞内物質のデータを、その後の説明変数の候補として採用しないようにする繰り返し制御手段を含むことを特徴とする。
【００２５】
第２３発明は、第１８発明ないし第２２発明のいずれか一つに記載の発明において、前記打切りサンプル仮設定手段が、打切りサンプルの仮イベント発生時間をイベントサンプル及び打切りサンプルの最大観測時間を考慮した基準観測時間に基づいて乱数を用いて仮設定することを特徴とする。
第２４発明は、第２３発明において、前記打切りサンプル仮設定手段が、打切りサンプルの仮イベント発生時間を、乱数を用いて、実際の打切り時間を下限とし、基準観測時間を上限とする所定時間に設定することを特徴とする。
【００２６】
第２５発明は、第１８発明ないし第２４発明のいずれか一つに記載の発明において、前記第１次生存時間解析モデルが部分最小自乗モデルであることを特徴とする。
第２６発明は、第１８発明ないし第２５発明のいずれか一つに記載の発明において、生体状態の変化データと、少なくとも選抜された第１次説明変数群に基づいて、第２次生存時間解析モデルの第２次選抜基準を算出しつつ、所定選抜方法によって第２次生存時間解析モデルを最適化する第２次生存時間解析モデル最適化手段とを含むことを特徴とする。
【００２７】
第２７発明は、第２６発明において、選抜された第１次説明変数群のみを第２次生存時間解析モデル構築の際の説明変数として採用することを特徴とする。
第２８発明は、第２６発明において、選抜された第１次説明変数群のみならず、他の説明変数を入れたものを第２次生存時間解析モデル構築の際の説明変数として採用することを特徴とする。
【００２８】
第２９発明は、第２６発明ないし第２８発明のいずれか一つに記載の発明において、第２次生存時間解析モデル最適化手段が、最適化された第１次生存時間解析モデルの潜在変数を用いて第２次生存時間解析モデルを最適化することを特徴とする。
第３０発明は、第２６発明ないし第２９発明のいずれか一つに記載の発明において、前記第２次生存時間解析モデルが回帰モデルであることを特徴とする。
【００２９】
第３１発明は、第３０発明において、前記回帰モデルが比例ハザードモデル又は加速モデルであることを特徴とする。
第３２発明は、第１８発明ないし第３１発明のいずれか一つに記載の発明において、選抜規準が残差平方和、相関係数、予測残差平方和、予測相関係数、赤池情報規準(AIC)、最小記述長(MDL)に基いて算出されることを特徴とする。
【００３０】
第３３発明は、予測対象サンプルの遺伝子発現のデータ、細胞内物質のデータの少なくとも一方のデータを入力する入力手段と、第１８発明ないし第２５発明のいずれか一つに記載の発明によって決定された第１次生存時間解析モデルとその選抜された第１次説明変数群に基づいて生存時間解析に関する予測出力を得る予測手段と、を有することを特徴とする。
第３４発明は、予測対象サンプルの遺伝子発現のデータ、細胞内物質のデータの少なくとも一方のデータを入力する入力手段と、第２６発明ないし第３２発明のいずれか一つに記載の発明によって決定された第２次生存時間解析モデルと、その選抜された第２次説明変数群に基づいて生存時間解析に関する予測出力を得る予測手段と、を有することを特徴とする。
【００３１】
以下、上記第１発明〜第３４発明について説明する。
「生存時間解析」とは、前述したように時間の経過とともに確率的に発生するイベントを扱う解析方法を言う。イベントは一般的には解析の対象となる事象を言い、本明細書においては生体状態の変化が該当する。イベントには、死亡以外にも、疾病の発症、再発、疾病のステージの進行など多くの生命現象が採用できる。生存時間解析の対象となる時間については時計やカレンダーで把握される実時間を測定単位とする場合が多いが、ＱＯＬ(Quality  of  Life)を考慮した時間を用いる場合や、時間経過と一定の相関関係を有する測定単位、例えば、一日一定量の投薬を行った場合に時間軸の代わりに投薬量を取る変形などを採用することも可能である。
【００３２】
「生存時間解析に使用する生命表」としては、例えば、カプラン・マイヤー法又はカトラー・エドラー法による生命表を適用することができる。
打切りサンプルがある生体状態の変化データについては、打切りデータの不適切な扱いによるバイアスを混入することがないように工夫することが大事である。このバイアスをなくす方法として、第１発明に係る一形態では、生存時間データの打切りサンプルに乱数などの確率的手法を用いてランダムに仮イベント発生時間を割り当てる。そして生体状態の変化に関するデータに生存時間解析に使用するカプラン・マイヤー法又はカトラー・エドラー法による生命表を適用して変化（イベント）の発生しない確率を算出して、そのイベント非発生確率に基づいて個々のサンプルの相対的な確率（例えば、生存率）を算出する。
【００３３】
この個々のサンプルの相対的な確率を算出する一方法として、生命表としてカプラン・マイヤー法を採用した場合は、カプラン・マイヤー法での生存率は、集団としての生存率を示すが、ここでは個人ｊの事象（イベントまたは打切り）発生時点での集団としての生存率をその個人ｊの相対的なイベント発生傾向を意味すると捉えなおすことによって個人の相対的な生存率を算出することができる。
なお、上記算出処理は、各種の生命表によって把握できる集団としての生存率を個々のサンプルの相対的生存率に変換する方法とも考えられ、Ｐ'i＝1−ＰKM(Ti)で示される後述の実施形態に記載した方法がその一例である。
【００３４】
なお、第１発明における「仮定した分布に基づいた変換又は仮定を前提としない変換」は、生体状態の変化データとしての前記相対的確率を説明変数の多項式で解析できるようにするために行われるものである。分布を仮定した場合は、確率を対数変換後に負の数にしたものを生体状態の変化を観測した時間で割るという変換、確率を対数変換後に負の数にしたものをさらに対数にしたものを生体状態の変化を観測した時間で割るという変換などがある。さらに、分布を仮定しない場合にはロジット変換などが採用できる。
【００３５】
第１発明において、第１次生存時間解析モデルに使用するモデルは限定されるものではないが、遺伝子発現量などの膨大な数を説明変数として扱え、オーバーフィットの問題を解決し、次元圧縮とモデルフィットを同時に行える相関モデルであることが望ましく、通常は部分最小自乗モデル(PLSモデル)が採用され、他には重回帰分析、主成分分析、ロジスティック回帰モデル、SIMCAモデル、正準判別分析モデルなどが採用できる可能性がある。
【００３６】
また、第１発明において「細胞内物質のデータ」としては、例えば細胞内物質の量が挙げられ、例えば、代謝中間体を含めた代謝産物全部であるメタボロームを意味する。
「遺伝子発現」の用語は、ｍＲＮＡ発現（トランスクリプトーム）やｍＲＮＡによる翻訳の結果として生じるタンパク質（プロテオーム）を含むものとして使用する。
【００３７】
また、遺伝子発現の量、細胞内物質の量は測定値に限定されるものではなく、加工計算された値、相対的な値、間接的に物質量を示す量でもよい。例えば、質量スペクトルでタンパク質の発現量を測定できることを応用して、生体状態の変化と質量スペクトルとを直接関係づける生存時間解析モデルを構築することもできる。ＤＮＡチップで各スポットの測定量を説明変数として採用することもできるし、そのスポットの測定量の加工値も利用することもできる。
【００３８】
また、第１次選抜基準としては交差検証成績だけでなく、残差平方和、相関係数、予測残差平方和、予測相関係数、赤池情報規準(AIC)、最小記述長(MDL) その他の統計量も使用することができる。所定選抜方法としては、遺伝的アルゴリズム(ＧＡ)、ステップワイズ（ＳＷ）、シミュレーテッドアニーリング(ＳＡ)、マルコフチェインモンテカルロ法(ＭＣＭＣ)などが例示できる。交差検証成績(cross  validation)の方法としては３分割法（3-fold），５分割法（5-fold），リーブ・ワン・アウト法（leave-one-out）などを用いることができる。
【００３９】
なお、本明細書において、「最適化」とは第１選抜基準、第２選抜基準がそのときの条件の範囲で改善がみられなくなるまで改良したことを意味しており、全ての説明変数の組み合わせの中で最適なものを見いだしたことを確認したことを意味するものではない。
【００４０】
第４発明において、「遺伝子発現のデータ、細胞内物質のデータの一部のみ使用する」場合としては、各一組において、遺伝子発現のデータ、細胞内物質のデータを採用するデータが全て異なっている場合（各一組が重複しない場合）と、遺伝子発現のデータ、細胞内物質のデータを採用するデータの一部が異なっている場合（各一組が重複する部分がある場合）とがある。
第７発明はともに実際の打切り時間以後のある時間に仮イベント発生時間を仮設定する場合（いわゆる右側打切り）を想定している。
【００４１】
本発明において、モデル構築処理と、診断予測処理とは必ずしも時間的、空間的に連続した処理である必要はない。例えば、胃癌、肺癌、肝臓癌、……と診断された多数の患者の生存時間解析データと各遺伝子発現データを使用して第１発明〜第１５発明までの処理により、各癌診断の相関モデルの構築をデータ解析センターにおいて集中的に行い、第１６発明、第１７発明に記載したように、そのデータ解析の成果として構築された癌診断用の相関モデル（選抜された説明変数群が重要であるという情報も含む）を世界各地の病院において診断に利用する形態などが考えられる。
なお、第１８発明〜第３４発明はデータ解析装置又はデータ解析プログラムとしても認識できるものである。
【００４２】
生存時間解析(比例ハザードモデルや加速モデル)に限らず、従来多くの統計的手法は膨大な数の説明変数を扱うことができない。例えば、重回帰分析ではサンプル数(Ｎｓ)が説明変数の数(Ｎｐ)より少ない時には計算不可能であり、通常はＮｓ／Ｎｐ＞３〜５が望ましいといわれている。そこで、本発明に係るデータ解析方法及びシステムでは、最初から膨大な数の説明変数全体を扱うのではなく、第１次生存時間解析モデルによって選抜した説明変数セットを用いて最終的な第２次生存時間解析モデルの構築を行なう。
【００４３】
第１次生存時間解析モデルを最適化するためには、第１次生存時間解析モデルの出来映えを第１次選抜基準で評価し、所定探索方法でモデルに使用する説明変数セットを最適化することができる。そして、その後の第２次生存時間解析モデルの探索において用いる説明変数をこれら第１次選抜説明変数群に限定する。あるいは第５発明などに記載したように、探索方法としてＧＡを採用し、その初期染色体に第１次選抜説明変数群を含めるのである。このような方法であれば、第２次生存時間解析モデルの探索範囲が限定されることになり、最終モデルの最適化が非常に効率的となる。
【発明の効果】
【００４４】
以下、各発明の効果などについて説明する。
なお、データ解析システムに係る発明の効果もデータ解析方法の効果に対応しているのでそれらの効果の記載は省略する。
第１発明であれば、打切りサンプル仮設定処理によって打切りサンプルの実際には観測されなかったイベント発生時間を統計的にランダムさを持って仮設定するので、生存時間解析においてバイアスが混入することを抑制しつつ、打切りサンプルの仮イベント発生時間を設定することができ、打切りサンプルデータを有効に生かしつつ、高精度なデータ解析を行える可能性を高めることができる。
【００４５】
第２発明であれば、前記繰返し処理を行うことによって、選抜される説明変数の全体としての容量を大きくでき、１回だけの打切りサンプル仮設定処理に比べて、仮設定された打切りサンプルの仮イベント発生時間のバイアスをより打ち消すように説明変数を選抜される可能性を高めることができる。
第３発明であれば、ランダムに設定された一組を複数回繰り返すことによって得られた複数の第１次説明変数群内において、重複して現れる説明変数は、より重要な説明変数である確率が高いので、そのような説明変数を選抜することにより重要な説明変数を漏らすことなく、第１次説明変数の選抜処理を行うことができる。
【００４６】
第４発明であれば、前記繰返し処理の際に、遺伝子発現のデータ、細胞内物質のデータの一部のみ使用するので、選抜される説明変数に多様性を確保することができる。
第５発明であれば、前記繰返し処理の際に、前に採用された遺伝子発現のデータ、細胞内物質のデータとは違うものを説明変数として採用することになるので、選抜される説明変数に多様性を確保することができる。
【００４７】
第６発明であれば、打切りサンプルの仮イベント発生時間をイベントサンプル及び打切りサンプルの最大観測時間を考慮した基準観測時間に基づいて設定するので、打切りサンプルの仮イベント発生時間をサンプルが有する時間分布に対応させて妥当性を高くして仮設定することができる。
【００４８】
第７発明であれば、打切りサンプルの仮イベント発生時間を、乱数を用いて、実際の打切り時間を下限とし、基準観測時間を上限とする所定時間に設定するので、バイアスのないランダムな仮イベント発生時間を設定できる可能性が高まる。
第８発明であれば、前記第１次生存時間解析モデルを部分最小自乗モデルとすることにより、説明変数の次元を圧縮しつつ、モデルフィットを同時に行うことができ、遺伝子発現量のような膨大な説明変数の処理に好適となる。
【００４９】
第９発明であれば、少なくとも選抜された第１次説明変数群に基づいて第２次生存時間解析モデル最適化処理を行うので、重要な説明変数を漏らすことなく、説明変数の選択の範囲を狭めることができ、最適化処理時間を短縮して遺伝子発現量などの膨大な説明変数の解析を効率的に行うことができる。
第１６発明及び第１７発明であれば、患者等の遺伝子発現量などを得ることで生存時間を予測することができるので、癌などの重大疾病の的確な医療診断、予後治療、及びＱＯＬの向上に役立てることができる。
【発明を実施するための最良の形態】
【００５０】
以下、本発明の実施の形態を図面に基づき説明する。
（第１実施形態）
図１は本発明の第１実施形態に係る基本的なフローチャートである。
図１に示すように、この第１実施形態に係るデータ解析方法は、患者の生体状態の変化データと遺伝子発現量などを示したサンプルのデータを、患者のカルテ（電子カルテを含む）、各種データファイルから、手入力、転送、コピーなどの方法によってコンピュータの記憶装置に入力（ステップ１）し、第１次生存時間解析モデルの第１次選抜基準を算出しつつ、所定選抜方法によって第１次生存時間解析モデルを最適化することによって重要な説明変数である第１次説明変数群を選抜し（ステップ２）、生体状態の変化データと、少なくとも絞り込まれた第１次説明変数群に基づいて、第２次生存時間解析モデルの第２次選抜基準を算出しつつ、所定選抜方法によって第２次生存時間解析モデルを最適化する（ステップ３）ことによって行われる。
【００５１】
ステップ１の入力データ（サンプル集合）の取得の方法は各種の形態が考えられるが、生命表の元データと説明変数となるｍ個（たとえば２０００個）の遺伝子発現量データの形で与えられることが多い。なお、この実施形態では図１に示すステップ１及びステップ２の処理が第１発明に対応する部分となり、ステップ３の処理が第９発明に対応する部分となる。
【００５２】
図２は第１次説明変数選抜処理（図１のステップ２）に係るフローチャートの一例である。
図２において、打切りサンプルの仮イベント発生時間を仮設定処理し（ステップ１１）、図１のステップ１において得られた生体状態の変化データから第１次生存時間解析モデルの目的関数を算出する（ステップ１２）。算出された目的変数と遺伝子発現量で表される説明変数に基づいて第１次生存時間解析モデル（ＰＬＳモデル）の第１次選抜基準を算出しつつ、所定選抜方法によってＰＬＳモデルを最適化し（ステップ１３）、第１次生存時間解析モデルによって選抜された説明変数群を第１次説明変数群として保存する（ステップ１４）。
【００５３】
次に前記ステップ１１〜ステップ１４の処理を設定されたｎ回繰り返した否かを判別し（ステップ１５）、繰り返していないと判別された場合は、所定回数になるまで繰り返し（ステップ１５のＮＯ）、ｎ回繰り返した後は（ステップ１５のＹＥＳ）、ｎ回の選抜された説明変数群に基づいて所定の選定方法に基づいて最も重要な説明変数を決定する（ステップ１６）。
なお、ステップ１１からステップ１５までの繰り返し処理により得られた複数回の選抜説明変数群に基づいて第１次説明変数を記憶装置の初期プールに記憶するまでの処理が第２発明に対応する。
【００５４】
図３は打切りサンプルのイベント発生時間を仮設定する処理（図２のステップ１１）の一例を示したフローチャートである。
図３において、入力した全てのサンプルから最大観測時間を設定し（ステップ２１）、この集団についての基準観測期間を設定し（ステップ２２）、各打切りサンプルの実際の打切り時間を抽出し（ステップ２３）、乱数を用いて各打切りサンプルの仮イベント発生時間を実際の打切り時間を下限とし、基準観測時間を上限となる所定の時間に設定する（ステップ２４）。
【００５５】
さらに詳しく説明する。
まず、患者iの症例について、イベントサンプルのイベント発生時間又は打切りサンプルの実際の打切り時間をtiとし、イベントか打切りかを示す変数をdiとする。
最大観測時間の設定については、例えば、tiの最大値を求めて最大観測時間(Tmax)に代入する。基準観測時間の設定については、例えば、Tmaxをa倍（ａは正の実数）したものを基準観測時間(Tk)とする。
次に各患者iにつき、di=censor(つまり打ち切り症例)なら乱数（random）を用いて、
【００５６】
【数１】
2005222422.tif 000003
【００５７】
を計算することにより、仮イベント発生時間をTiに代入する。
一方、di=event(つまりイベントの発生した症例)なら、
Ti←ti
として、イベント発生時間tiをそのままTiに代入する。
これを全ての患者iについて行うことにより、イベントの発生した患者のみならず、打切りの患者についても、仮イベント発生時間が設定されたことになる。
【００５８】
図３におけるステップ２１は、患者集団の生命表における観測時間が５ヶ月〜２４ヶ月の範囲ならば、例えば、tiの最大値を取ることによって最大観測時間Tmax＝２４ヶ月として設定する。また、ステップ２２の基準観測時間TkをTk＝ａ×Tmaxというように設定し、ａは１．０〜２．０程度の定数とし、例えば、ａ＝１．５とすれば、Tk＝１．５・Tmaxになる。ステップ２４では、患者iが打切り症例である場合には、上記(1)式に基づいて、ti＝１１ヶ月、Tk＝３６ヶ月とすると、Ti＝１１＋(３６−１１)×randomとなり、０＜random＜１として、仮イベント発生時間を実際の打切り時間の１１ヶ月以上で、基準観測時間Tk＝３６ヶ月の間の全くランダムな時間に設定することができる。
【００５９】
図４は図２においてステップ１２で示した目的変数の設定の一例を示す図である。
図４において、ステップ１１で得られたTiデータに基づいて、カプラン・マイヤー法によって生存率ＰKM(Ti)（イベント非発生率ＰKM(Ti)）を算出し（ステップ３１）、その生存率ＰKM(Ti)を用いて各患者ｉの相対的な生存率Ｐ'iを求め（ステップ３２）、その相対的生存率Ｐ'iのロジット変換を行って目的変数とする（ステップ３３）。
さらに詳しく説明する。
カプラン・マイヤー法によれば、各死亡（イベント）発生時間Tjの生存率ｐKM(Tj)は、di=eventならば、
【００６０】
【数２】
2005222422.tif 000004
【００６１】
によって算出される。
一方、di=censorならば、
ｐKM(Tj)＝1
となる。
ここでＮalive はTjの直後の生存観測症例数であり、Ｎalive ＋Ｎdead at the timeはTj直前までの生存観測症例数である。なお、Ｎdead at the timeは、イベント発生時間Tjにイベントが生じた症例数である。
また、カプラン・マイヤー法による生存率ＰKM(Ti)は、
【００６２】
【数３】
2005222422.tif 000005
【００６３】
となる。つまり、ｐKM(Tj)をTjの小さい順にTi まで乗算したものが生存率ＰKM(Ti)となるのである。
生存率ＰKM(Ti)は患者集団全体の生存曲線ＰKM(t)となるわけであるが、集団の生存率が高い時間にイベントが発生した症例は、集団の中では相対的に生存率の低い症例と見なすことができる。逆に集団の生存率が低い時間にイベントが発生した症例は、集団の中では相対的に生存率の高い症例と見なすことができる。
【００６４】
このことを用いてステップ３２に示す各患者iの相対的な生存率Ｐ'iを
Ｐ'i＝1−ＰKM(Ti)
によって算出する。
最後にその相対的な生存率Ｐ'iを Logit変換、即ち、
【００６５】
【数４】
2005222422.tif 000006
【００６６】
を行なって目的変数とする。
なお、生存率ＰKM(Ti)を「累積生存率」、生存率ｐKM(Tj)を「その時点の生存率」と定義されることもあるが、通例的に、単に生存率と呼ばれることも多い。ここでは単に生存率とした。
図５（Ａ）（Ｂ）はそれぞれ、多様性を持った説明変数選抜の一例を示す図である。
図５（Ａ）は図２に示したフローチャートのステップ１１とステップ１２の間にステップ１１ａを設け、そのステップ１１ａにおいて乱数によって幾つかの症例や説明変数を解析対象から除外するフローチャートを示したものである。
そのような例として生理化学的な事実や推測に基づいて互いに関連、近似すると思われる説明変数をコンピュータの判別処理で省く形態や、いくつかの説明変数間の相関係数を計算して、その中で特に相関係数の高い説明変数を除外する形態が考えられる。
【００６７】
図５（Ｂ）は、図２に示したフローチャートのステップ１４とステップ１５の間にステップ１４ａを設け、そのステップ１４ａにおいて既に選抜された説明変数を以後の選抜対象の候補から除外するフローチャートを示したものである。
【００６８】
図６は図２においてステップ１３で示した所定選抜基準に基づいて説明変数を選択して部分最小自乗法モデル（ＰＬＳモデル）を最適化する処理の一例を示すフローチャートである。
ここでは簡単に説明するため、第１次選抜基準は交差検証成績を独立変数として持つ関数として採用している。
【００６９】
図６において、まず、相関モデル作成用のデータが入力された状態から、まず、初期設定として、選択された説明変数の数Ｎｐを０とし、交差検証成績ＣＶの最良値ＣＶoを−∞とする（Ｓ１０１）。次に、説明変数の選択を行うために、まず、説明変数を指す番号ｉを１とし(Ｓ１０２)、第ｉ変数（遺伝子発現の量）を仮に採用して(Ｓ１０３)、部分最小自乗法を実行し、交差検証成績ＣＶを計算する(Ｓ１０４、図６参照）。
ここで、交差検証成績ＣＶの計算として、リーブ・ワン・アウト処理を用いる。これは、たとえば５０個のサンプルからなる訓練集合において、１番から５０番の全てを順次１個ずつ除いて残りの４９個のサンプルで構築したモデルで、その時除いた１個を予測した結果と実際の結果とを比較し、その誤差が大きい場合に、仮に選択した説明変数（第ｉ変数）が適していないと判断する手法である。
【００７０】
もし、得られた交差検証成績ＣＶが現在の最良値ＣＶoより最適化されれば（Ｓ１０５でＹＥＳ）、第ｉ変数を採用し、かつ、成績ＣＶを新しい最良値ＣＶoに更新する（Ｓ１０６）。しかし、得られた成績ＣＶが最良値ＣＶoより大きくなければ（Ｓ１０５でＮＯ）、第ｉ変数を採用しない(Ｓ１０７)。そして、ステップＳ１０２に戻り、同様の処理を繰り返す。この処理を交差検証成績ＣＶが改善されなくなる(Ｓ１０８でＮＯ）まで繰り返す。
【００７１】
ここで、相関モデルに採用する説明変数については１つずつ段階的に増加（追加）または減少（除外）して交差検証成績ＣＶを評価判定している。すなわち、全体としての合致度合いがよくなるように各説明変数を解析に加えるかどうかを逐次判定しながら、説明変数の取捨選択を行い、これを、全体としての合致度合いがよくならなくなるまで繰り返す。以上の処理で改善があると、ふたたびステップ１０２の初め(ｉ＝１)に戻り、それまでに選択されている説明変数を基に、さらに説明変数の選択を繰り返す。そして、最適化されたＰＬＳモデル（選抜説明変数も含む）を決定する（Ｓ１０９）。
【００７２】
なお、決定されたＰＬＳモデルの予測力を判断するために、訓練集合とテスト集合とに予め分割しておいたテスト集合を用いてデータ解析することもできる。このようなＰＬＳモデルの評価は必ずしも必要でないが、予測力を判断するには有効である。
【００７３】
図７はリーブ・ワン・アウト処理を含む交差検証成績ＣＶの計算（図６のＳ１０４の処理）のフローチャートの一例である。
具体的には交差検証成績を予測誤差自乗和(PRESS)に基づいて算出している。
選択された変数について交差検証成績が計算される。図７において、まず、PRESSの初期値を０とし（Ｓ１２１）、ｎ個の集合内のサンプルを指す番号ｊを１とし(Ｓ１２２)、第ｊサンプル以外のｎ−１個のサンプルで部分最小自乗法を実行し(Ｓ１２３)、第ｊサンプルの目的変数を予測する（Ｓ１２４）。次に予測と実測との差の自乗を計算してPRESSに加算する（Ｓ１２５、Ｓ１２６）。次に番号ｊを１増加し(Ｓ１２２)、同様の処理をおこなう。これを番号ｊ＝ｎまで各サンプルについて繰り返す。
【００７４】
得られたPRESSは、１個のサンプルを順次除外して計算した予測値と実測値との差の平方和であり、予測誤差を表わす量である。この予測残差自乗和PRESSの符号を変えたものを交差検証成績ＣＶとする(Ｓ１２７）。このフローチャートでは、交差検証法を用いて、入力変数（説明変数）を段階的に１つずつ追加・除外しながら、交差検証成績（ＣＶ=−PRESS）を最適化している。
【００７５】
図６、図７に示すような交差検証法を用いて、選抜基準の一例としての交差検証成績を最適化することによりＰＬＳモデルが構築され、図２に示すステップ１４において決定されたＰＬＳモデルによる説明変数を第１次説明変数群としてコンピュータの記憶装置に記憶する。同様にコンピュータの記憶装置にＰＬＳモデルのスコア行列、ローディング行列の値も記憶しておく。
【００７６】
図２に示すステップ１５において例えば繰り返し回数を複数の所定回数に設定して、ステップ１６において例えば重複する説明変数を最終的に選抜する説明変数として決定することにより、第１次説明変数群による初期プールの設定が完了する。
この初期プールの設定状態では図１のステップ２の処理が終わった状態であり、図１のステップ３に示すように第１次説明変数群に基づいて所定選抜方法によって最終モデルを構築する。ステップ２が終了した段階では、説明変数の個数は通常の統計的手法又は多変量解析が適用できる程度に説明変数の個数は少なく選抜されている場合が多いので、所定の生体状態の変化を目的変数として、選抜された説明変数、その潜在変数を用いて相関モデルを構築することにより、生体状態の変化に大きく寄与する遺伝子を特定することが可能になる。
【００７７】
この第２次生存時間解析モデルとしては、比例ハザードモデルや加速モデルが採用できる。比例ハザード法とはコックス回帰とも呼ばれ、寿命データ解析のためのノンパラメトリックな手法の一つであり、性別、年齢、遺伝子発現などの説明変数(Cox回帰では共変量ともいう、また生存時間解析という性格上、予後因子という表現をされることもある)を含めてモデルを構築できる利点がある。
また、第２次生存時間解析モデルの適応の度合いを数値化する方法としては、赤池情報規準(AIC)、最小記述長(MDL)などの選抜基準の採用が考えられる。又は交差検証成績によってモデルの適応の度合いを数値化することもできる。数量化としては第１生存時間解析モデルの最適化処理で説明したように交差検証成績の予測誤差自乗和(PRESS)などの代表値を用いることができる。
【００７８】
また、交差検証成績の方法としては、第１生存時間解析モデルと同様にサンプルを１個づつ逐次除外する（リーブ・ワン・アウト法）や、複数のサンプルを除外（リーブ・ｎ・アウト法）して交差検証成績を評価してもよい。さらに、３分割法（3-fold）、５分割法（5-fold）なども用いることができる。いずれにせよ、選抜基準としての適応度さえ求めることができれば、遺伝的アルゴリズム(ＧＡ)、ステップワイズ（ＳＷ）、シミュレーテッドアニーリング(ＳＡ)などの探索方法によってモデルの最適化が可能になる。
【００７９】
図１３はデータ解析システムの一例を示す概略的なブロック図である。
このデータ解析システムは、例えば、インターネット１１等のネットワーク１０に接続可能なコンピュータ１（例えば相関モデル構築サーバ）に設定される。
そのコンピュータ１は、前記打切りサンプル仮設定手段としての打切りサンプル仮設定部２と、生体状態の変化データと複数の遺伝子発現のデータ、細胞内物質のデータの少なくとも一方のデータを記憶するデータ記憶部３と、前記第１次生存時間解析モデルを最適化する手段としてのＰＬＳモデル最適化部４と、最適化によって絞り込まれた重要な説明変数群を記憶する第１次説明変数群記憶部５とを少なくとも有している。
なお、第１選抜基準算出部、説明変数を所定方法で選択する選択部は、ＰＬＳモデル最適化部４内に設けられている。
【００８０】
また、最終モデルを構築する場合は、例えば、第２６発明などに記載されたように、生体状態の変化と、少なくとも選抜された第１次説明変数群に基づいて、比例ハザードモデルの第２次選抜基準（ＡＩＣなど）を算出しつつ、所定選抜方法によって比例ハザードモデルを最適化する前記第２次生存時間解析モデル最適化手段としての比例ハザードモデル最適化部６を有している。
【００８１】
コンピュータ１に入力される上記各データは、大学、研究所が提供する生存時間解析に関する生体関係情報を提供するウエブサイト、疾病関係データベースなどの各種データベースから得ることも可能である。また、ＤＮＡチップなどを用いて患者の遺伝子発現量などのデータを取得するともに、カルテ・データべース１５から生体状態の変化データを病院コンピュータ１４に入力し、それらのデータをネットワーク１０を介してコンピュータ１へ入力することも可能である。
図示しない入力部、選択部、出力部などはコンピュータ１側に設ける場合が一般的であるが、記憶部、入力部、選択部、出力部などは、コンピュータ１側に設ける形態の他、ネットワーク上の記憶装置、機能実現プログラムで構成しても良い。
【００８２】
本データ解析システムは図１３に示したようなネットワーク型のコンピュータ、クライアントサーバー型のコンピュータのみならず、スタンドアロン型のコンピュータでも構成できる。
スタンドアロン型のコンピュータの場合は、主にソフトウエア処理で各処理部を構成し、記憶部を内部記憶装置や外部記憶装置で構成する。生体関係情報提供部は、光学的、磁気的記憶装置などの大容量記録媒体で構成されることが多い。
この場合も入力部、出力部はコンピュータに付随する入力装置、出力装置を利用するとともに、必要に応じて、入力、選択、出力を効率化する各種のソフトウエアで構成することもできる。選択部としては、内部記憶装置などに記憶されたサンプルの集合としてデータを選択する方法としては、データが記憶された少なくとも一つのファイルを指定する方法や、データが記憶されたファイル内の特定データを個別に指定する方法などが考えられる。
【００８３】
また、ネットワークに病院コンピュータ１４を接続し、患者の遺伝子発現量などのデータをＤＮＡチップなどを用いて取得し、病院コンピュータ１４に入力することで、コンピュータ１４内にある第３３発明、第３４発明に記載した予測手段を用いて生存時間等を予測することもできる。変形としては、予測手段は前記コンピュータ１内に設けることもでき、ネットワークを介して患者の遺伝子発現量のデータを送信し、予測出力を病院コンピュータ１４で受信するようにしてもよい。このような形態であると相関モデル構築サーバで所定期間毎に更新された最も新しく、予測精度の高い生存時間解析結果を世界各地の病院で利用することができる。
【実施例１】
【００８４】
以下、具体的な実施例について説明する。
上記本発明の方法によって相関モデルの探索範囲が効果的に限定されるとともに、打切りサンプルを適切に扱うモデル構築となっている実施例を示す。
データとしてはDLBCL(びまん性大細胞型Ｂリンパ種)の患者２４０人の予後データと7399遺伝子発現データ｛A Rosenwaldet al.. N. Engl.J.Med.346(25),1937-1947（2002）
(http://llmpp.nih.gov/DLBCL/)｝を解析対象とした。スポットパターンでχ１又はχ２が０となるものを除いた７３９９スポットについてｌｏｇ（χ１／χ２）を計算して説明変数とした。
【００８５】
（１：第１次説明変数群の選抜）  
実際には観測されなかった打切りサンプル（censorデータ）の仮イベント発生時間を、乱数を用いて打切り時間と最大観測時間Ｔｍａｘの１．５倍との間のある時間に仮設定した。
ここで、図４において説明したようにカプラン・マイヤー法による生命表に基づいて生存率ＰKM(Ti)を求め、その生存率ＰKM(Ti)を各サンプルの相対的確率Ｐ'iに変換し、その相対的確率Ｐ'iのlogit変換値、即ち、log{ P'i/ (1−P'i) }を目的変数としてステップワイズな説明変数増加によって相関係数を最適化しながらＰＬＳモデルを構築した。５０変数に至るまでの計算処理には約１５分要した。これを図２のステップ１５でｎ＝８とて８回繰り返し、合計４００個の説明変数が計算された。そして、図２のステップ１６に一例として示されるようにこの４００個に基づいて第１次選抜説明変数群を選んだ。具体的には、４００個の説明変数の中には重複して数えられた説明変数があるので、それら重複した個数を除くと、３３４個の説明変数を選抜できた。
【００８６】
（２：説明変数候補を限定したステップワイズな説明変数選抜）
上記３３４個の第１次選抜説明変数群にステップワイズによって更に選んだ説明変数を加えながら、比例ハザードモデルのＡＩＣ値を最適化(最小化)した。選抜説明変数の数(Ｎｐ)が２０に至るまでに約１０時間の処理を必要とした。
図８は実施例１の実験結果を示す図であり、図８に示すように横軸に選抜説明変数の数Ｎｐを取り、縦軸に適応度（ＡＩＣ／−２）を取ってプロットすると実線のようになった。
比較のため、７３９９個の説明変数から上記第１次説明変数選抜を経由することなく、比例ハザードモデルを構築し、プロットすると破線のようになった。この比較例の処理にはＮｐ＝２０に至るまで約１週間が必要であった。
【００８７】
図８にて実線と破線を比較すると、ほぼＮｐ＜１０の範囲では殆ど同じであり、またＮｐ＞１０の範囲では、同程度のＡＩＣ値で推移していることから、本発明の第１次選抜方法は重要な説明変数をもらすことなく説明変数候補を７３９９個から３３４個に絞り込むことに成功していることがわかる。また、約１０時間と約１週間とを比較すると、コンピュータの実行時間も大幅に短縮していることがわかる。
特に、Ｎｐ＞２０の範囲は７３９９個の説明変数を用いた場合には実行時間が非常に長くなり、実際上探索が不可能になるが、本実施例の方法によって３３４個に絞り込むことによって探索が可能になった。
【実施例２】
【００８８】
（説明変数候補を限定したＧＡによる説明変数選抜）
説明変数候補を第１次選抜説明変数に限定した場合(３３４個)と、限定しない場合(７３９９個)でＧＡの効率を比較した。 ＧＡの適応度はＡＩＣ／−２とし、プールサイズを１００とし、UX法、エリート戦略、トーナメント方式とした。
まず、遺伝子アルゴリズムについての予備的説明を行う。
遺伝的アルゴリズム(ＧＡ; Genetic Algorithm)は、進化論をヒントとした非線型最適化手法であり、ここ数年間、頻繁に用いられて定着してきた。遺伝子アルゴリズムを説明した文献としては、伊庭斉志;「遺伝的アルゴリズムの基礎」;オーム社;(1994)などがある。
ＧＡでは最適化対象がビット列などで表現される。このような表現の形式を遺伝型(genotype)と呼び、ひとつひとつの遺伝型の具体的な文字列を染色体(chromosome)と呼ぶ。最適化の指標である適応度(fitness)は遺伝型より一義的に決定される。
【００８９】
（処理１） 初期染色体プールの準備
ランダムに染色体の集合を準備する。染色体の数はプールサイズと呼ばれる。
（処理２）  交叉、突然変異、選択による最適化
[2-1] 交叉(crossover)
集合より2つの染色体を選抜し、ビット列のランダムな交叉によって新しい染色体を創生する。交叉を１個所に限る一点交叉(onepoint crossover)と全ビット間で交叉を行なう一様交叉(unform crossover,UX)などがあり、一般的にはUX法が優れているとされている。
[2-2] 突然変異
集合あるいは新しい染色体の任意のビットをランダムに選んで反転させる。
[2-3] 選択(selection)
各染色体の適応度を計算し、染色体プールを次世代のものに置き換える。適応度に比例した確率で染色体を選抜するルーレット方式や親世代と子世代との間で適応度の対決を行なうトーナメント方式などがある。これらいずれにおいても乱数が用いられることが多く、一度得られた最適染色体が失なわれることもある。このようなことが無いよう、プール内の最適染色体は必ず次世代に選抜されるとするエリート戦略なども用いられる。
【００９０】
（処理３）  終了
世代数がある回数に達した場合や、最適適応度に変化が見られなくなった場合に処理2は終了し、最適解は最適染色体の遺伝型およびその適応度として出力される。ＧＡは基本的には乱数を用いた方法であるため、処理１〜処理３が何回か試みられ(Run)、その中の最適染色体が用いられることが多い。
【００９１】
本実施例に用いたGAの具体的な設定を示すと、次のとおりである。
(1)genotype={b1,b2,b3,,,,,,bNg}
遺伝子発現量gjを説明変数に用いる場合にはbj=ON状態、
用いない場合にはbj=OFF状態とする。
Ngは説明変数候補の数
(2)fitness=-AIC
(3)プールサイズ＝100
(4)初期染色体は、平均でmin_of(Ni, Ng, 300)/2個のビットがON状態となるように乱数を用いて準備する。ここで300は実装の都合上設定された定数である。
(5)交叉は、2つの染色体を乱数にて抽出し、各ビットをそれぞれ1/2の確率でいずれかの親染色体から引き継いだ子染色体を作成する。
(6)突然変異は、ON状態ビット数が増加する変異と減少する変異とが同じ確率で発生するように乱数でビット反転を行なった。
(i)ON状態ビットの反転確率＝1.1/ON状態ビット数
(ii)OFF状態ビットの反転確率＝1.1/OFF状態ビット数
(7)選抜は、子染色体を作成する毎にトーナメント方式でプールを置換する。
(i)子染色体が新たな最適解となる場合は無条件で置換
(ii)トーナメント相手が最適解の場合には無置換
(iii)子染色体のfitness＜トーナメント相手のfitnessの場合には0.25の確率で置換
(iv)子染色体のfitness＞トーナメント相手のfitnessの場合には0.75の確率で置換
(8)終了は、交叉→突然変異→選抜の繰り返し回数の上限を最初10000とし、最適解が見出される毎に1000追加されるものとした。
【００９２】
この第２実施例の結果を図９に示す。図９においても実線が第１次選抜説明変数を３３４個に限定した場合であり、破線が７３９９個から始め、限定しない場合を示す。破線で示すように限定しない場合には探索範囲が広すぎてＧＡが有効に働いていないが、限定することでＧＡが働いていることがわかる。
【実施例３】
【００９３】
（初期染色体を設定したＧＡによる説明変数選抜）
第１次選抜説明変数群を初期染色体の一部に設定した場合と、設定しない場合でＧＡの効率を比較した。８回行なったＰＬＳモデル構築にて選抜された説明変数セットをそれぞれ染色体とし、合計８個の良好な初期染色体を準備した。プールサイズ１００の残る９２個の初期染色体は、実施例２と同様にランダムに設定した。その他のＧＡ条件も実施例２と同じである。図１０はその実験結果を示す図である。図１０では実線で示す初期染色体の一部に設定した場合は、破線で示す設定しない場合に比べて、良好な初期染色体の投入によって既に良好な比例ハザードモデルが得られており、ＧＡを続けるとさらに改善されることがわかる。
【００９４】
以上、実施例１〜実施例３によって、本実施形態に係る方法によって、第１次説明変数を選抜が、説明変数候補を限定したステップワイズな説明変数選抜、説明変数候補を限定したＧＡによる説明変数選抜、初期染色体を設定したＧＡによる説明変数選抜のいずれにおいても有効であることが確認できた。
【実施例４】
【００９５】
図１１は実施例２の方法に基づくＧＡ最適化の途中で得られたＡＩＣ／２＝５１８（ｎｐ＝７４）のモデルからステップワイズな変数増減法によって最適化を進めて得られた結果のモデル（比例ハザードモデル）であり、最適解(ＡＩＣ／２＝４４８、小数点以下四捨五入)のプロット図である。
図１１では横軸に生存時間を取り、縦軸にハザード値を取ってある。図１２は同じDLBCL(びまん性大細胞型Ｂリンパ種)の患者２４０人の予後データと７３９９遺伝子発現データを用いて前記特許文献２に記載された方法（特許文献２の図９のデータに対応する、ＡＩＣ／２＝６１９、小数点以下四捨五入）によって得られたプロット図である。なお、実施例４の方法では、ＡＩＣ／２＝４４８であることから、特許文献２の方法におけるＡＩＣ／２＝６１９に比べて数値的にも優れていることが分かる。
【００９６】
特許文献２の生存時間解析モデルでは、図１２に示すように×で示すイベント(dead)群と○で示す打切り (alive)群が明確には分離していないが、図１１の生存時間解析モデルでは綺麗に分離していることが分かる。
この結果によって本実施形態に係る方法は打切りデータにバイアスを入れることなく、適切に扱うことのできるデータ解析方法であることが分かる。
【図面の簡単な説明】
【００９７】
【図１】図１は本発明の第１実施形態に係る基本的なフローチャートである。
【図２】図２は第１次説明変数選抜処理に係るフローチャートの一例である。
【図３】図３は打切りサンプルの仮イベント発生時間の設定処理に係るフローチャートの一例である。
【図４】図４は目的変数の設定処理に係るフローチャートの一例である。
【図５】図５（Ａ）（Ｂ）はそれぞれ、多様性を持った説明変数の設定の一例を示す図である。
【図６】図６は部分最小自乗法モデル（ＰＬＳモデル）を最適化する処理の一例を示すフローチャートである。
【図７】図７は交差検証成績ＣＶの計算に係るフローチャートの一例である。
【図８】図８は実施例１の実験結果を示す図である。
【図９】図９は実施例２の実験結果を示す図である。
【図１０】図１０は実施例３の実験結果を示す図である。
【図１１】図１１は実施例１に係るデータを用いて本実施形態に係る方法によって得られたハザード値の分布を示した図である。
【図１２】図１２は実施例１に係るデータを用いて前記特許文献２に記載された方法によって得られたハザード値の分布を示した図である。
【図１３】図１３は本実施形態に係るデータ解析システムの一構成例を示す概略ブロック図である。
【符号の説明】
【００９８】
１…コンピュータ、２…打切りサンプル仮設定部、３…データ記憶部、４…ＰＬＳモデル最適化部、５…第１次説明変数群記憶部、６…比例ハザートモデル最適化部。
【図１】
2005222422.tif 000007
【図２】
2005222422.tif 000008
【図３】
2005222422.tif 000009
【図４】
2005222422.tif 000010
【図５】
2005222422.tif 000011
【図６】
2005222422.tif 000012
【図７】
2005222422.tif 000013
【図８】
2005222422.tif 000014
【図９】
2005222422.tif 000015
【図１０】
2005222422.tif 000016
【図１１】
2005222422.tif 000017
【図１２】
2005222422.tif 000018
【図１３】
2005222422.tif 000019