JP2006209276A - 生体試料分類システム及びその方法 - Google Patents

Abstract

【課題】本発明の課題は各生体試料から得た複数時点における遺伝子発現データを用いて、生体試料を分類する方法及び装置を提供するものである。
【解決手段】生体試料についての時系列からなる特徴量をまとめてパターンとして扱うことで、生体試料の特徴や生体反応の特徴をよりよく反映した生体試料の分類方法を実現する。
【選択図】図１

Description

本発明は、生体試料の分類システム及びその方法に関するものである。

ＤＮＡチップ（ＤＮＡマイクロアレイ）は数万種に上る遺伝子の発現量を一度に測定できる手法として生物学、医学薬学等の様々な分野で活用されている。しかし、網羅的な解析ができる反面、１解析あたり数万の数値データが得られるため、人の手だけでその結果を把握し解釈することが困難であり、計算機を用いて遺伝子発現データを分析し可視化するツールが広く用いられている。このようなツールの代表的なものに最も類似性の高いデータ同士を結合して１つのデータ（クラスタ）を作る手順を繰り返し行うことで階層型のクラスタを生成する階層型クラスタリング方法がある（非特許文献１）。また、予めユーザによってクラスタの数ｋを指定して（１）ｋ個のクラスタの中心を適当な方法で仮に決める、（２）どの中心に最も類似性が高いかを基準に全データをｋ個のクラスタに振り分ける、（３）振り分けられたデータに基づいて各クラスタの中心を計算しなおす、（４）（２）〜（３）の手順を中心が収束するまで繰り返す、という手順でｋ個のクラスタを生成するＫ−Ｍｅａｎｓクラスタリング等が知られている（非特許文献２）。他に、クラスタリング手法の改良を試みた例としては、ニューラルネットワークを用いて遺伝子をその時系列発現パターンに基づいてクラスタリングする方法がある（特許文献１）。

特開2002-175306

M.B. Eisen et al.（Proc. Natl. Acad. Sci. USA Vol.95, pp, 14863-14868, December 1998） Soukas, A. et al.（Genes & Development vol.14, pp, 963-980, 2000）

生体試料から得られる遺伝子発現量及び、臨床検査値等その他の生体試料の特徴量の性質として、非常に誤差が大きいことが挙げられる。この誤差は測定誤差だけでなく、生体試料にもともと存在する個体差等の実験誤差が大きな部分を占めている。

また、従来の分類手法の解析対象は、生体試料から得られた１実験条件のみでの遺伝子発現データであるか、または複数条件での遺伝子発現データを何の構造もなく並列に並べたデータの集合に限られていた。測定され得る特徴量に比べてその背景にある生理的現象は遥かに複雑であり、生体試料を特徴付けて背景にある生理現象を推し量るためにはある１時点における特徴量だけでは不十分である。例えば、日常的に測定される健康状態の指標である体温や血圧を例にとっても、ある１時点における測定ではあまり意味が無く、平常時の値や過去数日間の値の推移を持って健康状態が推定されている。
以上のような性質を有する生体試料を対象とした解析においては、様々な条件の変化に伴う特徴量の変動パターンを捕らえる必要がある。

生体試料間で変化させる条件の例としては、サンプリングの時点、化合物の投与量、処理や処置の時間等がある。
例えば、「医薬品の遺伝毒性試験に関するガイドライン」（厚生省医薬安全局、平成１１年）の中で「げっ歯類を用いる小核試験」においては、化合物投与の用量とサンプリングの時点を共に複数設定することが推奨されており、結果を解釈する際にも用量相関性関係に留意するよう指摘されている。
しかし、従来の分類手法はこれらの系列を持ったデータの系列を無視して別個のデータとして扱うか、系列の順序を無視して並列な集合として扱うことしかできなかった。

複数の生体試料についてのＧ個の特徴量からなるＴ時点の時系列データの入力に対して、各特徴量の時系列パターンを考慮した試料間の類似度を算出する方法を提供する。
図１は、本発明の構成を表している。
構成の全体は、クラスタリングに必要な情報を入力しクラスタリング結果を表示する、入出力部と、入出力部から送られた条件に基づいてクラスタリングの計算を行い入出力部へ計算結果を返すクラスタリング部からなっている。

