JP2005301789A - クラスタ解析装置、クラスタ解析方法、及びクラスタ解析プログラム - Google Patents

Abstract

【課題】 信頼性の高いクラスタリングを実現し、生命情報に関するデータの解析を行なう。
【解決手段】 ＣＰＵ３は、コリレーション族のカーネル関数を用いて、Ｍ次元のパターンベクトルに対し、確率密度関数を推定し、推定した確率密度関数を基に、密度分布モード探索アルゴリズムを用いてパターンベクトルのクラスタリングを行なう。以上を複数スケールのカーネルで行い、同時に図示する。
【選択図】 図１

Description

本発明は、生命情報に関する複数のデータをクラスタリングするクラスタ解析装置、クラスタ解析方法、クラスタ解析プログラムに関するものである。

遺伝子発現解析では、サンプル毎に複数の遺伝子の発現量が表された遺伝子発現行列に対し、樹状図を作成してクラスタリングを行なう階層化クラスタリングが広く用いられている。図１７は、階層化クラスタリングによる遺伝子発現解析を説明するための図面である。Ｈは、遺伝子発現行列を示し、各列は遺伝子の種類を示し、各行はサンプルの種類を示している。行列Ｈ内にある矩形状の領域は、行列Ｈの１成分を示し、あるサンプルにおける遺伝子の発現量が色を用いて表示されている。そして、階層化クラスタリングは、この行列Ｈにおいて、枝Ｖを用いて類似する遺伝子同士を繋ぎ合わせていき、トーナメント表のような樹状図を作成し、この樹状図を基に、クラスタリングを行っている。

しかしながら、階層化クラスタリングは、データ相互の関係によってボトムアップ的に樹状図を作成するものであるため、母集団中のデータが例えば１個出入りするだけで、全く異なる樹状図が作成されることが頻繁に発生する。そのため、信頼性の高いクラスタリングができないという問題があった。

本発明の目的は、母集団中に小数のデータの出入りがある場合であっても、クラスタリング結果が大きく変動せず、信頼性の高いクラスタリング結果が得られ、高精度に生命情報に関するデータをクラスタリングすることができるクラスタ解析装置、クラスタ解析方法、及びクラスタ解析プログラムを提供することである。

本発明にかかるクラスタ解析装置は、生命情報に関する複数のデータをクラスタリングするクラスタ解析装置であって、カーネル関数を用い、前記データの確率密度分布を推定する推定手段と、推定された確率密度分布を基に、密度分布モード探索アルゴリズムを用いて前記データのクラスタリングを行なうクラスタリング手段とを備えることを特徴とする。

また、前記生命情報に関するデータが、遺伝子発現プロファイルデータであることが好ましい。

また、上記クラスタ解析装置において、前記カーネル関数のスケールを設定するスケール設定手段を更に備え、前記推定手段は、前記スケール毎に、前記データの確率密度分布を推定し、前記クラスタリング手段は、前記スケール毎に前記データのクラスタリングを行なうことが好ましい。
前記クラスタリング手段によるクラスタリング結果を基に、各データが、スケール毎に、どのクラスタに属するかをクラスタ別に色分けして示す画像を作成する画像作成手段を更に備え、
また、上記クラスタ解析装置において、前記画像作成手段は、各クラスタがスケール全域に亘って連続して表示されるように前記画像を作成することが好ましい。

また、前記密度分布モード探索アルゴリズムは、ミーンシフトアルゴリズムであることが好ましい。

また、前記カーネル関数は、式（Ａ）で表され、前記密度分布モードアルゴリズムは、各データの座標ベクトルを始点ｖ₁として、式（Ｂ）の演算を繰り返し実行し、ｖ_i+1のｖ_iに対する変化が、ほぼなくなった点を収束点モードとし、前記収束点モードが同一となる始点ｖ₁を１つのクラスタとすることが好ましい。

但し、Ｋ（ｘ_n，ｘ）は、カーネル関数を示し、ｇはｃｏｓθ_nの増大に応じて出力値が増大する任意の関数を示し、ｇ´はｇの微分を示し、ｕ_n＝ｘ_n／｜ｘ_n｜を示し、ｖ₁＝ｘ／｜ｘ｜を示し、ｘ_nは生命情報に関するデータの座標ベクトルを示し、ｘは生命情報に関するデータｘ_nが配置される空間の任意の点の座標ベクトルを示し、θ_nはｘ_nとｘとがなす角度を示し、｜ｘ｜はｘの大きさを示し、｜ｘ_n｜はｘ_nの大きさを示し、ｎは１以上の整数を示す。

本発明にかかるクラスタ解析方法は、生命情報に関する複数のデータをコンピュータによりクラスタリングするクラスタ解析方法であって、前記コンピュータが、カーネル関数を用い、前記データの確率密度分布を推定するステップと、前記コンピュータが、推定された確率密度分布を基に、密度分布モード探索アルゴリズムを用いて、前記データのクラスタリングを行なうステップとを備えることを特徴とする。

本発明にかかるクラスタ解析プログラムは、生命情報に関する複数のデータをクラスタリングするクラスタ解析プログラムであって、カーネル関数を用い、前記データの確率密度分布を推定する推定手段と、推定された確率密度分布を基に、密度分布モード探索アルゴリズムを用いて、前記データのクラスタリングを行なうクラスタリング手段としてコンピュータを機能させることを特徴とする。

請求項１、７、及び８記載の発明によれば、カーネル関数を用いて、生命情報に関する複数のデータに対する確率密度分布が推定され、推定された確率密度分布のモード（極大値）を求める密度分布モード探索アルゴリズムによって、各データはクラスタリングされている。この確率密度分布は、カーネル関数を用いて算出されており、母集団中に小数のデータの出入りがあっても大きく変動するものではない。そして、この確率密度分布のモードを求めることでクラスタリングが行なわれているため、母集団中に小数のデータの出入りがある場合であっても、クラスタリング結果は大きく変動せず、信頼性の高いクラスタリングを行なうことができる。

