(57)【要約】【課題】遺伝子の経時的発現量データに基づく遺伝子クラスタリングを、先験的なデータ予測を行うことなく、高精度に行うことができる遺伝子クラスタリングツールの提供。【解決手段】（１）遺伝子の経時的な発現量変化を示すデータから、各データ間の類似度を反映した特徴量を算出するステップS1と、（２）全ての遺伝子間の組合せについて、算出された特徴量から類似度行列Ｍの固有ベクトルを算出するステップS2と、（３）類似度行列Ｍを、固有ベクトルの固有値を維持したまま、ブール行列Ｎに変換するステップS3と、（４）ブール行列Ｎに基づいて各データをクラスタリングするステップS4と、を少なくとも行う遺伝子クラスタリングプログラムを提供する。【選択図】図１
2010157214.tif 000002
【特許請求の範囲】
【請求項１】
（１）遺伝子の経時的な発現量変化を示すデータから、各データ間の類似度を反映した特徴量を算出するステップと、
（２）全ての遺伝子間の組合せについて、算出された特徴量から類似度行列Ｍの固有ベクトルを算出するステップと、
（３）類似度行列Ｍを、固有ベクトルの固有値を維持したまま、ブール行列Ｎに変換するステップと、
（４）ブール行列Ｎに基づいて各データをクラスタリングするステップと、
を少なくとも行う遺伝子クラスタリングプログラム。
【請求項２】
前記（１）のステップにおいて、線形回帰分析又はウェーブレット変換により、前記データから前記特徴量を算出する請求項１記載の遺伝子クラスタリングプログラム。
【請求項３】
前記（２）のステップにおいて、カーネル法又はコサイン類似度により、前記特徴量から前記固有ベクトルを算出する請求項２記載の遺伝子クラスタリングプログラム。
【請求項４】
前記（３）のステップにおいて、FSNN（Filter by Symmetric Nearest Neighbors）アルゴリズムにより、類似度行列Ｍをブール行列Ｎに変換する請求項３記載の遺伝子クラスタリングプログラム。
【請求項５】
さらに、前記（３）のステップにおいて、FSNNアルゴリズムによる処理後、グラフラプラシアン、マルコフ連鎖、DSA（Doubly-Stochastic Approximation）アルゴリズム又はDSS（Doubly-Stochastic Scaling）アルゴリズムのいずれかにより、行列を標準化する請求項４記載の遺伝子クラスタリングプログラム。
【請求項６】
前記（４）のステップにおいて、EM（Expectation Maximization）アルゴリズムとCP（Complete Positive factorization）アルゴリズムにより、ソフトクラスタリングを行う請求項５記載の遺伝子クラスタリングプログラム。
【請求項７】
さらに、前記（４）のステップにおいて、ソフトクラスタリング後、BAV（Bregman-Arthur-Vassilvitskiiinitialization）アルゴリズムにより、ハードクラスタリングを行う請求項６記載の遺伝子クラスタリングプログラム。
【請求項８】
請求項１記載の遺伝子クラスタリングプログラムをコンピュータが読み取り可能に記録した記録媒体。
【請求項９】
（１）遺伝子の経時的な発現量変化を示すデータから、各データ間の類似度を反映した特徴量を算出するステップと、
（２）算出された特徴量から、全ての遺伝子間の組合せについて類似度行列Ｍの固有ベクトルを算出するステップと、
（３）類似度行列Ｍを、固有ベクトルの固有値を維持したまま、ブール行列Ｎに変換するステップと、
（４）ブール行列Ｎに基づいて各データをクラスタリングするステップと、
を少なくとも行う遺伝子クラスタリング方法。
【請求項１０】
（１）遺伝子の経時的な発現量変化を示すデータから、各データ間の類似度を反映した特徴量を算出する手段と、
（２）算出された特徴量から、全ての遺伝子間の組合せについて類似度行列Ｍの固有ベクトルを算出する手段と、
（３）類似度行列Ｍを、固有ベクトルの固有値を維持したまま、ブール行列Ｎに変換する手段と、
（４）ブール行列Ｎに基づいて各データをクラスタリングする手段と、
を少なくとも備える遺伝子クラスター解析装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、遺伝子クラスタリングプログラム、遺伝子クラスタリング方法及び遺伝子クラスター解析装置に関する。より詳しくは、各遺伝子が示す経時的な発現量変化の類似度に基づいて、それぞれの遺伝子を所定のクラスターに分類する遺伝子クラスタリングプログラム等に関する。
【背景技術】
【０００２】
システム生物学の分野では、遺伝子の経時的な発現量変化や局在変化、活性変化の測定データに基づいて、遺伝子が細胞内で構築するシグナルネットワークを解明しようとする試みがなされている。
【０００３】
細胞内シグナルネットワークは、ダイナミックに変化する多階層のネットワーク構造によって構築されている。近年、細胞内シグナルネットワークを構成するひとつの基本的なネットワーク構造として、「ボウタイ（bow-tie）シグナルネットワーク」の存在が提唱されてきている（非特許文献１及び非特許文献2）。
【０００４】
ボウタイシグナルネットワーク（以下、単に「ボウタイネットワーク」という）は、そのネットワーク構造を蝶ネクタイの形に見立て、蝶ネクタイの結び目に、刺激に対する細胞の応答を規定する「分類器」として機能するコア分子の存在を仮定するものである。すなわち、ボウタイネットワークでは、細胞内及び細胞間のシグナル伝達における多種多様なインプットは、結び目に位置するコア分子に集約される。そして、コア分子がインプットに応じて細胞内濃度を変化させることにより、その濃度状態に応じてシグナル下流に位置する所定の遺伝子集団（クラスター）が活性化し、特定のアウトプットが現れる。
【０００５】
ボウタイネットワークは、免疫細胞間のシグナル伝達の他、代謝シグナル伝達（非特許文献１）やToll-likeレセプターシグナル伝達（非特許文献2）、上皮成長因子シグナル伝達（非特許文献3）への関与も報告されている。ボウタイネットワークは、強固でありながら進化のための柔軟性をも備えた優れたネットワーク構造であることが明らかにされつつある（非特許文献4及び非特許文献5）。
【０００６】
ボウタイネットワークにおいては、シグナル下流に位置する遺伝子は、所定のコア分子濃度のもとに一つの遺伝子集団（クラスター）としてクラスタリングされる。遺伝子の経時的な発現量変化や局在変化、活性変化の測定データから、各遺伝子が属するクラスターを明らかにし、ボウタイネットワークを解析するためには、ネットワーク全体の構成を把握して、各クラスター間の関係を予測するための優れた幾何学的ツールが必要となる。
【０００７】
これまでに、K平均（k-means）法（非特許文献6）や、階層的クラスタリング（非特許文献7）、自己組織化マップ（非特許文献8）に基づくツールが開発されている。
【０００８】
しかし、演算処理をひとつの工程のみで行うこれらのツールには、それぞれ欠点が存する。すなわち、階層的クラスタリングでは、クラスターを重ね合わせていくことにより、個々のデータ要素を階層化するため、柔軟性のない系統樹しか作成できない。また、階層的クラスタリングでは、遺伝子を一対一の類似性によってクラスター化していくため、最終的に一つのクラスター内に位置づけられた遺伝子間に生物学的な関連性がみられないという事態が生じ得る。
【０００９】
自己組織化マップ（Self Organizing Map; SOM）に基づくツール（例えば、「GENECLUSTER」）は、特にデータの予備解析に優れた能力を発揮するが、クラスター数の予測初期値のためのグリッドサイズを予め設定することが必要となる。従来のK平均法でも、初めにクラスター数の設定が必要で、クラスタリング結果がこの設定数に依存するため、生物学的に意味のない結果が得られてしまうことがある。
【００１０】
「GENEPattern」（非特許文献9）は、これらの従来ツールを水平統合して得られた、現状入手可能な最も有用なツールである。しかし、遺伝子の経時的発現量データ等から正確に各遺伝子をクラスタリングして、ボウタイネットワークを解明するためには、いまだ十分な能力を有しているとはいえない。
【先行技術文献】
【非特許文献】
【００１１】
【非特許文献１】“The Edinburghhuman metabolic network reconstruction and its functional analysis.” Molecular System Biology, 2007;3:135
【非特許文献２】“A comprehensive map of the toll-like receptor signaling network.” Molecular System Biology, 2006;2:2006.0015
【非特許文献３】“A comprehensive pathway map of epidermal growth factor receptor signaling.” Molecular System Biology, 2005;1:2005.0010
【非特許文献４】“Bow ties, metabolism and disease.” Trends in Biotechnology, 2004;22(9): 446-50
【非特許文献５】“Biological robustness.” Nature Reviews Genetics, 2004;5(11):826-37
【非特許文献６】“Systematic determination of genetic network architecture.” Nature Genetics, 1999;22(3):281-285
【非特許文献７】“Cluster analysis and display of genome-wide expression patterns.” Proceeding of National Academyof Sciences, 1998;95(25):14863-14868
【非特許文献８】“Interpreting patterns of gene expression with self-organizing maps: Methods and application to hematopoietic differentiation.” Proceeding of National Academy of Sciences, 1999;96(6):2907-2912
【非特許文献９】“GenePattern2.0.” Nature Genetics, 2006;38:500-501
【発明の概要】
【発明が解決しようとする課題】
【００１２】
そこで、本発明は、遺伝子の経時的発現量データに基づく遺伝子クラスタリングを、先験的なデータ予測を行うことなく、高精度に行うことができる遺伝子クラスタリングツールを提供することを主な目的とする。
【課題を解決するための手段】
【００１３】
上記課題解決のため、本発明は、（１）遺伝子の経時的な発現量変化を示すデータから、各データ間の類似度を反映した特徴量を算出するステップと、（２）全ての遺伝子間の組合せについて、算出された特徴量から類似度行列Ｍの固有ベクトルを算出するステップと、（３）類似度行列Ｍを、固有ベクトルの固有値を維持したまま、ブール行列Ｎに変換するステップと、（４）ブール行列Ｎに基づいて各データをクラスタリングするステップと、を少なくとも行う遺伝子クラスタリングプログラムを提供する。
この遺伝子クラスタリングプログラムは、前記（１）のステップにおいて、線形回帰分析又はウェーブレット変換により、前記データから前記特徴量を算出する。
前記（２）のステップにおいては、カーネル法又はコサイン類似度により、前記特徴量から前記固有ベクトルを算出する。
また、前記（３）のステップにおいては、FSNN（Filter by Symmetric Nearest Neighbors）アルゴリズムにより、類似度行列Ｍをブール行列Ｎに変換する。
さらに、前記（３）のステップでは、FSNNアルゴリズムによる処理後、グラフラプラシアン、マルコフ連鎖、DSA（Doubly-Stochastic Approximation）アルゴリズム又はDSS（Doubly-Stochastic Scaling）アルゴリズムのいずれかにより、行列を標準化する。
この遺伝子クラスタリングプログラムは、前記（４）のステップにおいて、EM（Expectation Maximization）アルゴリズムとCP（Complete Positive factorization）アルゴリズムにより、ソフトクラスタリングを行う。
さらに、前記（４）のステップでは、ソフトクラスタリング後、BAV（Bregman-Arthur-Vassilvitskiiinitialization）アルゴリズムにより、ハードクラスタリングを行う。
【００１４】
また、本発明は、上記遺伝子クラスタリングプログラムをコンピュータが読み取り可能に記録した記録媒体を提供する。
併せて、本発明は、（１）遺伝子の経時的な発現量変化を示すデータから、各データ間の類似度を反映した特徴量を算出するステップと、（２）算出された特徴量から、全ての遺伝子間の組合せについて類似度行列Ｍの固有ベクトルを算出するステップと、（３）類似度行列Ｍを、固有ベクトルの固有値を維持したまま、ブール行列Ｎに変換するステップと、（４）ブール行列Ｎに基づいて各データをクラスタリングするステップと、を少なくとも行う遺伝子クラスタリング方法を提供する。
さらに、（１）遺伝子の経時的な発現量変化を示すデータから、各データ間の類似度を反映した特徴量を算出する手段と、（２）算出された特徴量から、全ての遺伝子間の組合せについて類似度行列Ｍの固有ベクトルを算出する手段と、（３）類似度行列Ｍを、固有ベクトルの固有値を維持したまま、ブール行列Ｎに変換する手段と、（４）ブール行列Ｎに基づいて各データをクラスタリングする手段と、を少なくとも備える遺伝子クラスター解析装置をも提供する。
【発明の効果】
【００１５】
本発明により、遺伝子の経時的発現量データに基づく遺伝子クラスタリングを、先験的なデータ予測を行うことなく、高精度に行うことができる遺伝子クラスタリングツールが提供される。
【図面の簡単な説明】
【００１６】
【図１】本発明に係る遺伝子クラスタリングプログラムの処理ステップを説明するフローチャートである。
【図２】本発明に係る遺伝子クラスタリングプログラムで処理される遺伝子の経時的な発現量変化を示すデータの一例を説明する図である。
【図３】ウェーブレット変換によるデータ処理を説明する概念図である。
【図４】遺伝子の経時的発現量変化のヒストグラム化を説明する概念図である。
【図５】特徴量の算出ステップの前後でのデータ次元の変化を説明する概念図である。
【図６】遺伝子 iの対称最近傍を説明する概念図である。
【図７】類似度行列Ｍからブール行列Ｎへの変換処理を説明する概念図である。
【図８】ブール行列とDSS行列を説明する概念図である。
【図９】本発明に係る遺伝子クラスタリングプログラムにおいて、最終的なクラスタリング結果が得られるまでのデータ処理を説明する概念図である。
【図１０】本発明に係る遺伝子クラスター解析装置の構成例を示すブロック図である。
【発明を実施するための形態】
【００１７】
本発明に係る遺伝子クラスタリングプログラムは、（１）遺伝子の経時的な発現量変化を示すデータから、各データ間の類似度を反映した特徴量を算出するステップと、（２）全ての遺伝子間の組合せについて、算出された特徴量から類似度行列Ｍの固有ベクトルを算出するステップと、（３）類似度行列Ｍを、固有ベクトルの固有値を維持したまま、ブール行列Ｎに変換するステップと、（４）ブール行列Ｎに基づいて各データをクラスタリングするステップと、を少なくとも行うものである。以下、これらのステップについて順に説明する。
【００１８】
１．特徴量の算出
本ステップは、上記（１）の「遺伝子の経時的な発現量変化を示すデータから、各データ間の類似度を反映した特徴量を算出するステップ」に対応するステップである（図1中、S1参照）。
【００１９】
まず、遺伝子の経時的な発現量変化を示すデータから、スケーリング関数係数D4-20を用いて、線形回帰分析又はウェーブレット変換（Haar wavelet transform又はDaubechies wavelet transform）により、各データ間の類似度を反映した特徴量を算出する。図２に、遺伝子の経時的な発現量変化を示すデータの一例を示す。ここでは、３つの遺伝子（gene）a, b, cについて、４つのタイムポイント（time）1, 2, 3, 4で発現量を測定して得たデータを示す。
【００２０】
線形回帰分析は、発現量の変化を示す変動曲線を比較するための簡便な方法である。他方、ウェーブレット変換は、変動曲線の全ての経時的情報を集約することができる。従って、ウェーブレット変換では、従来の解析手法では不完全な測定データとして解析対象から排除されてしまうような１つのタイムポイントでしか発現データが得られていない遺伝子であっても解析対象とすることが可能である。
【００２１】
図３に、ウェーブレット変換（Harr wavelet transform）によるデータ処理の概念図を示す。
【００２２】
ウェーブレット変換では、遺伝子の経時的発現量変化データ（ここでは、経時的に9, 7, 3, 5と変化するデータ）を、変動曲線に代えてヒストグラムとして処理し、このヒストグラムを、例えば４つで一組のHarrウェーブレットコンポーネントに分解する（図３（Ａ）参照）。
【００２３】
このデータは、４次元では平均［9, 7, 3, 5］、２次元では平均［8, 4］、係数［1, -1］と表され、１次元では平均［6］、係数［2］と表される（図３（Ｂ）参照）。従って、このデータは、１次元ウェーブレット変換により、［6（basis）, 2, 1, -1（係数）］と処理される（図３（Ｃ）参照）。ウェーブレット変換では、このように遺伝子の経時的発現量変化データをヒストグラムとして処理することで、変動曲線として処理する場合に比べて、顕著に少ない係数による最適なフィッティングが可能となる。
【００２４】
図４に、図１に示したような遺伝子の経時的な発現量変化を示すデータを、ヒストグラム化する概念図を示す。図４（Ａ）中、実線あるいは点線で示す発現量変化は、（Ｂ）に示すようなヒストグラムに変換できる。
【００２５】
本ステップでは、遺伝子の経時的発現量変化データをこのように変換したヒストグラムとして処理し、特徴量を上記のような係数の集まりとして算出することにより、データ次元を減少させる。図５に、本ステップの前後でのデータ次元の変化を概念図として示す。
【００２６】
２．類似度行列（Similarity matrix）の固有ベクトルの算出
次に、全ての遺伝子間の組み合せについて、算出された特徴量からカーネル（Heat kernel）法又はコサイン類似度により、類似度行列M（semi-definite positive matrix M）の固有ベクトルを算出する。以下、類似度行列Ｍを、単に「行列Ｍ」と称する。
【００２７】
このステップは、本発明に係る遺伝子クラスタリングプログラムの（２）のステップ「全ての遺伝子間の組合せについて、算出された特徴量から類似度行列Ｍの固有ベクトルを算出するステップ」に対応するステップ（図1中、S2参照）である。
【００２８】
（2-1）カーネル法による行列M
２つの遺伝子をi及びj（i及びjは１以上の整数）とすると、カーネル法による行列Mの横列（row）i、縦列（column）jの入力は、以下の式（１）によって定義される。この入力は、遺伝子i及び遺伝子jの類似度を示す。
【００２９】
【数１】
2010157214.tif 000003
【００３０】
（2-2）コサイン類似度による行列M
また、コサインシミラリティ法による行列Mでは、以下の式（２）によって定義される。
【００３１】
【数２】
2010157214.tif 000004
【００３２】
３．標準化
（3-1）FSNNによる変換
この行列Mを、LLE（Local Linear Embedding）アルゴリズムに類似したFSNN（Filter by Symmetric Nearest Neighbors）アルゴリズムによって処理する。これにより、LLEアルゴリズムに比べ、より処理が容易なブール行列（Boolean similarity matrix）を得ることができる。FSNNアルゴリズム及びLLEアルゴリズムについては、それぞれ「A simple locally adaptive nearest neighbor rule with application to pollution forecasting. International Journal on Pattern Recognition and Artificial Intelligence, 17: 1-14, 2003」及び「Nonlinear dimensionality reduction by locally linear embedding. Science, 290: 2323-2326, 2000」を参照。
【００３３】
このステップは、本発明に係る遺伝子クラスタリングプログラムの（３）のステップ「類似度行列Ｍを、固有ベクトルの固有値を維持したまま、ブール行列Ｎに変換するステップ」に対応するステップである（図1中、S3参照）。
【００３４】
FSNNアルゴリズムでは、まず、1以上の整数qを指定する。そして、行列Mの各横列iについて、遺伝子iのｑ最近傍を定義する縦列を与えるｑ最大入力を求める。
【００３５】
次に、行列Mの各横列iについて、遺伝子iの最近傍及び遺伝子iが最近傍となる遺伝子を集約することにより、各遺伝子iの対称最近傍（symmetric nearest neighbors）を算出する（「A simple locally adaptive nearest neighbor rule with application to pollution forecasting. International Journal on Pattern Recognition and Artificial Intelligence, 17: 1-14, 2003」参照）。図６に、算出された遺伝子iの対称最近傍の概念図を示す。
【００３６】
最後に、行列Mを、対称最近傍を表示するブール行列Nに置換する。以下、ブール行列Nを、単に「行列N」と称する。行列Nでは、遺伝子i, jが対称最近傍である場合mijに対して1を付与し、対称最近傍でない場合には0を付与する。これにより、対称性の行列Mを、非対称性でもよい行列Nに、固有値を維持させたまま標準化することができる。図７に、類似度行列Ｍからブール行列Ｎへの変換処理の概念図を示す。
【００３７】
LLEアルゴリズムは標準的なq最近傍法に基づいて、各遺伝子の近傍を制約された条件下で再構築するものである。すなわち、LLEアルゴリズムでは、最終的に得られる行列が対称性でなければならないという制約がある。これに対して、FSNNアルゴリズムではこのような制約がないため、簡便で高速な処理を行うことができる。
【００３８】
（3-2）DSS等による行列の標準化
FSNNアルゴリズムによる標準化に続けて、さらにグラフラプラシアン、マルコフ連鎖、DSA（Doubly-Stochastic Approximation）アルゴリズムやDSS（Doubly-Stochastic Scaling）アルゴリズムによる標準化のいずれかを行うことにより、固有値の摂動（perturbation）を減少させることができる。なお、このうち、DSSアルゴリズムは、本発明に係る遺伝子クラスタリングプログラムのため、本発明者らによって新たに作製された新規なアルゴリズムである。
【００３９】
グラフラプラシアンを用いた標準化では、行列Nは以下の式（３）によって定義される。グラフラプラシアンについては、「On spectral Clustering: Analysis and an Algorithm. Neural Information Processing Systems, 2001」を参照。
【００４０】
【数３】
2010157214.tif 000005
【００４１】
マルコフチェーンを用いた標準化では、行列Nは以下の式（４）によって定義される。マルコフチェーンについては、「Soft Membership for spectral clustering, with application to permeable language distinction. Pattern Recognition, 2008」を参照。 
【００４２】
【数４】
2010157214.tif 000006
【００４３】
DSAアルゴリズムを用いた標準化では、行列Nは以下の式（５）に対する解によって定義される。DSAアルゴリズムについては、「Doubly Stochastic Normalization for Spectral Clustering. Neural Information and Processing Systems, 2006」を参照。
【００４４】
【数５】
2010157214.tif 000007
【００４５】
より詳細には、DSAアルゴリズムを用いた標準化は、Xの初期値X0をMとして初期化を行った後、以下の式（６）を解くステップと、この式（６）を満たさない入力Xtを０に置換するステップと、を繰り返すことにより行われる。
【００４６】
【数６】
2010157214.tif 000008
【００４７】
このDSAアルゴリズムを用いた標準化の問題点として、行列Mにおける固有ベクトルの偏差に比して、行列Nでの偏差が大きくなるという点がある。このDSAアルゴリズムにおける偏差拡大を抑制するため、本発明に係る遺伝子クラスタリングプログラムにおいては、DSAアルゴリズムを用いた標準化に先立ってDSSアルゴリズムによる標準化を行うことが好ましい。これにより、偏差のシフトを抑制してマニホールド(manifold)を保持することができる。
【００４８】
DSSアルゴリズムを用いた標準化では、以下の式（７）に基づいて２値検索を行い、行列Mを行列κＭに置換する。得られた行列κＭは、二重確率行列に類似したものであり、当初の固有ベクトルを維持したものとなっている。図８に、DSSアルゴリズムにより得られる、新しい遺伝子座標による行列κＭ（Ｂ）と、ブール行列Ｎ（Ａ）の概念図を示す。
【００４９】
【数７】
2010157214.tif 000009
【００５０】
より詳細には、mi.、m.j、m..を、mi.＝ΣjMij、m.j＝Σimij、m..＝Σi,jmijにそれぞれ固定する。そして、式（７）を変化させてその最小値（左のパラメータがブレグマン発散の凸にある）を求め、以下の式（８）を満たすκを求める。
【００５１】
【数８】
2010157214.tif 000010
【００５２】
これにより、κは以下の式（９）によって示すことができ、簡便な２値検索によって任意の確度でκの近似値を得ることができる。２値検索によればデータ処理量を顕著に抑制することができるため、迅速な処理が可能となる。
【００５３】
【数９】
2010157214.tif 000011
【００５４】
先に説明したDSAアルゴリズムは、当初のマニホールド(manifold)と、次に説明するクラスタリング予測との間で、中間的な段階の行列Mを得るものである。前者の多様体の幾何学的な表示を、可能な限り、後者のクラスタリング予測の表示に近似させるため、DSAアルゴリズムの反復回数はできるだけ少ないことが望ましい。述べたように、DSSアルゴリズムは幾何学的な表示を変化させないが、DSAアルゴリズムを用いた標準化に先立ってDSSアルゴリズムによる標準化を行うことで、DSAアルゴリズムの反復回数を顕著に少なくすることができる。
【００５５】
４．クラスタリング
（4-1）ソフトクラスタリング
ソフトメンバーシップクラスタリングのためのアルゴリズムには、EM（Expectation Maximization）アルゴリズムとCP（Complete Positive factorization）アルゴリズムを用いて行うことができる。EMアルゴリズムはパラメータ密度を最適化し、CPアルゴリズムは情報量のロスを最小限に抑えるために寄与する。
【００５６】
このステップは、本発明に係る遺伝子クラスタリングプログラムの（４）のステップ「ブール行列Ｎに基づいて各データをクラスタリングするステップ」に対応するステップである（図1中、S4参照）。
【００５７】
まず、標準化のステップで得られた行列Nにおいて、クラスターナンバーkを固定し、以下の式（１０）の解を求める。行列Gは、各遺伝子がクラスターkのソフトメンバーシップである確率を与える。
【００５８】
【数１０】
2010157214.tif 000012
【００５９】
ソフトメンバーシップでは１つの遺伝子が２以上のクラスターに分類され得る。このソフトメンバーシップは、１つの遺伝子が１つのクラスターのみに属するハードメンバーシップとして変換することができる。
【００６０】
（4-2）ハードクラスタリング
ハードメンバーシップクラスタリングのためのアルゴリズムには、ブレグマンK平均（Bregman k-means: BkM）法と、階層型クラスタリング（Hierarchical clustering：HC）アルゴリズムと、AP（Affinity Propagation）アルゴリズムと、を用いる。BkM法は、K平均法を一般化して、指数型分布族のメンバーシップへ適用可能としたものである。また、HCは、集積的なクラスタリングアルゴリズムである。さらに、APアルゴリズムは、クラスターナンバーを予め設定することなく、類似度行列に基づいて、空間中に格ノードが中心的なノード候補（exemplar）を選ぶ確率（responsibility：当選確率）と、中心ノードが他のノードを自分のグループに属させる確率（availability：属性確率）と、を元にクラスターへの収束をさせるアルゴリズムである。以下、「確率」を「確率メッセージ」という場合がある（「Clustering by Passing Messages Between Data Points. Science 315: 972-976, 2007」参照）。
【００６１】
EMとBkMとに用いられる有力な初期化手法を統合したBAV（Bregman-Arthur-Vassilvitskiiinitialization）アルゴリズムにより、最適条件に近いハードクラスタリングを可能にすることができる（「k-Means++: the advantages of careful seeding. ACM-SIAM Symposium on Discrete Algorithms, 2007」参照）。
【００６２】
BAVアルゴリズムにおいては、まず、クラスターの中心をランダムに選択し、プロトタイプを作成する。そして、他の中心候補に対しては、確率p(x)で求める中心xに存する偏向分布を用いる。この確率p(x)は、以下の式（１１）及び式（１２）によって表される。ここで、プロトタイプとは、幾何学的構造のために各遺伝子が最も意味のある座標軸をとる状態を意味する。
【００６３】
【数１１】
2010157214.tif 000013
【００６４】
【数１２】
2010157214.tif 000014
【００６５】
この式（１２）は、クラスターの中心Cとして選択された遺伝子の１つと遺伝子Xとの間の最小発散を表す。ベルグマン発散のジェネレーターは任意とすることができる。特に、ブレグマン発散として、２乗ユークリッド距離を選択した場合は、Arthur-Vassilvitskii初期化に回帰される。
【００６６】
その後、通常のk-Means法と同様にBkMによる処理を実行し、各遺伝子をクラスターの中心に割り振って新たな中心の計算を行う処理を反復する。新たな中心の計算は通常のk-Means法のアルゴリズム（クラスターメンバーの相加平均を新しい中心とする）に従えばよいが、各遺伝子の再割り振りはより一般的な規則に従って行われる。すなわち、遺伝子Xが下記式（１３）を満たすCの解である設定された中心cに関連しているものとして行われる。
【００６７】
【数１３】
2010157214.tif 000015
【００６８】
HCアルゴリズムは、距離の二乗和の中の増分（Ward linkage distance）と最短距離（Single linkage distance）とを両方使って実行される。
【００６９】
APアルゴリズムにおいて、アルゴリズムのインプットは遺伝子間の類似度行列Mである。この類似度行列Mには、負の数を含んでいてもよい。行列Mの横列（row）i、縦列（column）jの入力は、以下の式（１４）によって定義される。
【００７０】
【数１４】
2010157214.tif 000016
【００７１】
ここで、相似性の代わりにベルグマン分散を用いる点については、上記の「Clustering by Passing Messages Between Data Points. Science 315: 972-976, 2007」も言及していな点に注意を要する。
【００７２】
APアルゴリズムにおいて、各クラスターの最も平均的なプロファイルを示す遺伝子群をモデルとして定義される。APアルゴリズムは、遺伝子間での属性確率メッセージと当選確率メッセージの交換を反復することにより実行される。遺伝子iからモデル候補となる遺伝子jへ送られる属性確率メッセージは下記式（１５）によって定義できる。
【００７３】
【数１５】
2010157214.tif 000017
【００７４】
初回の反復処理では属性確率メッセージは０に設定され、当選確率メッセージが更新されると属性確率メッセージは次の２つの規則に従って更新される。
【００７５】
まず、jと異なるiについては、モデル候補となる遺伝子jから遺伝子iへ送られる属性確率メッセージは次の式（１６）によって示される。
【００７６】
【数１６】
2010157214.tif 000018
【００７７】
また、遺伝子jの自己属性確率メッセージは、次の式（１７）によって算出される。
【００７８】
【数１７】
2010157214.tif 000019
【００７９】
アルゴリズムの収束を促進するため、各メッセージは、新しい値がそれ以前の値のλ倍に更新値の（1−λ）倍を加えたものに等しくなるように抑制される。アルゴリズムは、m回の反復処理の間にモデルが変化しなくなった時点で終了する。
【００８０】
クラスターの計算は、各遺伝子ペアのaij＋rijに基づいて行われる。遺伝子iについて、この合計を最大にするjの値は、i＝jである場合においてiをモデルとして特定するか、又は遺伝子iのモデルとなる遺伝子を特定する。各クラスターは、モデルの遺伝子と、これをモデルとする遺伝子群に対応する。
【００８１】
本発明に係る遺伝子クラスタリングプログラムにおいて、最終的なクラスタリング結果が得られるまでのデータ処理を、図９に概念図としてまとめた。
【００８２】
以上に説明した各種の処理を行う遺伝子クラスタリングプログラムをユーザに提供する提供媒体としては、磁気ディスク、CD-ROM、固体メモリなどの記録媒体の他、ネットワーク、衛星などの通信媒体を利用することができる。
【００８３】
本発明に係る遺伝子クラスター解析装置は、（１）遺伝子の経時的な発現量変化を示すデータから、各データ間の類似度を反映した特徴量を算出する手段と、（２）算出された特徴量から、全ての遺伝子間の組合せについて類似度行列Ｍの固有ベクトルを算出する手段と、（３）類似度行列Ｍを、固有ベクトルの固有値を維持したまま、ブール行列Ｎに変換する手段と、（４）ブール行列Ｎに基づいて各データをクラスタリングする手段と、を少なくとも備えるものである。この遺伝子クラスター解析装置は、上記の遺伝子クラスタリングプログラムを通常のコンピュータにインストールして構成することができる。
【００８４】
図１０は、本発明に係る遺伝子クラスター解析装置の構成例を示すブロック図である。遺伝子クラスター解析装置１において、内部バス１０は、例えばPCI（Peripheral Component Interconnect）またはローカルバス等により構成され、CPU１１、ROM１２、RAM１３、およびインタフェース１４を相互に接続している。各部は、この内部バス１０を介してデータの授受を行う。CPU１１は、ROM１２に記憶されている遺伝子クラスタリングプログラムに従って処理を実行する。RAM１３には、CPU１１が各種の処理を実行する上において必要なデータやプログラム等が適宜記憶される。インタフェース１４には、キーボード１５とマウス１６が接続されており、ユーザは、これらを用いてパラメータ等の設定を行うことができる。インタフェース１４は、これらから出力された操作信号をCPU１１に出力する。また、インタフェース１４には、モニタ１７とハードディスク１８が接続されている。モニタ１７は、CPU１１に制御され、所定の画像を表示する。CPU１１は、ハードディスク１８に対して、インタフェース１４を介してデータまたはプログラム等の記録または読み出しを行うことができる。
【００８５】
本発明に係る遺伝子クラスタリング方法は、上述の遺伝子クラスタリングプログラムや遺伝子クラスター解析装置によって実施が可能であり、（１）遺伝子の経時的な発現量変化を示すデータから、各データ間の類似度を反映した特徴量を算出するステップと、（２）算出された特徴量から、全ての遺伝子間の組合せについて類似度行列Ｍの固有ベクトルを算出するステップと、（３）類似度行列Ｍを、固有ベクトルの固有値を維持したまま、ブール行列Ｎに変換するステップと、（４）ブール行列Ｎに基づいて各データをクラスタリングするステップと、を少なくとも行うものである。
【００８６】
この遺伝子クラスタリング方法によれば、最適化された複数のアルゴリズムを用いることで、多次元のデータを、遺伝子発現パターンの類似性に関する生物学的情報を維持したまま、低次元の多様体に変換することができる。従って、先験的なデータ予測を行うことなく、高度なクラスタリング解析を迅速に行って、正確なクラスタリング構造を検出することが可能である。
【産業上の利用可能性】
【００８７】
本発明に係る遺伝子クラスタリングプログラムは、遺伝子の経時的発現量データに基づく遺伝子クラスタリングを、先験的なデータ予測を行うことなく、高精度に行うことができるため、ボウタイネットワーク等の細胞内シグナルネットワーク解析のために有効に用いることができる。
【図１】
2010157214.tif 000020
【図２】
2010157214.tif 000021
【図３】
2010157214.tif 000022
【図４】
2010157214.tif 000023
【図５】
2010157214.tif 000024
【図６】
2010157214.tif 000025
【図７】
2010157214.tif 000026
【図８】
2010157214.tif 000027
【図９】
2010157214.tif 000028
【図１０】
2010157214.tif 000029