JPWO2003085548A1 - データ解析装置および方法 - Google Patents
データ解析装置および方法 Download PDFInfo
- Publication number
- JPWO2003085548A1 JPWO2003085548A1 JP2003582665A JP2003582665A JPWO2003085548A1 JP WO2003085548 A1 JPWO2003085548 A1 JP WO2003085548A1 JP 2003582665 A JP2003582665 A JP 2003582665A JP 2003582665 A JP2003582665 A JP 2003582665A JP WO2003085548 A1 JPWO2003085548 A1 JP WO2003085548A1
- Authority
- JP
- Japan
- Prior art keywords
- variable
- state
- data analysis
- cross
- living body
- 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
-
- 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
- G16B40/00—ICT specially adapted for biostatistics; ICT specially adapted for bioinformatics-related machine learning or data mining, e.g. knowledge discovery or pattern finding
- G16B40/10—Signal processing, e.g. from mass spectrometry [MS] or from PCR
-
- 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
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
-
- 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
-
- 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
- G16B40/00—ICT specially adapted for biostatistics; ICT specially adapted for bioinformatics-related machine learning or data mining, e.g. knowledge discovery or pattern finding
- G16B40/20—Supervised data analysis
-
- 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
- G16B40/00—ICT specially adapted for biostatistics; ICT specially adapted for bioinformatics-related machine learning or data mining, e.g. knowledge discovery or pattern finding
Abstract
生体の状態と複数の遺伝子発現の量および/または細胞内物質の量との相関モデルを決定するデータ解析において、生体の状態または時間とともに確率的に発生する生体の状態の変化を目的変数とし、複数の遺伝子発現の量および/または細胞内物質の量を説明変数とするデータの集合において、データに含まれる説明変数を選択し、選択された説明変数と目的変数とを含む相関モデルについて交差検証成績を計算し、その結果を評価判定する。ここで、交差検証成績が改善しなくなるまで、説明変数の選択、交差検証成績の計算、その結果の評価判定を行い、部分最小自乗法モデルを決定する。これにより、多変量の遺伝子発現情報の効果的な情報処理を提供する。
Description
技術分野
本発明は、生体の状態と遺伝子発現の量および/または細胞内物質の量との多変量解析処理並びそれを基に可能となる測定機材、検定方法などに関するものである。
背景技術
2000年6月のヒトゲノムの解読宣言以降、ゲノムに書かれた遺伝情報がどのように発現して機能しているかのを解明するポストゲノム時代に突入したと言われている。ヒトゲノム計画の進展の中で、ゲノム発現状態を測定する方法論も進展してきた。トランスクリプトーム(mRNA)測定手段としてオリゴヌクレオチドアレイやマイクロチップが知られている。またプロテオーム(蛋白質)測定手段として、以前からある2次元電気泳動に加えて、最近では質量分析の方法が進歩してきた。また抗体チップなどの先進の技術も注目されている。これらの測定技術は、生体の状態パラメータを短時間に一挙に測定できることがそれまでの技術と比較して画期的であるといえる。
遺伝子発現状態を効率的に測定する技術として次のものがあげられる。トランスクリプトーム(mRNAの総体)を特定するものとして、基盤に複数種のDNAを担持し、それに相補的なmRNAを検出するDNAチップが知られている。代表的なDNAチップには、遺伝子チップやDNAマイクロアレイがある。また、プロテオーム(蛋白質の総体)を特定するものには、2次元電気泳動、抗体チップ、質量スペクトルを用いるものがある。またメタボローム(代謝中間体を含めた代謝産物の総体)を測定する手法も質量分析などによって試みられており、進展が見られる。
生体内の細胞の状態は遺伝子産物の発現によってよく記述されるため、従来の診断マーカーでは情報が不足している場面でも、精度のより高い診断が可能になるという期待も出てきている。たとえば、次のような研究があげられる。
P.O.Brownらは、DNAチップによってリンパ腫患者の細胞のトランスクリプトームを測定し、クラスター解析によって悪性と良性のリンパ腫(DLBCL)を別クラスターに分離した(Nature 403(3),503−11(2000))。しかし、これは因果関係(相関関係)のモデルを得る方法ではなく、どの遺伝子がどの程度重要かを判断できない。
A.Alaiyaらは、2次元電気泳動によって子宮がん患者40人の細胞のプロテオームを測定し、うち22人のデータから部分最小自乗法診断モデルを構築し、悪性度を説明した(Int.J.Cancer,86,731−36(2000);Electrophoresis,21,1210−17(2000);国際公開WO 00/70340)。その際、全変数モデルにおいて1553変数からloadingの大きな170変数に限定することによって交差検証成績がよくなり(Q2=0.84)、残り18患者の深刻度(3段階)を11/18の比率で正答した。交差検証法がモデル構築の際の指標になるという考えが表明されている。しかし、この方法では、loadingを得る際にまず全変数モデルが成立しなければならない。また、それ以外の変数選択手法が考案されていない。
J.Khanらは、DNAチップによって小児がん患者の細胞を測定し、ニューラルネットワークによって悪性度を説明した(Nature Medicine,7(6),673−79(2001))。小児がん(SRBCT)患者88人のトランスクリプトーム(6567遺伝子)を測定し、うち63人のデータから主成分分析によって10次元に圧縮し、次に、人工ニューラルネットワーク診断モデルを構築した。ここで、影響力のある上位遺伝子を交差検証法によって絞り込み、96遺伝子で最良の成績(100%)を得た。このモデルで残り25人を予測し、93〜100%の結果を得た。しかし、この方法でも、影響力を得る際にまず全変数モデルが成立しなければならない。またそれ以外の変数選択手法が考案されていない。10次元のような少ない変数の場合を扱えるが、変数の数が膨大な場合には適用できない。
また、最近になってDNAチップの解析に部分最小自乗法を用いる研究がD.M.RockeとD.V.Nguyenによって報告されるに至った(国際公開WO 02/25405;Bioinformatics 18(1),39−50(2002);Bioinformatics 18(9),1216−26(2002);Bioinformatics 18(12),1625−32(2002))。部分最小自乗法の潜在変数を線型判別分析などの多変量解析の説明変数として用いた場合に良好な結果が得られることが報告されている。これは部分最小自乗法が次元圧縮とモデルフィットを同時に行なうことのできる方法であるために可能となったものである。報告に示された実施例では部分最小自乗法がDNAチップ情報のモデル構築方法として優れたものであることが示されている。しかし報告においては重要な遺伝子発現量を選抜する手段としての最小自乗法の適用については触れられておらず、事前の前処理によって選択された説明変数を全て用いて解析が行なわれているという点において上述のA.Alaiyaらの研究と同様の課題を含んでいる。
従来の診断マーカーでは情報が不足している場面でも、遺伝子発現情報を活用することで、より精度(解像度)の高い診断が可能になるという期待も出てきている。遺伝子発現状態の測定結果は、膨大な情報量が得られることが従来にはなかった特徴であり、逆に情報量が多いために、効果的なデータ処理なくしてデータの活用はありえない。したがって、有用な知識を獲得するためには効果的な情報処理が欠かせない。前に説明したように、現状ではクラスター解析を中心とする方法が用いられているが、主成分分析などの方法も採用されている。クラスター解析や主成分分析は、教師付学習方法ではないため、病状の因果関係(相関関係)のモデルを得ることはできない。すなわち、どの遺伝子がどの程度重要かを解析結果から得ることができないのが難点である。一方、部分最小自乗法は次元圧縮とモデルフィットを同時に行なう強力な多変量解析手法であるが、変数の数が膨大になった場合にしばしば有意な結果が得られない事態に直面する。したがって、膨大な遺伝子発現情報などから有用な知識を獲得できるような効果的な情報処理が望まれている。また、そのような情報処理の結果を基にした効率的な測定機材、検定処理などが期待されている。
発明の開示
(発明が解決しようとする技術的課題)
この発明の目的は、多変量の遺伝子発現情報、細胞内物質情報の効果的な情報処理を提供することである。
また、この発明の目的は、効率的な検定処理を提供することである。
(その解決方法)
本発明に係るデータ解析装置は、生体の状態または時間とともに確率的に発生する生体の状態の変化を目的変数とし、複数の遺伝子発現の量および/または細胞内物質の量を説明変数とする相関モデルを決定するデータ解析装置であって、生体の状態或いはそれを導出するデータまたは時間とともに確率的に発生する生体の状態の変化に関するデータと、複数の遺伝子発現の量および/または細胞内物質の量からなるサンプルの集合を入力する入力手段と、(1)説明変数を選択する選択手段と、(2)部分最小自乗法を実行して交差検証成績を計算する計算手段または上記生体の状態の変化に関するデータにカプラン・マイヤー法又はカトラー・エデラー法による生命表を適用して変化の発生しなかったものの確率を計算して得られた確率を、仮定した分布に基づいた変換または仮定を前提としない変換をし、該変換結果を目的変数とする部分最小自乗法を実行して交差検証成績を計算する計算手段と、(3)上記(2)の計算手段の結果を評価し、説明変数の採用、不採用を判定する評価判定手段とを有し、(4)上記(1)の選択手段と上記(2)の計算手段と上記(3)の評価判定手段とを実行して部分最小自乗法モデルの少なくとも交差検証成績を独立変数として持つ関数を改善し続けて部分最小自乗法モデルを決定する決定手段とからなる。選択手段は、たとえば、説明変数を逐次取捨選択したり、遺伝的アルゴリズムを用いて説明変数を選択する。計算手段は、たとえば、1個のサンプルを逐次除外したり、複数のサンプルを逐次除外して部分最小自乗法を実行して交差検証成績を計算する。評価判定手段は、たとえば、計算手段の結果から、各計算において除外したサンプルの遺伝子発現から予測される生体の状態を示す目的変数値と、前記除外したサンプルの生体の状態を示す目的変数値との誤差の代表値を求め、当該誤差の代表値が小さくなった場合に、その交差検証成績が改善されたと判定し、説明変数を取捨選択しながら交差検証成績の評価判定を繰り返す。あるいは交差検証成績ではなく、少なくとも部分最小自乗法モデルの交差検証成績を独立変数として持つ関数が改善するかどうかを評価判定の基準として用いることもできる。決定手段は、たとえば、選択手段と計算手段と評価判定手段とを繰り返し実行して部分最小自乗法モデルの交差検証成績を改善し続けて部分最小自乗法モデルを決定する。また、選択手段と計算手段とを複数のコンピュータで実行させることもできる。こうして、相関モデルを構成するとき、交差検証成績を基準に最適化させることにより説明変数を取捨選択し、説明変数の次元を減らして良好なモデルを得る。
上述の、仮定した分布に基づいた変換または仮定を前提としない変換は、生体の状態の変化の確率が説明変数の多項式で解析できるようにするために行なうものである。分布を仮定した場合には、確率を対数変換後に負の数にしたものを状態の変化を観測した時間で割るという変換、確率を対数変換後に負の数にしたものをさらに対数にしたものを状態の変化を観測した時間で割るという変換、または確率を1より減じたものをプロビット変換したものを計算して状態の変化を観測した時間で割るという変換などが考えられる。一方、分布を仮定しない場合にはロジット変換といった方法が考えられる。変換の方法は分布にどのような仮定が成り立つかどうかあるいはなりたたないかどうかを判断することにより、それぞれの場合に応じて適切に選ぶことができる。少なくとも部分最小自乗法モデルの交差検証成績を独立変数として持つ関数としては、たとえば、前記誤差の代表値と選抜された説明変数の数の関数が考えられ、あるいはその他の独立変数を含むものであってよい。望ましくは、関数は誤差の代表値の単調減少関数であり、説明変数の数の単調減少関数である。計算量を増やさないためには簡単に計算できる関数が望ましい。具体的には−PRESS×alphaNPという関数が考えられる。ここでPRESSは予測残差自乗和であり、NPは採用された説明変数の数であり、alphaは1または1より大きい実数である。また、−PRESS×(NP+beta)gammaや−PRESS×(beta−NP)−gammaなる関数も考えられる。ここで、gammaは正の実数である。
説明変数の個数を少なくすると、通常の統計的手法または多変量解析手法が適用可能になる。本発明では部分最小自乗法を用いて選抜された説明変数を統計手法又は多変量解析手法の説明変数として、より良好なモデルを得る。或いは選抜された説明変数を用いた部分最小自乗法モデルの潜在変数を統計手法又は多変量解析手法の説明変数として、より良好なモデルを得る。ここで潜在変数とは、部分最小自乗法において通常用いられているものであって、目的変数(Yil)と説明変数(Xij)の背後に共通する次元数の少ない潜在変数(Tik)を抽出することが部分最小自乗法の次元圧縮であり、モデルフィットである。
Yil=Σ Qkl×Tik+Fil
Xij=Σ Pkj×Tik+Eij
(iはサンプル番号、lは目的変数番号、jは説明変数番号、kは潜在変数番号、F,Eは残差)
また、統計的手法又は多変量解析手法としては、重回帰分析法、線型判別分析法、適応最小自乗法、ロジスティック回帰分析法、比例ハザード解析法、マハラノビス距離を用いる判別分析法、kNN法、人工ニューラルネットワークなどが挙げられる。
本発明者等は、また、Q2やPRESS値などの交差検証成績に加えて、説明変数の個数を第2の独立変数として含む関数を最適化することで選抜される説明変数を任意に絞り込むことができることを新たに見出した。通常の統計的手法や多変量解析手法では、抽出される説明変数の個数NPの望ましい範囲がサンプル数との兼ね合いで決まっている場合がある。そのような場合、関数を、目的とする選抜数によって任意に変更できる。関数形をたとえば−PRESS×alphaNPとした場合、説明変数の個数を数個から数十個に絞り込むためには通常は定数alphaとして1.0〜3.0の値が望ましい。より望ましくは、alphaは1.0〜2.0の値となる。他の関数形f(PRESS,NP)であっても、実際に選択される説明変数の数MPおよびその時のPRESS値PRESS_MPの周辺で、f(PRESS_MP÷alpha,MP+1)≒f(PRESS_MP,MP)となるような関数は、変数選択という点では同様の効果を持つ場合がある。こうして、適当な関数形を用いることにより、望ましい範囲の個数NPの説明変数を選抜できる。このようにして、交差検証成績を用いて決定されたモデルに採用されている説明変数をさらに絞り込むと、統計的手法又は多変量解析手法によるモデルを構築できる。したがって、その性質が十分解明されている統計的手法又は多変量解析手法を採用して解析を加えることができる。
また、目的変数として、時間とともに確率的に発生する生体の状態の変化から導出された量を用いて、時間とともに確率的に発生する生体の状態の変化と複数の遺伝子発現の量および/または細胞内物質の量との相関モデルを決定できる。「時間とともに確率的に発生する生体の状態の変化」とはたとえば生存時間である。ここで、前述の部分最小自乗法に、カプラン・マイヤー法又はカトラー・エデラー法と、ロジット(logit)変換とを組み合わせる。部分最小自乗法での目的変数は、時間とともに確率的に発生する生体の状態の変化に関するデータにカプラン・マイヤー法又はカトラー・エデラー法による生命表を適用して変化の発生しなかったものの確率を計算し、これをロジット変換した値である。ロジット(logit)値とは、分類分けされたデータの、ある分類の割合(確率)Pを基に、次式logit=log{P/(1−P)}にて計算される値である。ロジット値を目的変数とする部分最小自乗法を実行して交差検証成績を計算する。こうして、先に説明したのと同様に、部分最小自乗法の交差検証成績を考慮した説明変数の抽出を行って、生存時間解析を行える。
説明変数の個数を少なくすると、通常の統計的手法または多変量解析手法が適用可能になる。そこで、決定されたモデルに採用されている説明変数又はその潜在変数を用い、時間とともに確率的に発生する生体の状態の変化を説明する統計的手法又は多変量解析手法によるモデルを構築する。たとえば、ロジット値を目的変数として求めた説明変数を用いて、他の統計的手法又は多変量解析手法(たとえば比例ハザード法や、パラメトリックな分布にあてはめた回帰分析法)を行なうことによって、より良好なモデルを得ることができる。比例ハザード法とは、Coxによって考案された方法であり、生存率の解析に時間を考慮し、かつ、多変量を扱える。比例ハザード法では、観測されている個々ごとにハザード値と呼ばれる生存率を左右する値があり、それを導く関数がある(モデルが仮定されている)として解析される。カプラン−マイヤー法は、集団全体または群ごとの生存率の推移を示す。また、パラメトリックな分布とは、ガウスが提案した正規分布から計算された確率分布のことであり、生存時間解析では指数分布、ワイブル分析、対数正規分布が用いられる。指数分布などへの当て嵌めで、数式中に多項式があり、前述の部分最小自乗法の交差検証成績を考慮した説明変数の抽出が適用される。
入力手段で説明変数として入力される複数の遺伝子の発現量および/または細胞内物質の量とは、必ずしも物質の絶対的な濃度の測定値に限定されるものではなく、加工計算された値、相対的な値、間接的に物質量を表す量などでもよい。たとえば、質量スペクトルで蛋白質の発現量を測定することができることを応用して、生体の状態を表わす目的変数と、質量スペクトルとを直接関係づける相関モデルを構築することができる。またAffymetrix社タイプのDNAチップ(ジーンチップ)では、単一のスポットが単一の遺伝子発現を特定するとは限らず、複数個のスポットが集まってはじめて単一の遺伝子発現を特定することもある。ここでもまた、各スポットの測定量を説明変数として、直接、生体の状態を説明する相関モデルを得ることができる。更には、タンパク質の電気泳動パターンの各ピークは単一のタンパク質に帰属できず、複数個のタンパク質の重ねあわせであることも多い。このような場合にも生体の状態を説明する説明変数として各ピーク強度を用いることができる。このことは、上述のAlaiyaらは子宮癌の診断の説明変数として電気泳動パターンのピーク強度を採用していることから明らかである。前述のようにポストシークエンス時代のトランスクリプトーム解析、プロテオーム解析、メタボローム解析という研究分野では、生体(細胞)内の物質を総体として把握することから出発することを特徴とする実験的アプローチが注目されている。ひとつひとつの物質の絶対的定量は必須事項ではなく、これらの実験方法によって定量される物質の量を直接、間接に表現する測定値やその加工計算値が、生体の状態を説明する説明変数と成り得る。また以上の物質量を表現する説明変数以外に、場合によっては問診データなどの他の説明変数を追加すると、さらに有効な解析結果が得られる場合もある。
本発明に係るデータ解析方法は、生体の状態または時間とともに確率的に発生する生体の状態の変化を目的変数とし、複数の遺伝子発現の量および/または細胞内物質の量を説明変数とする相関モデルを決定するデータ解析方法であって、生体の状態或いはそれを導出するデータまたは時間とともに確率的に発生する生体の状態の変化に関するデータと、複数の遺伝子発現の量および/または細胞内物質の量からなるサンプルの集合を入力する入力ステップと、(1)説明変数を選択する選択ステップと、(2)部分最小自乗法を実行して交差検証成績を計算する計算ステップまたは前記生体の状態の変化に関するデータにカプラン・マイヤー法又はカトラー・エデラー法による生命表を適用して変化の発生しなかったものの確率を計算して得られた確率を、仮定した分布に基づいた変換または仮定を前提としない変換をし、該変換結果を目的変数とする部分最小自乗法を実行して交差検証成績を計算する計算ステップと、(3)前記(2)の計算ステップの結果を評価し、説明変数の採用、不採用を判定する評価判定ステップとを有し、(4)前記(1)の選択ステップと前記(2)の計算ステップと前記(3)の評価判定ステップとを実行して部分最小自乗法モデルの少なくとも交差検証成績を独立変数として持つ関数を改善し続けて部分最小自乗法モデルを決定する決定ステップとからなる。
このデータ解析方法において、選択ステップは、たとえば、説明変数を逐次取捨選択したり、遺伝的アルゴリズムを用いて説明変数を選択する。計算ステップは、たとえば、1個のサンプルを逐次除外したり、複数のサンプルを逐次除外して部分最小自乗法を実行して交差検証成績を計算する。評価判定ステップは、たとえば、計算ステップの結果から、各計算において除外したサンプルの遺伝子発現から予測される生体の状態を示す目的変数値と、前記除外したサンプルの生体の状態を示す目的変数値との誤差の代表値を求め、当該誤差の代表値が小さくなった場合に、その交差検証成績が改善されたと判定し、説明変数を取捨選択しながら交差検証成績の評価判定を繰り返す。決定ステップは、たとえば、選択ステップと計算ステップと評価判定ステップとを繰り返し実行して部分最小自乗法モデルの交差検証成績を改善し続けて部分最小自乗法モデルを決定する。また、選択ステップと計算ステップとを複数のコンピュータで実行させることもできる。
本発明に係るデータ解析プログラムは、生体の状態または時間とともに確率的に発生する生体の状態の変化を目的変数とし、複数の遺伝子発現の量および/または細胞内物質の量を説明変数とする相関モデルを決定する、コンピュータにより実行されるデータ解析プログラムであって、生体の状態或いはそれを導出するデータまたは時間とともに確率的に発生する生体の状態の変化に関するデータと、複数の遺伝子発現の量および/または細胞内物質の量からなるサンプルの集合を入力する入力ステップと、(1)説明変数を選択する選択ステップと、(2)部分最小自乗法を実行して交差検証成績を計算する計算ステップまたは前記生体の状態の変化に関するデータにカプラン・マイヤー法又はカトラー・エデラー法による生命表を適用して変化の発生しなかったものの確率を計算して得られた確率を、仮定した分布に基づいた変換または仮定を前提としない変換をし、該変換結果を目的変数とする部分最小自乗法を実行して交差検証成績を計算する計算ステップと、(3)前記(2)の計算ステップの結果を評価し、説明変数の採用、不採用を判定する評価判定ステップとを有し、(4)前記(1)の選択ステップと前記(2)の計算ステップと前記(3)の評価判定ステップとを実行して部分最小自乗法モデルの少なくとも交差検証成績を独立変数として持つ関数を改善し続けて部分最小自乗法モデルを決定する決定ステップとからなる。
このデータ解析プログラムにおいて、選択ステップは、たとえば、説明変数を逐次取捨選択したり、遺伝的アルゴリズムを用いて説明変数を選択する。計算ステップは、たとえば、1個のサンプルを逐次除外したり、複数のサンプルを逐次除外して部分最小自乗法を実行して交差検証成績を計算する。評価判定ステップは、たとえば、計算ステップの結果から、各計算において除外したサンプルの遺伝子発現から予測される生体の状態を示す目的変数値と、前記除外したサンプルの生体の状態を示す目的変数値との誤差の代表値を求め、少なくとも当該誤差の代表値を独立変数として持つ関数である当該誤差の代表値の単調減少関数の値が小さくなった場合に、その交差検証成績が改善されたと判定し、説明変数を取捨選択しながら交差検証成績の評価判定を繰り返す。決定ステップは、たとえば、選択ステップと計算ステップと評価判定ステップとを繰り返し実行して少なくとも部分最小自乗法モデルの交差検証成績を独立変数として持つ関数を改善し続けて部分最小自乗法モデルを決定する。また、選択ステップと計算ステップとを複数のコンピュータで実行させることもできる。さらには、前記の説明変数の選択において、たとえば、初期状態では説明変数を全く含まないか、或いは、初期状態では全説明変数を含むこともできる。
前記のデータ解析プログラムにおいて、上記の生体の状態は、たとえば病気のタイプをあらわす測定値、病気の重篤度をあらわす測定値、病気のタイプをあらわす医療診断の結果、病気の重篤度をあらわす医療診断の結果、あるいはそれらを2次加工した数値である。例えば後の実施例で示すように、患者の生存時間を予測することは、QOL(quality of life:生活の質)を含めた治療計画や人生設計などを判断する上で重要な情報をもたらすものであり、社会的に価値のある診断モデルを提供することができる。また癌の再発可能性を予測することは、QOLを考慮した治療計画を立案し、医師または当の患者が選択の判断をするうえで、貴重な情報をもたらすものである。
また、本発明は、決定された前記相関モデル及び予測対象のサンプルについて当該モデルにおいて採用された説明変数を入力する入力手段と、入力された該説明変数に基づいて該サンプルの生体の状態を予測判定する予測判定手段からなるデータ解析装置、前記で決定された相関モデル及び予測対象のサンプルについて当該モデルにおいて採用された説明変数を入力する入力ステップと、入力された該説明変数に基づいて該サンプルの生体の状態を予測判定する予測判定ステップからなるデータ解析方法及び前記で決定された相関モデル及び予測対象のサンプルについて当該モデルにおいて採用された説明変数を入力する入力ステップと、入力された該説明変数に基づいて該サンプルの生体の状態を予測判定する予測判定ステップからなるデータ解析プログラムも包含する。
本発明に係るコンピュータにより読取可能な記録媒体は、上記のいずれかのプログラムを記録する。
本発明に係るびまん性大細胞型Bリンパ腫の重篤度検定用の細胞内物質測定機材および測定方法並びにびまん性大細胞型Bリンパ腫の重篤度検定方法は、実質的にジーンバンクアクセッション番号がU15085、M23452、X52479、U70426、H57330及びS69790からなる遺伝子群の発現を検出する。さらに、ジーンバンクアクセッション番号がU03398、M65066、AK001546、BC003536、X00437、U12979、H96306、AA830781及びAA804793からなる群から選択される少なくとも一つの遺伝子の発現を検出してもよい。
また、本発明に係る乳癌の重篤度検定用の細胞内物質測定機材および測定方法並びに乳癌の重篤度検定方法は、実質的にジーンバンクアクセッション番号がAA598572、AA703058及びAA453345からなる遺伝子産物を含む細胞内物質を検出する。さらに、ジーンバンクアクセッション番号がAA406242、H73335、W84753、N71160、AA054669、N32820及びR05667からなる群から選択される少なくとも一つの遺伝子産物を含む細胞内物質を検出してもよい。
また、本発明に係る乳癌の再発性検定用の細胞内物質測定機材および測定方法並びに乳癌の再発性検定方法は、実質的にジーンバンクアクセッション番号がW84753、H08581、AA045730及びAI250654からなる遺伝子産物を含む細胞内物質を検出する。さらに、ジーンバンクアクセッション番号がAA448641、R78516、R05934、AA629838及びH53037からなる群から選択される少なくとも一つの遺伝子産物を含む細胞内物質を検出してもよい。
また、本発明に係る乳癌の再発性検定用の細胞内物質測定機材および測定方法並びに乳癌の再発性検定方法は、実質的にジーンバンクアクセッション番号がAA434397、T83209、N53427、N29639、AA485739、AA425861、H84871、T64312、T59518及びAA037488からなる遺伝子産物を含む細胞内物質を検出する。さらに、ジーンバンクアクセッション番号がAA406231の遺伝子産物を含む細胞内物質を検出してもよい。
また、本発明に係る乳癌の再発性検定用の細胞内物質測定機材および測定方法並びに乳癌の再発性検定方法は、実質的にジーンバンクアクセッション番号がH11482、T64312及びAA045340からなる遺伝子産物を含む細胞内物質を検出する。
細胞内物質測定機材としては、DNAマイクロアレイ、ジーンチップ、オリゴDNA型のDNAチップ、電気化学DNAチップ(ECAチップ)、繊維型DNAチップ、磁性ビーズDNAチップ(PSS)、糸巻きDNAチップ(PSS)、などのDNAチップ、マクロアレイ、抗体チップ、測定用試薬キットなどが挙げられる。また、上記の機材を適宜組み込んだ測定機械であってもよい。
発明を実施するための最良の形態
以下、添付の図面を参照して本発明の実施の形態を説明する。
以下に、選択された生体の状態と遺伝子発現の量および/または細胞内物質の量との相関モデルの決定について説明する。ここで、遺伝子発現の用語は、mRNA発現(トランスクリプトーム)や、mRNAによる翻訳の結果として生じる蛋白質(プロテオーム)を含むものとして用いる。また、細胞内物質の量とはここではたとえば、代謝中間体を含めた代謝産物全部であるメタボロームを意味する。たとえば、トランスクリプトーム(mRNA)やプロテオーム(蛋白質)の解析において、各サンプルデータは、生体の状態と遺伝子発現の量などからなる。各サンプルはたとえば1000個以上の膨大な遺伝子発現の量を含む。生体の状態は、たとえば病気のタイプまたは病気の診断指標であるが、より一般的には生体情報であればよい。「病気の診断指標」には、病気の進行度合いのほか、病気のタイプ、重篤度、深刻度などの表現で表わされるものも含む。ここで、遺伝子発現の量などの測定データは膨大な情報量からなるので、コンピュータを用いた効率的な多変量解析が必要である。
データ収集において、予めいくつかのサンプルについて生体の状態(たとえば診断指標)を判定し、また、そのサンプルされたものから細胞液を獲得し、その細胞液中の多くの遺伝子産物の発現の量などを測定する。本発明の実施の形態のデータ解析では、こうして得られた遺伝子産物の発現の量などと生体の状態(たとえば診断指標)を入力し、相関モデル(たとえば部分最小自乗法モデル)を得る。ここで、コンピュータによる多変量解析プログラムを用いて、診断指標を目的変数とし、遺伝子発現の量および/または細胞内物質の量を説明変数とする因果関係型の解析を行なって、各説明変数の重要性や影響度に関する情報を得る。また、前記目的変数は、必ずしも測定値そのものである必要はなく、ロジット変換を行なった値や群を表す離散値を用いても良く、その場合、より有意な解析結果を得ることもできる。
本発明者らは、遺伝子発現による医療診断という分野において、データ解析における交差検証(cross validation)の成績を少なくとも独立変数のひとつとして持つ関数を最適化するように変数を選択することによって良好な相関モデル(たとえば部分最小自乗法モデル)が得られることを見出した。交差検証法では、手持ちのデータを複数群に分割し、その一部のデータ群(訓練集合)だけを使ってフィットしたモデルを用いて残る別のデータ群(テスト集合)を予測することによって、モデルの予測力を試す。通常の部分最小自乗法(PLS)においては潜在変数の次元選択に交差検証法が用いられているが、ここでは、部分最小自乗法において、潜在変数を1次元に固定し、1以上の入力変数(説明変数)を逐次取捨選択しながら、交差検証成績(たとえば平方和の予測誤差)を少なくとも独立変数のひとつとして持つ関数を最適化した。ただし本発明の効果は潜在変数の次元を1に限定するものではない。その結果、全変数を採用した場合には有意な相関モデルを得られなかった場合にも、良好でかつ予測力のある相関モデルが得られることが判明したのである。この交差検証法を用いた変数選択の逐次取捨選択により、安定な相関モデルが得られる。また本発明者らは、関数形を適切に設定することによって説明変数を絞り込むことにより、部分最小自乗法以外の統計学又は多変量解析の良好な相関モデルを得ることが可能となり、それぞれ生体の状態を記述する目的変数にふさわしい相関モデルを得ることができることを見出した。なお、ここでいう「最適化」とは、交差検証成績が、説明変数を取捨選択するための、そのときの解析条件の範囲で、改善がみられなくなるまで改良したことを意味しており、交差検証成績がすべての説明変数の組合せの中で最適なものを見出したという意味ではない。この変数選択手法を用いると、病状を決定する因子を少数に特定し、廉価な診断用材料(DNAチップ、抗体チップ、DNA含有ベクターなど)を設計でき、それ自体独自の価値を持つものである。また、この変数選択手法は、予め設定される各種の変数選択条件と共に運用することが可能である。
上に述べたように、説明変数は、交差検証成績を基準に逐次取捨選択される。ここで、取捨選択のため、交差検証成績を少なくとも独立変数のひとつとして持つ関数を用いる。説明変数を追加する場合は、その説明変数について、前記関数が改善されなかったと判定された場合には当該説明変数を除外し、改善されたと判定された場合には当該説明変数を追加する。また、説明変数を除外する場合は、その説明変数について、前記関数が改善されなかったと判定された場合には当該説明変数を除外せず、改善されたと判定された場合には当該説明変数を除外する。ここで、1以上の説明変数を選択した場合に、交差検証成績評価は次のように進める。n個のサンプルからいくつかのサンプルを逐次除外して部分最小自乗法モデルを求め、各モデルにおいて除外したサンプルの遺伝子発現の量から予測される生体の状態を示す目的変数と、除外したサンプルの生体の状態を示す目的変数との各々の誤差の代表値を求める。「代表値」とは、和、平均、最大値、中位値、最頻値などのデータを特徴づける値をいう。そして、当該誤差の代表値を少なくともひとつの独立変数とする関数が小さくなった場合に、交差検証成績が改善されたと判定し、当該説明変数を追加または削除する。この交差検証成績評価を、説明変数を取捨選択しながら逐次繰り返して、前記関数を改善し続ける。改善されなくなれば交差検証成績を最適化したとして説明変数の取捨選択を終了する。その結果、取捨選択により絞り込んだ数の説明変数からなる最適な部分最小自乗法モデルが得られる。具体的には、計算手段において計算される交差検証成績の数値指標として予想残差自乗和(PRESS)を採用し、評価判定手段において予想残差自乗和の値が説明変数あたり一定の閾値以下の比率で小さくなる場合に、その説明変数を採用すると判定することにより、上記の処理は実行可能である。
因果関係型の解析手法においてはオーバーフィット(over fitting)を避けるための工夫が必要となる。ここでいうオーバーフィットとは、説明変数が多すぎるためにたまたま予測結果と実績とが一致するものの、本当の相関関係をとらえ損なっているため、モデルフィットに用いたデータ以外に予測能力を持たないことをいう。ここでは、相関モデルとして部分最小自乗法を用いるが、部分最小自乗法は次元圧縮とモデルフィットを同時に行なう強力な多変量解析手法であり、オーバーフィットの問題に比較的強いとされている。しかし遺伝子発現状態解析のように膨大な変数を扱う場合には、有意な結果が得られない事態に直面する。従来技術として説明したAlaiyaやKhanの手法は全変数モデルが有意に成立することを前提としているので、変数の絞込みには一般的には適用できない。これに対し、本発明では、交差検証予測結果を最適にするように変数を絞り込むことにより、オーバーフィットを減らすことができた。また、本発明は、前記Khanの手法とは異なり、主成分分析などの前処理を介さない方法である。従来技術では、説明変数が膨大な場合には、有意なモデルを得ることができないことから、予め、全説明変数を基にたとえば、主成分分析などで次元圧縮する前処理をし、これによって得られた説明変数によって解析する方法が用いられる。しかし、この方法では、構成したモデルで予測を行なうためには、モデル構成の基となった全説明変数が必ず必要となり、たとえば、説明変数が遺伝子発現の量であれば、診断用遺伝子チップに担持する遺伝子としては、モデル構成に用いた遺伝子の全てが必要となるか、または別の手法を用いて変数選択することが必要となる。一方、本発明においては、説明変数の選択によって説明変数を絞り込んでいるので、たとえば、説明変数が遺伝子発現の量であれば、診断用遺伝子チップに担持する遺伝子は、選択された説明変数に相当する遺伝子を担持すれば良いことになる。
なお、Todeschiniらは、有機化合物の大気中の分解を予測するため、遺伝的アルゴリズムによって交差検証成績を最適化するように変数選択を行ない、重回帰モデルを得ている(P.Gramatics,V.Consonni & R.Todeschini,Chemosphere 38(5),1371−78(1999))。53化合物と175記述子でモデル構築を行ない(Q2=0.79)、7変数が選択され、98化合物の予測を行なった(Q2=0.75)。交差検証成績を最適化するように変数選択を行なっている点では、本実施形態と同様の手法である。しかし、重回帰モデルを採用しているために、説明変数の選択過程を通じて選択される変数は少数個にとどまらざるを得ず、複数の遺伝子発現の量および/または細胞内物質の量の解析には適用できない。本発明者らの調査した範囲では、Q2やPRESS値を最適化する方法では、選抜される説明変数は百程度から数百程度にわたり、重回帰モデルでは解析が不能となる。またTodeschiniらは、説明変数を絞り込むための有効な方法について言及していない。これは、もともとの説明変数の候補がたかだか175個であり、説明変数を絞り込むために特別の工夫をする必要がないからである。遺伝子発現解析の分野はこれとは全く異なり、数十から数百のサンプル数に対して、数百から数千、数万の説明変数候補が存在する。したがってこれまでとは異なる工夫が必要となる。
本実施形態では、生体の状態と複数の遺伝子発現の量および/または細胞内物質の量との相関モデルを決定するとき、交差検証成績を少なくとも独立変数のひとつとして持つ関数を最適化させるように説明変数を逐次追加・除外することによって、説明変数を選抜して、良好な相関モデルを得る。このようなアプローチの優位性は、下記の実施例から推測されるように、次のとおりである。
1)病気や生体現象の背後で働いている重要な遺伝子やメカニズムを推定/特定でき、理解が深まる。
2)重要な遺伝子産物や細胞内物質だけに絞った廉価な診断用材料(DNAチップ、抗体チップなど)の設計が可能になる。
本実施形態では、交差検証成績を少なくとも独立変数のひとつとして持つ関数を最適化するように説明変数を段階的に取捨選択するが、たとえば具体的には、ステップワイズ(step wise)法に代表される説明変数を選択する選択手段と、リーブ・ワン・アウト(leave−one−out)法に代表される交差検証法に部分最小自乗法を適用して計算する計算手段と、前記計算手段の結果を評価し、説明変数の採用、不採用を判定する評価判定手段とを組合せて用いる。すなわち、m個の説明変数の中から1以上の説明変数を選択し、次いで、部分最小自乗法を実行して交差検証成績を計算し、さらに、該計算結果を評価して、選択した説明変数の採用、不採用を判定する。この評価判定では、計算手段の結果から、各計算において除外したサンプルの遺伝子発現から予測される生体の状態を示す目的変数値と、前記除外したサンプルの生体の状態を示す目的変数値との誤差の代表値を求め、少なくとも当該誤差の代表値を独立変数として持つ関数である当該誤差の代表値の単調減少関数の値が小さくなった場合に説明変数の取捨選択を判定する。このように、選択手段と計算手段と評価判定手段とを用いて、少なくとも部分最小自乗法モデルの交差検証成績を独立変数として持つ関数を改善し続けて、その改善がみられなくなるまで改良し、部分最小自乗法モデルを決定する。なお、本実施形態では、サンプルを1個づつ逐次除外している(リーブ・ワン・アウト法)が、その代わりに、複数のサンプルを除外して交差検証成績を評価してもよい(リーブ・n・アウト法)し、また、Khan et al.により用いられた3分割法(three−fold)等の他の方法を用いることもできる。3分割法では、説明変数をランダムにシャッフルして3つのグループに分ける。その中の2つのグループを用いてモデルを構成し、残りの1つのグループでモデルを評価する。また、説明変数の選択方法としてはステップワイズ法、非線形アルゴリズム(たとえば遺伝的アルゴリズムなど)を用いてもよく、変数選択に関して予め何らかの条件が分っていれば、それに応じて探索範囲を限定できる。
次に、データの収集と解析について具体的に説明する。図1は、遺伝子発現解析システムを示す。データ収集のため、予めいくつかのサンプルについて診断指標(たとえば病気のタイプないし進行度合いを含む)を判定し、また、そのサンプルされたものから細胞液を獲得し、DNAチップを用いてその細胞液中の多くの遺伝子産物の発現の量を測定する。測定には、共焦点型レーザスキャナ(たとえばAffymetrix社、428アレイスキャナ)10を用いる。吸光度によりmRNAの量が測定される。このデータ収集は公知の方法である。測定データは、コンピュータ12に送られ解析される。コンピュータ12は、CPU14を備えた通常の構成のコンピュータであり、それに接続される記憶装置(たとえばハードディスク装置)16の記録媒体(たとえばハードディスク)には、測定データ18や解析ソフト20が格納される。この解析ソフト20を用いてデータ18が解析され、生体の状態と遺伝子発現の量などとの相関モデルが決定される。
なお、説明変数の選択と、交差検証法に部分最小自乗法を適用する計算とを複数のコンピュータで実行させてもよい。交差検証予測の計算を複数個のコンピュータに分散させることで計算を加速することができる。
図2は、コンピュータ12により実行される、生体の状態と遺伝子発現の量などとの相関モデルを得るためのデータ解析ソフト20のフローチャートを示す。ここでは簡単に説明するため、少なくとも部分最小自乗法モデルの交差検証成績を独立変数として持つ関数として−PRESSを採用しているが、発明の範囲を限定するものでなく、実施例2〜5においては別の関数を採用している。まず、相関モデル作成用のデータを入力する(S10)。データはたとえばDNAチップを用いて収集したものである。入力データ(サンプル集合)は、それぞれ目的変数(たとえば診断指標)とm個(たとえば2000個)の説明変数(たとえば遺伝子発現の量)からなる。また、場合によっては、上述のデータ(訓練集合)以外に、テスト集合のデータを入力する。ここでテスト集合とは交差検証の評価のためのデータ群を意味するのではなく、モデル決定が終了した後にモデルの予測力をテストするためのデータ群である。
まず、初期設定として、選択された説明変数の数を0とし、交差検証成績CVの最良値CV0を−∞とする(S12)。次に、説明変数の選択を行う。まず、説明変数を指す番号iを1とし(S14)、第i変数(遺伝子発現の量)を仮に採用して(S16)、部分最小自乗法を実行し、交差検証成績CVを計算する(S18、図3参照)。ここで、リーブ・ワン・アウト処理を用いる。これは、たとえば50個のサンプルからなる訓練集合において、1番から50番の全てを順次1個づつ除いて残りの49個のサンプルで予測した結果と、その時除いた1個の結果とを比較し、その誤差が大きい場合に、仮に選択した説明変数(第i変数)が適していないと判断する手法である。もし、得られた成績CVが現在の最良値CV0より最適化されれば(S20でYES)、第i変数を採用し、かつ、成績CVを新らしい最良値CV0に更新する(S22)。しかし、得られた成績CVが最良値CV0より大きくなければ(S20でNO)、第i変数を採用しない(S24)。そして、ステップS14に戻り、同様の処理を繰り返す。この処理を交差検証成績CVが改善されなくなる(S26でNO)まで繰り返す。ここで、相関モデルに採用する説明変数については1つづつ段階的に増加(追加)または減少(除外)して成績CVを評価判定している。すなわち、全体としての合致度合いがよくなるように各説明変数を解析に加えるかどうかを逐次判定しながら、説明変数の取捨選択を行い、これを、全体としての合致度合いがよくならなくなるまで繰り返す。以上の処理で改善があると、ふたたびステップS14の初め(i=1)に戻り、それまでに選択されている説明変数を基に、さらに説明変数の選択を繰り返す。なお、ここではモデルの予測力を判断するために、訓練集合とテスト集合とに予め分割しておいたデータ集合を用いてデータ解析しており、上述の解析は、訓練集合を用いて行なった結果であるので、この結果からテスト集合について予測を行い、実測データとの合致度を評価(S28)している。このような評価は必ずしも必要でないが、予測力を判断するには有効である。
図3は、リーブ・ワン・アウト処理を含む交差検証成績CVの計算(図2、S18)のフローチャートを示す。ここで、選択された変数について交差検証成績が計算される。まず、PRESSの初期値を0とする(S180)。次に、n個の集合内のサンプルを指す番号jを1とし(S182)、第jサンプル以外のn−1個のサンプルで部分最小自乗法を実行し(S184)、第jサンプルの目的変数を予測する(S186)。差の自乗を計算してPRESSに加算する(S190)。次に番号jを1増加し(S182)、同様の処理をおこなう。これを番号j=nまで各サンプルについて繰り返す。得られたPRESSは、1個のサンプルを順次除外して計算した予測値と実測値との差の平方和であり、予測誤差を表わす量である。この予測残差自乗和PRESSの符号を変えたものを交差検証成績CVとする(S192)。
本実施形態では、交差検証法を用いて、入力変数(説明変数)を段階的に1つづつ追加・除外しながら、交差検証成績(CV=−PRESS)を最適化する。ここで、説明変数の段階的な追加・除外の内容を理解しやすくするため、以下で、さらに具体的に5つのモデル構築手法について説明する。これらは、説明変数の逐次的な選択の手順が異なる。
図4は、第1のモデル構築手法を示す。データ集合においてどの説明変数も選択されていない状態を初期状態とする(S112)。次に、1番目の説明変数から最後(m番目)の説明変数までの未だ選択されていない説明変数ごとに逐次、その説明変数を選択した場合に交差検証成績が改善するかどうかを、リーブ・ワン・アウト処理を用いた交差検証成績評価ステップ(S118)を繰り返しながら判定(S120)し、改善する場合にはその説明変数を追加する(S114〜S124)。そのような改善と追加がなくなる(S126でNO)まで、1番目の説明変数から上記逐次判定操作を繰り返す。
さらに詳しく説明すると、まず、初期設定として、選択された説明変数の数NPを0とし、交差検証成績CVの最良値CV0を−∞とする(S112)。次に、説明変数の選択を行う。まず、変数iを1とし(S114)、第i変数を仮に採用する(S116)。ただし、第i変数がすでに採用されていれば(S115でYES)、ステップS114に戻る。次に、部分最小自乗法を実行し、交差検証成績CVを計算する(S118)。ここで、リーブ・ワン・アウト処理を用いる。もし、得られた成績CVが現在の最良値CV0より最適化されれば(S120でYES)、第i変数を採用し、かつ、成績CVを新らしい最良値CV0に更新する(S122)。しかし、得られた成績CVが最良値CV0より大きくなければ(S120でNO)、第i変数を採用しない(S124)。そして、ステップS114に戻り、同様の処理を繰り返す。この処理を交差検証成績CVが改善されなくなる(S126でNO)まで繰り返す。以上の処理で改善があると、ふたたびステップS114に戻り、新しいループを開始する。ここで、それまでに選択されている変数を基に、さらに変数の選択を繰り返す。こうして、データ集合を用いて選択された変数を用いた相関モデルが得られる。
図5は、第2のモデル構築手法を示す。この手法では、全ての説明変数が選択されている状態を初期状態とする(S212)。次に、1番目の説明変数から最後(m番目)の説明変数までの選択されている説明変数ごとに逐次、その説明変数を除外した場合に交差検証成績が改善するかどうかを、リーブ・ワン・アウト処理を用いた交差検証成績評価ステップ(S218)を繰り返しながら判定(S220)し、改善する場合にはその説明変数を除外する(S214〜S224)。そのような改善と除外がなくなる(S226でNO)まで、1番目の説明変数から上記逐次判定操作を繰り返す。
さらに詳しく説明すると、まず、初期設定として、選択された説明変数の数NPをmとし、交差検証成績CVの最良値CV0を−∞とする(S212)。すなわち、すべての説明変数を選択する。次に、説明変数の選択を行う。まず、変数iを1とし(S214)、第i変数を仮に除外する(S216)。ただし、第i変数がすでに除外されていれば(S215でYES)、ステップS214に戻る。部分最小自乗法を実行し、交差検証成績CVを計算する(S218)。ここで、リーブ・ワン・アウト処理を用いる。もし、得られた成績CVが現在の最良値CV0より最適化されれば(S220でYES)、第i変数を除外し、かつ、成績CVを新らしい最良値CV0に更新する(S222)。しかし、得られた成績CVが最良値CV0より大きくなければ(S220でNO)、第i変数を除外しない(S224)。そして、ステップS214に戻り、同様の処理を繰り返す。この処理を交差検証成績CVが改善されなくなる(S226でNO)まで繰り返す。以上の処理で改善があると、ふたたびステップS214に戻り、新しいループを開始する。ここで、それまでに選択されている変数を基に、さらに変数の選択を繰り返す。こうして、データ集合を用いて選択された変数を用いた相関モデルが得られる。
図6は、第3のモデル構成手法を示す。この手法は、第1と第2の手法の直列的な組合せである。まず、どの説明変数も選択されていない状態を初期状態とする(S112)。次に、1番目の説明変数から最後(m番目)の説明変数までの未だ選択されていない説明変数ごとに逐次、その説明変数を選択した場合に交差検証成績が改善するかどうかを、リーブ・ワン・アウト処理を用いた交差検証成績評価ステップを繰り返しながら判定し、改善する場合にはその説明変数を追加選択し、そのような改善と追加がなくなるまで1番目の説明変数から上記逐次判定操作を繰り返す(S114〜S126)。次に、1番目の説明変数から最後(m番目)の説明変数までの選択されている説明変数ごとに逐次、その説明変数を除外した場合に交差検証成績が改善するかどうかを、リーブ・ワン・アウト処理を用いた交差検証成績評価ステップを繰り返しながら判定し、改善する場合にはその説明変数を除外し、そのような改善と除外がなくなるまで、1番目の説明変数から上記逐次判定操作を繰り返す(S214〜S226)。
図7は、第4のモデル構築手法を示す。この手法は、第3の手法の変形である。まず、どの説明変数も選択されていない状態を初期状態とする(S112)。次に、1番目の説明変数から最後(m番目)の説明変数までの未だ選択されていない説明変数ごとに逐次、その説明変数を選択した場合に交差検証成績が改善するかどうかを、リーブ・ワン・アウト処理を用いた交差検証成績評価ステップ(S118)を繰り返しながら判定(S120)し、改善する場合にはその説明変数を追加選択する(S114〜S124)。そのような改善と追加がなくなる(S126でNO)まで、1番目の説明変数から上記逐次判定操作を繰り返す。次に、1番目の説明変数から最後(m番目)の説明変数までの選択されている説明変数ごとに逐次、その説明変数を除外した場合に交差検証成績が改善するかどうかを、リーブ・ワン・アウト処理を用いた交差検証成績評価ステップ(S218)を繰り返しながら判定(S220)し、改善する場合にはその説明変数を除外する(S214〜224)。そのような改善と除外がなくなる(S226でNO)まで、1番目の説明変数から上記逐次判定操作を繰り返す。上記逐次判定追加改善ステップまたは上記逐次判定除外改善ステップで少なくとも一度改善があれば(S227でYES)、ステップS112に戻り、上記操作(S112〜S227)を繰り返す。これを改善がなくなる(S227でNO)までおこなう。
図8は、第5のモデル構築手法を示す。この手法は、第1と第2のスキームの並列的な組合せである。どの説明変数も選択されていない状態を初期状態とする(S112)。次に、1番目の説明変数から最後(m番目)の説明変数までの説明変数ごとに逐次、その説明変数が選択されていない場合にはその説明変数を選択した場合に交差検証成績が改善するかどうかを、リーブ・ワン・アウト処理を用いた交差検証成績評価ステップ(S118)を繰り返しながら判定(S120)し、改善する場合にはその説明変数を追加する(S114〜S124)。また、選択する説明変数ごとに、その説明変数がすでに選択されている場合には、その説明変数を除外した場合に交差検証成績が改善するかどうかを、リーブ・ワン・アウト処理を用いた交差検証成績評価ステップ(S218)を繰り返しながら判定(S220)し、改善する場合にはその説明変数を除外する(S216〜S224)。そのような改善と追加または除外がなくなる(S126でNO)まで、1番目の説明変数から上記逐次判定操作を繰り返す。
次に、第4のモデル構築手法(図7)を適用した場合を、表1のデータ集合を例として説明する。このデータ集合に対して、部分最小自乗法による解析を用いて相関モデルを求める。表1のデータでは、サンプルの数nは10であり、また、説明を容易にするため、説明変数の数mは19と少なくしている。表1において、p1は目的変数を表わし、p2〜p20は説明変数を表わす。(ただし表1では、表示の便宜のため、p16以降のデータを省略している。)第4手法(図7)のステップS114、S214とは異なり、説明変数を表わすiはp20からp2まで逆に逐次処理することとした。CV評価値としてここでは予測残差自乗和(PRESS)を採用した。PRESSが小さいほど、CV評価値はよい。初期状態では、採用された説明変数の数NPは0であり、PRESS=∞(CV0=−∞)である。
先に述べたように、変数はp20からp2まで逆の順で処理する。表2は、表1のサンプルについて、左端の数字は、変数の取捨選択で改善がみられた10の段階を示す。なお、0は初期状態を意味する。次の列の「追加」と「除外」は、追加のループと除外のループの処理であることを意味する。次の列の変数は、追加または除外された変数を示す。次の列は、交差検証成績(PRESSをサンプル数で割ったもの)を示す。右端の列は、その段階で選択されている変数を示す。
初期状態では、変数は全くない状態であり、PRESSは∞である。表2に示すように、最初、p20を説明変数として採用すると、PRESS=0.111となり、初期値に比べて改善されるので、説明変数p20の追加を実施する。次に、変数p19を加えてp19とp20の2つを説明変数とすると、PRESS=0.129となり改善をもたらさないので、p19は追加しない。次に、説明変数p18を加えるとPRESS=0.090となり、改善するので、p18を追加し、p18とp20を説明変数とする。以下同様に表2に示すように続く。(ここで、p10を追加採用するのは、小数点以下4桁目で改善されているためである。)説明変数p20〜p2の1回目のループを終了した時点で、説明変数がp3、p6、p10、p16、p18およびp20となり、PRESS=0.60となる。2回目のループでは、説明変数p12が追加され、PRESS=0.55となる。3回目のループでは追加による改善がなく、ひとまずS114〜S126の追加処理を終了し、S214に移る。この時点での部分最小自乗法のフィットならびにリーブ・ワン・アウト予測状況は表3のとおりである。
表3は、10のサンプルについて、表2の7で示す段階まで処理が進んだ時点での部分最小自乗法のフィットならびにリーブ・ワン・アウト予測状況を示す。ここで、モデル予測とリーブ・ワン・アウト予測のそれぞれにおいて、計算値と実測値との誤差を示す。さらに、その下側に、誤差の自乗平均、相関係数Rの自乗および予測相関係数Qの自乗を示す。
次に、ステップS214から始まる除外処理の1回目のループにおいて、説明変数p10とp20を除外することが改善をもたらした。2回目のループでは改善がなく、ステップS214〜S226を終了するが、ステップS227の判断により再度S112に戻る。次に、追加処理の1回目のループにおいて、P13の追加だけが改善をもたらしたが、続く除外処理の1回目のループでは、改善がなかった。もう一度ステップS112に戻り、ステップS114〜S126およびステップS214〜S226では改善がなくなったのを確認して、処理を終了した。こうして選択された説明変数は、p3、p6、p12、p13、p16およびp18の5個であり、PRESS=0.048となった。詳細は表4のとおりである。
表4は、表2の段階10まで処理が進んだ時点での部分最小自乗法のフィットならびにリーブ・ワン・アウト予測状況を示す。
なお、説明変数の数が多い時に強いとされる部分最小自乗法であるが、p20〜p2の全てを説明変数として採用した場合には、表5に示すように、PRESS=0.124となった。すなわち、リーブ・ワン・アウト処理は、平均値からの誤差(0.093)よりも悪い成績をもたらす。
実施例.
次に、実施例を挙げて本発明をさらに詳細に説明するが、本発明はこれらの例によって何ら限定されるものではない。
実施例1: 部分最小自乗法の交差検証成績を考慮した特徴抽出によるDLBCL患者のデータ解析.
P.O.Brownらのホームページ(http://llmpp.nih.gov/lymphoma/)より入手した28名のDLBCL(リンパ腫)患者のデータを、20名のデータからなる訓練集合と8名のデータからなるテスト集合に分けた。目的変数に生存月数を採用し、説明変数には18432スポットのうち、28データにおいてch1、ch2ともに正の数となる12832スポットのlog(ch1/ch2)値を採用した。
訓練集合において部分最小自乗法(PLS)のモデル決定を試みた。12832変数全てを用いて部分最小自乗法の解析をしたところ、リーブ・ワン・アウト予測は有意(Q2>0.5)にはならなかった。次にリーブ・ワン・アウト予測誤差が最小になるように説明変数を段階的に1つづつ増減した。モデル構成手法としては前述の第3のモデル構成手法において説明変数の追加及び除外の順番並びにリーブ・ワン・アウト処理におけるサンプルの除外の順番が異なるほかは同様な方法を用いた。すなわち、どの説明変数も選択されていない状態を初期状態とする(S112)。次に、最後(m番目)の説明変数から最初(1番目)の説明変数までの未だ選択されていない説明変数ごとに逐次、その説明変数を選択した場合に交差検証成績が改善するかどうかを、リーブ・ワン・アウト処理(ここでは、最後(n番目)のサンプルから最初(1番目)のサンプルを逐次除外した)を用いた交差検証成績評価ステップを繰り返しながら判定し、改善する場合にはその説明変数を追加選択し、そのような改善と追加がなくなるまでm番目の説明変数から上記逐次判定操作を繰り返す(S114〜S126)。次に、最後(m番目)の説明変数から最初(1番目)の説明変数までの選択されている説明変数ごとに逐次、その説明変数を除外した場合に交差検証成績が改善するかどうかを、リーブ・ワン・アウト処理{ここでも最後(n番目)のサンプルから逐次除外した}を用いた交差検証成績評価ステップを繰り返しながら判定し、改善する場合にはその説明変数を除外し、そのような改善と除外がなくなるまで、最後(m番目)の説明変数から上記逐次判定操作を繰り返す(S214〜S226)。その結果、有意なモデル(R2=0.988、Q2=0.895、NP=342)を得た。図9は、このデータについての最小自乗法成績を示す。図9において、ひし形(fit)は訓練集合のデータ(20人)を示し、三角(cv)は、それらについての交差検証成績のデータを示す。また、四角(test)はテスト集合のデータ(8人)を示す。得られた部分最小自乗法モデルは、テスト集合のうち、4/8をきわめて良好に、また1/8を良好に予測するものであった。
なお、上述の多変量解析によるデータ解析では、扱ったサンプルはDNAチップを用いて得たデータであった。しかし、このデータ解析は、DNAチップを用いて得たデータに限定されるものではなく、蛋白質発現量、細胞内物質の量などのデータに対しても有用であろうことは容易に推測されることである。
以下の実施例2〜7では、部分最小自乗法を用いて選抜した少ない個数の説明変数について、通常の統計的手法または多変量解析手法(比例ハザード法、重回帰分析、適応最小自乗法、ロジスティック回帰分析法、線型判別分析法など)を適用する。
実施例2: 部分最小自乗法の交差検証成績を考慮した特徴抽出と比例ハザード解析による240名のDLBCL患者の生存時間解析.
RosenwaldらがWeb上(http://llmpp.nih.gov/DLBCL/)で公開している240名のDLBCL(びまん性大細胞型Bリンパ腫)のデータセットをダウンロードして用いた。全データを訓練集合として利用した。スポットパターンでχ1またはχ2が0となるものを除いた7399スポットについてlog(χ1/χ2)を計算して説明変数とした。本実施例では実施例1と異なり、生存時間として観測打切り時間と死亡時間とが混在していることを考慮してカプラン・マイヤー(Kaplan−Meier)法による生命表を適用して事象発生時点での生存確率(PKM)を求め、ロジット変換(log(PKM/1−PKM))した値を目的変数とした。カプラン・マイヤー法による生命表は集団としての生存確率を示すが、ここでは、個人jを含む集団としての事象発生時点での残存確率(変化の発生しなかったものが残存する確率)を個人jの事象発生時点での残存時間に読み代えるという新規な考え方を用いている。また、この確率をロジット変換して、変化の発生傾向を表現するロジット値に変換して、目的変数とした。訓練集合内の交差検証はリーブ・ワン・アウト法によって行ない、PRESS×1.02NPが小さくなるようにパラメータを逐次取捨選択して部分最小自乗法モデルを得た。ここで、交差検証成績(CV=−PRESS)の代わりに、少なくとも交差検証成績を独立変数として持つ関数の1つである関数−PRESS×1.02NPを改善して部分最小自乗法モデルを得た。ここでPRESSはリーブ・ワン・アウト予測の残差自乗和であり、NPは、選択された説明変数の数である。
図7のフロー中の交差検証成績CVを−PRESS×1.02NPと読み換えて、処理を実行することにより、下記の19個の遺伝子の発現が説明変数として選抜された。ここでdata IDはWebデータ元でのID番号を示す。またACCESSIONはGenBankのアクセション番号であり、アクセション番号の無い行はデータ元でのみ明らかとなっている遺伝子(Unknown)ないしESTであり、論文記載の方法によって入手することができる。
これらの遺伝子の発現を説明変数の候補として比例ハザード(hazard)解析を試みた。比例ハザード法とは、生存率の解析に時間を考慮した統計的手法である。解析の実行はプログラムパッケージJMP(JMP Sales SAS Campus Drive Cary,NC 27513 USA)を用いて行なった。変数削除基準としてP≧0.05を採用した変数減少法によって更に絞り込んだ結果、14遺伝子の発現からなる以下の比例ハザード式が得られた。ここでGenbank(ジーンバンク)のアクセション番号ないしdata IDで示される各項は、各遺伝子のlog(χ1/x2)値であり、またPは統計的な有意性が成り立たない危険率である。この式の右辺から求められるハザード値(hazard)が大きいほど、死亡傾向が大きい。
Rosenwaldらは、単相関の比例ハザード解析を行なって、5群(17遺伝子)の診断指標を選抜している。図10に、本実施例で得られたハザード値(Hazard、図中Hazard(pls(14))と示した)とRosenwaldらの診断指標がどの程度、生存時間を説明できているかを比較した。Rosenwaldらの5群のパラメータを同時に用いた比例ハザード式ではProlifirationパラメータがP>0.05で統計的に有意でないなどの問題を有していため、これを除く4群のパラメータを同時に含めたハザード値も比較のために掲載した(図中Hazard(Rosenwald/4para)と示した)。ここで、菱形は死亡した人または打ち切った人のデータを示し、四角は生存している人のデータを示す。
これらの診断指標のうち、本実施例で求めたハザード値と生存時間との相関は際立って明白である。即ちハザード値は生存時間につれて減衰しており、大きなハザード値の患者は長く生きることが出来ないことが示されている。一方、Rosenwaldらの指標はいずれも生存時間を診断するには不十分なものである。数百、数千という数のパラメータの中から効率的に最適のパラメータセットを見出すことは比例ハザード解析だけではできないことである。しかし以上のようにカプラン・マイヤー法、ロジット変換、部分最小自乗法の交差検証成績を考慮した特徴抽出、比例ハザード解析を組み合わせることで、従来に無い、有効な診断指標を得ることができた。統計学的に異質なモデルをこのように組み合わせることによってこのような良好な結果が得られたことは意外でもあり、興味深いことであった。患者の生存時間を予測することは、QOLを含めた治療計画や人生設計などを判断する上で重要な情報をもたらすものであり、本実施例で求められた診断モデルは社会的に価値のあるものである。
また、変数削除基準としてP≧0.001を採用した変数減少法によって更に絞り込むと、6遺伝子の発現からなる以下の比例ハザード式が得られた。このように、変数削除基準を変えることにより、選択される説明変数の数を制御できる。
図11は、右辺を計算して求められるハザード値を縦軸とし、生存時間を横軸としたプロットを示す。図10と同様に、図11において、菱形は死亡した人または打ち切った人のデータを示し、四角は生存している人のデータを示す。
実施例3: 部分最小自乗法の交差検証成績を考慮した特徴抽出と比例ハザード解析による40名の乳癌患者の生存時間解析.
SorleらがWeb上(http://genome−www.stanford.edu/breast_cacer/mopo_clinical/)で公開している乳癌患者のデータセットをダウンロードして用いた。全データを訓練集合として利用した。データセットの大部分は、タイプA,Bという2種類のDNAチップで測定されたそれぞれ40名、24名の患者よりなるが、ここではタイプAのデータを用いた。生存時間データより実施例2と同様にロジット値を求め、目的変数とした。説明変数としては、データに欠測のある遺伝子を除いた6891件のLOG_RAT2N_MEAN値を採用した。そして、少なくとも交差検証成績を独立変数として持つ関数の1つである、交差検証成績と説明変数NPの関数PRESS×1.13NPが小さくなるようにパラメータを逐次取捨選択して部分最小自乗法モデルを得た。図7のフロー中の交差検証成績CVを−PRESS×1.13NPと読み換えて、処理を実行することにより、下記の10個の遺伝子の発現が説明変数として選抜された。
これらを説明変数の候補として、比例ハザード解析において変数削除基準としてP≧0.05を採用した変数減少法を試み、7遺伝子の発現からなる以下の比例ハザード式が得られた。ここでアクセッション番号で示される各項はそれぞれの遺伝子のLOG_RAT2N_MEANである。
図12に、右辺を計算して求められるハザード値を縦軸とし、生存時間を横軸としたプロットを示す。ここでもハザード値が優れた診断指標となることが示されている。図12において、菱形は死亡した人または打ち切った人のデータを示し、四角は生存している人のデータを示す。
変数削除基準としてP≧0.001を採用した変数減少法によって更に絞り込んだ。これにより、3遺伝子の発現からなる以下の比例ハザード式が得られた。このように、変数削除基準を変えることにより、説明変数の数を制御できた。
図13は、右辺を計算して求められるハザード値を縦軸とし、生存時間を横軸としたプロットを示す。ここで、菱形は死亡した人のデータを示し、四角は生存している人のデータを示す。
実施例4: 部分最小自乗法の交差検証成績を考慮した特徴抽出と重回帰分析による40名の乳癌患者の再発予測解析.
SorleらのDNAチップAで6891遺伝子の発現が測定された40名の患者をデータセットとして用いた。再発の有無を目的変数として、PRESS×1.10NPが小さくなるようにパラメータを逐次取捨選択して11遺伝子の発現からなる部分最小自乗法モデルを得た。
次に、選抜された遺伝子発現を説明変数とし、再発の有無を目的変数として、通常の多変数解析法の一つである重回帰分析によって判別分析を実行した。解析の実行はプログラムパッケージJMPを用いて行なった。変数削除基準としてP≧0.15を採用した変数減少法によってさらに絞り込んだ結果、10遺伝子の発現からなる以下の重回帰式が得られた。この式で計算されるOLS値が正の時は再発の可能性が高く、負の時は低い。
上式に含まれる各パラメータをそれぞれ1つ用いて判別分析式を作成した場合のP値及び決定係数を以下の表6に示す。
単独では有意とはならない(P>0.05)パラメータが3つ存在し、また、どのパラメータも決定係数が小さい。従って、パラメータを1つずつ吟味するだけでは、上式のような良好な判別式は得られなかった。また数百、数千という数のパラメータの中から効率的に最適のパラメータセットを見出すことは重回帰分析だけではできないことである。しかし、以上のように、部分最小自乗法の交差検証成績を考慮して特徴抽出することにより、従来に無い、有効な診断指標を得ることができた。乳癌の再発可能性を予測することは、QOLを考慮した治療計画を立案し判断するうえで、社会的に求められているところのものである。
実施例5: 部分最小自乗法の交差検証成績を考慮した特徴抽出と適応最小自乗法による40+24名の乳癌患者の再発予測解析.
DNAチップのタイプA(40名)とタイプB(24名)に共通する3448遺伝子に限って解析を試みた。PRESS×1.17NPが小さくなるようにパラメータを逐次取捨選択して部分最小自乗法モデルを得た。選抜された遺伝子発現を説明変数とし、適応最小自乗法によって判別分析を実行した結果、次式が得られた。次式で計算されるALS値が0.5より大きいと再発の危険性が存在する。
下記の表7にみるように、H11482は単相関では有意ではなく、他の変数と同時に用いることで初めて把握できたパラメータである。また、表8は、上式を用いてタイプBの患者を予測した結果である。本判別式の感度=81.8%、特異度=53.8%となり、χ2=3.233(5%<P<10%)、予測判別正解率=66.7%、という統計的に有意な結果を得た。タイプA、BはDNAチップの構成の相違に基づく測定誤差が存在すると思われるデータであるにもかかわらず、タイプAで訓練したモデルでタイプBの予測に危険率10%以下で成功したことは勇気付けられる結果である。
また、PRESS×1.12NPが小さくなるように選んだ場合には、以下の遺伝子の発現を説明変数とする部分最小自乗法モデルを得た。
これらを説明変数の候補として、リーブ・ワン・アウトを指標にして、さらに絞り込んだ結果、次の判別式を得た。
パラメータを1つずつ吟味するだけでは、上式のような良好な判別式は得られなかった。また数百、数千という数のパラメータの中から効率的に最適のパラメータセットを見出すことは、適応最小自乗法、ロジスティック回帰分析、その他の判別分析手法だけではできないことである。しかし、以上のように、部分最小自乗法の交差検証成績を考慮して特徴抽出することにより、従来に無い、有効な診断指標を得ることができた。
実施例6: 部分最小自乗法の交差検証成績を考慮した特徴抽出とロジスティック回帰分析法または線型判別分析法による40+24名の乳癌患者の再発予測解析.
実施例5での1つめの適応最小自乗法による解析をロジスティック回帰分析法に置き換えた場合、次の判別式が得られた。
右辺で求められるLORA値が正の場合には再発の危険性が存在する。係数の比率や相関係数は実施例5の適応最小自乗法の場合と異なるものの、各患者の識別結果は全く同一であった。またタイプBの患者を予測した結果も表7と同じになった。
次に、実施例5での適応最小自乗法による解析を線型判別分析に置き換えて解析して、次の判別式が得られた。
右辺で求められるLDA値が正の場合には再発の危険性が存在する。係数の比率や相関係数は、実施例5の適応最小自乗法の場合と異なり、各患者の識別結果も若干異なったが、概ね同一であった。また、タイプBの患者を予測した結果も表7と同じになった。
以上の実施例4,5,6では、乳癌の再発の有無を目的変数としている。したがって、部分最小自乗法の交差検証成績を考慮して特徴抽出する方法が、目的変数が名義尺度や順序尺度などのデータである場合にも有効であることが示された。なお、名義尺度とは、対象(サンプル)をある分類に属するかどうかを測り分けるときの分類で、分類の間に大小や順序はない。また、順序尺度とは、対象の特定の分類について測り分けるときの分類であり、分類の間に大小、高低といった順序がある。
実施例7: 部分最小自乗法の交差検証成績を考慮した特徴抽出と比例ハザード解析による40名の乳癌患者の再発時間解析.
実施例4と同じデータを用いて、再発の時系列データを基に実施例2と同様の方法で求めたロジット値を目的変数として、PRESS×1.15NPが小さくなるようにパラメータを逐次取捨選択して9遺伝子の発現からなる部分最小自乗法モデルを得た。これらの遺伝子発現の測定値を説明変数として比例ハザード解析において変数削除基準としてP≧0.05を採用した変数減少法を試み、8遺伝子からなる、以下の比例ハザード式が得られた。
図14は、右辺を計算して求められるハザード値を縦軸とし、再発時間を横軸としたプロットを示す。ここで、菱形は再発しない人のデータを示し、四角は再発している人のデータを示す。ここでもハザード値が優れた診断指標となっており、生存時間に限らず、時間とともに確率的に発生する生体の状態の変化を解析する手法として、本発明の手法が有効であることが示されている。
変数削除基準としてP≧0.005を採用した変数減少法によって更に絞り込んだ場合には、4遺伝子の発現からなる以下の比例ハザード式が得られた。
図15は、右辺を計算して求められるハザード値を縦軸とし、再発時間を横軸としたプロットを示す。ここで、菱形は再発しない人のデータを示し、四角は再発している人のデータを示す。
実施例8: Genbankアクセッション番号H11482、T64312、AA045340を含む乳癌再発性診断用DNAチップの作成と測定.
実験医学別冊「ゲノム機能研究プロトコール」(ISBN4−89706−932−7 C3047)p34−38記載の関直彦、永杉友美、東孝典、吉川勉、鈴木収、村松正明らの方法に準じてDNAチップの作成と測定を行なった。Genbankアクセッション番号H11482、T64312、AA045340のcDNAを用いた。
プローブ用の各PCR産物をエタノール(和光純薬,Cat#057−00456)で沈殿させ、2μg/μlとなるようにDDWで調整する。ニトロセルロース(GibcoBRL Cat#41051−012)4mg/mlのDMSO溶液を等量加え、よく混和させて100℃で5分間熱変性を行ない、氷上で急冷する。次いで室温に戻し、DNAスポッターSPBIO2000(日立ソフトエンジニアリング)を用いてカルボジイミドスライドガラス(日清紡)へのスポッティングを速やかに行なう。スポットの乾燥を確認し、Ultraviolet crosslinker(アマシャムファルマシアバイオテック社)を用いて60mJ/cm2で紫外クロスリンク処理を行ない、ガラスラックに立てて室温保存する。
3%BSA、0.2M NaCl、0.1M Tris(PH7.5)、0.05% Triton X−100よりなるブロッキング液に上記マイクロアレイを浸け、約30分間放置する。次いで、ガラスに付着している溶液をよく切り、37℃で乾燥させる。TEバッファー(PH8.0,ニッポンジーンCat #316−90025)で3回軽く洗い、プレートホルダーに入れて軽く遠心(1000rpm,1分間)して余分な水分を除去する。
次に、乳腺正常株SV−40及び乳癌細胞株MCF−7、MDA−MB−468又はT−47−Dの各細胞液より、TRIZOL(GibcoBRL,Cat#15596−018)、Oligotex dT30〈Super〉(TaKaRa,Cat#W9021A)を用いてマニュアルに従って、mRNAを精製する。2μgのmRNAをDEPC処理した6.4μlのDDWに溶かし、Oligo dTプライマー9μl、5×SuperScript IIバッファー(GibcoBRL,Cat#18089−011)6μl、DTT(SuperScriptの付属)3μl、50×dNTP 0.6μl、Cy3−dUTP(アマシャムファルマシアバイオテク Cat# PA53022)又はCy5−dUTP(アマシャムファルマシアバイオテク Cat# PA55022)3μl、SuperScript II 2μlよりなる溶液を加え、42℃で2時間反応させる。途中1時間経過時点で、SuperScript IIを1μlを追加する。1.5μlアルカリバッファー(1N NaOH/20nM EDTA)を加え、65℃で10分間反応させ、TEバッファーを270μl、1N HClを1.5μl加えて、Cy3,Cy5ラベルの反応液を2つまとめて1本のMicrocon−YM−30(Millipore/Amicon,Cat#42410)に移す。10,000rpmで上のカップに残る液量が約10μlになるまで遠心を続け、カップを通りぬける液を別のチューブに移し替え、その後、上のカップにTEバッファー500μl、Human Cot−1 DNA(GibcoBRL Cat#15279−011)20μgを加え、再び液量が10μl以下になるまで遠心を続ける。3,000rpmで3分間遠心し、蛍光標識したDNAを回収する。DDWとyeast RNA(Sigma,Cat#R7125)50μg、poly(A)(ロッシュダイアグノスティクス,Cat#108 626)50μgを加えて20μlにし、PCR用のチューブに移し換え、さらに4.25μl 20×SSC(GibcoBRL,Cat#15553−035)と0.75μl 10% SDS(GibcoBRL,Cat#15553−035)を加え、PCR用の機器で100℃、1分間熱変性させ、次いで、室温で30分間放置して、ゆっくり冷却する。
蛍光標識したDNAの全量をカバーガラスにのせ、泡が入らないように注意しながら前記マイクロアレイにかぶせ、水で濡らしたキムタオルを底に敷いたハイブリダイゼーションチェンバーに入れて密閉する。毎分2〜4サイクルで軽く振とうさせながら、65℃で一晩ハイブリダイズする。ハイブリダイゼーションチェンバーからマイクロアレイを取り出し、カバーガラスが載ったままの状態で静かに2×SSC/0.1% SDS溶液中に入れ、5分間シェイキングし、カバーガラスが自然にはがれるのを待つ。カバーガラスがはがれたところでマイクロアレイをスライドガラスラックに入れ、もう一度2×SSC/0.1% SDS溶液中で5分間軽く振とうして洗う。さらに0.2×SSC/0.1% SDS40℃で5分間2回洗い、0.2×SSCでリンスする。マイクロアレイを別の乾いたプレパラートケースに移し、マイクロタイタープレート用の遠心機で軽く遠心して(1000rpm,1分室温)マイクロアレイ上の水分を除く。そして、ScanArray4000(GSI luminonics社)でシグナルを読み込み、解析ソフトにはQuant Array(GSI luminonics社)およびChip Space(日立ソフトウェアエンジニアリング)を用いる。
実施例9: 遺伝的アルゴリズムによる部分最小自乗法モデルの最適化.
実施例4で用いたSorleらのDNAチップAで6891遺伝子の発現が測定された40名の患者をデータセットとして用いた。遺伝的アルゴリズムは、たとえば、伊庭斉志;「遺伝的アルゴリズムの基礎」(オーム社(1994))に説明されている。前記データを用い、遺伝的アルゴリズムによる説明変数選択を行なった。以下において「」で区切られた用語は遺伝的アルゴリズムで通常用いられる専門用語であり、特に必要な場合には解説を加えている。「適合度」(fitness)には−PRESS×1.01Npを採用した。各「個体」の「遺伝型」は説明変数を採用する場合には1、採用しない場合には0をとる数列{b1,b2,b3,...}とした。
個体集合のサイズを100個とし、初期の個体の「遺伝型」(GTYPE)は、平均でmin_of(Ns,Ng,300)/2個の説明変数が採用となるように乱数を用いて準備した。ここでNsはサンプル数(患者数)、Ngは説明変数の候補の数、300は実装の都合上設定された定数である。
集合よりランダムに2つの個体を選抜し、「遺伝型」の「一様交叉」を行なったものの一方を新しい「個体」とした。即ち、「各遺伝子座」ごとに1/2の確率でいずれかの「親個体」の数列値(0または1)を選びそれを代入したものを新しい「個体」とした。続いて新しい「個体」の「各遺伝子座」毎に、1の場合(説明変数が採用されている場合)には1.1/採用された説明変数の数の確率で、0の場合(採用されていない場合)には1.1/採用されていない説明変数候補の数の確率で、0←→1を反転させた。
上述の「交叉・突然変異オペレーション」によって準備された新しい「個体」の「適合度」と、ランダムに選抜された「トーナメント相手」となる集合中の「個体」の「適合度」とを比較し、新しい「個体」の適合度が勝った場合には0.75の確率で、劣った場合には0.25の確率で「個体」の置き換えを行なった。ただし、「トーナメント相手」が集合中の最適解のものである場合には置き換えを禁止するという「エリート戦略」を採用した。
以上の「交叉」→「突然変異」→「選抜」サイクルを繰り返して最適化を行なった。ここではサイクル数を集合サイズで割ったものを「世代数」とする。最大「世代数」の初期値を100とし、新しい最適解が見出されるたびに最大「世代数」を10増加させながら、実行「世代数」が最大「世代数」に至るまでサイクルを繰り返した。
以上の初期集合の準備〜最適化の繰り返しおよび終了にいたる一連の処理を一回のラン(run)とし、15回のランを行なった。図16は、15回のランにおける最適化の様子をまとめている。最良の結果は25個の説明変数を用いたものである。
実施例10: 階層型人工ニューラルネットワーク(MLP)によるモデル構築.
実施例5の乳癌患者の再発性判別解析において、DNAチップtype A(40名)とtype B(24名)に共通する3448遺伝子より、PRESS×1.17Npが小さくなるようにしてPLS−CVで特徴抽出された3つの説明変数を用いた。
解析方法について説明すると、MLPは3層とし、中間層(tk)において一度だけシグモイド変換を行なう構造とし、図17の4つのトポロジーを試みた。ネットワークの重みの学習はBack propagation(逆伝播)アルゴリズムによって行なった。中間層(tk)において一度だけシグモイド変換を行なう3層MLPを用いた。
ネットワークトポロジーIおよびトポロジーIIbの結果は以下のとおりであった。なお、トポロジーIIa及びトポロジーIIcは、トポロジーIIbに劣るものであった。
実施例11: 潜在変数を用いた比例ハザードモデルの構築.
実施例3のPLS−CV法で選抜された10遺伝子の発現量を説明変数とし、目的変数として生存確率のlogit値を用いてPLSの解析過程で作成される潜在変数を1個抽出した。その抽出した潜在変数を説明変数にして比例ハザードモデルによる解析を試みた結果、作成された式はP≦0.0001で有意となった。図18に右辺を計算して得られるハザード値を縦軸とし、生存時間を横軸にしたプロットを示す。
本技術で得られたハザード式の予測の性能を評価するために、用いた40例の中から1例を除外し、残りの39例のデータを用いてハザード式を作成し、除外した1例のハザード値を予測した。39例からのハザード式によって予測した値と40例からのハザード式からの計算値をプロットした図19より、本技術はハザード値の予測において良好な成績を示した。
発明の効果について以下に説明すると、生体の状態と複数の遺伝子発現の量および/または細胞内物質の量との相関モデルを決定するとき、説明変数の選択と交差検証法とを用いて変数を絞り込むことができる。これにより、良好でかつ予測力のある多変量解析モデル(相関モデル)が得られる。特に遺伝子発現の量のように、説明変数の数がたとえば1000以上と膨大な場合に有用である。変数の数を少なくすることにより、病気や生体現象の背後で働いている重要な遺伝子やメカニズムを推定/特定でき、理解が深まる。また、重要な遺伝子産物や細胞内物質だけに絞った廉価な診断用材料(DNAチップ、DNA含有ベクター、抗体チップなど)を設計し、提供できる。
また、時間とともに確率的に発生する生体の状態の変化から導出された量を目的変数として用いて、時間とともに確率的に発生する生体の状態の変化と複数の遺伝子発現の量および/または細胞内物質の量との相関モデルを決定できる。
また、部分最小自乗法を用いて説明変数の個数を少なくすると、通常の統計的手法または多変量解析手法が適用可能になる。
【図面の簡単な説明】
図1は、遺伝子発現解析システムのブロック図である。
図2は、解析ソフトのフローチャートである。
図3は、交差検証成績CVの計算のフローチャートである。
図4は、変数選択の第1モデル構築手法のフローチャートである。
図5は、変数選択の第2モデル構築手法のフローチャートである。
図6は、変数選択の第3モデル構築手法のフローチャートである。
図7は、変数選択の第4モデル構築手法のフローチャートである。
図8は、変数選択の第5モデル構築手法のフローチャートである。
図9は、最小自乗法モデルの成績を示すグラフである。
図10は、DLBCL患者の生存時間と診断指標のプロット各種比較の図である。
図11は、実施例2のDLBCL患者の生存時間診断指標のプロットの図である。
図12は、実施例3の乳癌患者の生存時間診断指標のプロットの図である。
図13は、実施例3の乳癌患者の変数削除基準としてP≧0.0005を採用したときの生存時間診断指標のプロットの図である。
図14は、実施例7の乳癌患者の再発時間診断指標のプロットの図である。
図15は、実施例7の乳癌患者の変数削除基準としてP≧0.025を採用したときの再発時間診断指標のプロットの図である。
図16は、実施例9の遺伝的アルゴリズムによる部分最小自乗法モデルの最適化の様子を示す図である。
図17は、実施例10の階層型人工ニューラルネットワークにおける4つのトポロジーを示す図である。
図18は、実施例11の潜在変数を用いた比例ハザードモデルの乳癌患者の生存時間診断指標のグラフである。
図19は、実施例11の潜在変数を用いた比例ハザードモデルの乳癌患者の生存時間診断指標の予測値と計算値のグラフである。
本発明は、生体の状態と遺伝子発現の量および/または細胞内物質の量との多変量解析処理並びそれを基に可能となる測定機材、検定方法などに関するものである。
背景技術
2000年6月のヒトゲノムの解読宣言以降、ゲノムに書かれた遺伝情報がどのように発現して機能しているかのを解明するポストゲノム時代に突入したと言われている。ヒトゲノム計画の進展の中で、ゲノム発現状態を測定する方法論も進展してきた。トランスクリプトーム(mRNA)測定手段としてオリゴヌクレオチドアレイやマイクロチップが知られている。またプロテオーム(蛋白質)測定手段として、以前からある2次元電気泳動に加えて、最近では質量分析の方法が進歩してきた。また抗体チップなどの先進の技術も注目されている。これらの測定技術は、生体の状態パラメータを短時間に一挙に測定できることがそれまでの技術と比較して画期的であるといえる。
遺伝子発現状態を効率的に測定する技術として次のものがあげられる。トランスクリプトーム(mRNAの総体)を特定するものとして、基盤に複数種のDNAを担持し、それに相補的なmRNAを検出するDNAチップが知られている。代表的なDNAチップには、遺伝子チップやDNAマイクロアレイがある。また、プロテオーム(蛋白質の総体)を特定するものには、2次元電気泳動、抗体チップ、質量スペクトルを用いるものがある。またメタボローム(代謝中間体を含めた代謝産物の総体)を測定する手法も質量分析などによって試みられており、進展が見られる。
生体内の細胞の状態は遺伝子産物の発現によってよく記述されるため、従来の診断マーカーでは情報が不足している場面でも、精度のより高い診断が可能になるという期待も出てきている。たとえば、次のような研究があげられる。
P.O.Brownらは、DNAチップによってリンパ腫患者の細胞のトランスクリプトームを測定し、クラスター解析によって悪性と良性のリンパ腫(DLBCL)を別クラスターに分離した(Nature 403(3),503−11(2000))。しかし、これは因果関係(相関関係)のモデルを得る方法ではなく、どの遺伝子がどの程度重要かを判断できない。
A.Alaiyaらは、2次元電気泳動によって子宮がん患者40人の細胞のプロテオームを測定し、うち22人のデータから部分最小自乗法診断モデルを構築し、悪性度を説明した(Int.J.Cancer,86,731−36(2000);Electrophoresis,21,1210−17(2000);国際公開WO 00/70340)。その際、全変数モデルにおいて1553変数からloadingの大きな170変数に限定することによって交差検証成績がよくなり(Q2=0.84)、残り18患者の深刻度(3段階)を11/18の比率で正答した。交差検証法がモデル構築の際の指標になるという考えが表明されている。しかし、この方法では、loadingを得る際にまず全変数モデルが成立しなければならない。また、それ以外の変数選択手法が考案されていない。
J.Khanらは、DNAチップによって小児がん患者の細胞を測定し、ニューラルネットワークによって悪性度を説明した(Nature Medicine,7(6),673−79(2001))。小児がん(SRBCT)患者88人のトランスクリプトーム(6567遺伝子)を測定し、うち63人のデータから主成分分析によって10次元に圧縮し、次に、人工ニューラルネットワーク診断モデルを構築した。ここで、影響力のある上位遺伝子を交差検証法によって絞り込み、96遺伝子で最良の成績(100%)を得た。このモデルで残り25人を予測し、93〜100%の結果を得た。しかし、この方法でも、影響力を得る際にまず全変数モデルが成立しなければならない。またそれ以外の変数選択手法が考案されていない。10次元のような少ない変数の場合を扱えるが、変数の数が膨大な場合には適用できない。
また、最近になってDNAチップの解析に部分最小自乗法を用いる研究がD.M.RockeとD.V.Nguyenによって報告されるに至った(国際公開WO 02/25405;Bioinformatics 18(1),39−50(2002);Bioinformatics 18(9),1216−26(2002);Bioinformatics 18(12),1625−32(2002))。部分最小自乗法の潜在変数を線型判別分析などの多変量解析の説明変数として用いた場合に良好な結果が得られることが報告されている。これは部分最小自乗法が次元圧縮とモデルフィットを同時に行なうことのできる方法であるために可能となったものである。報告に示された実施例では部分最小自乗法がDNAチップ情報のモデル構築方法として優れたものであることが示されている。しかし報告においては重要な遺伝子発現量を選抜する手段としての最小自乗法の適用については触れられておらず、事前の前処理によって選択された説明変数を全て用いて解析が行なわれているという点において上述のA.Alaiyaらの研究と同様の課題を含んでいる。
従来の診断マーカーでは情報が不足している場面でも、遺伝子発現情報を活用することで、より精度(解像度)の高い診断が可能になるという期待も出てきている。遺伝子発現状態の測定結果は、膨大な情報量が得られることが従来にはなかった特徴であり、逆に情報量が多いために、効果的なデータ処理なくしてデータの活用はありえない。したがって、有用な知識を獲得するためには効果的な情報処理が欠かせない。前に説明したように、現状ではクラスター解析を中心とする方法が用いられているが、主成分分析などの方法も採用されている。クラスター解析や主成分分析は、教師付学習方法ではないため、病状の因果関係(相関関係)のモデルを得ることはできない。すなわち、どの遺伝子がどの程度重要かを解析結果から得ることができないのが難点である。一方、部分最小自乗法は次元圧縮とモデルフィットを同時に行なう強力な多変量解析手法であるが、変数の数が膨大になった場合にしばしば有意な結果が得られない事態に直面する。したがって、膨大な遺伝子発現情報などから有用な知識を獲得できるような効果的な情報処理が望まれている。また、そのような情報処理の結果を基にした効率的な測定機材、検定処理などが期待されている。
発明の開示
(発明が解決しようとする技術的課題)
この発明の目的は、多変量の遺伝子発現情報、細胞内物質情報の効果的な情報処理を提供することである。
また、この発明の目的は、効率的な検定処理を提供することである。
(その解決方法)
本発明に係るデータ解析装置は、生体の状態または時間とともに確率的に発生する生体の状態の変化を目的変数とし、複数の遺伝子発現の量および/または細胞内物質の量を説明変数とする相関モデルを決定するデータ解析装置であって、生体の状態或いはそれを導出するデータまたは時間とともに確率的に発生する生体の状態の変化に関するデータと、複数の遺伝子発現の量および/または細胞内物質の量からなるサンプルの集合を入力する入力手段と、(1)説明変数を選択する選択手段と、(2)部分最小自乗法を実行して交差検証成績を計算する計算手段または上記生体の状態の変化に関するデータにカプラン・マイヤー法又はカトラー・エデラー法による生命表を適用して変化の発生しなかったものの確率を計算して得られた確率を、仮定した分布に基づいた変換または仮定を前提としない変換をし、該変換結果を目的変数とする部分最小自乗法を実行して交差検証成績を計算する計算手段と、(3)上記(2)の計算手段の結果を評価し、説明変数の採用、不採用を判定する評価判定手段とを有し、(4)上記(1)の選択手段と上記(2)の計算手段と上記(3)の評価判定手段とを実行して部分最小自乗法モデルの少なくとも交差検証成績を独立変数として持つ関数を改善し続けて部分最小自乗法モデルを決定する決定手段とからなる。選択手段は、たとえば、説明変数を逐次取捨選択したり、遺伝的アルゴリズムを用いて説明変数を選択する。計算手段は、たとえば、1個のサンプルを逐次除外したり、複数のサンプルを逐次除外して部分最小自乗法を実行して交差検証成績を計算する。評価判定手段は、たとえば、計算手段の結果から、各計算において除外したサンプルの遺伝子発現から予測される生体の状態を示す目的変数値と、前記除外したサンプルの生体の状態を示す目的変数値との誤差の代表値を求め、当該誤差の代表値が小さくなった場合に、その交差検証成績が改善されたと判定し、説明変数を取捨選択しながら交差検証成績の評価判定を繰り返す。あるいは交差検証成績ではなく、少なくとも部分最小自乗法モデルの交差検証成績を独立変数として持つ関数が改善するかどうかを評価判定の基準として用いることもできる。決定手段は、たとえば、選択手段と計算手段と評価判定手段とを繰り返し実行して部分最小自乗法モデルの交差検証成績を改善し続けて部分最小自乗法モデルを決定する。また、選択手段と計算手段とを複数のコンピュータで実行させることもできる。こうして、相関モデルを構成するとき、交差検証成績を基準に最適化させることにより説明変数を取捨選択し、説明変数の次元を減らして良好なモデルを得る。
上述の、仮定した分布に基づいた変換または仮定を前提としない変換は、生体の状態の変化の確率が説明変数の多項式で解析できるようにするために行なうものである。分布を仮定した場合には、確率を対数変換後に負の数にしたものを状態の変化を観測した時間で割るという変換、確率を対数変換後に負の数にしたものをさらに対数にしたものを状態の変化を観測した時間で割るという変換、または確率を1より減じたものをプロビット変換したものを計算して状態の変化を観測した時間で割るという変換などが考えられる。一方、分布を仮定しない場合にはロジット変換といった方法が考えられる。変換の方法は分布にどのような仮定が成り立つかどうかあるいはなりたたないかどうかを判断することにより、それぞれの場合に応じて適切に選ぶことができる。少なくとも部分最小自乗法モデルの交差検証成績を独立変数として持つ関数としては、たとえば、前記誤差の代表値と選抜された説明変数の数の関数が考えられ、あるいはその他の独立変数を含むものであってよい。望ましくは、関数は誤差の代表値の単調減少関数であり、説明変数の数の単調減少関数である。計算量を増やさないためには簡単に計算できる関数が望ましい。具体的には−PRESS×alphaNPという関数が考えられる。ここでPRESSは予測残差自乗和であり、NPは採用された説明変数の数であり、alphaは1または1より大きい実数である。また、−PRESS×(NP+beta)gammaや−PRESS×(beta−NP)−gammaなる関数も考えられる。ここで、gammaは正の実数である。
説明変数の個数を少なくすると、通常の統計的手法または多変量解析手法が適用可能になる。本発明では部分最小自乗法を用いて選抜された説明変数を統計手法又は多変量解析手法の説明変数として、より良好なモデルを得る。或いは選抜された説明変数を用いた部分最小自乗法モデルの潜在変数を統計手法又は多変量解析手法の説明変数として、より良好なモデルを得る。ここで潜在変数とは、部分最小自乗法において通常用いられているものであって、目的変数(Yil)と説明変数(Xij)の背後に共通する次元数の少ない潜在変数(Tik)を抽出することが部分最小自乗法の次元圧縮であり、モデルフィットである。
Yil=Σ Qkl×Tik+Fil
Xij=Σ Pkj×Tik+Eij
(iはサンプル番号、lは目的変数番号、jは説明変数番号、kは潜在変数番号、F,Eは残差)
また、統計的手法又は多変量解析手法としては、重回帰分析法、線型判別分析法、適応最小自乗法、ロジスティック回帰分析法、比例ハザード解析法、マハラノビス距離を用いる判別分析法、kNN法、人工ニューラルネットワークなどが挙げられる。
本発明者等は、また、Q2やPRESS値などの交差検証成績に加えて、説明変数の個数を第2の独立変数として含む関数を最適化することで選抜される説明変数を任意に絞り込むことができることを新たに見出した。通常の統計的手法や多変量解析手法では、抽出される説明変数の個数NPの望ましい範囲がサンプル数との兼ね合いで決まっている場合がある。そのような場合、関数を、目的とする選抜数によって任意に変更できる。関数形をたとえば−PRESS×alphaNPとした場合、説明変数の個数を数個から数十個に絞り込むためには通常は定数alphaとして1.0〜3.0の値が望ましい。より望ましくは、alphaは1.0〜2.0の値となる。他の関数形f(PRESS,NP)であっても、実際に選択される説明変数の数MPおよびその時のPRESS値PRESS_MPの周辺で、f(PRESS_MP÷alpha,MP+1)≒f(PRESS_MP,MP)となるような関数は、変数選択という点では同様の効果を持つ場合がある。こうして、適当な関数形を用いることにより、望ましい範囲の個数NPの説明変数を選抜できる。このようにして、交差検証成績を用いて決定されたモデルに採用されている説明変数をさらに絞り込むと、統計的手法又は多変量解析手法によるモデルを構築できる。したがって、その性質が十分解明されている統計的手法又は多変量解析手法を採用して解析を加えることができる。
また、目的変数として、時間とともに確率的に発生する生体の状態の変化から導出された量を用いて、時間とともに確率的に発生する生体の状態の変化と複数の遺伝子発現の量および/または細胞内物質の量との相関モデルを決定できる。「時間とともに確率的に発生する生体の状態の変化」とはたとえば生存時間である。ここで、前述の部分最小自乗法に、カプラン・マイヤー法又はカトラー・エデラー法と、ロジット(logit)変換とを組み合わせる。部分最小自乗法での目的変数は、時間とともに確率的に発生する生体の状態の変化に関するデータにカプラン・マイヤー法又はカトラー・エデラー法による生命表を適用して変化の発生しなかったものの確率を計算し、これをロジット変換した値である。ロジット(logit)値とは、分類分けされたデータの、ある分類の割合(確率)Pを基に、次式logit=log{P/(1−P)}にて計算される値である。ロジット値を目的変数とする部分最小自乗法を実行して交差検証成績を計算する。こうして、先に説明したのと同様に、部分最小自乗法の交差検証成績を考慮した説明変数の抽出を行って、生存時間解析を行える。
説明変数の個数を少なくすると、通常の統計的手法または多変量解析手法が適用可能になる。そこで、決定されたモデルに採用されている説明変数又はその潜在変数を用い、時間とともに確率的に発生する生体の状態の変化を説明する統計的手法又は多変量解析手法によるモデルを構築する。たとえば、ロジット値を目的変数として求めた説明変数を用いて、他の統計的手法又は多変量解析手法(たとえば比例ハザード法や、パラメトリックな分布にあてはめた回帰分析法)を行なうことによって、より良好なモデルを得ることができる。比例ハザード法とは、Coxによって考案された方法であり、生存率の解析に時間を考慮し、かつ、多変量を扱える。比例ハザード法では、観測されている個々ごとにハザード値と呼ばれる生存率を左右する値があり、それを導く関数がある(モデルが仮定されている)として解析される。カプラン−マイヤー法は、集団全体または群ごとの生存率の推移を示す。また、パラメトリックな分布とは、ガウスが提案した正規分布から計算された確率分布のことであり、生存時間解析では指数分布、ワイブル分析、対数正規分布が用いられる。指数分布などへの当て嵌めで、数式中に多項式があり、前述の部分最小自乗法の交差検証成績を考慮した説明変数の抽出が適用される。
入力手段で説明変数として入力される複数の遺伝子の発現量および/または細胞内物質の量とは、必ずしも物質の絶対的な濃度の測定値に限定されるものではなく、加工計算された値、相対的な値、間接的に物質量を表す量などでもよい。たとえば、質量スペクトルで蛋白質の発現量を測定することができることを応用して、生体の状態を表わす目的変数と、質量スペクトルとを直接関係づける相関モデルを構築することができる。またAffymetrix社タイプのDNAチップ(ジーンチップ)では、単一のスポットが単一の遺伝子発現を特定するとは限らず、複数個のスポットが集まってはじめて単一の遺伝子発現を特定することもある。ここでもまた、各スポットの測定量を説明変数として、直接、生体の状態を説明する相関モデルを得ることができる。更には、タンパク質の電気泳動パターンの各ピークは単一のタンパク質に帰属できず、複数個のタンパク質の重ねあわせであることも多い。このような場合にも生体の状態を説明する説明変数として各ピーク強度を用いることができる。このことは、上述のAlaiyaらは子宮癌の診断の説明変数として電気泳動パターンのピーク強度を採用していることから明らかである。前述のようにポストシークエンス時代のトランスクリプトーム解析、プロテオーム解析、メタボローム解析という研究分野では、生体(細胞)内の物質を総体として把握することから出発することを特徴とする実験的アプローチが注目されている。ひとつひとつの物質の絶対的定量は必須事項ではなく、これらの実験方法によって定量される物質の量を直接、間接に表現する測定値やその加工計算値が、生体の状態を説明する説明変数と成り得る。また以上の物質量を表現する説明変数以外に、場合によっては問診データなどの他の説明変数を追加すると、さらに有効な解析結果が得られる場合もある。
本発明に係るデータ解析方法は、生体の状態または時間とともに確率的に発生する生体の状態の変化を目的変数とし、複数の遺伝子発現の量および/または細胞内物質の量を説明変数とする相関モデルを決定するデータ解析方法であって、生体の状態或いはそれを導出するデータまたは時間とともに確率的に発生する生体の状態の変化に関するデータと、複数の遺伝子発現の量および/または細胞内物質の量からなるサンプルの集合を入力する入力ステップと、(1)説明変数を選択する選択ステップと、(2)部分最小自乗法を実行して交差検証成績を計算する計算ステップまたは前記生体の状態の変化に関するデータにカプラン・マイヤー法又はカトラー・エデラー法による生命表を適用して変化の発生しなかったものの確率を計算して得られた確率を、仮定した分布に基づいた変換または仮定を前提としない変換をし、該変換結果を目的変数とする部分最小自乗法を実行して交差検証成績を計算する計算ステップと、(3)前記(2)の計算ステップの結果を評価し、説明変数の採用、不採用を判定する評価判定ステップとを有し、(4)前記(1)の選択ステップと前記(2)の計算ステップと前記(3)の評価判定ステップとを実行して部分最小自乗法モデルの少なくとも交差検証成績を独立変数として持つ関数を改善し続けて部分最小自乗法モデルを決定する決定ステップとからなる。
このデータ解析方法において、選択ステップは、たとえば、説明変数を逐次取捨選択したり、遺伝的アルゴリズムを用いて説明変数を選択する。計算ステップは、たとえば、1個のサンプルを逐次除外したり、複数のサンプルを逐次除外して部分最小自乗法を実行して交差検証成績を計算する。評価判定ステップは、たとえば、計算ステップの結果から、各計算において除外したサンプルの遺伝子発現から予測される生体の状態を示す目的変数値と、前記除外したサンプルの生体の状態を示す目的変数値との誤差の代表値を求め、当該誤差の代表値が小さくなった場合に、その交差検証成績が改善されたと判定し、説明変数を取捨選択しながら交差検証成績の評価判定を繰り返す。決定ステップは、たとえば、選択ステップと計算ステップと評価判定ステップとを繰り返し実行して部分最小自乗法モデルの交差検証成績を改善し続けて部分最小自乗法モデルを決定する。また、選択ステップと計算ステップとを複数のコンピュータで実行させることもできる。
本発明に係るデータ解析プログラムは、生体の状態または時間とともに確率的に発生する生体の状態の変化を目的変数とし、複数の遺伝子発現の量および/または細胞内物質の量を説明変数とする相関モデルを決定する、コンピュータにより実行されるデータ解析プログラムであって、生体の状態或いはそれを導出するデータまたは時間とともに確率的に発生する生体の状態の変化に関するデータと、複数の遺伝子発現の量および/または細胞内物質の量からなるサンプルの集合を入力する入力ステップと、(1)説明変数を選択する選択ステップと、(2)部分最小自乗法を実行して交差検証成績を計算する計算ステップまたは前記生体の状態の変化に関するデータにカプラン・マイヤー法又はカトラー・エデラー法による生命表を適用して変化の発生しなかったものの確率を計算して得られた確率を、仮定した分布に基づいた変換または仮定を前提としない変換をし、該変換結果を目的変数とする部分最小自乗法を実行して交差検証成績を計算する計算ステップと、(3)前記(2)の計算ステップの結果を評価し、説明変数の採用、不採用を判定する評価判定ステップとを有し、(4)前記(1)の選択ステップと前記(2)の計算ステップと前記(3)の評価判定ステップとを実行して部分最小自乗法モデルの少なくとも交差検証成績を独立変数として持つ関数を改善し続けて部分最小自乗法モデルを決定する決定ステップとからなる。
このデータ解析プログラムにおいて、選択ステップは、たとえば、説明変数を逐次取捨選択したり、遺伝的アルゴリズムを用いて説明変数を選択する。計算ステップは、たとえば、1個のサンプルを逐次除外したり、複数のサンプルを逐次除外して部分最小自乗法を実行して交差検証成績を計算する。評価判定ステップは、たとえば、計算ステップの結果から、各計算において除外したサンプルの遺伝子発現から予測される生体の状態を示す目的変数値と、前記除外したサンプルの生体の状態を示す目的変数値との誤差の代表値を求め、少なくとも当該誤差の代表値を独立変数として持つ関数である当該誤差の代表値の単調減少関数の値が小さくなった場合に、その交差検証成績が改善されたと判定し、説明変数を取捨選択しながら交差検証成績の評価判定を繰り返す。決定ステップは、たとえば、選択ステップと計算ステップと評価判定ステップとを繰り返し実行して少なくとも部分最小自乗法モデルの交差検証成績を独立変数として持つ関数を改善し続けて部分最小自乗法モデルを決定する。また、選択ステップと計算ステップとを複数のコンピュータで実行させることもできる。さらには、前記の説明変数の選択において、たとえば、初期状態では説明変数を全く含まないか、或いは、初期状態では全説明変数を含むこともできる。
前記のデータ解析プログラムにおいて、上記の生体の状態は、たとえば病気のタイプをあらわす測定値、病気の重篤度をあらわす測定値、病気のタイプをあらわす医療診断の結果、病気の重篤度をあらわす医療診断の結果、あるいはそれらを2次加工した数値である。例えば後の実施例で示すように、患者の生存時間を予測することは、QOL(quality of life:生活の質)を含めた治療計画や人生設計などを判断する上で重要な情報をもたらすものであり、社会的に価値のある診断モデルを提供することができる。また癌の再発可能性を予測することは、QOLを考慮した治療計画を立案し、医師または当の患者が選択の判断をするうえで、貴重な情報をもたらすものである。
また、本発明は、決定された前記相関モデル及び予測対象のサンプルについて当該モデルにおいて採用された説明変数を入力する入力手段と、入力された該説明変数に基づいて該サンプルの生体の状態を予測判定する予測判定手段からなるデータ解析装置、前記で決定された相関モデル及び予測対象のサンプルについて当該モデルにおいて採用された説明変数を入力する入力ステップと、入力された該説明変数に基づいて該サンプルの生体の状態を予測判定する予測判定ステップからなるデータ解析方法及び前記で決定された相関モデル及び予測対象のサンプルについて当該モデルにおいて採用された説明変数を入力する入力ステップと、入力された該説明変数に基づいて該サンプルの生体の状態を予測判定する予測判定ステップからなるデータ解析プログラムも包含する。
本発明に係るコンピュータにより読取可能な記録媒体は、上記のいずれかのプログラムを記録する。
本発明に係るびまん性大細胞型Bリンパ腫の重篤度検定用の細胞内物質測定機材および測定方法並びにびまん性大細胞型Bリンパ腫の重篤度検定方法は、実質的にジーンバンクアクセッション番号がU15085、M23452、X52479、U70426、H57330及びS69790からなる遺伝子群の発現を検出する。さらに、ジーンバンクアクセッション番号がU03398、M65066、AK001546、BC003536、X00437、U12979、H96306、AA830781及びAA804793からなる群から選択される少なくとも一つの遺伝子の発現を検出してもよい。
また、本発明に係る乳癌の重篤度検定用の細胞内物質測定機材および測定方法並びに乳癌の重篤度検定方法は、実質的にジーンバンクアクセッション番号がAA598572、AA703058及びAA453345からなる遺伝子産物を含む細胞内物質を検出する。さらに、ジーンバンクアクセッション番号がAA406242、H73335、W84753、N71160、AA054669、N32820及びR05667からなる群から選択される少なくとも一つの遺伝子産物を含む細胞内物質を検出してもよい。
また、本発明に係る乳癌の再発性検定用の細胞内物質測定機材および測定方法並びに乳癌の再発性検定方法は、実質的にジーンバンクアクセッション番号がW84753、H08581、AA045730及びAI250654からなる遺伝子産物を含む細胞内物質を検出する。さらに、ジーンバンクアクセッション番号がAA448641、R78516、R05934、AA629838及びH53037からなる群から選択される少なくとも一つの遺伝子産物を含む細胞内物質を検出してもよい。
また、本発明に係る乳癌の再発性検定用の細胞内物質測定機材および測定方法並びに乳癌の再発性検定方法は、実質的にジーンバンクアクセッション番号がAA434397、T83209、N53427、N29639、AA485739、AA425861、H84871、T64312、T59518及びAA037488からなる遺伝子産物を含む細胞内物質を検出する。さらに、ジーンバンクアクセッション番号がAA406231の遺伝子産物を含む細胞内物質を検出してもよい。
また、本発明に係る乳癌の再発性検定用の細胞内物質測定機材および測定方法並びに乳癌の再発性検定方法は、実質的にジーンバンクアクセッション番号がH11482、T64312及びAA045340からなる遺伝子産物を含む細胞内物質を検出する。
細胞内物質測定機材としては、DNAマイクロアレイ、ジーンチップ、オリゴDNA型のDNAチップ、電気化学DNAチップ(ECAチップ)、繊維型DNAチップ、磁性ビーズDNAチップ(PSS)、糸巻きDNAチップ(PSS)、などのDNAチップ、マクロアレイ、抗体チップ、測定用試薬キットなどが挙げられる。また、上記の機材を適宜組み込んだ測定機械であってもよい。
発明を実施するための最良の形態
以下、添付の図面を参照して本発明の実施の形態を説明する。
以下に、選択された生体の状態と遺伝子発現の量および/または細胞内物質の量との相関モデルの決定について説明する。ここで、遺伝子発現の用語は、mRNA発現(トランスクリプトーム)や、mRNAによる翻訳の結果として生じる蛋白質(プロテオーム)を含むものとして用いる。また、細胞内物質の量とはここではたとえば、代謝中間体を含めた代謝産物全部であるメタボロームを意味する。たとえば、トランスクリプトーム(mRNA)やプロテオーム(蛋白質)の解析において、各サンプルデータは、生体の状態と遺伝子発現の量などからなる。各サンプルはたとえば1000個以上の膨大な遺伝子発現の量を含む。生体の状態は、たとえば病気のタイプまたは病気の診断指標であるが、より一般的には生体情報であればよい。「病気の診断指標」には、病気の進行度合いのほか、病気のタイプ、重篤度、深刻度などの表現で表わされるものも含む。ここで、遺伝子発現の量などの測定データは膨大な情報量からなるので、コンピュータを用いた効率的な多変量解析が必要である。
データ収集において、予めいくつかのサンプルについて生体の状態(たとえば診断指標)を判定し、また、そのサンプルされたものから細胞液を獲得し、その細胞液中の多くの遺伝子産物の発現の量などを測定する。本発明の実施の形態のデータ解析では、こうして得られた遺伝子産物の発現の量などと生体の状態(たとえば診断指標)を入力し、相関モデル(たとえば部分最小自乗法モデル)を得る。ここで、コンピュータによる多変量解析プログラムを用いて、診断指標を目的変数とし、遺伝子発現の量および/または細胞内物質の量を説明変数とする因果関係型の解析を行なって、各説明変数の重要性や影響度に関する情報を得る。また、前記目的変数は、必ずしも測定値そのものである必要はなく、ロジット変換を行なった値や群を表す離散値を用いても良く、その場合、より有意な解析結果を得ることもできる。
本発明者らは、遺伝子発現による医療診断という分野において、データ解析における交差検証(cross validation)の成績を少なくとも独立変数のひとつとして持つ関数を最適化するように変数を選択することによって良好な相関モデル(たとえば部分最小自乗法モデル)が得られることを見出した。交差検証法では、手持ちのデータを複数群に分割し、その一部のデータ群(訓練集合)だけを使ってフィットしたモデルを用いて残る別のデータ群(テスト集合)を予測することによって、モデルの予測力を試す。通常の部分最小自乗法(PLS)においては潜在変数の次元選択に交差検証法が用いられているが、ここでは、部分最小自乗法において、潜在変数を1次元に固定し、1以上の入力変数(説明変数)を逐次取捨選択しながら、交差検証成績(たとえば平方和の予測誤差)を少なくとも独立変数のひとつとして持つ関数を最適化した。ただし本発明の効果は潜在変数の次元を1に限定するものではない。その結果、全変数を採用した場合には有意な相関モデルを得られなかった場合にも、良好でかつ予測力のある相関モデルが得られることが判明したのである。この交差検証法を用いた変数選択の逐次取捨選択により、安定な相関モデルが得られる。また本発明者らは、関数形を適切に設定することによって説明変数を絞り込むことにより、部分最小自乗法以外の統計学又は多変量解析の良好な相関モデルを得ることが可能となり、それぞれ生体の状態を記述する目的変数にふさわしい相関モデルを得ることができることを見出した。なお、ここでいう「最適化」とは、交差検証成績が、説明変数を取捨選択するための、そのときの解析条件の範囲で、改善がみられなくなるまで改良したことを意味しており、交差検証成績がすべての説明変数の組合せの中で最適なものを見出したという意味ではない。この変数選択手法を用いると、病状を決定する因子を少数に特定し、廉価な診断用材料(DNAチップ、抗体チップ、DNA含有ベクターなど)を設計でき、それ自体独自の価値を持つものである。また、この変数選択手法は、予め設定される各種の変数選択条件と共に運用することが可能である。
上に述べたように、説明変数は、交差検証成績を基準に逐次取捨選択される。ここで、取捨選択のため、交差検証成績を少なくとも独立変数のひとつとして持つ関数を用いる。説明変数を追加する場合は、その説明変数について、前記関数が改善されなかったと判定された場合には当該説明変数を除外し、改善されたと判定された場合には当該説明変数を追加する。また、説明変数を除外する場合は、その説明変数について、前記関数が改善されなかったと判定された場合には当該説明変数を除外せず、改善されたと判定された場合には当該説明変数を除外する。ここで、1以上の説明変数を選択した場合に、交差検証成績評価は次のように進める。n個のサンプルからいくつかのサンプルを逐次除外して部分最小自乗法モデルを求め、各モデルにおいて除外したサンプルの遺伝子発現の量から予測される生体の状態を示す目的変数と、除外したサンプルの生体の状態を示す目的変数との各々の誤差の代表値を求める。「代表値」とは、和、平均、最大値、中位値、最頻値などのデータを特徴づける値をいう。そして、当該誤差の代表値を少なくともひとつの独立変数とする関数が小さくなった場合に、交差検証成績が改善されたと判定し、当該説明変数を追加または削除する。この交差検証成績評価を、説明変数を取捨選択しながら逐次繰り返して、前記関数を改善し続ける。改善されなくなれば交差検証成績を最適化したとして説明変数の取捨選択を終了する。その結果、取捨選択により絞り込んだ数の説明変数からなる最適な部分最小自乗法モデルが得られる。具体的には、計算手段において計算される交差検証成績の数値指標として予想残差自乗和(PRESS)を採用し、評価判定手段において予想残差自乗和の値が説明変数あたり一定の閾値以下の比率で小さくなる場合に、その説明変数を採用すると判定することにより、上記の処理は実行可能である。
因果関係型の解析手法においてはオーバーフィット(over fitting)を避けるための工夫が必要となる。ここでいうオーバーフィットとは、説明変数が多すぎるためにたまたま予測結果と実績とが一致するものの、本当の相関関係をとらえ損なっているため、モデルフィットに用いたデータ以外に予測能力を持たないことをいう。ここでは、相関モデルとして部分最小自乗法を用いるが、部分最小自乗法は次元圧縮とモデルフィットを同時に行なう強力な多変量解析手法であり、オーバーフィットの問題に比較的強いとされている。しかし遺伝子発現状態解析のように膨大な変数を扱う場合には、有意な結果が得られない事態に直面する。従来技術として説明したAlaiyaやKhanの手法は全変数モデルが有意に成立することを前提としているので、変数の絞込みには一般的には適用できない。これに対し、本発明では、交差検証予測結果を最適にするように変数を絞り込むことにより、オーバーフィットを減らすことができた。また、本発明は、前記Khanの手法とは異なり、主成分分析などの前処理を介さない方法である。従来技術では、説明変数が膨大な場合には、有意なモデルを得ることができないことから、予め、全説明変数を基にたとえば、主成分分析などで次元圧縮する前処理をし、これによって得られた説明変数によって解析する方法が用いられる。しかし、この方法では、構成したモデルで予測を行なうためには、モデル構成の基となった全説明変数が必ず必要となり、たとえば、説明変数が遺伝子発現の量であれば、診断用遺伝子チップに担持する遺伝子としては、モデル構成に用いた遺伝子の全てが必要となるか、または別の手法を用いて変数選択することが必要となる。一方、本発明においては、説明変数の選択によって説明変数を絞り込んでいるので、たとえば、説明変数が遺伝子発現の量であれば、診断用遺伝子チップに担持する遺伝子は、選択された説明変数に相当する遺伝子を担持すれば良いことになる。
なお、Todeschiniらは、有機化合物の大気中の分解を予測するため、遺伝的アルゴリズムによって交差検証成績を最適化するように変数選択を行ない、重回帰モデルを得ている(P.Gramatics,V.Consonni & R.Todeschini,Chemosphere 38(5),1371−78(1999))。53化合物と175記述子でモデル構築を行ない(Q2=0.79)、7変数が選択され、98化合物の予測を行なった(Q2=0.75)。交差検証成績を最適化するように変数選択を行なっている点では、本実施形態と同様の手法である。しかし、重回帰モデルを採用しているために、説明変数の選択過程を通じて選択される変数は少数個にとどまらざるを得ず、複数の遺伝子発現の量および/または細胞内物質の量の解析には適用できない。本発明者らの調査した範囲では、Q2やPRESS値を最適化する方法では、選抜される説明変数は百程度から数百程度にわたり、重回帰モデルでは解析が不能となる。またTodeschiniらは、説明変数を絞り込むための有効な方法について言及していない。これは、もともとの説明変数の候補がたかだか175個であり、説明変数を絞り込むために特別の工夫をする必要がないからである。遺伝子発現解析の分野はこれとは全く異なり、数十から数百のサンプル数に対して、数百から数千、数万の説明変数候補が存在する。したがってこれまでとは異なる工夫が必要となる。
本実施形態では、生体の状態と複数の遺伝子発現の量および/または細胞内物質の量との相関モデルを決定するとき、交差検証成績を少なくとも独立変数のひとつとして持つ関数を最適化させるように説明変数を逐次追加・除外することによって、説明変数を選抜して、良好な相関モデルを得る。このようなアプローチの優位性は、下記の実施例から推測されるように、次のとおりである。
1)病気や生体現象の背後で働いている重要な遺伝子やメカニズムを推定/特定でき、理解が深まる。
2)重要な遺伝子産物や細胞内物質だけに絞った廉価な診断用材料(DNAチップ、抗体チップなど)の設計が可能になる。
本実施形態では、交差検証成績を少なくとも独立変数のひとつとして持つ関数を最適化するように説明変数を段階的に取捨選択するが、たとえば具体的には、ステップワイズ(step wise)法に代表される説明変数を選択する選択手段と、リーブ・ワン・アウト(leave−one−out)法に代表される交差検証法に部分最小自乗法を適用して計算する計算手段と、前記計算手段の結果を評価し、説明変数の採用、不採用を判定する評価判定手段とを組合せて用いる。すなわち、m個の説明変数の中から1以上の説明変数を選択し、次いで、部分最小自乗法を実行して交差検証成績を計算し、さらに、該計算結果を評価して、選択した説明変数の採用、不採用を判定する。この評価判定では、計算手段の結果から、各計算において除外したサンプルの遺伝子発現から予測される生体の状態を示す目的変数値と、前記除外したサンプルの生体の状態を示す目的変数値との誤差の代表値を求め、少なくとも当該誤差の代表値を独立変数として持つ関数である当該誤差の代表値の単調減少関数の値が小さくなった場合に説明変数の取捨選択を判定する。このように、選択手段と計算手段と評価判定手段とを用いて、少なくとも部分最小自乗法モデルの交差検証成績を独立変数として持つ関数を改善し続けて、その改善がみられなくなるまで改良し、部分最小自乗法モデルを決定する。なお、本実施形態では、サンプルを1個づつ逐次除外している(リーブ・ワン・アウト法)が、その代わりに、複数のサンプルを除外して交差検証成績を評価してもよい(リーブ・n・アウト法)し、また、Khan et al.により用いられた3分割法(three−fold)等の他の方法を用いることもできる。3分割法では、説明変数をランダムにシャッフルして3つのグループに分ける。その中の2つのグループを用いてモデルを構成し、残りの1つのグループでモデルを評価する。また、説明変数の選択方法としてはステップワイズ法、非線形アルゴリズム(たとえば遺伝的アルゴリズムなど)を用いてもよく、変数選択に関して予め何らかの条件が分っていれば、それに応じて探索範囲を限定できる。
次に、データの収集と解析について具体的に説明する。図1は、遺伝子発現解析システムを示す。データ収集のため、予めいくつかのサンプルについて診断指標(たとえば病気のタイプないし進行度合いを含む)を判定し、また、そのサンプルされたものから細胞液を獲得し、DNAチップを用いてその細胞液中の多くの遺伝子産物の発現の量を測定する。測定には、共焦点型レーザスキャナ(たとえばAffymetrix社、428アレイスキャナ)10を用いる。吸光度によりmRNAの量が測定される。このデータ収集は公知の方法である。測定データは、コンピュータ12に送られ解析される。コンピュータ12は、CPU14を備えた通常の構成のコンピュータであり、それに接続される記憶装置(たとえばハードディスク装置)16の記録媒体(たとえばハードディスク)には、測定データ18や解析ソフト20が格納される。この解析ソフト20を用いてデータ18が解析され、生体の状態と遺伝子発現の量などとの相関モデルが決定される。
なお、説明変数の選択と、交差検証法に部分最小自乗法を適用する計算とを複数のコンピュータで実行させてもよい。交差検証予測の計算を複数個のコンピュータに分散させることで計算を加速することができる。
図2は、コンピュータ12により実行される、生体の状態と遺伝子発現の量などとの相関モデルを得るためのデータ解析ソフト20のフローチャートを示す。ここでは簡単に説明するため、少なくとも部分最小自乗法モデルの交差検証成績を独立変数として持つ関数として−PRESSを採用しているが、発明の範囲を限定するものでなく、実施例2〜5においては別の関数を採用している。まず、相関モデル作成用のデータを入力する(S10)。データはたとえばDNAチップを用いて収集したものである。入力データ(サンプル集合)は、それぞれ目的変数(たとえば診断指標)とm個(たとえば2000個)の説明変数(たとえば遺伝子発現の量)からなる。また、場合によっては、上述のデータ(訓練集合)以外に、テスト集合のデータを入力する。ここでテスト集合とは交差検証の評価のためのデータ群を意味するのではなく、モデル決定が終了した後にモデルの予測力をテストするためのデータ群である。
まず、初期設定として、選択された説明変数の数を0とし、交差検証成績CVの最良値CV0を−∞とする(S12)。次に、説明変数の選択を行う。まず、説明変数を指す番号iを1とし(S14)、第i変数(遺伝子発現の量)を仮に採用して(S16)、部分最小自乗法を実行し、交差検証成績CVを計算する(S18、図3参照)。ここで、リーブ・ワン・アウト処理を用いる。これは、たとえば50個のサンプルからなる訓練集合において、1番から50番の全てを順次1個づつ除いて残りの49個のサンプルで予測した結果と、その時除いた1個の結果とを比較し、その誤差が大きい場合に、仮に選択した説明変数(第i変数)が適していないと判断する手法である。もし、得られた成績CVが現在の最良値CV0より最適化されれば(S20でYES)、第i変数を採用し、かつ、成績CVを新らしい最良値CV0に更新する(S22)。しかし、得られた成績CVが最良値CV0より大きくなければ(S20でNO)、第i変数を採用しない(S24)。そして、ステップS14に戻り、同様の処理を繰り返す。この処理を交差検証成績CVが改善されなくなる(S26でNO)まで繰り返す。ここで、相関モデルに採用する説明変数については1つづつ段階的に増加(追加)または減少(除外)して成績CVを評価判定している。すなわち、全体としての合致度合いがよくなるように各説明変数を解析に加えるかどうかを逐次判定しながら、説明変数の取捨選択を行い、これを、全体としての合致度合いがよくならなくなるまで繰り返す。以上の処理で改善があると、ふたたびステップS14の初め(i=1)に戻り、それまでに選択されている説明変数を基に、さらに説明変数の選択を繰り返す。なお、ここではモデルの予測力を判断するために、訓練集合とテスト集合とに予め分割しておいたデータ集合を用いてデータ解析しており、上述の解析は、訓練集合を用いて行なった結果であるので、この結果からテスト集合について予測を行い、実測データとの合致度を評価(S28)している。このような評価は必ずしも必要でないが、予測力を判断するには有効である。
図3は、リーブ・ワン・アウト処理を含む交差検証成績CVの計算(図2、S18)のフローチャートを示す。ここで、選択された変数について交差検証成績が計算される。まず、PRESSの初期値を0とする(S180)。次に、n個の集合内のサンプルを指す番号jを1とし(S182)、第jサンプル以外のn−1個のサンプルで部分最小自乗法を実行し(S184)、第jサンプルの目的変数を予測する(S186)。差の自乗を計算してPRESSに加算する(S190)。次に番号jを1増加し(S182)、同様の処理をおこなう。これを番号j=nまで各サンプルについて繰り返す。得られたPRESSは、1個のサンプルを順次除外して計算した予測値と実測値との差の平方和であり、予測誤差を表わす量である。この予測残差自乗和PRESSの符号を変えたものを交差検証成績CVとする(S192)。
本実施形態では、交差検証法を用いて、入力変数(説明変数)を段階的に1つづつ追加・除外しながら、交差検証成績(CV=−PRESS)を最適化する。ここで、説明変数の段階的な追加・除外の内容を理解しやすくするため、以下で、さらに具体的に5つのモデル構築手法について説明する。これらは、説明変数の逐次的な選択の手順が異なる。
図4は、第1のモデル構築手法を示す。データ集合においてどの説明変数も選択されていない状態を初期状態とする(S112)。次に、1番目の説明変数から最後(m番目)の説明変数までの未だ選択されていない説明変数ごとに逐次、その説明変数を選択した場合に交差検証成績が改善するかどうかを、リーブ・ワン・アウト処理を用いた交差検証成績評価ステップ(S118)を繰り返しながら判定(S120)し、改善する場合にはその説明変数を追加する(S114〜S124)。そのような改善と追加がなくなる(S126でNO)まで、1番目の説明変数から上記逐次判定操作を繰り返す。
さらに詳しく説明すると、まず、初期設定として、選択された説明変数の数NPを0とし、交差検証成績CVの最良値CV0を−∞とする(S112)。次に、説明変数の選択を行う。まず、変数iを1とし(S114)、第i変数を仮に採用する(S116)。ただし、第i変数がすでに採用されていれば(S115でYES)、ステップS114に戻る。次に、部分最小自乗法を実行し、交差検証成績CVを計算する(S118)。ここで、リーブ・ワン・アウト処理を用いる。もし、得られた成績CVが現在の最良値CV0より最適化されれば(S120でYES)、第i変数を採用し、かつ、成績CVを新らしい最良値CV0に更新する(S122)。しかし、得られた成績CVが最良値CV0より大きくなければ(S120でNO)、第i変数を採用しない(S124)。そして、ステップS114に戻り、同様の処理を繰り返す。この処理を交差検証成績CVが改善されなくなる(S126でNO)まで繰り返す。以上の処理で改善があると、ふたたびステップS114に戻り、新しいループを開始する。ここで、それまでに選択されている変数を基に、さらに変数の選択を繰り返す。こうして、データ集合を用いて選択された変数を用いた相関モデルが得られる。
図5は、第2のモデル構築手法を示す。この手法では、全ての説明変数が選択されている状態を初期状態とする(S212)。次に、1番目の説明変数から最後(m番目)の説明変数までの選択されている説明変数ごとに逐次、その説明変数を除外した場合に交差検証成績が改善するかどうかを、リーブ・ワン・アウト処理を用いた交差検証成績評価ステップ(S218)を繰り返しながら判定(S220)し、改善する場合にはその説明変数を除外する(S214〜S224)。そのような改善と除外がなくなる(S226でNO)まで、1番目の説明変数から上記逐次判定操作を繰り返す。
さらに詳しく説明すると、まず、初期設定として、選択された説明変数の数NPをmとし、交差検証成績CVの最良値CV0を−∞とする(S212)。すなわち、すべての説明変数を選択する。次に、説明変数の選択を行う。まず、変数iを1とし(S214)、第i変数を仮に除外する(S216)。ただし、第i変数がすでに除外されていれば(S215でYES)、ステップS214に戻る。部分最小自乗法を実行し、交差検証成績CVを計算する(S218)。ここで、リーブ・ワン・アウト処理を用いる。もし、得られた成績CVが現在の最良値CV0より最適化されれば(S220でYES)、第i変数を除外し、かつ、成績CVを新らしい最良値CV0に更新する(S222)。しかし、得られた成績CVが最良値CV0より大きくなければ(S220でNO)、第i変数を除外しない(S224)。そして、ステップS214に戻り、同様の処理を繰り返す。この処理を交差検証成績CVが改善されなくなる(S226でNO)まで繰り返す。以上の処理で改善があると、ふたたびステップS214に戻り、新しいループを開始する。ここで、それまでに選択されている変数を基に、さらに変数の選択を繰り返す。こうして、データ集合を用いて選択された変数を用いた相関モデルが得られる。
図6は、第3のモデル構成手法を示す。この手法は、第1と第2の手法の直列的な組合せである。まず、どの説明変数も選択されていない状態を初期状態とする(S112)。次に、1番目の説明変数から最後(m番目)の説明変数までの未だ選択されていない説明変数ごとに逐次、その説明変数を選択した場合に交差検証成績が改善するかどうかを、リーブ・ワン・アウト処理を用いた交差検証成績評価ステップを繰り返しながら判定し、改善する場合にはその説明変数を追加選択し、そのような改善と追加がなくなるまで1番目の説明変数から上記逐次判定操作を繰り返す(S114〜S126)。次に、1番目の説明変数から最後(m番目)の説明変数までの選択されている説明変数ごとに逐次、その説明変数を除外した場合に交差検証成績が改善するかどうかを、リーブ・ワン・アウト処理を用いた交差検証成績評価ステップを繰り返しながら判定し、改善する場合にはその説明変数を除外し、そのような改善と除外がなくなるまで、1番目の説明変数から上記逐次判定操作を繰り返す(S214〜S226)。
図7は、第4のモデル構築手法を示す。この手法は、第3の手法の変形である。まず、どの説明変数も選択されていない状態を初期状態とする(S112)。次に、1番目の説明変数から最後(m番目)の説明変数までの未だ選択されていない説明変数ごとに逐次、その説明変数を選択した場合に交差検証成績が改善するかどうかを、リーブ・ワン・アウト処理を用いた交差検証成績評価ステップ(S118)を繰り返しながら判定(S120)し、改善する場合にはその説明変数を追加選択する(S114〜S124)。そのような改善と追加がなくなる(S126でNO)まで、1番目の説明変数から上記逐次判定操作を繰り返す。次に、1番目の説明変数から最後(m番目)の説明変数までの選択されている説明変数ごとに逐次、その説明変数を除外した場合に交差検証成績が改善するかどうかを、リーブ・ワン・アウト処理を用いた交差検証成績評価ステップ(S218)を繰り返しながら判定(S220)し、改善する場合にはその説明変数を除外する(S214〜224)。そのような改善と除外がなくなる(S226でNO)まで、1番目の説明変数から上記逐次判定操作を繰り返す。上記逐次判定追加改善ステップまたは上記逐次判定除外改善ステップで少なくとも一度改善があれば(S227でYES)、ステップS112に戻り、上記操作(S112〜S227)を繰り返す。これを改善がなくなる(S227でNO)までおこなう。
図8は、第5のモデル構築手法を示す。この手法は、第1と第2のスキームの並列的な組合せである。どの説明変数も選択されていない状態を初期状態とする(S112)。次に、1番目の説明変数から最後(m番目)の説明変数までの説明変数ごとに逐次、その説明変数が選択されていない場合にはその説明変数を選択した場合に交差検証成績が改善するかどうかを、リーブ・ワン・アウト処理を用いた交差検証成績評価ステップ(S118)を繰り返しながら判定(S120)し、改善する場合にはその説明変数を追加する(S114〜S124)。また、選択する説明変数ごとに、その説明変数がすでに選択されている場合には、その説明変数を除外した場合に交差検証成績が改善するかどうかを、リーブ・ワン・アウト処理を用いた交差検証成績評価ステップ(S218)を繰り返しながら判定(S220)し、改善する場合にはその説明変数を除外する(S216〜S224)。そのような改善と追加または除外がなくなる(S126でNO)まで、1番目の説明変数から上記逐次判定操作を繰り返す。
次に、第4のモデル構築手法(図7)を適用した場合を、表1のデータ集合を例として説明する。このデータ集合に対して、部分最小自乗法による解析を用いて相関モデルを求める。表1のデータでは、サンプルの数nは10であり、また、説明を容易にするため、説明変数の数mは19と少なくしている。表1において、p1は目的変数を表わし、p2〜p20は説明変数を表わす。(ただし表1では、表示の便宜のため、p16以降のデータを省略している。)第4手法(図7)のステップS114、S214とは異なり、説明変数を表わすiはp20からp2まで逆に逐次処理することとした。CV評価値としてここでは予測残差自乗和(PRESS)を採用した。PRESSが小さいほど、CV評価値はよい。初期状態では、採用された説明変数の数NPは0であり、PRESS=∞(CV0=−∞)である。
先に述べたように、変数はp20からp2まで逆の順で処理する。表2は、表1のサンプルについて、左端の数字は、変数の取捨選択で改善がみられた10の段階を示す。なお、0は初期状態を意味する。次の列の「追加」と「除外」は、追加のループと除外のループの処理であることを意味する。次の列の変数は、追加または除外された変数を示す。次の列は、交差検証成績(PRESSをサンプル数で割ったもの)を示す。右端の列は、その段階で選択されている変数を示す。
初期状態では、変数は全くない状態であり、PRESSは∞である。表2に示すように、最初、p20を説明変数として採用すると、PRESS=0.111となり、初期値に比べて改善されるので、説明変数p20の追加を実施する。次に、変数p19を加えてp19とp20の2つを説明変数とすると、PRESS=0.129となり改善をもたらさないので、p19は追加しない。次に、説明変数p18を加えるとPRESS=0.090となり、改善するので、p18を追加し、p18とp20を説明変数とする。以下同様に表2に示すように続く。(ここで、p10を追加採用するのは、小数点以下4桁目で改善されているためである。)説明変数p20〜p2の1回目のループを終了した時点で、説明変数がp3、p6、p10、p16、p18およびp20となり、PRESS=0.60となる。2回目のループでは、説明変数p12が追加され、PRESS=0.55となる。3回目のループでは追加による改善がなく、ひとまずS114〜S126の追加処理を終了し、S214に移る。この時点での部分最小自乗法のフィットならびにリーブ・ワン・アウト予測状況は表3のとおりである。
表3は、10のサンプルについて、表2の7で示す段階まで処理が進んだ時点での部分最小自乗法のフィットならびにリーブ・ワン・アウト予測状況を示す。ここで、モデル予測とリーブ・ワン・アウト予測のそれぞれにおいて、計算値と実測値との誤差を示す。さらに、その下側に、誤差の自乗平均、相関係数Rの自乗および予測相関係数Qの自乗を示す。
次に、ステップS214から始まる除外処理の1回目のループにおいて、説明変数p10とp20を除外することが改善をもたらした。2回目のループでは改善がなく、ステップS214〜S226を終了するが、ステップS227の判断により再度S112に戻る。次に、追加処理の1回目のループにおいて、P13の追加だけが改善をもたらしたが、続く除外処理の1回目のループでは、改善がなかった。もう一度ステップS112に戻り、ステップS114〜S126およびステップS214〜S226では改善がなくなったのを確認して、処理を終了した。こうして選択された説明変数は、p3、p6、p12、p13、p16およびp18の5個であり、PRESS=0.048となった。詳細は表4のとおりである。
表4は、表2の段階10まで処理が進んだ時点での部分最小自乗法のフィットならびにリーブ・ワン・アウト予測状況を示す。
なお、説明変数の数が多い時に強いとされる部分最小自乗法であるが、p20〜p2の全てを説明変数として採用した場合には、表5に示すように、PRESS=0.124となった。すなわち、リーブ・ワン・アウト処理は、平均値からの誤差(0.093)よりも悪い成績をもたらす。
実施例.
次に、実施例を挙げて本発明をさらに詳細に説明するが、本発明はこれらの例によって何ら限定されるものではない。
実施例1: 部分最小自乗法の交差検証成績を考慮した特徴抽出によるDLBCL患者のデータ解析.
P.O.Brownらのホームページ(http://llmpp.nih.gov/lymphoma/)より入手した28名のDLBCL(リンパ腫)患者のデータを、20名のデータからなる訓練集合と8名のデータからなるテスト集合に分けた。目的変数に生存月数を採用し、説明変数には18432スポットのうち、28データにおいてch1、ch2ともに正の数となる12832スポットのlog(ch1/ch2)値を採用した。
訓練集合において部分最小自乗法(PLS)のモデル決定を試みた。12832変数全てを用いて部分最小自乗法の解析をしたところ、リーブ・ワン・アウト予測は有意(Q2>0.5)にはならなかった。次にリーブ・ワン・アウト予測誤差が最小になるように説明変数を段階的に1つづつ増減した。モデル構成手法としては前述の第3のモデル構成手法において説明変数の追加及び除外の順番並びにリーブ・ワン・アウト処理におけるサンプルの除外の順番が異なるほかは同様な方法を用いた。すなわち、どの説明変数も選択されていない状態を初期状態とする(S112)。次に、最後(m番目)の説明変数から最初(1番目)の説明変数までの未だ選択されていない説明変数ごとに逐次、その説明変数を選択した場合に交差検証成績が改善するかどうかを、リーブ・ワン・アウト処理(ここでは、最後(n番目)のサンプルから最初(1番目)のサンプルを逐次除外した)を用いた交差検証成績評価ステップを繰り返しながら判定し、改善する場合にはその説明変数を追加選択し、そのような改善と追加がなくなるまでm番目の説明変数から上記逐次判定操作を繰り返す(S114〜S126)。次に、最後(m番目)の説明変数から最初(1番目)の説明変数までの選択されている説明変数ごとに逐次、その説明変数を除外した場合に交差検証成績が改善するかどうかを、リーブ・ワン・アウト処理{ここでも最後(n番目)のサンプルから逐次除外した}を用いた交差検証成績評価ステップを繰り返しながら判定し、改善する場合にはその説明変数を除外し、そのような改善と除外がなくなるまで、最後(m番目)の説明変数から上記逐次判定操作を繰り返す(S214〜S226)。その結果、有意なモデル(R2=0.988、Q2=0.895、NP=342)を得た。図9は、このデータについての最小自乗法成績を示す。図9において、ひし形(fit)は訓練集合のデータ(20人)を示し、三角(cv)は、それらについての交差検証成績のデータを示す。また、四角(test)はテスト集合のデータ(8人)を示す。得られた部分最小自乗法モデルは、テスト集合のうち、4/8をきわめて良好に、また1/8を良好に予測するものであった。
なお、上述の多変量解析によるデータ解析では、扱ったサンプルはDNAチップを用いて得たデータであった。しかし、このデータ解析は、DNAチップを用いて得たデータに限定されるものではなく、蛋白質発現量、細胞内物質の量などのデータに対しても有用であろうことは容易に推測されることである。
以下の実施例2〜7では、部分最小自乗法を用いて選抜した少ない個数の説明変数について、通常の統計的手法または多変量解析手法(比例ハザード法、重回帰分析、適応最小自乗法、ロジスティック回帰分析法、線型判別分析法など)を適用する。
実施例2: 部分最小自乗法の交差検証成績を考慮した特徴抽出と比例ハザード解析による240名のDLBCL患者の生存時間解析.
RosenwaldらがWeb上(http://llmpp.nih.gov/DLBCL/)で公開している240名のDLBCL(びまん性大細胞型Bリンパ腫)のデータセットをダウンロードして用いた。全データを訓練集合として利用した。スポットパターンでχ1またはχ2が0となるものを除いた7399スポットについてlog(χ1/χ2)を計算して説明変数とした。本実施例では実施例1と異なり、生存時間として観測打切り時間と死亡時間とが混在していることを考慮してカプラン・マイヤー(Kaplan−Meier)法による生命表を適用して事象発生時点での生存確率(PKM)を求め、ロジット変換(log(PKM/1−PKM))した値を目的変数とした。カプラン・マイヤー法による生命表は集団としての生存確率を示すが、ここでは、個人jを含む集団としての事象発生時点での残存確率(変化の発生しなかったものが残存する確率)を個人jの事象発生時点での残存時間に読み代えるという新規な考え方を用いている。また、この確率をロジット変換して、変化の発生傾向を表現するロジット値に変換して、目的変数とした。訓練集合内の交差検証はリーブ・ワン・アウト法によって行ない、PRESS×1.02NPが小さくなるようにパラメータを逐次取捨選択して部分最小自乗法モデルを得た。ここで、交差検証成績(CV=−PRESS)の代わりに、少なくとも交差検証成績を独立変数として持つ関数の1つである関数−PRESS×1.02NPを改善して部分最小自乗法モデルを得た。ここでPRESSはリーブ・ワン・アウト予測の残差自乗和であり、NPは、選択された説明変数の数である。
図7のフロー中の交差検証成績CVを−PRESS×1.02NPと読み換えて、処理を実行することにより、下記の19個の遺伝子の発現が説明変数として選抜された。ここでdata IDはWebデータ元でのID番号を示す。またACCESSIONはGenBankのアクセション番号であり、アクセション番号の無い行はデータ元でのみ明らかとなっている遺伝子(Unknown)ないしESTであり、論文記載の方法によって入手することができる。
これらの遺伝子の発現を説明変数の候補として比例ハザード(hazard)解析を試みた。比例ハザード法とは、生存率の解析に時間を考慮した統計的手法である。解析の実行はプログラムパッケージJMP(JMP Sales SAS Campus Drive Cary,NC 27513 USA)を用いて行なった。変数削除基準としてP≧0.05を採用した変数減少法によって更に絞り込んだ結果、14遺伝子の発現からなる以下の比例ハザード式が得られた。ここでGenbank(ジーンバンク)のアクセション番号ないしdata IDで示される各項は、各遺伝子のlog(χ1/x2)値であり、またPは統計的な有意性が成り立たない危険率である。この式の右辺から求められるハザード値(hazard)が大きいほど、死亡傾向が大きい。
Rosenwaldらは、単相関の比例ハザード解析を行なって、5群(17遺伝子)の診断指標を選抜している。図10に、本実施例で得られたハザード値(Hazard、図中Hazard(pls(14))と示した)とRosenwaldらの診断指標がどの程度、生存時間を説明できているかを比較した。Rosenwaldらの5群のパラメータを同時に用いた比例ハザード式ではProlifirationパラメータがP>0.05で統計的に有意でないなどの問題を有していため、これを除く4群のパラメータを同時に含めたハザード値も比較のために掲載した(図中Hazard(Rosenwald/4para)と示した)。ここで、菱形は死亡した人または打ち切った人のデータを示し、四角は生存している人のデータを示す。
これらの診断指標のうち、本実施例で求めたハザード値と生存時間との相関は際立って明白である。即ちハザード値は生存時間につれて減衰しており、大きなハザード値の患者は長く生きることが出来ないことが示されている。一方、Rosenwaldらの指標はいずれも生存時間を診断するには不十分なものである。数百、数千という数のパラメータの中から効率的に最適のパラメータセットを見出すことは比例ハザード解析だけではできないことである。しかし以上のようにカプラン・マイヤー法、ロジット変換、部分最小自乗法の交差検証成績を考慮した特徴抽出、比例ハザード解析を組み合わせることで、従来に無い、有効な診断指標を得ることができた。統計学的に異質なモデルをこのように組み合わせることによってこのような良好な結果が得られたことは意外でもあり、興味深いことであった。患者の生存時間を予測することは、QOLを含めた治療計画や人生設計などを判断する上で重要な情報をもたらすものであり、本実施例で求められた診断モデルは社会的に価値のあるものである。
また、変数削除基準としてP≧0.001を採用した変数減少法によって更に絞り込むと、6遺伝子の発現からなる以下の比例ハザード式が得られた。このように、変数削除基準を変えることにより、選択される説明変数の数を制御できる。
図11は、右辺を計算して求められるハザード値を縦軸とし、生存時間を横軸としたプロットを示す。図10と同様に、図11において、菱形は死亡した人または打ち切った人のデータを示し、四角は生存している人のデータを示す。
実施例3: 部分最小自乗法の交差検証成績を考慮した特徴抽出と比例ハザード解析による40名の乳癌患者の生存時間解析.
SorleらがWeb上(http://genome−www.stanford.edu/breast_cacer/mopo_clinical/)で公開している乳癌患者のデータセットをダウンロードして用いた。全データを訓練集合として利用した。データセットの大部分は、タイプA,Bという2種類のDNAチップで測定されたそれぞれ40名、24名の患者よりなるが、ここではタイプAのデータを用いた。生存時間データより実施例2と同様にロジット値を求め、目的変数とした。説明変数としては、データに欠測のある遺伝子を除いた6891件のLOG_RAT2N_MEAN値を採用した。そして、少なくとも交差検証成績を独立変数として持つ関数の1つである、交差検証成績と説明変数NPの関数PRESS×1.13NPが小さくなるようにパラメータを逐次取捨選択して部分最小自乗法モデルを得た。図7のフロー中の交差検証成績CVを−PRESS×1.13NPと読み換えて、処理を実行することにより、下記の10個の遺伝子の発現が説明変数として選抜された。
これらを説明変数の候補として、比例ハザード解析において変数削除基準としてP≧0.05を採用した変数減少法を試み、7遺伝子の発現からなる以下の比例ハザード式が得られた。ここでアクセッション番号で示される各項はそれぞれの遺伝子のLOG_RAT2N_MEANである。
図12に、右辺を計算して求められるハザード値を縦軸とし、生存時間を横軸としたプロットを示す。ここでもハザード値が優れた診断指標となることが示されている。図12において、菱形は死亡した人または打ち切った人のデータを示し、四角は生存している人のデータを示す。
変数削除基準としてP≧0.001を採用した変数減少法によって更に絞り込んだ。これにより、3遺伝子の発現からなる以下の比例ハザード式が得られた。このように、変数削除基準を変えることにより、説明変数の数を制御できた。
図13は、右辺を計算して求められるハザード値を縦軸とし、生存時間を横軸としたプロットを示す。ここで、菱形は死亡した人のデータを示し、四角は生存している人のデータを示す。
実施例4: 部分最小自乗法の交差検証成績を考慮した特徴抽出と重回帰分析による40名の乳癌患者の再発予測解析.
SorleらのDNAチップAで6891遺伝子の発現が測定された40名の患者をデータセットとして用いた。再発の有無を目的変数として、PRESS×1.10NPが小さくなるようにパラメータを逐次取捨選択して11遺伝子の発現からなる部分最小自乗法モデルを得た。
次に、選抜された遺伝子発現を説明変数とし、再発の有無を目的変数として、通常の多変数解析法の一つである重回帰分析によって判別分析を実行した。解析の実行はプログラムパッケージJMPを用いて行なった。変数削除基準としてP≧0.15を採用した変数減少法によってさらに絞り込んだ結果、10遺伝子の発現からなる以下の重回帰式が得られた。この式で計算されるOLS値が正の時は再発の可能性が高く、負の時は低い。
上式に含まれる各パラメータをそれぞれ1つ用いて判別分析式を作成した場合のP値及び決定係数を以下の表6に示す。
単独では有意とはならない(P>0.05)パラメータが3つ存在し、また、どのパラメータも決定係数が小さい。従って、パラメータを1つずつ吟味するだけでは、上式のような良好な判別式は得られなかった。また数百、数千という数のパラメータの中から効率的に最適のパラメータセットを見出すことは重回帰分析だけではできないことである。しかし、以上のように、部分最小自乗法の交差検証成績を考慮して特徴抽出することにより、従来に無い、有効な診断指標を得ることができた。乳癌の再発可能性を予測することは、QOLを考慮した治療計画を立案し判断するうえで、社会的に求められているところのものである。
実施例5: 部分最小自乗法の交差検証成績を考慮した特徴抽出と適応最小自乗法による40+24名の乳癌患者の再発予測解析.
DNAチップのタイプA(40名)とタイプB(24名)に共通する3448遺伝子に限って解析を試みた。PRESS×1.17NPが小さくなるようにパラメータを逐次取捨選択して部分最小自乗法モデルを得た。選抜された遺伝子発現を説明変数とし、適応最小自乗法によって判別分析を実行した結果、次式が得られた。次式で計算されるALS値が0.5より大きいと再発の危険性が存在する。
下記の表7にみるように、H11482は単相関では有意ではなく、他の変数と同時に用いることで初めて把握できたパラメータである。また、表8は、上式を用いてタイプBの患者を予測した結果である。本判別式の感度=81.8%、特異度=53.8%となり、χ2=3.233(5%<P<10%)、予測判別正解率=66.7%、という統計的に有意な結果を得た。タイプA、BはDNAチップの構成の相違に基づく測定誤差が存在すると思われるデータであるにもかかわらず、タイプAで訓練したモデルでタイプBの予測に危険率10%以下で成功したことは勇気付けられる結果である。
また、PRESS×1.12NPが小さくなるように選んだ場合には、以下の遺伝子の発現を説明変数とする部分最小自乗法モデルを得た。
これらを説明変数の候補として、リーブ・ワン・アウトを指標にして、さらに絞り込んだ結果、次の判別式を得た。
パラメータを1つずつ吟味するだけでは、上式のような良好な判別式は得られなかった。また数百、数千という数のパラメータの中から効率的に最適のパラメータセットを見出すことは、適応最小自乗法、ロジスティック回帰分析、その他の判別分析手法だけではできないことである。しかし、以上のように、部分最小自乗法の交差検証成績を考慮して特徴抽出することにより、従来に無い、有効な診断指標を得ることができた。
実施例6: 部分最小自乗法の交差検証成績を考慮した特徴抽出とロジスティック回帰分析法または線型判別分析法による40+24名の乳癌患者の再発予測解析.
実施例5での1つめの適応最小自乗法による解析をロジスティック回帰分析法に置き換えた場合、次の判別式が得られた。
右辺で求められるLORA値が正の場合には再発の危険性が存在する。係数の比率や相関係数は実施例5の適応最小自乗法の場合と異なるものの、各患者の識別結果は全く同一であった。またタイプBの患者を予測した結果も表7と同じになった。
次に、実施例5での適応最小自乗法による解析を線型判別分析に置き換えて解析して、次の判別式が得られた。
右辺で求められるLDA値が正の場合には再発の危険性が存在する。係数の比率や相関係数は、実施例5の適応最小自乗法の場合と異なり、各患者の識別結果も若干異なったが、概ね同一であった。また、タイプBの患者を予測した結果も表7と同じになった。
以上の実施例4,5,6では、乳癌の再発の有無を目的変数としている。したがって、部分最小自乗法の交差検証成績を考慮して特徴抽出する方法が、目的変数が名義尺度や順序尺度などのデータである場合にも有効であることが示された。なお、名義尺度とは、対象(サンプル)をある分類に属するかどうかを測り分けるときの分類で、分類の間に大小や順序はない。また、順序尺度とは、対象の特定の分類について測り分けるときの分類であり、分類の間に大小、高低といった順序がある。
実施例7: 部分最小自乗法の交差検証成績を考慮した特徴抽出と比例ハザード解析による40名の乳癌患者の再発時間解析.
実施例4と同じデータを用いて、再発の時系列データを基に実施例2と同様の方法で求めたロジット値を目的変数として、PRESS×1.15NPが小さくなるようにパラメータを逐次取捨選択して9遺伝子の発現からなる部分最小自乗法モデルを得た。これらの遺伝子発現の測定値を説明変数として比例ハザード解析において変数削除基準としてP≧0.05を採用した変数減少法を試み、8遺伝子からなる、以下の比例ハザード式が得られた。
図14は、右辺を計算して求められるハザード値を縦軸とし、再発時間を横軸としたプロットを示す。ここで、菱形は再発しない人のデータを示し、四角は再発している人のデータを示す。ここでもハザード値が優れた診断指標となっており、生存時間に限らず、時間とともに確率的に発生する生体の状態の変化を解析する手法として、本発明の手法が有効であることが示されている。
変数削除基準としてP≧0.005を採用した変数減少法によって更に絞り込んだ場合には、4遺伝子の発現からなる以下の比例ハザード式が得られた。
図15は、右辺を計算して求められるハザード値を縦軸とし、再発時間を横軸としたプロットを示す。ここで、菱形は再発しない人のデータを示し、四角は再発している人のデータを示す。
実施例8: Genbankアクセッション番号H11482、T64312、AA045340を含む乳癌再発性診断用DNAチップの作成と測定.
実験医学別冊「ゲノム機能研究プロトコール」(ISBN4−89706−932−7 C3047)p34−38記載の関直彦、永杉友美、東孝典、吉川勉、鈴木収、村松正明らの方法に準じてDNAチップの作成と測定を行なった。Genbankアクセッション番号H11482、T64312、AA045340のcDNAを用いた。
プローブ用の各PCR産物をエタノール(和光純薬,Cat#057−00456)で沈殿させ、2μg/μlとなるようにDDWで調整する。ニトロセルロース(GibcoBRL Cat#41051−012)4mg/mlのDMSO溶液を等量加え、よく混和させて100℃で5分間熱変性を行ない、氷上で急冷する。次いで室温に戻し、DNAスポッターSPBIO2000(日立ソフトエンジニアリング)を用いてカルボジイミドスライドガラス(日清紡)へのスポッティングを速やかに行なう。スポットの乾燥を確認し、Ultraviolet crosslinker(アマシャムファルマシアバイオテック社)を用いて60mJ/cm2で紫外クロスリンク処理を行ない、ガラスラックに立てて室温保存する。
3%BSA、0.2M NaCl、0.1M Tris(PH7.5)、0.05% Triton X−100よりなるブロッキング液に上記マイクロアレイを浸け、約30分間放置する。次いで、ガラスに付着している溶液をよく切り、37℃で乾燥させる。TEバッファー(PH8.0,ニッポンジーンCat #316−90025)で3回軽く洗い、プレートホルダーに入れて軽く遠心(1000rpm,1分間)して余分な水分を除去する。
次に、乳腺正常株SV−40及び乳癌細胞株MCF−7、MDA−MB−468又はT−47−Dの各細胞液より、TRIZOL(GibcoBRL,Cat#15596−018)、Oligotex dT30〈Super〉(TaKaRa,Cat#W9021A)を用いてマニュアルに従って、mRNAを精製する。2μgのmRNAをDEPC処理した6.4μlのDDWに溶かし、Oligo dTプライマー9μl、5×SuperScript IIバッファー(GibcoBRL,Cat#18089−011)6μl、DTT(SuperScriptの付属)3μl、50×dNTP 0.6μl、Cy3−dUTP(アマシャムファルマシアバイオテク Cat# PA53022)又はCy5−dUTP(アマシャムファルマシアバイオテク Cat# PA55022)3μl、SuperScript II 2μlよりなる溶液を加え、42℃で2時間反応させる。途中1時間経過時点で、SuperScript IIを1μlを追加する。1.5μlアルカリバッファー(1N NaOH/20nM EDTA)を加え、65℃で10分間反応させ、TEバッファーを270μl、1N HClを1.5μl加えて、Cy3,Cy5ラベルの反応液を2つまとめて1本のMicrocon−YM−30(Millipore/Amicon,Cat#42410)に移す。10,000rpmで上のカップに残る液量が約10μlになるまで遠心を続け、カップを通りぬける液を別のチューブに移し替え、その後、上のカップにTEバッファー500μl、Human Cot−1 DNA(GibcoBRL Cat#15279−011)20μgを加え、再び液量が10μl以下になるまで遠心を続ける。3,000rpmで3分間遠心し、蛍光標識したDNAを回収する。DDWとyeast RNA(Sigma,Cat#R7125)50μg、poly(A)(ロッシュダイアグノスティクス,Cat#108 626)50μgを加えて20μlにし、PCR用のチューブに移し換え、さらに4.25μl 20×SSC(GibcoBRL,Cat#15553−035)と0.75μl 10% SDS(GibcoBRL,Cat#15553−035)を加え、PCR用の機器で100℃、1分間熱変性させ、次いで、室温で30分間放置して、ゆっくり冷却する。
蛍光標識したDNAの全量をカバーガラスにのせ、泡が入らないように注意しながら前記マイクロアレイにかぶせ、水で濡らしたキムタオルを底に敷いたハイブリダイゼーションチェンバーに入れて密閉する。毎分2〜4サイクルで軽く振とうさせながら、65℃で一晩ハイブリダイズする。ハイブリダイゼーションチェンバーからマイクロアレイを取り出し、カバーガラスが載ったままの状態で静かに2×SSC/0.1% SDS溶液中に入れ、5分間シェイキングし、カバーガラスが自然にはがれるのを待つ。カバーガラスがはがれたところでマイクロアレイをスライドガラスラックに入れ、もう一度2×SSC/0.1% SDS溶液中で5分間軽く振とうして洗う。さらに0.2×SSC/0.1% SDS40℃で5分間2回洗い、0.2×SSCでリンスする。マイクロアレイを別の乾いたプレパラートケースに移し、マイクロタイタープレート用の遠心機で軽く遠心して(1000rpm,1分室温)マイクロアレイ上の水分を除く。そして、ScanArray4000(GSI luminonics社)でシグナルを読み込み、解析ソフトにはQuant Array(GSI luminonics社)およびChip Space(日立ソフトウェアエンジニアリング)を用いる。
実施例9: 遺伝的アルゴリズムによる部分最小自乗法モデルの最適化.
実施例4で用いたSorleらのDNAチップAで6891遺伝子の発現が測定された40名の患者をデータセットとして用いた。遺伝的アルゴリズムは、たとえば、伊庭斉志;「遺伝的アルゴリズムの基礎」(オーム社(1994))に説明されている。前記データを用い、遺伝的アルゴリズムによる説明変数選択を行なった。以下において「」で区切られた用語は遺伝的アルゴリズムで通常用いられる専門用語であり、特に必要な場合には解説を加えている。「適合度」(fitness)には−PRESS×1.01Npを採用した。各「個体」の「遺伝型」は説明変数を採用する場合には1、採用しない場合には0をとる数列{b1,b2,b3,...}とした。
個体集合のサイズを100個とし、初期の個体の「遺伝型」(GTYPE)は、平均でmin_of(Ns,Ng,300)/2個の説明変数が採用となるように乱数を用いて準備した。ここでNsはサンプル数(患者数)、Ngは説明変数の候補の数、300は実装の都合上設定された定数である。
集合よりランダムに2つの個体を選抜し、「遺伝型」の「一様交叉」を行なったものの一方を新しい「個体」とした。即ち、「各遺伝子座」ごとに1/2の確率でいずれかの「親個体」の数列値(0または1)を選びそれを代入したものを新しい「個体」とした。続いて新しい「個体」の「各遺伝子座」毎に、1の場合(説明変数が採用されている場合)には1.1/採用された説明変数の数の確率で、0の場合(採用されていない場合)には1.1/採用されていない説明変数候補の数の確率で、0←→1を反転させた。
上述の「交叉・突然変異オペレーション」によって準備された新しい「個体」の「適合度」と、ランダムに選抜された「トーナメント相手」となる集合中の「個体」の「適合度」とを比較し、新しい「個体」の適合度が勝った場合には0.75の確率で、劣った場合には0.25の確率で「個体」の置き換えを行なった。ただし、「トーナメント相手」が集合中の最適解のものである場合には置き換えを禁止するという「エリート戦略」を採用した。
以上の「交叉」→「突然変異」→「選抜」サイクルを繰り返して最適化を行なった。ここではサイクル数を集合サイズで割ったものを「世代数」とする。最大「世代数」の初期値を100とし、新しい最適解が見出されるたびに最大「世代数」を10増加させながら、実行「世代数」が最大「世代数」に至るまでサイクルを繰り返した。
以上の初期集合の準備〜最適化の繰り返しおよび終了にいたる一連の処理を一回のラン(run)とし、15回のランを行なった。図16は、15回のランにおける最適化の様子をまとめている。最良の結果は25個の説明変数を用いたものである。
実施例10: 階層型人工ニューラルネットワーク(MLP)によるモデル構築.
実施例5の乳癌患者の再発性判別解析において、DNAチップtype A(40名)とtype B(24名)に共通する3448遺伝子より、PRESS×1.17Npが小さくなるようにしてPLS−CVで特徴抽出された3つの説明変数を用いた。
解析方法について説明すると、MLPは3層とし、中間層(tk)において一度だけシグモイド変換を行なう構造とし、図17の4つのトポロジーを試みた。ネットワークの重みの学習はBack propagation(逆伝播)アルゴリズムによって行なった。中間層(tk)において一度だけシグモイド変換を行なう3層MLPを用いた。
ネットワークトポロジーIおよびトポロジーIIbの結果は以下のとおりであった。なお、トポロジーIIa及びトポロジーIIcは、トポロジーIIbに劣るものであった。
実施例11: 潜在変数を用いた比例ハザードモデルの構築.
実施例3のPLS−CV法で選抜された10遺伝子の発現量を説明変数とし、目的変数として生存確率のlogit値を用いてPLSの解析過程で作成される潜在変数を1個抽出した。その抽出した潜在変数を説明変数にして比例ハザードモデルによる解析を試みた結果、作成された式はP≦0.0001で有意となった。図18に右辺を計算して得られるハザード値を縦軸とし、生存時間を横軸にしたプロットを示す。
本技術で得られたハザード式の予測の性能を評価するために、用いた40例の中から1例を除外し、残りの39例のデータを用いてハザード式を作成し、除外した1例のハザード値を予測した。39例からのハザード式によって予測した値と40例からのハザード式からの計算値をプロットした図19より、本技術はハザード値の予測において良好な成績を示した。
発明の効果について以下に説明すると、生体の状態と複数の遺伝子発現の量および/または細胞内物質の量との相関モデルを決定するとき、説明変数の選択と交差検証法とを用いて変数を絞り込むことができる。これにより、良好でかつ予測力のある多変量解析モデル(相関モデル)が得られる。特に遺伝子発現の量のように、説明変数の数がたとえば1000以上と膨大な場合に有用である。変数の数を少なくすることにより、病気や生体現象の背後で働いている重要な遺伝子やメカニズムを推定/特定でき、理解が深まる。また、重要な遺伝子産物や細胞内物質だけに絞った廉価な診断用材料(DNAチップ、DNA含有ベクター、抗体チップなど)を設計し、提供できる。
また、時間とともに確率的に発生する生体の状態の変化から導出された量を目的変数として用いて、時間とともに確率的に発生する生体の状態の変化と複数の遺伝子発現の量および/または細胞内物質の量との相関モデルを決定できる。
また、部分最小自乗法を用いて説明変数の個数を少なくすると、通常の統計的手法または多変量解析手法が適用可能になる。
【図面の簡単な説明】
図1は、遺伝子発現解析システムのブロック図である。
図2は、解析ソフトのフローチャートである。
図3は、交差検証成績CVの計算のフローチャートである。
図4は、変数選択の第1モデル構築手法のフローチャートである。
図5は、変数選択の第2モデル構築手法のフローチャートである。
図6は、変数選択の第3モデル構築手法のフローチャートである。
図7は、変数選択の第4モデル構築手法のフローチャートである。
図8は、変数選択の第5モデル構築手法のフローチャートである。
図9は、最小自乗法モデルの成績を示すグラフである。
図10は、DLBCL患者の生存時間と診断指標のプロット各種比較の図である。
図11は、実施例2のDLBCL患者の生存時間診断指標のプロットの図である。
図12は、実施例3の乳癌患者の生存時間診断指標のプロットの図である。
図13は、実施例3の乳癌患者の変数削除基準としてP≧0.0005を採用したときの生存時間診断指標のプロットの図である。
図14は、実施例7の乳癌患者の再発時間診断指標のプロットの図である。
図15は、実施例7の乳癌患者の変数削除基準としてP≧0.025を採用したときの再発時間診断指標のプロットの図である。
図16は、実施例9の遺伝的アルゴリズムによる部分最小自乗法モデルの最適化の様子を示す図である。
図17は、実施例10の階層型人工ニューラルネットワークにおける4つのトポロジーを示す図である。
図18は、実施例11の潜在変数を用いた比例ハザードモデルの乳癌患者の生存時間診断指標のグラフである。
図19は、実施例11の潜在変数を用いた比例ハザードモデルの乳癌患者の生存時間診断指標の予測値と計算値のグラフである。
Claims (61)
- 生体の状態または時間とともに確率的に発生する生体の状態の変化を目的変数とし、複数の遺伝子発現の量および/または細胞内物質の量を説明変数とする相関モデルを決定するデータ解析装置であって、
生体の状態或いはそれを導出するデータまたは時間とともに確率的に発生する生体の状態の変化に関するデータと、複数の遺伝子発現の量および/または細胞内物質の量からなるサンプルの集合を入力する入力手段と、
(1)説明変数を選択する選択手段と、
(2)部分最小自乗法を実行して交差検証成績を計算する計算手段または前記生体の状態の変化に関するデータにカプラン・マイヤー法又はカトラー・エデラー法による生命表を適用して変化の発生しなかったものの確率を計算して得られた確率を、仮定した分布に基づいた変換または仮定を前提としない変換をし、該変換結果を目的変数とする部分最小自乗法を実行して交差検証成績を計算する計算手段と、
(3)前記(2)の計算手段の結果を評価し、説明変数の採用、不採用を判定する評価判定手段とを有し、
(4)前記(1)の選択手段と前記(2)の計算手段と前記(3)の評価判定手段とを実行して部分最小自乗法モデルの少なくとも交差検証成績を独立変数として持つ関数を改善し続けて部分最小自乗法モデルを決定する決定手段とからなることを特徴とするデータ解析装置。 - 目的変数が生体の状態であって、前記入力手段で入力するデータが生体の状態或いはそれを導出するデータであって、前記(2)の計算手段が部分最小自乗法を実行して交差検証成績を計算する計算手段であることを特徴とする請求項1に記載のデータ解析装置。
- 目的変数が時間とともに確率的に発生する生体の状態の変化であって、前記入力手段で入力するデータが時間とともに確率的に発生する生体の状態の変化に関するデータであって、前記(2)の計算手段が前記生体の状態の変化に関するデータにカプラン・マイヤー法又はカトラー・エデラー法による生命表を適用して変化の発生しなかったものの確率を計算して得られた確率を、仮定した分布に基づいた変換または仮定を前提としない変換をし、該変換結果を目的変数とする部分最小自乗法を実行して交差検証成績を計算する計算手段であることを特徴とする請求項1に記載のデータ解析装置。
- さらに、前記の決定手段にて決定された部分最小自乗法モデルに採用されている説明変数又は該モデルの潜在変数を用い、統計的手法又は多変量解析手法によるモデルを構築する最終モデル決定手段を備えることを特徴とする請求項1、2又は3に記載のデータ解析装置。
- 前記の選択手段において、説明変数を逐次取捨選択することを特徴とする請求項1〜4のいずれかに記載のデータ解析装置。
- 前記の選択手段において、遺伝的アルゴリズムを用いて説明変数を選択することを特徴とする請求項1〜4のいずれかに記載のデータ解析装置。
- 前記の計算手段において、1個のサンプルを逐次除外して部分最小自乗法を実行して交差検証成績を計算することを特徴とする請求項1〜6のいずれかに記載のデータ解析装置。
- 前記の計算手段において、複数のサンプルを逐次除外して部分最小自乗法を実行して交差検証成績を計算することを特徴とする請求項1〜6のいずれかに記載のデータ解析装置。
- 前記計算手段において、各計算において除外したサンプルの遺伝子発現から予測される生体の状態を示す目的変数値と、前記除外したサンプルの生体の状態を示す目的変数値との誤差の代表値を求め、交差検証成績の指標として当該誤差を用いることを特徴とする請求項7又は8に記載のデータ解析装置。
- 前記関数が交差検証成績であることを特徴とする請求項1〜9のいずれかに記載のデータ解析装置。
- 前記関数が交差検証成績と選択された説明変数の個数との関数であることを特徴とする請求項1〜9のいずれかに記載のデータ解析装置。
- 前記の決定手段において、少なくとも交差検証成績を独立変数として持つ関数を改善しながら評価判定を繰り返すことを特徴とする請求項5に記載のデータ解析装置。
- 前記(1)の選択手段と前記(2)の計算手段とを複数のコンピュータで実行させることを特徴とする請求項1〜12のいずれかに記載のデータ解析装置。
- 請求項1、2、3又は4で決定された相関モデル及び予測対象のサンプルについて当該モデルにおいて採用された説明変数を入力する入力手段と、入力された該説明変数に基づいて該サンプルの生体の状態を予測判定する予測判定手段からなることを特徴とするデータ解析装置。
- 生体の状態を名義尺度、順序尺度或いは連続量で表現する目的変数とする請求項2に記載のデータ解析装置。
- 最終モデル決定手段が用いる前記の統計的手法又は多変量解析手法が、比例ハザード法又はパラメトリックな分布にあてはめた回帰分析法であることを特徴とする請求項2又は4に記載のデータ解析装置。
- 生体の状態または時間とともに確率的に発生する生体の状態の変化を目的変数とし、複数の遺伝子発現の量および/または細胞内物質の量を説明変数とする相関モデルを決定するデータ解析方法であって、
生体の状態或いはそれを導出するデータまたは時間とともに確率的に発生する生体の状態の変化に関するデータと、複数の遺伝子発現の量および/または細胞内物質の量からなるサンプルの集合を入力する入力ステップと、
(1)説明変数を選択する選択ステップと、
(2)部分最小自乗法を実行して交差検証成績を計算する計算ステップまたは前記生体の状態の変化に関するデータにカプラン・マイヤー法又はカトラー・エデラー法による生命表を適用して変化の発生しなかったものの確率を計算して得られた確率を、仮定した分布に基づいた変換または仮定を前提としない変換をし、該変換結果を目的変数とする部分最小自乗法を実行して交差検証成績を計算する計算ステップと、
(3)前記(2)の計算ステップの結果を評価し、説明変数の採用、不採用を判定する評価判定ステップとを有し、
(4)前記(1)の選択ステップと前記(2)の計算ステップと前記(3)の評価判定ステップとを実行して部分最小自乗法モデルの少なくとも交差検証成績を独立変数として持つ関数を改善し続けて部分最小自乗法モデルを決定する決定ステップとからなることを特徴とするデータ解析方法。 - 目的変数が生体の状態であって、前記入力ステップで入力するデータが生体の状態或いはそれを導出するデータであって、前記(2)の計算ステップが部分最小自乗法を実行して交差検証成績を計算する計算ステップであることを特徴とする請求項17に記載のデータ解析方法。
- 目的変数が時間とともに確率的に発生する生体の状態の変化であって、前記入力ステップで入力するデータが時間とともに確率的に発生する生体の状態の変化に関するデータであって、前記(2)の計算ステップが前記生体の状態の変化に関するデータにカプラン・マイヤー法又はカトラー・エデラー法による生命表を適用して変化の発生しなかったものの確率を計算して得られた確率を、仮定した分布に基づいた変換または仮定を前提としない変換をし、該変換結果を目的変数とする部分最小自乗法を実行して交差検証成績を計算する計算ステップであることを特徴とする請求項17に記載のデータ解析方法。
- さらに、前記の決定ステップにて決定された部分最小自乗法モデルに採用されている説明変数又は該モデルの潜在変数を用い、統計的手法又は多変量解析手法によるモデルを構築する最終モデル決定ステップを備えることを特徴とする請求項17、18又は19に記載のデータ解析方法。
- 前記の選択ステップにおいて、説明変数を逐次取捨選択することを特徴とする請求項17〜20のいずれかに記載のデータ解析方法。
- 前記の選択ステップにおいて、遺伝的アルゴリズムを用いて説明変数を選択することを特徴とする請求項17〜20のいずれかに記載のデータ解析方法。
- 前記の計算ステップにおいて、1個のサンプルを逐次除外して部分最小自乗法を実行して交差検証成績を計算することを特徴とする請求項17〜22のいずれかに記載のデータ解析方法。
- 前記の計算ステップにおいて、複数のサンプルを逐次除外して部分最小自乗法を実行して交差検証成績を計算することを特徴とする請求項17〜22のいずれかに記載のデータ解析方法。
- 前記計算ステップにおいて、各計算において除外したサンプルの遺伝子発現から予測される生体の状態を示す目的変数値と、前記除外したサンプルの生体の状態を示す目的変数値との誤差の代表値を求め、交差検証成績の指標として当該誤差を用いることを特徴とする請求項23又は24に記載のデータ解析方法。
- 前記関数が交差検証成績であることを特徴とする請求項17〜25のいずれかに記載のデータ解析方法。
- 前記関数が交差検証成績と選択された説明変数の個数との関数であることを特徴とする請求項17〜25のいずれかに記載のデータ解析方法。
- 前記決定ステップにおいて、少なくとも交差検証成績を独立変数として持つ関数を改善しながら評価判定を繰り返すことを特徴とする請求項21に記載のデータ解析方法。
- 前記(1)の選択ステップと前記(2)の計算ステップとを複数のコンピュータで実行させることを特徴とする請求項17〜28のいずれかに記載のデータ解析方法。
- 請求項17、18、19又は20で決定された相関モデル及び予測対象のサンプルについて当該モデルにおいて採用された説明変数を入力する入力ステップと、入力された該説明変数に基づいて該サンプルの生体の状態を予測判定する予測判定ステップからなることを特徴とするデータ解析方法。
- 生体の状態を名義尺度、順序尺度或いは連続量で表現する目的変数とする請求項18に記載のデータ解析方法。
- 前記の統計的手法又は多変量解析手法が、比例ハザード法又はパラメトリックな分布にあてはめた回帰分析法によるモデルを構築する最終モデル決定ステップとからなることを特徴とする請求項18又は20に記載のデータ解析方法。
- 生体の状態または時間とともに確率的に発生する生体の状態の変化を目的変数とし、複数の遺伝子発現の量および/または細胞内物質の量を説明変数とする相関モデルを決定する、コンピュータにより実行されるデータ解析プログラムであって、
生体の状態或いはそれを導出するデータまたは時間とともに確率的に発生する生体の状態の変化に関するデータと、複数の遺伝子発現の量および/または細胞内物質の量からなるサンプルの集合を入力する入力ステップと、
(1)説明変数を選択する選択ステップと、
(2)部分最小自乗法を実行して交差検証成績を計算する計算ステップまたは前記生体の状態の変化に関するデータにカプラン・マイヤー法又はカトラー・エデラー法による生命表を適用して変化の発生しなかったものの確率を計算して得られた確率を、仮定した分布に基づいた変換または仮定を前提としない変換をし、該変換結果を目的変数とする部分最小自乗法を実行して交差検証成績を計算する計算ステップと、
(3)前記(2)の計算ステップの結果を評価し、説明変数の採用、不採用を判定する評価判定ステップとを有し、
(4)前記(1)の選択ステップと前記(2)の計算ステップと前記(3)の評価判定ステップとを実行して部分最小自乗法モデルの少なくとも交差検証成績を独立変数として持つ関数を改善し続けて部分最小自乗法モデルを決定する決定ステップとからなることを特徴とするデータ解析プログラム。 - 目的変数が生体の状態であって、前記入力ステップで入力するデータが生体の状態或いはそれを導出するデータであって、前記(2)の計算ステップが部分最小自乗法を実行して交差検証成績を計算する計算ステップであることを特徴とする請求項33に記載のデータ解析プログラム。
- 目的変数が時間とともに確率的に発生する生体の状態の変化であって、前記入力ステップで入力するデータが時間とともに確率的に発生する生体の状態の変化に関するデータであって、前記(2)の計算ステップが前記生体の状態の変化に関するデータにカプラン・マイヤー法又はカトラー・エデラー法による生命表を適用して変化の発生しなかったものの確率を計算して得られた確率を、仮定した分布に基づいた変換または仮定を前提としない変換をし、該変換結果を目的変数とする部分最小自乗法を実行して交差検証成績を計算する計算ステップであることを特徴とする請求項33に記載のデータ解析プログラム。
- さらに、前記の決定ステップにて決定された部分最小自乗法モデルに採用されている説明変数又は該モデルの潜在変数を用い、統計的手法又は多変量解析手法によるモデルを構築する最終モデル決定ステップを備えることを特徴とする請求項33、34又は35に記載のデータ解析プログラム。
- 前記の選択ステップにおいて、説明変数を逐次取捨選択することを特徴とする請求項33〜36のいずれかに記載のデータ解析プログラム。
- 前記の選択ステップにおいて、遺伝的アルゴリズムを用いて説明変数を選択することを特徴とする請求項33〜36のいずれかに記載のデータ解析プログラム。
- 前記の計算ステップにおいて、1個のサンプルを逐次除外して部分最小自乗法を実行して交差検証成績を計算することを特徴とする請求項33〜38のいずれかに記載のデータ解析プログラム。
- 前記の計算ステップにおいて、複数のサンプルを逐次除外して部分最小自乗法を実行して交差検証成績を計算することを特徴とする請求項33〜38のいずれかに記載のデータ解析プログラム。
- 前記計算ステップにおいて、各計算において除外したサンプルの遺伝子発現から予測される生体の状態を示す目的変数値と、前記除外したサンプルの生体の状態を示す目的変数値との誤差の代表値を求め、交差検証成績の指標として当該誤差を用いることを特徴とする請求項39又は40に記載のデータ解析プログラム。
- 前記関数が交差検証成績であることを特徴とする請求項33〜41のいずれかに記載のデータ解析プログラム。
- 前記関数が交差検証成績と選択された説明変数の個数との関数であることを特徴とする請求項33〜41のいずれかに記載のデータ解析プログラム。
- 前記決定ステップにおいて、少なくとも交差検証成績を独立変数として持つ関数を改善しながら評価判定を繰り返すことを特徴とする請求項37に記載のデータ解析プログラム。
- 前記(1)の選択ステップと前記(2)の計算ステップとを複数のコンピュータで実行させることを特徴とする請求項33〜44のいずれかに記載のデータ解析プログラム。
- 請求項33、34、35又は36で決定された相関モデル及び予測対象のサンプルについて当該モデルにおいて採用された説明変数を入力する入力ステップと、入力された該説明変数に基づいて該サンプルの生体の状態を予測判定する予測判定ステップからなることを特徴とするデータ解析プログラム。
- 生体の状態を名義尺度、順序尺度或いは連続量で表現する目的変数とする請求項34に記載のデータ解析プログラム。
- 前記の統計的手法又は多変量解析手法が、比例ハザード法又はパラメトリックな分布にあてはめた回帰分析法によるモデルを構築する最終モデル決定ステップとからなることを特徴とする請求項34又は36に記載のデータ解析プログラム。
- 前記の説明変数の選択において、初期状態では説明変数を全く含まないことを特徴とする請求項37に記載のプログラム。
- 前記の説明変数の選択において、初期状態では全説明変数を含むことを特徴とする請求項37に記載のプログラム。
- 前記の生体の状態が病気のタイプをあらわす測定値、病気の重篤度をあらわす測定値、病気のタイプをあらわす医療診断の結果、病気の重篤度をあらわす医療診断の結果、あるいはそれらを2次加工した数値であることを特徴とする請求項37〜50のいずれかに記載のプログラム。
- 請求項33〜請求項48のいずれかに記載されたプログラムを記録した、コンピュータにより読み取り可能な記録媒体。
- 実質的にジーンバンクアクセッション番号がU15085、M23452、X52479、U70426、H57330及びS69790からなる遺伝子群の発現を検出することを特徴とするびまん性大細胞型Bリンパ腫の重篤度検定用の細胞内物質測定機材および測定方法並びにびまん性大細胞型Bリンパ腫の重篤度検定方法。
- さらにジーンバンクアクセッション番号がU03398、M65066、AK001546、BC003536、X00437、U12979、H96306、AA830781及びAA804793からなる群から選択される少なくとも一つの遺伝子の発現を検出することを特徴とする請求項53に記載のびまん性大細胞型Bリンパ腫の重篤度検定用の細胞内物質測定機材および測定方法並びにびまん性大細胞型Bリンパ腫の重篤度検定方法。
- 実質的にジーンバンクアクセッション番号がAA598572、AA703058及びAA453345からなる遺伝子産物を含む細胞内物質を検出することを特徴とする乳癌の重篤度検定用の細胞内物質測定機材および測定方法並びに乳癌の重篤度検定方法。
- さらにジーンバンクアクセッション番号がAA406242、H73335、W84753、N71160、AA054669、N32820及びR05667からなる群から選択される少なくとも一つの遺伝子産物を含む細胞内物質を検出することを特徴とする請求項55に記載の乳癌の重篤度検定用の細胞内物質測定機材および測定方法並びに乳癌の重篤度検定方法。
- 実質的にジーンバンクアクセッション番号がW84753、H08581、AA045730及びAI250654からなる遺伝子産物を含む細胞内物質を検出することを特徴とする乳癌の再発性検定用の細胞内物質測定機材および測定方法並びに乳癌の再発性検定方法。
- さらにジーンバンクアクセッション番号がAA448641、R78516、R05934、AA629838及びH53037からなる群から選択される少なくとも一つの遺伝子産物を含む細胞内物質を検出することを特徴とする請求項57に記載の乳癌の再発性検定用の細胞内物質測定機材および測定方法並びに乳癌の再発性検定方法。
- 実質的にジーンバンクアクセッション番号がAA434397、T83209、N53427、N29639、AA485739、AA425861、H84871、T64312、T59518及びAA037488からなる遺伝子産物を含む細胞内物質を検出することを特徴とする乳癌の再発性検定用の細胞内物質測定機材および測定方法並びに乳癌の再発性検定方法。
- さらにジーンバンクアクセッション番号がAA406231の遺伝子産物を含む細胞内物質を検出することを特徴とする請求項59に記載の乳癌の再発性検定用の細胞内物質測定機材および測定方法並びに乳癌の再発性検定方法。
- 実質的にジーンバンクアクセッション番号がH11482、T64312及びAA045340からなる遺伝子産物を含む細胞内物質を検出することを特徴とする乳癌の再発性検定用の細胞内物質測定機材および測定方法並びに乳癌の再発性検定方法。
Applications Claiming Priority (5)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
JP2002102743 | 2002-04-04 | ||
JP2002102743 | 2002-04-04 | ||
JP2002352645 | 2002-12-04 | ||
JP2002352645 | 2002-12-04 | ||
PCT/JP2003/004059 WO2003085548A1 (fr) | 2002-04-04 | 2003-03-31 | Dispositif et procede d'analyse de donnees |
Publications (1)
Publication Number | Publication Date |
---|---|
JPWO2003085548A1 true JPWO2003085548A1 (ja) | 2005-08-11 |
Family
ID=28793526
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
JP2003582665A Pending JPWO2003085548A1 (ja) | 2002-04-04 | 2003-03-31 | データ解析装置および方法 |
Country Status (8)
Country | Link |
---|---|
US (1) | US20050159896A1 (ja) |
EP (1) | EP1498825A1 (ja) |
JP (1) | JPWO2003085548A1 (ja) |
KR (1) | KR20040111456A (ja) |
CN (1) | CN1647067A (ja) |
AU (1) | AU2003220998A1 (ja) |
CA (1) | CA2481485A1 (ja) |
WO (1) | WO2003085548A1 (ja) |
Families Citing this family (28)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US7747391B2 (en) | 2002-03-01 | 2010-06-29 | Maxygen, Inc. | Methods, systems, and software for identifying functional biomolecules |
US20050084907A1 (en) | 2002-03-01 | 2005-04-21 | Maxygen, Inc. | Methods, systems, and software for identifying functional biomolecules |
US20070224126A1 (en) * | 2004-06-01 | 2007-09-27 | Therese Dufresne | Index and Method of use of Adapted Food Compositions for Dysphagic Persons |
JP5002811B2 (ja) * | 2004-10-26 | 2012-08-15 | 国立大学法人横浜国立大学 | 多変数モデル解析システム、方法、プログラム、およびプログラム媒体 |
JPWO2006088208A1 (ja) * | 2005-02-21 | 2008-07-10 | 大日本住友製薬株式会社 | 生体の生理変化の予測方法および装置 |
JPWO2006098192A1 (ja) * | 2005-03-16 | 2008-08-21 | 味の素株式会社 | 生体状態評価装置、生体状態評価方法、生体状態評価システム、生体状態評価プログラム、評価関数作成装置、評価関数作成方法、評価関数作成プログラムおよび記録媒体 |
JP4714869B2 (ja) * | 2005-12-02 | 2011-06-29 | 国立大学法人山口大学 | 有効因子抽出システム |
EP1804172B1 (en) * | 2005-12-20 | 2021-08-11 | Roche Diagnostics GmbH | PCR elbow determination using curvature analysis of a double sigmoid |
JP5011830B2 (ja) * | 2006-06-09 | 2012-08-29 | 富士通セミコンダクター株式会社 | データ処理方法、データ処理プログラム、該プログラムを記録した記録媒体およびデータ処理装置 |
US9500656B2 (en) * | 2006-08-10 | 2016-11-22 | Millennium Pharmaceuticals, Inc. | Methods for the identification, assessment, and treatment of patients with cancer therapy |
JP5307996B2 (ja) * | 2006-09-06 | 2013-10-02 | 株式会社Dnaチップ研究所 | 判別因子セットを特定する方法、システム及びコンピュータソフトウェアプログラム |
EP2069990B1 (en) * | 2006-09-20 | 2019-03-13 | Koninklijke Philips N.V. | A molecular diagnostics decision support system |
US8374795B2 (en) * | 2008-05-13 | 2013-02-12 | Roche Molecular Systems, Inc. | Systems and methods for step discontinuity removal in real-time PCR fluorescence data |
JP2012256182A (ja) * | 2011-06-08 | 2012-12-27 | Sharp Corp | データ解析装置、データ解析方法およびデータ解析プログラム |
CN102539326B (zh) * | 2012-01-13 | 2014-03-12 | 江苏大学 | 茶叶汤色品质的量化评价方法 |
JP5794160B2 (ja) * | 2012-01-26 | 2015-10-14 | 富士通株式会社 | 説明変数の決定のための情報処理装置、情報処理方法及びプログラム |
JP2014100249A (ja) * | 2012-11-19 | 2014-06-05 | Toshiba Corp | 血管解析装置、医用画像診断装置、血管解析方法、及び血管解析プログラム |
SG11201505977RA (en) | 2013-01-31 | 2015-08-28 | Codexis Inc | Methods, systems, and software for identifying bio-molecules using models of multiplicative form |
CN103324866A (zh) * | 2013-03-26 | 2013-09-25 | 张弘 | Ripple系统 |
JP6059122B2 (ja) * | 2013-10-11 | 2017-01-11 | カルチュア・コンビニエンス・クラブ株式会社 | 顧客データ解析システム |
US9928516B2 (en) * | 2013-12-30 | 2018-03-27 | Nice Ltd. | System and method for automated analysis of data to populate natural language description of data relationships |
EP3155592B1 (en) * | 2014-06-10 | 2019-09-11 | Leland Stanford Junior University | Predicting breast cancer recurrence directly from image features computed from digitized immunohistopathology tissue slides |
EP3316159A4 (en) | 2015-06-25 | 2019-08-14 | Advanced Telecommunications Research Institute International | PREDICTIVE EQUIPMENT BASED ON A SYSTEM ASSOCIATED WITH SEVERAL INSTITUTIONS AND PREDICTION PROGRAM |
WO2017170803A1 (ja) | 2016-03-29 | 2017-10-05 | 株式会社国際電気通信基礎技術研究所 | 医薬組成物又は食品組成物、及び有効成分の体内での効果の評価方法 |
EP3640946A1 (en) * | 2018-10-15 | 2020-04-22 | Sartorius Stedim Data Analytics AB | Multivariate approach for biological cell selection |
US11410064B2 (en) * | 2020-01-14 | 2022-08-09 | International Business Machines Corporation | Automated determination of explanatory variables |
JP7214672B2 (ja) * | 2020-03-13 | 2023-01-30 | 株式会社東芝 | 情報処理装置、情報処理方法及びコンピュータプログラム |
CN111652302B (zh) * | 2020-05-28 | 2023-05-23 | 泰康保险集团股份有限公司 | 一种解释保险核保分类结果的方法、装置 |
Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JPH05233011A (ja) * | 1992-02-18 | 1993-09-10 | Nippon Telegr & Teleph Corp <Ntt> | 独立要因抽出法 |
JPH09167152A (ja) * | 1995-12-19 | 1997-06-24 | Hitachi Ltd | 対話的モデル作成方法 |
WO2000070340A2 (en) * | 1999-05-14 | 2000-11-23 | Karolinska Innovations Ab | Materials and methods relating to disease diagnosis |
WO2002025405A2 (en) * | 2000-09-19 | 2002-03-28 | The Regents Of The University Of California | Methods for classifying high-dimensional biological data |
-
2003
- 2003-03-31 CA CA002481485A patent/CA2481485A1/en not_active Abandoned
- 2003-03-31 KR KR10-2004-7015515A patent/KR20040111456A/ko not_active Application Discontinuation
- 2003-03-31 JP JP2003582665A patent/JPWO2003085548A1/ja active Pending
- 2003-03-31 US US10/509,886 patent/US20050159896A1/en not_active Abandoned
- 2003-03-31 CN CNA038075237A patent/CN1647067A/zh active Pending
- 2003-03-31 WO PCT/JP2003/004059 patent/WO2003085548A1/ja active Application Filing
- 2003-03-31 EP EP03715637A patent/EP1498825A1/en not_active Withdrawn
- 2003-03-31 AU AU2003220998A patent/AU2003220998A1/en not_active Abandoned
Patent Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JPH05233011A (ja) * | 1992-02-18 | 1993-09-10 | Nippon Telegr & Teleph Corp <Ntt> | 独立要因抽出法 |
JPH09167152A (ja) * | 1995-12-19 | 1997-06-24 | Hitachi Ltd | 対話的モデル作成方法 |
WO2000070340A2 (en) * | 1999-05-14 | 2000-11-23 | Karolinska Innovations Ab | Materials and methods relating to disease diagnosis |
WO2002025405A2 (en) * | 2000-09-19 | 2002-03-28 | The Regents Of The University Of California | Methods for classifying high-dimensional biological data |
Non-Patent Citations (3)
Title |
---|
ALIZADEH A A: "Distinct types of diffuse large B-cell lymphoma identified by gene expression profiling", NATURE, vol. V403, JPN5003008150, 3 February 2000 (2000-02-03), GB, pages 503 - 512, ISSN: 0001246421 * |
J.KHAN: "Classification and diagnostic prediction of cancers using gene expression profiling and artificial n", NATURE MEDICINE, vol. V7 N6, JPN5004013405, 2001, pages 673 - 679, ISSN: 0001246423 * |
P.GRAMATICA: "QSAR STUDY ON THE TROPOSPHERIC DEGRADATION OF ORGANIC COMPOUNDS", CHEMOSPHERE, vol. V38 N6, JPN5004013404, 1999, pages 1371 - 1378, ISSN: 0001246422 * |
Also Published As
Publication number | Publication date |
---|---|
CN1647067A (zh) | 2005-07-27 |
EP1498825A1 (en) | 2005-01-19 |
KR20040111456A (ko) | 2004-12-31 |
CA2481485A1 (en) | 2003-10-16 |
US20050159896A1 (en) | 2005-07-21 |
WO2003085548A1 (fr) | 2003-10-16 |
AU2003220998A1 (en) | 2003-10-20 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
JPWO2003085548A1 (ja) | データ解析装置および方法 | |
JP6888123B2 (ja) | 深層畳み込みニューラルネットワークを事前訓練するための深層学習ベースの技術 | |
US20020095260A1 (en) | Methods for efficiently mining broad data sets for biological markers | |
JP2008511058A (ja) | コンピュータシステムを用いるデータ品質および/または部分異数染色体の決定 | |
JP2023507252A (ja) | パッチ畳み込みニューラルネットワークを用いる癌分類 | |
JP2021505977A (ja) | 体細胞突然変異のクローン性を決定するための方法及びシステム | |
JP2003021630A (ja) | 臨床診断サービスを提供するための方法 | |
JP2005524124A (ja) | システムの診断構成要素を識別するための方法および装置 | |
Win et al. | Cancer recurrence prediction using machine learning | |
Jiang et al. | MHAttnSurv: Multi-head attention for survival prediction using whole-slide pathology images | |
Al-Yousef et al. | A novel computational approach for biomarker detection for gene expression-based computer-aided diagnostic systems for breast cancer | |
Nayak et al. | Deep learning approaches for high dimension cancer microarray data feature prediction: A review | |
Sha et al. | Feature selection for polygenic risk scores using genetic algorithm and network science | |
Jaber et al. | The importance of data classification using machine learning methods in microarray data | |
Abdullahi et al. | Pretrained convolutional neural networks for cancer genome classification | |
US20240047081A1 (en) | Designing Chemical or Genetic Perturbations using Artificial Intelligence | |
Chebli et al. | Unlocking the Potential of DNA Microarray for Accurate Cancer Diagnosis with Deep Learning | |
Praveena et al. | The Potential Uses of Data Science and Deep Learning Techniques in Mining Biological Data: A Comprehensive Analysis | |
Gulla | An integrated systems biology approach to investigate transcriptomic data of thyroid carcinoma | |
Hardin | Microarray data from a statistician’s point of view | |
JP2010514001A (ja) | 特徴の順位付け | |
Kaulgud et al. | A review on detection of prostate cancer techniques | |
WO2023129621A1 (en) | Rare variant polygenic risk scores | |
Elyasigomari | Analysis of microarray and next generation sequencing data for classification and biomarker discovery in relation to complex diseases | |
Murthy et al. | Complexity-Reduced Tumor Classification System using Microarray Gene Expression Dataset |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
A621 | Written request for application examination |
Free format text: JAPANESE INTERMEDIATE CODE: A621 Effective date: 20060328 |
|
RD03 | Notification of appointment of power of attorney |
Free format text: JAPANESE INTERMEDIATE CODE: A7423 Effective date: 20080124 |
|
A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20090210 |
|
A02 | Decision of refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A02 Effective date: 20090616 |