請求項２記載の発明によれば、生命情報に関するデータを遺伝子発現プロファイルデータとしたため、信頼性の高い遺伝子発現解析を行なうことができる。

請求項３記載の発明によれば、カーネル関数のスケールが設定され、設定されたスケール毎にデータのクラスタリングが行なわれているため、スケール毎のクラスタリング結果を一度に得ることができる。カーネル関数を用いてクラスタリングを行なった場合、スケールに応じてクラスタの数が変動する。遺伝子発現解析の分野では、スケールの値を変更し、解析対象となる遺伝子の種類等に応じて好ましいクラスタリング結果を得たいという要望がある。本発明では、スケール毎のクラスタリング結果が一度に算出されているため、ユーザは、この結果を基に、所望するスケールの値を容易に決定することができる。

請求項４記載の発明によれば、同一のクラスタが全スケールに亘って、連続して表示されるような、スケール毎に各クラスタが色分けして表示された画像が作成されるため、この画面を見たユーザは、スケールの決定をより直感的に行なうことができる。

請求項５記載の発明によれば、ミーンシフトアルゴリズムによりクラスタリングが行なわれているため、旧来の勾配法アルゴリズムに比較してモード検出の計算効率を向上させることができる。

請求項６記載の発明によれば、各データｘ_nが配置される空間（データｘの要素数がＭの場合、Ｍ次元の空間）において、原点を中心とし、半径を１とするＭ次元の球の表面上に、データｘ_nが投影され、投影されたデータｖ₁と、球表面上のある点ｕ_nとの原点から見た角度の大きさのみに基づいて、データｘ_nの確率密度分布が球表面上において推定される。

そして、投影点ｖ₁に対し、式（Ｂ）の演算が繰り返し実行されることにより、ｖ₁が、球表面上の確率密度分布のモードに向けて移動されていき、ｖ_i+1が、ｖ_iに対して、ほぼ変化しなくなった点が収束点とされ、収束点を同一とするｖ₁が１つのクラスタとされ、クラスタリングが行なわれている。そのため、公知のユークリッド距離に基づくカーネル関数を用いてクラスタリングを行なう場合に比べ、高次元のデータに対するクラスタリングの信頼性を高めることができる。

以下、本発明の一実施の形態によるクラスタ解析装置について図面を参照しながら説明する。図１は、本発明の一実施の形態によるクラスタ解析装置の構成を示すブロック図である。ここでは、例として、マイクロアレイなどによる遺伝子発現プロファイルデータに対する解析について説明する。

図１に示すクラスタ解析装置は、通常のコンピュータ等から構成され、入力装置１、ＲＯＭ（リードオンリメモリ）２、ＣＰＵ（中央演算処理装置）３、ＲＡＭ（ランダムアクセスメモリ）４、外部記憶装置５、表示装置６及び記録媒体駆動装置７を備える。各ブロックは内部のバスに接続され、このバスを介して種々のデータ等が入出力され、ＣＰＵ３の制御の下、種々の処理が実行される。

入力装置１は、キーボード、マウス等から構成され、操作者が種々のデータ及び操作指令等を入力するために使用される。例えば、入力装置１は、操作者が入力した遺伝子発現プロファイルデータを取得し、ＣＰＵ３の制御の下、ＲＡＭ４又は外部記憶装置５に出力する。

ＲＯＭ２には、ＢＩＯＳ（Basic Input/Output System）等のシステムプログラム等が記憶される。外部記憶装置５は、ハードディスクドライブ等から構成され、外部記憶装置５には所定のＯＳ（Operating System）及びクラスタ解析プログラム等が記憶される。ＣＰＵ３は、外部記憶装置５からクラスタ解析プログラム等を読み出し、各ブロックの動作を制御する。ＲＡＭ４は、ＣＰＵ３の作業領域等として用いられる。

表示装置６は、液晶表示装置等から構成され、ＣＰＵ３の制御の下に種々の操作画面及び推定結果画面等を表示する。また、必要に応じて推定結果等を印字する印刷装置を付加してもよい。

記録媒体駆動装置７は、ＣＤ−ＲＯＭドライブ、フレキシブルディスクドライブ等から構成される。なお、クラスタ解析プログラムを、ＣＤ−ＲＯＭ、フレキシブルディスク等のコンピュータ読み取り可能な記録媒体８に記録し、記録媒体駆動装置７により記録媒体８からクラスタ解析プログラムを読み出して外部記憶装置５にインストールして実行するようにしてもよい。また、図１に示すクラスタ解析装置が通信装置等を備え、クラスタ解析プログラムが所定のネットワークを介して図１に示すクラスタ解析装置に接続された他のコンピュータ等に記憶されている場合、当該コンピュータからネットワークを介してクラスタ解析プログラムをダウンロードして実行するようにしてもよい。

本実施の形態では、ＣＰＵ３がクラスタリング手段、推定手段、スケール設定手段、画像作成手段の一例に相当し、遺伝子発現量を示すパターンベクトル及びサンプルが生命情報に関するデータに相当する。

図２は、クラスタ解析装置の動作を示すフローチャートである。まず、ステップＳ１において、ＣＰＵ３は、遺伝子発現プロファイルデータを取得する。図３は、遺伝子発現プロファイルデータのデータ構造を示す模式図である。遺伝子発現プロファイルデータは、（Ｎ行×Ｍ列）行列Ｈの形で得られる。この行列Ｈを発現量行列と呼ぶ。Ｎは遺伝子の種類の個数を示し、Ｍはサンプル数を示す。行列Ｈの第（ｎ，ｊ）成分ｈ_njは、第ｎ遺伝子の第ｊサンプルにおける発現量を表す数値である。ここで、第ｎ横ベクトルを、パターンベクトルｘ_nと呼ぶ。パターンベクトルｘ_nは、１〜Ｍのサンプルに対するＭ個の発現量を表すベクトルである。そのため、１つのパターンベクトルｘ_nは、Ｍ次元空間内において、１つの点として表すことができる。本クラスタ解析装置は、Ｍ次元空間内にＮ個のパターンベクトルｘ_nを配置し、これらのパターンベクトルｘ_nをクラスタリングする。