入出力部は、情報入力部、条件指定部、表示部を有している。
情報入力部は、複数試料の特徴量データ、例えばＤＮＡチップの発現データを入力するための手段である。条件指定部は前記特徴量データのうちクラスタリングに用いる試料のリストと特徴量項目、時系列特徴量の類似度、特徴量データ変換条件、クラスタリングアルゴリズムといったクラスタリングに必要な条件を指定するための手段である。時系列特徴量の類似度については後で述べる。また、特徴量データ変換条件とは入力された特徴量データの数値をそのまま用いる、対数変換を行う、各試料の中央値または平均値で割り算する、各遺伝子の中央値または平均値で割り算する、何らかの処置を施した試料の特徴量データを未処置の試料の特徴量データで割り算するといった数値変換処理の条件を指す。クラスタリングアルゴリズムには、例えば階層型クラスタリングやＫ−Ｍｅａｎｓクラスタリング等を指定できる。情報入力部と条件指定部で入力された、特徴量データ及び指定された条件はクラスタリング部へ送られる。

次にクラスタリング部の内部の構成について述べる。
入出力部より送られた特徴量データ及び指定された条件は、まずクラスタリング前処理部に送られる。そこでは、前記特徴量データから条件指定部で指定された特徴量データ変換条件に従い特徴量データを変換し、指定された試料と遺伝子に該当する部分の発現データを抜き出してクラスタリング制御部へ送るといった処理がなされる。

クラスタリング制御部は、試料間類似度行列生成部とクラスタ生成部を有している。試料間類似度行列生成部は、適宜類似度算出部を呼び出しながらクラスタリングに必要な試料間類似度行列を生成するための手段である。また、クラスタ生成部は、クラスタリングアルゴリズム、例えば階層型クラスタリングによって試料間類似度行列に基づいてクラスタを生成する手段である。

類似度算出部は、特徴量選択部、類似度ベクトル算出部、類似度ベクトル変換部を有しており、試料間類似度行列生成部が指定した複数試料についての試料間類似度を、条件指定部で指定された類似度を用いて算出する。特徴量選択部は、試料間類似度行列生成部から送られてきた複数の生体試料の時系列特徴量データから、特徴量ベクトルを順次選択する。類似度ベクトル算出部は、特徴量ベクトルから類似度ベクトルを算出する。類似度ベクトル変換部は、類似度ベクトルをスカラー量で表される類似度に変換する。類似度算出部で行われる計算の詳細については後で述べる。図２は、試料間類似度行列生成部で行われる類似度行列算出のフローを表している。クラスタリング前処理部より送られてきた複数試料の特徴量データから、２つの試料の特徴量データを抜き出して類似度算出部に送り、類似度算出部から試料間類似度を得る。このこの動作を必要な２試料の組合せ全てについて繰り返し、試料間類似度行列を生成する。生成した試料間類似度行列はクラスタ生成部へ送られる。

図３は、類似度算出部で行われる特徴量データ間の類似度算出フローを表している。このフローは図２で示した２試料間の時系列特徴量類似度を算出する方法に相当する。
以下に、図３中の２つの生体試料間の類似度を算出する手順を、図４〜６を用いながら示す。

（１）試料間類似度行列生成部から送られてきた２つの生体試料の時系列特徴量データから２つの生体試料C_１、C_２に共通の特徴量ｇ_ｋを選択し（図４）、２つの時系列特徴量ベクトルｘ_１ｋ、ｘ_２ｋを得る（図５）。