ステップＳ２において、ＣＰＵ３は、入力された遺伝子発現プロファイルデータからｎ個のパターンベクトルｘ_nを生成し、Ｍ次元空間に配置する。図４は、２次元空間内に設定されたパターンベクトルｘ_nの確率密度分布を示した図面であり、横軸はサンプル１の発現量を示し、縦軸はサンプル２の発現量を示している。図４に示された複数の点Ｐ１，Ｐ２が、パターンベクトルｘ_nである。例えば、サンプル１の発現量が３、サンプル２の発現量が５のパターンベクトルｘ_nがあったとすると、このパターンベクトルｘ_nは、横軸が３，縦軸が２の座標にプロットされる。

ステップＳ３において、ＣＰＵ３は、カーネル関数のスケールに初期値を設定する。そして、スケールが設定されたカーネル関数を用いて、Ｍ次元空間に配置されたパターンベクトルｘ_nの確率密度分布を推定する（ステップＳ４）。カーネル関数をＫ（ｘ_n，ｘ）とすると、確率密度分布を示す確率密度関数ｆ（ｘ）は、式（１）で表される。ｘ_nはパターンベクトルを示し、Ｎはパターンベクトルｘ_nの数を示し、ｘは空間内の任意の点を示す。また、Ｋ（ｘ_n，ｘ）≧０、∫Ｋ（ｘ_n，ｘ）ｄｘ_n＝１である。

カーネル関数Ｋ（ｘ_n，ｘ）としては、ＲＢＦ族のカーネル関数が、公知である。式（２）はＲＢＦ族のカーネル関数の一般形を示している。ｇはｘとｘ_nとのユークリッド距離が小さいほど出力値が大きくなる任意の関数である。

ここで、ＲＢＦ（radial basis function:動径基底関数）族の１次元空間を例に挙げ、スケールと確率密度分布との関係について説明する。図５は、カーネル関数と確率密度分布との関係を説明するための図面であり、（ａ）は半径がｒ１のカーネル関数Ｋ（ｘ_n，ｘ）を示し、（ｂ）は半径がｒ２（ｒ２＜ｒ１）のカーネル関数Ｋ（ｘ_n，ｘ）を示し、（ｂ）は（ａ）のカーネル関数を用いて算出された確率密度関数ｆ（ｘ）を示し、（ｄ）は（ｃ）のカーネル関数を用いて算出された確率密度関数ｆ（ｘ）を示している。（ｂ），（ｄ）のｘ軸上に描かれた複数の単線ＳＬはパターンベクトルｘ_nが存在する座標を示している。このカーネル関数の半径がスケールである。

（ａ）及び（ｂ）において、ｘ₀はカーネル関数の中心を示し、ｘ_nはｘ軸上に存在するあるパターンベクトルｘ_nを示している。ＲＢＦ族のカーネル関数は、中心ｘ₀とパターンベクトルｘ_nとの距離が小さいほど、頂点の高さＫＨが高くなり、逆の場合、頂点の高さＫＨが低くなる。また、ＲＢＦ族のカーネル関数は、ｘ＝ｘ₀に対して対象な釣鐘状の形状を有している。そして、地点ｘにカーネル関数の中心ｘ₀を当てはめて、全パターンベクトルｘ_nに対するｘ軸の確率密度分布を求め、さらに、カーネル関数をｘ軸上にずらしていき、各地点ｘにカーネル関数の中心ｘ₀を当てはめたときのｘ軸の確率密度分布を算出し、算出した各地点ｘのｘ軸の確率密度分布の総和が、確率密度関数ｆ（ｘ）となる。

（ｂ）及び（ｄ）を比較すれば分かるように、カーネル関数のスケールの値が大きい場合、確率密度関数ｆ（ｘ）の形状はなだらかであるが、カーネル関数のスケールの値を小さくすると、確率密度関数ｆ（ｘ）の形状はシャープとなり、その結果、極大値（モード）の数が変動する。つまり、スケールを減少させると、確率密度関数ｆ（ｘ）のモードの数は増大し、一方、スケールを増大させると、確率密度関数ｆ（ｘ）のモードの数は減少する。本実施の形態では、スケールを大きい値から小さい値へ減少させるものとする。したがって、ステップＳ３においては、初期値としては、スケールの変動範囲の最大値が設定される。

なお、式（２）に示すＲＢＦ族のカーネル関数のより具体的な関数として、式（３）で示すガウシアン（Gaussian）カーネル関数が知られており、図４で示す確率密度分布は、式（３）を用いて算出されたものである。式（３）において、ｈはスケールを示している。図４は、ｈ＝０．５のときの確率密度分布である。

図４に示された複数の閉曲線は、確率密度分布の等高線ＣＬを示しており、内側にある閉曲線ほど確率密度が高い。また、Ｍ１〜Ｍ４は確率密度分布のモードを示している。

ステップＳ５において、ＣＰＵ３は、密度分布モード探索アルゴリズムの一例であるミーンシフトアルゴリズムを用いてクラスタリングを行なう。ミーンシフトアルゴリズムは、空間上の任意の点を始点として確率密度分布の勾配に沿って移動していくような点列を生成し、いずれかのモードに到達させ、同一のモードに到達した点を１つのクラスタとして纏めることによりクラスタリングを行なう手法である。図４に示す点Ｐ１を確率密度分布の勾配に沿って移動させていくと、モードＭ１に到達するため、点Ｐ１はモードＭ１に対応するクラスタである。また、点Ｐ２を確率密度分布の勾配に沿って移動させていくと、モードＭ２に到達するため、点Ｐ２はモードＭ２に対応するクラスタとなる。図４に示す点線によって区分けされた４つの領域Ｄ１〜Ｄ４は、それぞれモードＭ１〜Ｍ４に対応する領域を示しており、領域Ｄ１〜Ｄ４に属する各点に対してミーンシフトアルゴリズムを適用した場合、これら各点はそれぞれモードＭ１〜Ｍ４に到達する。そして、パターンベクトルｘ_nに対してミーンシフトアルゴリズムが適用され、どのモードに到達したかによって、パターンベクトルｘ_nはクラスタリングされる。

ステップＳ６において、カーネル関数に設定されたスケールが設定範囲の最終値でない場合（ステップＳ６でＮＯ）、ＣＰＵ３は、スケールの再設定を行い（ステップＳ７）、処理をステップＳ４に戻す。一方、スケールが設定範囲の最終値である場合（ステップＳ６でＹＥＳ）、処理が終了される。ＣＰＵ３は、カーネル関数として、式（３）に示すガウシアンカーネル関数を用いた場合、例えばスケールｈの値を初期値＝１、最終値＝０として０．１の刻み幅で変化させる。なお、後述するコリレーション族のカーネル関数であるｖＭＦカーネル関数を用いた場合は、スケールの増大に伴って、確率密度関数ｆ（ｘ）の形状はシャープになるため、初期値から最終値に至るまで、スケールの値を増大させる。この場合、スケールの刻み幅をスケールの値に応じて適宜変更してもよい。

図６は、ガウシアンカーネル関数を用いた場合の、２次元のパターンベクトルｘ_nに対する確率密度分布を示しており、（ａ）はスケール１、（ｂ）はスケール０．７、（ｃ）はスケール０．５、（ｄ）はスケール０．２の場合を示している。

図６（ａ）の場合、モードの数は１個であり、（ｂ）の場合、モードの数は２個となり、（ｃ）の場合、モードの数は４個となり、（ｄ）に至ってはモード数が９個となっている。そのため、スケールの値が減少するにつれて、等高線ＣＬの数が増大し、確率密度分布の勾配が急峻となり、それに応じてモードの数が増大していることが分かる。

図２に示すステップＳ４では、ＲＢＦ族のカーネル関数を用いたものを示したが、これに変えて本発明者によって開発されたコリレーション（correlation）族のカーネル関数を用いてもよい。式（４）はコリレーション族のカーネル関数の一般形を示している。

式（４）に示すθ_nは、パターンベクトルｘ_nが配置された空間のある地点ｘの座標ベクトルと、パターンベクトルｘ_nとのなす角度を示している。ｇはｃｏｓθ_nの値が大きくなるにつれて（θ_nの値が小さいほど）値が小さくなる任意の関数を示している。
式（４）に示すように、ｇの引数はｃｏｓθ_nであり、ｃｏｓθ_nはｘ_n・ｘ／｜ｘ_n｜｜ｘ｜で表され、ｘ_n／｜ｘ_n｜と、ｘ／｜ｘ｜とはそれぞれ単位ベクトルであるため、ｇは、パターンベクトルｘ_nと地点ｘとのなす角度のみに依存する関数となる。その結果、パターンベクトルｘ_nと地点ｘとは、原点Ｏを中心とするＭ次元の球の表面上に配置された点と考えることができる。

２次元空間においては、パターンベクトルｘ_nと地点ｘとは、図７に示すように、原点をＯを中心とする単位円Ｃの円周上に配置された点と考えることができる。空間内において、地点ｘの近傍に位置するパターンベクトルｘ₁が地点ｘとなす角度θ₁は、地点ｘの遠方に位置するパターンベクトルｘ₂が地点ｘとなす角度θ₂よりも小くなることは明白である。また、ｇは上述したように、θ_nの値が小さいほど、値が大きくなる関数である。そのため、ｇの値は、パターンベクトルｘ_nの近傍の地点ｘの方が、パターンベクトルｘ_nの遠方の地点ｘよりも大きくなる。そして、式（１）に示すように地点ｘの確率密度は、パターンベクトルｘ_n毎に算出されたｇの総和となるため、パターンベクトルｘ_nが密集する場所の地点ｘほど、確率密度の値が大きくなる。

このように、コリレーション族のカーネル関数を用いた確率密度分布推定は、距離に依存せず角度のみに依存するものであるため、ユークリッド距離に基づくＲＢＦ族のカーネル関数を用いた確率密度分布推定と比べ、高次元空間におけるクラスタリングの信頼性を高めることができる。

さらに、本発明者は、式（４）に示すカーネル関数の一般形に対するより具体的な関数として、式（５）に示すカーネル関数を開発した。

以下、この関数をｖＭＦ(von Mises-Fisher)関数と呼ぶ。κはスケールを示している。ｖＭＦ関数の本質的部分はｃｏｓθを含むｅｘｐ（κｕ）にあり、Ｚ（Ｍ，κ）^-1は、ｅｘｐ（κｕ）を正規化するための係数である。Ｚ（Ｍ，κ）は、式（６）によって表される。ここで、Ｂνは、ν（＝Ｍ／２−１）次の変形ベッセル関数であり、式（７）によって表される。ｖＭＦカーネル関数は、方向ベクトルデータを扱う際にしばしば用いられるｖＭＦ分布に基づいて算出されたものである。また、ｖＭＦ関数は、ｃｏｓθの増大に伴って単調増大するため、ｖＭＦは球面上の自然な窓関数としてふるまう。

図８は、２次元のｖＭＦカーネル関数を示したグラフである。図８から分かるように、ｖＭＦカーネル関数は、ほぼ楕円形状を有し、スケールκが増大するにつれて、楕円の長軸が大きくなり、かつ、楕円の中心が縦方向にずれていく関数である。

次に、本発明者が開発したコリレーション族のカーネル関数に基づくミーンシフトアルゴリズムについて説明する。コリレーション族のカーネル関数による確率密度関数ｆ（ｘ）の傾きは、式（８）に示すように、確率密度関数ｆ（ｘ）をｘに関して偏微分することよって得られる。そして、ｆ（ｘ）の傾きが０を満たすｘがモードとなる。したがって、式（８）に対するモードは、式（８）を０としたときの演算結果から式（９）によって表される。

但し、ｕ＝ｘ／｜ｘ｜、ｕ_n＝ｘ_n／｜ｘ_n｜、ｃｏｓθ_n＝ｕ・ｕ_nである。