（２）２つの生体試料におけるこの１つの特徴量ｇ_ｋを時系列の順序に並べた２つのベクトル同士の類似度ｓ_ｋ（ｋ＝１~Ｇ、ただしＧは特徴量の個数）を、２つの時系列データｘ_１ｋ、ｘ_２ｋの関数として求める（図５）。この関数の例としては、ピアソン相関係数やコサイン係数、スピアマン相関係数等に代表される各種相関係数、ユークリッド平方距離、標準化ユークリッド距離、ミンコフスキー距離、マハラノビス汎距離を用いてよい。または、特徴量の平均、中央値、最大値、最小値、最大値と最小値の差分、最大値または最小値を取る時点、特定の２時点間の差分または傾き、回帰直線の傾きまたは切片、回帰曲線の平均的な傾きまたは曲率、極値の大きさ、極値を取る時点、極値の数、傾きが初めて閾値を超えるまたは下回る時点、各種の統計的検定によって有意となる時点の数または初めて有意となる時点または初めて有意でなくなる時点等を特徴量の要約値として算出し、２つの特徴量ベクトルについての該要約値の差異を基準に類似度を与える関数を用いてもよい。特徴量の要約値として最大値を取る時点を用いて、特徴量ｇの試料C_ａと試料C_ｂの最大値を取る時点をそれぞれｔ_ａ，ｔ_ｂとし、時点間の距離をｄ_ｔ（ｔ_ａ，ｔ_ｂ）と表すことにすると、類似度ｓを例えば（総時点数）−ｄ_ｔ（ｔ_ａ，ｔ_ｂ）と定義することができる。図６のように総時点数＝４での時、特徴量ｇが図６の(a)のような値をとる場合には試料C_１と試料C_２の最大値を取る時点は共にｔ_３で、類似度ｓ＝４−ｄ_ｔ（ｔ_３，ｔ_３）＝４−０＝４。図６の(b)のような値をとる場合には試料C_１と試料C_２の最大値を取る時点はそれぞれｔ_３、ｔ_２で、類似度ｓ＝４−ｄ_ｔ（ｔ_３，ｔ_２）＝４−１＝３となる。

（３）全ての特徴量について前記（１）（２）の操作を繰り返すことで得られる類似度ｓ_ｋを成分に持つベクトルＳ＝（ｓ_１，・・・，ｓ_Ｇ）を、２つの生体試料の類似度を表す類似度ベクトルとする。

（４）類似度ベクトルの関数としてスカラー値で表現される２生体試料間の類似度を求める。この関数として例えば、ベクトルの要素の総和、m乗和（ｍは０でない）、または要素の最大値、最小値、中央値、または特定の範囲に含まれる要素の総和、m乗和（ｍは０でない）、または各要素を入力とするロジスティック関数の出力の総和、m乗和（ｍは０でない）等を用いてよい。

前記２生体試料間の類似度、または類似度の逆数や類似度の符号を変えた数等を用いて算出した非類似度を、公知の教師信号なしまたは教師信号ありのクラスタリングに適用することで、複数生体試料を各特徴量の時系列パターンを考慮して分類することが可能となる。

本発明は時系列データ以外にも、用量、線量など順序が規定できるデータセット一般に適用可能である。
また、生体試料の特徴量としては、遺伝子の発現量以外にも、タンパク質量や修飾済みタンパク質の量または非修飾タンパク質に対する存在比、特定ゲノム領域におけるゲノムDNAのコピー数、特定ゲノム領域におけるゲノムDNAまたはヌクレオソームの修飾の程度、代謝産物量、体重、身長、年齢、血中イオン濃度、血中の逸脱酵素濃度、各種血球数、その他の臨床検査値等、多数の生物学的特徴量の時系列データを扱うことができる。

各特徴量の時系列パターンに基づいた類似度は、より生理的現象に即した尺度であると考えられる。そのため、該類似度を用いることで生理的現象をより正確に捉えた生物学的分析が可能となる。

化合物の毒性を評価するため、ラットにＣ_１からＣ_５の５種類の化合物を投与し、時系列的に５点（ｔ＝｛０，６，１２，２４，４８｝，単位：時間、ｔ＝０は投与前を意味する。）で屠殺して肝臓のｔｏｔａｌＲＮＡを抽出した。オリゴ（ｄＴ）２４プライマーをアニールさせ、ｔｏｔａｌＲＮＡを鋳型とするｃＤＮＡを合成した後、ｔ＝０の試料はＣｙ３で、それ以外の試料はＣｙ５で蛍光標識した。化合物投与後（ｔ＝｛６，１２，２４，４８｝）のＣｙ５標識ｃＤＮＡの１つと投与前（ｔ＝０）のＣｙ３標識ｃＤＮＡとを組み合わせて、それぞれをＤＮＡチップ（日立製作所社製薬物応答解析用ラットＤＮＡチップ）上にかけ、競合ハイブリダイゼーションを６２℃のもと１２時間行った。洗浄後スキャナー（GSI-Lumonics社製ScanArray 5000）により各スポットの蛍光強度を測定し、数値化ソフトウエア（GSI-Lumonics社製QuantArray）を用いて各遺伝子におけるＣｙ５（投与後試料）とＣｙ３（投与前試料）との蛍光強度を求め、Ｃｙ５／Ｃｙ３の蛍光強度の比を、該スポットに対応する遺伝子の発現データとした（表１）。