そして、本発明者が開発したコリレーション族のカーネル関数に対するミーンシフトアルゴリズムは、式（１０）によって表される。

但し、ｉは１以上の整数を示し、ｖ₁は、Ｍ次元の球の表面に設定された始点である。式（１０）に示す演算を繰り返し実行し、ｖ_i+1の値がｖ_iに対してほとんど変化しなくなった時の空間上の地点を当該始点の収束点とする。そして、空間上の各点を始点とし、式（１０）に示す演算を繰り返し実行し、収束点がほぼ同じ点となる始点の集合が一つのクラスタとなる。したがって、ＣＰＵ３は、パターンベクトルｘ_nに対し、式（１０）の演算を繰り返し実行することにより、パターンベクトルｘ_nをクラスタリングする。

図９は、ＣＰＵ３が、２次元のパターンベクトルｘ_nに対し、コリレーション族のカーネル関数によるミーンシフトアルゴリズムを用いてクラスタリングを行なったときの、クラスタリング結果を示した図面であり、（ａ）はスケールκ＝０．１、（ｂ）はスケールκ＝１、（ｃ）はスケールκ＝１０、（ｄ）はスケールκ＝１００の場合を示している。なお、（ａ）〜（ｄ）において、同一のパターンベクトルｘ_nが用いられている。点Ｐはパターンベクトルｘ_nを示している。なお、カーネル関数としてはｖＭＦカーネル関数を用いている。

これらのグラフは、原点Ｏから見て、０〜３６０度のどの方向にパターンベクトルｘ_nが多く分布しているかを示しており、グラフＧが基準線ＤＬに対して外側にふくらんでいる場合、その方向へのパターンベクトルｘ_nの分布は密となり、グラフＧが基準線に対して内側にへこんでいる場合、その方向へのパターンベクトルｘ_nの分布は疎となる。パターンベクトルｘ_nが０〜３６０度の範囲で均一に存在する場合、グラフＧは基準線ＤＬに一致する。（ａ）及び（ｂ）に示すように、スケールκが小さい場合、グラフＧは、基準線ＤＬとほぼ同じ円形となっており、パターンベクトルｘ_nの分布状態を細かく検出することができていない。一方、（ｃ）及び（ｄ）に示すように、スケールκの値を大きくすると、パターンベクトルｘ_nの分布状態がより細かく検出することができる。例えば（ｃ）の場合、グラフＧには４個の山が形成されているため、パターンベクトルｘ_nを少なくとも４つのクラスタに分けることができる。（ｄ）の場合、グラフＧは８個の山が形成されているため、パターンベクトルｘ_nを少なくとも８個のクラスタに分けることができる。

このように、本実施の形態では、パターンベクトルｘ_nの確率密度分布を算出し、この確率密度分布を基に、クラスタリングを行なっているため、母集団中に小数のデータの出入りがある場合であっても、クラスタリング結果が大きく変動することなく、信頼性の高いクラスタリングを行なうことができる。

また、コリレーション族のカーネル関数を用いたミーンシフトアルゴリズムによるクラスタリングを行なっているため、Ｍの値が５０００を超えるような高次元のパターンベクトルｘ_nを取り扱う場合であっても、より信頼性の高いクラスタリングを行なうことができる。

次に、各スケールに対するクラスタリング結果が描画されたブリックビュー（Brick view）について説明する。図１０は、ブリックビューの一例を示した図面である。ブリックビューは縦軸にサンプル、横軸にスケールが設定されている。同一クラスタには、同じ色が付されている。どのクラスタにも属さない孤立点は、例えば黒色（図面では白色）で表示される。このブリックビューによれば、スケール毎にクラスタがどのように分布しているかを一目で理解することができる。また、このブリックビューは、従来の階層化クラスタリングによるクラスタリング結果から作成することもできる。

このブリックビューは、スケール毎に算出されたクラスタリング結果に対し、ＣＰＵ３が、クラスタリナンバリング処理と、サンプル配列処理とを実行することによって作成される。まず、クラスタリナンバリング処理を説明する。図１１は、クラスタリナンバリング処理を示したフローチャートである。図１２（ａ）は、クラスタリナンバリング処理を説明するための図面である。図１１及び図１２（ａ）において、ｉはスケールを示す数値であり、その値は、実際にカーネル関数に代入されるスケール値を大きい順または、小さい順に並べたときの、順序を示す整数値である。Ｌは、ｉ−１番目のスケールのＬ番目のクラスタを示す整数値である。ｋは、ｉ番目のスケールのｋ番目のクラスタを示す整数値である。また、サンプルは、上記パターンベクトルのことを示す。

ＣＰＵ３は、ｉに対し初期値である２を設定し（ステップＳ１１）、Ｌに対し初期値である１を設定し（ステップＳ１２）、ｋに対し初期値である１を設定する（ステップＳ１３）。

ステップＳ１４において、ＣＰＵ３は、ｉ−１番目のスケールのＬ番目のクラスタと、ｉ番目のスケールのｋ番目のクラスタとのサンプルのオーバラップ数を算出する。この場合、図１２（ａ）に示すように、クラスタＣ２１は、サンプルＳ１，Ｓ３，Ｓ４，Ｓ５，Ｓ６，Ｓ７を含み、クラスタＣ３１は、サンプルＳ１，Ｓ３，Ｓ４，Ｓ５を含み、両クラスタ共、サンプルＳ１，Ｓ３，Ｓ４，Ｓ５を含んでいるため、クラスタＣ２１とクラスタＣ３１とのオーバラップ数は４と算出される。また、クラスタＣ３２には、サンプルＳ２，Ｓ５，Ｓ６が含まれているため、クラスタＣ２１とクラスタＣ３２とのオーバラップ数は２と算出される。また、クラスタ３３には、サンプルＳ９，Ｓ１０が含まれているため、クラスタＣ２１と、クラスタＣ３３とのオーバラップ数は０と算出される。