次にＤＮＡチップにより測定された遺伝子の発現データに基づいた投与化合物の分類を試みた。
化合物ｃ投与後の時点ｔ（ｔ＝｛６，１２，２４，４８｝）における遺伝子ｇの発現データをＸ_ｇｔｃとする。

また、化合物ｃ投与後の遺伝子ｇの時系列発現データをＸ_ｇｃとする。（ｇ＝１~Ｇ、ただしＧは遺伝子の数でこの例ではＧ＝３６。）
２種の化合物Ｃ_ｉとＣ_ｊ（ｉ，ｊ＝｛１，２，３，４，５｝）をそれぞれ投与したラットの発現データを用いて、化合物Ｃ_ｉとＣ_ｊの類似度を、ｓ_ｉｊ＝Σ_ｇ（ｆ_ｓ（Ｘ_ｇｉ，Ｘ_ｇｊ））で定義する。

ここではピアソン相関係数を用て、ｆ_ｓ（Ｘ_ｇｉ，Ｘ_ｇｊ）＝（Σ_ｔ（Ｘ_ｇｔｉ−ｍ（Ｘ_ｇｉ）（Ｘ_ｇｔｊ−ｍ（Ｘ_ｇｊ））／（Σ_ｔ（Ｘ_ｇｔｉ−ｍ（Ｘ_ｇｉ））^２・Σ_ｔ（Ｘ_ｇｔｊ−ｍ（Ｘ_ｇｊ））^２）^１／２（ただし、ｍ（Ｘ_ｇｉ）は、Ｘ_{ｇｔｉ（ｔ＝｛６，１２，２４，４８｝）}の平均値）と定義した。ピアソン相関係数やコサイン係数、スピアマン相関係数等に代表される各種相関係数をもちいることで、発現データの値そのものではなく相対的な増減のパターンを反映した類似度とすることができる。

ピアソン相関係数以外の類似度、例えばユークリッド平方距離を用いる場合には、ｆ_ｓ（Ｘ_ｇｉ，Ｘ_ｇｊ）＝（Σ_ｔ（Ｘ_ｇｔｉ−Ｘ_ｇｔｊ）^２）^１／２とする。ユークリッド平方距離は発現データの値自体の大小を反映した類似度であり、仮に時系列で発現量の相対的な変化がない遺伝子でも、高い発現量で変化がないのかより低い発現量で変化がないのかを区別することができる。ユークリッド平方距離の他、標準化ユークリッド距離、ミンコフスキー距離、マハラノビス汎距離等を用いて同様の効果を得ることができる。

また、各試料における特徴量の要約値を求め、要約値の差異を基準とした類似度を用いても良い。例えば時系列データの総時点数ｎ_ｔ、時系列発現データＸ_ｇｃの最大値を取る時点を与える関数ｔ_ｍａｘ（Ｘ_ｇｃ）、時点ｔ_ａとｔ_ｂの距離を与える関数ｄ_ｔ（ｔ_ａ，ｔ_ｂ）とを用いてｆ_ｓ（Ｘ_ｇｉ，Ｘ_ｇｊ）＝ｎ_ｔ−ｄ_ｔ（ｔ_ｍａｘ（Ｘ_ｇｉ）−ｔ_ｍａｘ（Ｘ_ｇｊ））と定義することができる。この類似度は、相関係数やユークリッド平方距離のように全体の傾向を捉えるのではなくより目的に特化した類似度となる。最大値をとる時点の他にも極値を取る時点、極値の数、傾きが初めて閾値を超えるまたは下回る時点、各種の統計的検定によって有意となる時点の数または初めて有意となる時点または初めて有意でなくなる時点等を用いて、同様に目的に特化した類似度を得ることがることができる。