ステップＳ１５において、ｋが最終値でない場合（ステップＳ１５でＮＯ）、ＣＰＵ３は、ｋに１を加算し（ステップＳ１６）、処理をステップＳ１４に戻す。そして、ステップＳ１４〜Ｓ１６の処理が繰り返され、ｉ―１番目のスケールのＬ番目のクラスタと、ｉ番目のスケールの全クラスタとのサンプルのオーバラップ数が算出されていく。すなわち、図１２（ａ）に示すクラスタＣ１１に対し、クラスタＣ２１及びＣ２２のそれぞれのオーバラップ数が順次算出され、次に、クラスタＣ２１に対し、クラスタＣ３１，Ｃ３２，及びＣ３３それぞれのオーバラップ数が順次算出され、次に、クラスタＣ２２に対し、クラスタＣ３１，Ｃ３２，及びＣ３３のオーバラップ数が順次算出されるというようにして、クラスタ間のオーバラップ数が算出されていく。

一方、ステップＳ１５において、ｋが最終値である場合（ステップＳ１５でＹＥＳ）処理がステップＳ１７に進められる。ステップＳ１７において、ＣＰＵ３は、オーバラップ数が最大のクラスタにＬ番目のクラスタのクラスタラベルを承継する。この場合、図１２（ａ）に示すように、クラスタＣ１１と、クラスタＣ２１とのオーバラップ数は６であり、クラスタＣ１１とクラスタＣ２２とのオーバラップ数は４であるため、オーバラップ数が大きなクラスタＣ２１に、クラスタＣ１１のクラスタラベル１が承継され、クラスタＣ２２には、新たなクラスタラベル２が与えられる。

また、クラスタＣ２１とクラスタＣ３１とのオーバラップ数は４であり、クラスタＣ２１とクラスタＣ３２とのオーバラップ数は２であり、クラスタＣ２１とクラスタＣ３３とのオーバラップ数は０であるため、オーバラップ数が最大のクラスタであるクラスタＣ３１にクラスタＣ２１のクラスタラベルである１が承継される。

さらに、クラスタＣ２２とクラスタＣ３１とのオーバラップ数は０であり、クラスタＣ２２とクラスタＣ３２とのオーバラップ数は１であり、クラスタＣ２２とクラスタＣ３３とのオーバラップ数は、２であるため、オーバラップ数が最大のクラスタであるクラスタＣ３３にクラスタＣ２２のクラスタラベル２が承継される。そして、クラスタＣ３２には、新たにクラスタラベル３が付与される。

ステップＳ１８において、Ｌが最終値でない場合（ステップＳ１８でＮＯ）、Ｌに１が加算され（ステップＳ２２）、処理がステップＳ１３に戻される。一方、Ｌが最終値である場合（ステップＳ１８でＹＥＳ）、処理がステップＳ１９に進められる。ステップＳ１９において、ｉ番目のスケール中、クラスタラベルが承継されていないクラスタに対し、新たなクラスタラベルが付与される。この場合、クラスタラベルは、同一スケール内において、数が連続し、かつ、各クラスタ固有のクラスタラベルが付与される。

ステップＳ２０において、ｉが最終値でない場合（ステップＳ２０でＮＯ）、ｉに１が加算され（ステップＳ２１）、処理がステップＳ１２に戻される。一方、ｉが最終値の場合（ステップＳ２０でＹＥＳ）、処理が終了される。以上によって、図１２（ａ）に示すように、各クラスタにクラスタラベルが付与される。そして、各クラスタラベルに固有の色を割り当てることにより、各クラスタが色分けされる。すなわち、クラスタリナンバリング処理によって、ブリックビューを表示するにあたり、各サンプルに対して、スケール毎に付すべき色が決定されることとなる。

クラスタリナンバリング処理によって、前後するスケール間において、オーバラップするサンプル数が大きいクラスタには、同一の色が付され、この状態でスケール毎のクラスタリング結果を表示すれば、スケール毎のクラスタリング結果を理解することも可能である。しかしながら、このままでは、スケール方向においてクラスタが点在して表示されてしまう。例えば、図１２（ａ）において、クラスタＣ４５に、クラスタＣ３１のクラスタラベル１が承継されたとすると、スケール方向において、同一クラスタが連続して表示されないこととなってしまう。そこで、スケール方向において、同一クラスタが可能な限り連続して表示されるように、各サンプルを配列するサンプル配列処理を実行し、スケール毎のクラスタリング結果を分かりやすく示すブリックビューを作成する。以下、サンプル配列処理について説明する。

図１３は、サンプル配列処理を示すフローチャートである。図１４〜図１６は、Ａ〜Ｊのサンプル名を有する１０個のサンプルに対して、サンプル配列処理を施した時の処理結果を示した表である。以下、図１３〜図１６を用いて、サンプル配列処理を説明する。図１３〜図１６において、ｉはスケールの値を示し、ｆ（ｉ）は、ｉスケールにける各サンプルが有する得点を示し、その値は実数値をとる。δ（ｉ）は、ｆ（ｉ）に基づくサンプルの順位を示し、得点の高いサンプルほど高い順位が与えられ、その値は、１からサンプルの個数までの連続する整数となる。Δ（ｉ）は、ｉスケールにおいて、各サンプルがどのクラスタに属するかを示す数値であり、その値は、１から始まりクラスタ数分の整数値をとる。また、サンプルは、上述のパターンベクトルｘ_nを示す。

まず、ＣＰＵ３は、ｉに１を設定し（ステップＳ１０１）、各サンプルに固有の得点を付与し（ステップＳ１０２）、各サンプルを得点の大きい順に並べ替える（ステップＳ１０３）。この場合、図１４（ａ）に示すように、サンプルＡ〜Ｊに対し、０．８１〜０．０２まで、各サンプル固有の得点ｆ（１）が付与され、得点の高い順にサンプルＡ〜Ｊが配列される。また、図１４（ｂ）のδ（１）の列に示すように、サンプルＡ〜Ｊに対し、得点ｆ（１）の高い順に、１０〜１の順位が与えられる。ここで、得点ｆ（ｉ）は、例えば式（１）及び式（３）又は式（１）及び式（５）によって各スケールに対する確率密度関数を用いて算出される。

ステップＳ１０４において、ｉが最終値でない場合（ステップＳ１０４でＮＯ）、ＣＰＵ３は、ｉに１を加算する（ステップＳ１０５）。

ステップＳ１０６において、ＣＰＵ３は、クラスタ毎に、一つ前のスケールの各サンプルの順位δ（ｉ−１）の平均値を算出し、各クラスタに得点Ａ_c（ｉ）を付与する。ｃは、クラスタを識別するために便宜上付した数値（仮のクラスタ番号）を示し、以下に示すクラスタの順位Δ（ｉ）とは相違する。図１４（ｂ）の場合、クラスタ２の得点Ａ₂（２）は、クラスタ２を構成するサンプルのスケール１における順位δ（１）を用いて、Ａ₂（２）＝（１０＋９＋８＋７＋５）／５＝７．８と算出され、クラスタ１の得点Ａ₁（２）は、クラスタ１を構成するスケール１におけるサンプルの順位δ（１）を用いて、Ａ₁（２）＝（６＋４＋３＋２＋１）／５＝３．２と算出される。

ステップＳ１０７において、ＣＰＵ３は、クラスタの得点Ａｃ（ｉ）の大きい順に、スケールｉにおけるクラスタの順位Δ（ｉ）を決定する。図１４（ｂ）例では、Δ（２）の列に示すように、サンプルＡ，Ｈ，Ｂ，Ｃ，Ｅからなるクラスタ２（２は便宜上付した数値）は、Ａ₂（２）＝７．８であり、サンプルＤ，Ｆ，Ｇ，Ｉ，Ｊからなるクラスタ１（１は便宜上付した数値）、Ａ₁（２）＝３．２であり、Ａ₂（２）＞Ａ₁（２）であるため、クラスタ２はクラスタ順位Δ（２）＝２が与えられ、クラスタ１はクラスタ順位Δ（２）＝１が与えられる。すなわち、クラスタ２は、「２」のクラスタ番号が正式に与えられ、クラスタ１は「１」のクラスタ番号が正式に与えられることとなる。

ステップＳ１０８において、ＣＰＵ３は、式（１）を用いて各サンプルに対する得点ｆ（ｉ）を新たに付与する（ステップＳ１０８）。この場合、図１４（ｂ）のｆ（２）の列、あるいは、図１５（ａ）のｆ（３）の列に示すように、サンプルＡ〜Ｊに対し、新たな得点が付与される。

ステップＳ１０９において、ＣＰＵ３は、同一クラスタに属する各サンプルが連続して配置され、且つ、得点ｆ（ｉ）が高い順となるように各サンプルを並べ替える。この場合、図１４（ｂ）及び（ｃ）を比較すれば分かるように、クラスタ１の中に孤立して配置されたクラスタ２に属するサンプルＥは、得点ｆ（２）が０．４であるため、得点ｆ（２）が０．５２であるサンプルＣの下の位置に配置されている。また、クラスタ２の中に孤立して配置されたクラスタ１に属するサンプルＤは、得点ｆ（２）が０．５であるため、得点ｆ（２）が０．２１であるサンプルＦの上の位置に配置されている。これによって、サンプルＡ〜Ｊは、図１４（ｃ）に示すように並べ替えられる。

さらに、図１５（ａ）及び（ｂ）のΔ（３）の列を比較すれば分かるように、クラスタ２の中に孤立して配置されたクラスタ３に属するサンプルＥは、得点ｆ（３）が０．４であるため、得点ｆ（３）が０．９１であるクラスタ３に属するサンプルＢの下の位置に配置される。さらに、クラスタ２に属するサンプルＣ及びＤは、各々、得点（２）が０．５２、０．５であるため、サンプルＣの下にサンプルＤが配置される。

さらに、クラスタ３に属するサンプルＡ及びＨは、得点ｆ（３）が、各々、１．２１及び１．４２であるため、図１５（ａ）では、サンプルＡの下に位置していたサンプルＨが、図１５（ｂ）では、サンプルＡの上に配置される。

さらに、図１６（ａ）及び（ｂ）のΔ（４）の列を比較すれば分かるように、クラスタ４に属するサンプルＢ及びＥは、得点ｆ（４）が、各々、１．２及び１．４であるため、（ａ）では、サンプルＢの下に位置していたサンプルＥが、（ｂ）では、サンプルＢの下に配置される。

ステップＳ１１０において、ＣＰＵ３は、ステップＳ１０９の処理によって並べ替えられた各サンプルに対し、サンプルの順位δ（ｉ）を算出する。この場合、サンプルＡ〜Ｊは、例えば図１４（ｃ）のδ（２）の列で示される順位を示す数値が与えられ、処理がステップＳ１０４に戻される。一方、ステップＳ１０４において、ｉが最終値となった場合（ステップＳ１０４でＹＥＳ）、処理が終了される。

図１２（ｂ）を用い、このサンプル配列処理の説明をまとめると、ｉ＝１において、各サンプルの得点ｆ（１）が算出され、得点の高い順に、各サンプルが並べ替えられ、各サンプルの順位を示すδ（１）が算出される。次に、ｉ＝２において、クラスタを分類するための順位を示すΔ（２）が決定され、各サンプルの得点ｆ（２）が新たに算出され、得点ｆ（２）を基に、同一クラスタに属するサンプルが連続して配置されるように、各サンプルが並べ替えられ、各サンプルの順位δ（２）が決定される。このような処理が全スケールに亘って行なわれ、最終的に得られたδ（ｎ）によって、各サンプルの配列位置が決定される。すなわち、図１０に示す縦軸に対する各サンプルの座標が決定される。

そして、サンプル配列処理によって決定された順序で、各サンプルを図１０に示す縦軸上に配列し、配列したサンプルに対し、クラスタリナンバリング処理によって、スケール毎に決定された色を用い、各サンプルをスケール毎に色分けして表示すると、図１０に示すブリックビューが得られる。