最大値と最小値の差分、最大値または最小値を取る時点、特定の２時点間の差分または傾き、回帰直線の傾きまたは切片、回帰曲線の平均的な傾きまたは曲率、極値の大きさを利用した類似度を用いることで相関係数を用いた場合と同様に時系列における変化の傾向を捉える効果を得ることができると考えられる。
特徴量の要約値として特徴量の平均、中央値、最大値、最小値を用いた場合はユークリッド平方距離と似た効果を得ることができると考えられる。

表２に類似度の計算過程を示した。
化合物C_１とC_２の遺伝子G_１の時系列発現データはそれぞれ、Ｘ_{ｇ＝１，ｃ＝１}＝（１．０７，０．６２，０．６７，１．２１）、Ｘ_{ｇ＝１，ｃ＝２}＝（２．００，１．０６，０．６４，２．２１）で、Ｘ_{ｇ＝１，ｃ＝１}とＸ_{ｇ＝１，ｃ＝２}のピアソン係数は０．９５である。同様に他の遺伝子についてもピアソン相関係数を順次計算し、その総和が５．０１であるので、化合物C_１とC_２の類似度ｓ_１，２＝５．０１となる。

この手順を化合物の組合せについて繰り返すことで、類似度ｓ_ｉｊを成分に持つ、化合物Ｃ_１からＣ_５の類似度行列Ｓを得た（表３）。

この類似度行列Ｓと前出の遺伝子数Gを用いて、ｄ_ｉｊ＝Ｇ−ｓ_ｉｊを成分に持つ、化合物Ｃ_１からＣ_５の非類似度行列Ｄ＝ｄ_ｉｊ定義し、表４を得た。

この非類似度行列Ｄを用いて、化合物Ｃ_１からＣ_５についての階層型クラスタリングを行った。階層型クラスタリングを実施するに当たっては統計解析パッケージＲ（The R Project for Statistical Computingによる。http://www.r-project.org/ を参照されたい）のclusterパッケージに含まれるagnes関数を用いた。agnes関数は凝集法階層型クラスタリングのUPGMA法を実装したものである。agnes関数の詳細についてはKaufman, L. and Rousseeuw, P.J. (1990). Finding Groups in Data: An Introduction to Cluster Analysis. Wiley, New York.を参照されたい。

階層型クラスタリングの結果得られた化合物Ｃ_１からＣ_５の樹状図を図７に示した。各化合物の既知の肝毒性を照会したところ、化合物Ｃ_２とＣ_５は共に肝細胞の壊死を示すことがが報告されており、本手法によりラット肝臓における化合物投与時の生理的反応を反映した化合物の分類がなされていることを確認することができた。

図８は本実施例における装置の構成を示すものである。
発現データベース１とデータ解析コンピュータ２はネットワークで接続されている。発現データベース１には、例えば、予め測定された各化合物を投与した際の時系列遺伝子発現データを格納しておく。データ解析コンピュータ２には、各種のクラスタリング条件を指定するためのインタフェースが、ウェブブラウザ上に表示される。このインタフェースは、例えば、化合物のリストから取捨選択して選択された化合物を投与した際の発現データを解析に用いるためのチェックボックスと、データ変換条件をデータベース中の発現値を対数変換するか否かを選択するラジオボタンと、時系列発現データの類似度の算出方法を実施例１に示したピアソン相関係数の総和とユークリッド距離の総和のいずれかから選択するラジオボタンと、クラスタリングのアルゴリズムを階層型クラスタリングとＫ−Ｍｅａｎｓから選択するラジオボタンを含む。

利用者はデータ解析コンピュータ２を用いて、ウェブブラウザ上のインタフェースを通じてクラスタリング条件を指定し計算の実行を指示する。クラスタリングの計算は実施例１に示された方法に則ってデータ解析コンピュータ２にて行われ、クラスタリングの結果はデータ解析コンピュータ２のディスプレイに表示される。
クラスタリングの結果は、画像ファイルとしてデータ解析コンピュータ２の外部記憶装置に保存したり、紙面に印刷することも可能である。