本発明の一実施の形態によるクラスタ解析装置の構成を示すブロック図である。 クラスタ解析装置の動作を示すフローチャートである。 遺伝子発現プロファイルデータのデータ構造を示す模式図である。 ２次元空間内に設定されたパターンベクトルｘ_nの確率密度分布を示した図面であり、横軸はサンプル１の発現量を示し、縦軸はサンプル２の発現量を示している。 カーネル関数と確率密度分布との関係を説明するための図面であり、（ａ）は半径がｒ１のカーネル関数Ｋ（ｘ_n，ｘ）を示し、（ｂ）は半径がｒ２（ｒ２＜ｒ１）のカーネル関数（ｘ_n，ｘ）を示し、（ｂ）は（ａ）のカーネル関数を用いて算出された確率密度関数ｆ（ｘ）を示し、（ｄ）は（ｃ）のカーネル関数を用いて算出された確率密度関数ｆ（ｘ）を示している。 ガウシアンカーネル関数を用いて算出された、２次元のパターンベクトルｘ_nに対する確率密度分布を示しており、（ａ）はスケールｈ＝１、（ｂ）はスケール０．７、（ｃ）はスケール０．５、（ｄ）はスケール０．２の場合を示している。 コリレーション族のカーネル関数を説明するための図面である。 ２次元のｖＭＦカーネル関数を示したグラフである。 ２次元のパターンベクトルｘ_nに対し、ｖＭＦカーネル関数によるミーンシフトアルゴリズムを用いてクラスタリングを行なったときの、クラスタリング結果を示した図面であり、（ａ）はスケールκ＝０．１、（ｂ）はスケールκ＝１、（ｃ）はスケールκ＝１０、（ｄ）はスケールκ＝１００の場合を示している。 ブリックビューの一例を示した図面である。 クラスタリナンバリング処理を示したフローチャートである。 （ａ）は、クラスタリナンバリング処理を説明するための図面であり、（ｂ）は、サンプル配列処理を説明するための図面である。 サンプル配列処理を示すフローチャートである。 １０個のサンプルに対し、サンプル配列処理を施した時の処理結果を示した表である。 １０個のサンプルに対し、サンプル配列処理を施した時の処理結果を示した表である。 １０個のサンプルに対し、サンプル配列処理を施した時の処理結果を示した表である。 階層化クラスタリングを説明するための従来図である。

符号の説明

１ 入力装置
２ ＲＯＭ
３ ＣＰＵ
４ ＲＡＭ
５ 外部記憶装置
６ 表示装置
７ 記録媒体駆動装置
８ 記録媒体

Claims (8)

1. 生命情報に関する複数のデータをクラスタリングするクラスタ解析装置であって、
   カーネル関数を用い、前記データの確率密度分布を推定する推定手段と、
   推定された確率密度分布を基に、密度分布モード探索アルゴリズムを用いて前記データのクラスタリングを行なうクラスタリング手段とを備えることを特徴とするクラスタ解析装置。
2. 前記生命情報に関するデータが、遺伝子発現プロファイルデータであることを特徴とする請求項１記載のクラスタ解析装置。
3. 前記カーネル関数のスケールを設定するスケール設定手段を更に備え、
   前記推定手段は、前記スケール毎に、前記データの確率密度分布を推定し、
   前記クラスタリング手段は、前記スケール毎に前記データのクラスタリングを行なうことを特徴とする請求項１又は２記載のクラスタ解析装置。
4. 前記クラスタリング手段によるクラスタリング結果を基に、各データが、スケール毎に、どのクラスタに属するかをクラスタ別に色分けして示す画像を作成する画像作成手段を更に備え、
   前記画像作成手段は、各クラスタがスケール全域に亘って連続して表示されるように前記画像を作成することを特徴とする請求項３記載のクラスタ解析装置。
5. 前記密度分布モード探索アルゴリズムは、ミーンシフトアルゴリズムであることを特徴とする請求項１〜４のいずれかに記載のクラスタ解析装置。
6. 前記カーネル関数は、式（Ａ）で表され、
   前記密度分布モードアルゴリズムは、各データの座標ベクトルを始点ｖ₁として、式（Ｂ）の演算を繰り返し実行し、ｖ_i+1のｖ_iに対する変化が、ほぼなくなった点を収束点モードとし、前記収束点モードが同一となる始点ｖ₁を１つのクラスタとすることを特徴とする請求項１〜５のいずれかに記載のクラスタ解析装置。
   

   

   但し、Ｋ（ｘ_n，ｘ）は、カーネル関数を示し、ｇはｃｏｓθ_nの増大に応じて出力値が増大する任意の関数を示し、ｇ´はｇの微分を示し、ｕ_n＝ｘ_n／｜ｘ_n｜を示し、ｖ₁＝ｘ／｜ｘ｜を示し、ｘ_nは生命情報に関するデータの座標ベクトルを示し、ｘは前記生命情報に関するデータｘ_nが配置される空間の任意の点の座標ベクトルを示し、θ_nはｘ_nとｘとがなす角度を示し、｜ｘ｜はｘの大きさを示し、｜ｘ_n｜はｘ_nの大きさを示し、ｎは１以上の整数を示す。
7. 生命情報に関する複数のデータをコンピュータによりクラスタリングするクラスタ解析方法であって、
   前記コンピュータが、カーネル関数を用い、前記データの確率密度分布を推定するステップと、
   前記コンピュータが、推定された確率密度分布を基に、密度分布モード探索アルゴリズムを用いて、前記データのクラスタリングを行なうステップとを備えることを特徴とするクラスタ解析方法。
8. 生命情報に関する複数のデータをクラスタリングするクラスタ解析プログラムであって、
   カーネル関数を用い、前記データの確率密度分布を推定する推定手段と、
   推定された確率密度分布を基に、密度分布モード探索アルゴリズムを用いて、前記データのクラスタリングを行なうクラスタリング手段としてコンピュータを機能させることを特徴とするクラスタ解析プログラム。

Images