本発明の実施形態例を示す構成図である。 本発明の実施形態例に係る特徴量時系列データ間の類似度算出方法の一例を示すフローチャートである。 本発明の実施形態例に係る複数生体試料間の類似度行列算出方法の一例を示すフローチャートである。 本発明の実施形態例に係る生体試料特徴量時系列データの構造の一例を示す説明図である。 特徴量時系列データ間の類似度算出方法の一例を示す模式図である。 特徴量時系列データ間の類似度算出方法の一例を示す 階層型クラスタリングの結果得られた化合物Ｃ_１からＣ_５の樹状図を示した。 実施例２における装置の構成図を示す。

符号の説明

１・・・時系列遺伝子発現データを格納する発現データベース
２・・・各種クラスタリング条件を入力し、クラスタリング結果を出力するデータ解析コンピュータ
３・・・出力されたクラスタリング結果。

Claims (8)

1. 複数の生体試料の各々について、複数特徴量の組が順序をもって与えられたデータセットを入力させる入力手段と、
   前記特徴量のデータの変換条件を指定させる条件指定手段と、
   前記データセットから少なくとも２つの生体試料の組み合わせを選択するクラスタリング制御手段と、
   選択された前記生体試料間の類似度を算出する類似度算出手段と、
   算出された前記生体試料間の類似度に基づきクラスタを生成するクラスタ生成手段と、
   生成されたクラスタを表示する表示手段とを備えることを特徴とする生体試料分類システム。
2. 請求項１に記載の生体試料分類システムにおいて、前記類似度算出手段は、前記クラスタリング制御手段により選択された２つの生体試料の各々について、特徴量を順序をもって並べたベクトルから、前記複数の特徴量の各特徴量における前記２つの生体試料間の類似度を算出する特徴量パターン類似度算出手段と、算出された前記各特徴量における２つの生体試料間の類似度から、前記２つの生体試料間の類似度を算出する生体試料間類似度算出手段を有することを特徴とする生体試料分類システム。
3. 請求項２に記載の生体試料分類システムにおいて、前記特徴量パターン類似度算出手段は、前記２つの生体試料についての特徴量ベクトル間でピアソン相関係数を用いて類似度を算出することを特徴とする生体試料分類システム。
4. 請求項２に記載の生体試料分類システムにおいて、前記生体試料間類似度算出手段は、特徴量ｇについての特徴量パターン類似度ｓ_ｇを用いて、ｓ＝Σ_ｇ（ｓ_ｇ）^ｍ（ｍは０でない）で表される関数によって類似度ｓを算出することを特徴とする生体試料分類システム。
5. 複数の生体試料の各々について、複数特徴量の組が順序をもって与えられたデータセットを入力させる工程と、
   前記特徴量のデータの変換条件を指定させる工程と、
   前記データセットから少なくとも２つの生体試料の組み合わせを選択させる工程と、
   選択された前記生体試料間の類似度を算出する類似度算出工程と、
   算出された前記生体試料間の類似度に基づきクラスタを生成するクラスタ生成工程と、
   生成されたクラスタを表示する表示する工程とを有すること特徴とする生体試料分類方法。
6. 請求項５に記載の生体試料分類方法において、前記類似度算出工程は、前記選択された２つの生体試料の各々について、特徴量を順序をもって並べたベクトルから、前記複数の特徴量の各特徴量における前記２つの生体試料間の類似度を算出する工程と、算出された前記各特徴量における２つの生体試料間の類似度から、前記２つの生体試料間の類似度を算出する工程とを有することを特徴とする生体試料分類方法。
7. 請求項６に記載の生体試料分類方法において、前記各特徴量における２つの生体試料間の類似度は、前記２つの生体試料についての特徴量ベクトル間でピアソン相関係数を用いて類似度を算出することを特徴とする生体試料分類システム。
8. 請求項６に記載の生体試料分類方法において、前記２つの生体試料間の類似度は、特徴量ｇについての特徴量パターン類似度ｓ_ｇを用いて、ｓ＝Σ_ｇ（ｓ_ｇ）^ｍ（ｍは０でない）で表される関数によって類似度ｓを算出することを特徴とする生体試料分類システム。

Images