JP6029683B2 - データ解析装置、データ解析プログラム - Google Patents
データ解析装置、データ解析プログラム Download PDFInfo
- Publication number
- JP6029683B2 JP6029683B2 JP2014548349A JP2014548349A JP6029683B2 JP 6029683 B2 JP6029683 B2 JP 6029683B2 JP 2014548349 A JP2014548349 A JP 2014548349A JP 2014548349 A JP2014548349 A JP 2014548349A JP 6029683 B2 JP6029683 B2 JP 6029683B2
- Authority
- JP
- Japan
- Prior art keywords
- cluster
- data
- clustering
- sample data
- likelihood
- 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.)
- Expired - Fee Related
Links
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F16/00—Information retrieval; Database structures therefor; File system structures therefor
- G06F16/20—Information retrieval; Database structures therefor; File system structures therefor of structured data, e.g. relational data
- G06F16/23—Updating
- G06F16/2365—Ensuring data consistency and integrity
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F16/00—Information retrieval; Database structures therefor; File system structures therefor
- G06F16/20—Information retrieval; Database structures therefor; File system structures therefor of structured data, e.g. relational data
- G06F16/28—Databases characterised by their database models, e.g. relational or object models
- G06F16/284—Relational databases
- G06F16/285—Clustering or classification
-
- 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
-
- 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/30—Unsupervised 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
- G16B50/00—ICT programming tools or database systems specially adapted for bioinformatics
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Databases & Information Systems (AREA)
- Theoretical Computer Science (AREA)
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Data Mining & Analysis (AREA)
- Medical Informatics (AREA)
- Bioethics (AREA)
- Bioinformatics & Computational Biology (AREA)
- Spectroscopy & Molecular Physics (AREA)
- General Health & Medical Sciences (AREA)
- Biophysics (AREA)
- Evolutionary Biology (AREA)
- Biotechnology (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Public Health (AREA)
- Software Systems (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Evolutionary Computation (AREA)
- Epidemiology (AREA)
- Artificial Intelligence (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Computer Security & Cryptography (AREA)
- Information Retrieval, Db Structures And Fs Structures Therefor (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Description
本発明は、サンプルデータをクラスタリング解析する技術に関するものである。
再生医療、遺伝子を用いた診断、あるいは生命現象の基本的な理解のため、組織の平均としての遺伝子の発現量ばかりでなく、組織を構成する1つ1つの細胞の中味を定量的に分析することが重要視され始めている。このような個々の細胞の特性を1つずつを解析することを単一細胞解析と呼んでいる。
単一細胞中の生体分子の量が極めて微量であるため、単一細胞解析は、これまで細胞膜上の蛋白質など一部の生体分子を対象とする解析にしか用いられてこなかった。しかし近年では、技術の発展により単一細胞中の微量な生体分子を定量評価できるようになりつつある。
下記非特許文献1には、単一細胞中の遺伝子発現解析に関して、定量PCR装置を用いた方法によって、十分な精度で特定の遺伝子発現量を計測できる方法が開示されている。下記非特許文献2には、同様に単一細胞中に遺伝子発現解析に関して、大規模DNAシーケンサ(次世代シーケンサ)を用いて、ほぼすべての遺伝子の発現を定量化する方法が開示されている。非特許文献2には、細胞の種類を特定するためのデータ解析方法も開示されている。今後、ゲノム配列、細胞内の蛋白質、さらに細胞内の様々な生体分子が単一細胞レベルで同定されるようになることが予想される。
Nature Method,Vol.6,No.7(2009),pp.503
Genome Research,Vol.21,No.7(2011),pp.1088
上記のような単一細胞解析技術の進歩によって、これまで均一であると仮定され解析されてきた細胞組織は、単一細胞解析で得られるデータを用いて、これまで知られているよりも詳細なグループ、下位組織を構成していることが明らかとなるであろう。その結果、非常に多数の細胞から構成されている人体などの個体中の複雑な生命現象は、細胞のデータによって分類された細胞のグループで構成され、このグループ間を様々な生化学的なシグナルがやりとりされるネットワークとして生命が把握され、生命科学分野、特に医療あるいは創薬分野において大きなインパクトを与えると考えられる。
たとえば、これまで均一であると考えられてきたがん組織をグループ分けし、各グループに対して遺伝子変異を解析することによって、より適切な分子診断薬を選択することができるようになる可能性がある。また、血中の免疫細胞中の遺伝子発現量を解析することにより、様々な病気を診断することができる可能性が示唆されているが、免疫細胞の詳細な分類によって、より精度の高い診断ができるようになると考えられる。
しかし、データのみを用いて細胞を分類するアルゴリズムおよびこれを応用した解析/診断装置は、細胞を分類して医療向けの診断において用いるために必要な特性が必ずしも十分ではない。ここでいう必要な特性の例としては、例えば細胞のグループ化(以下ではクラスタリングと呼ぶ)において、最適なグループ数(以下ではクラスタ数と呼ぶ)が事前に分かっていない場合であっても、適切なクラスタリングを実施することなどが考えられる。特に、データ数の少ない例外的なクラスタを独立のクラスタとするのか、他のデータ数の多いクラスタの一部とするのかを判定することは、従来の解析/診断装置においては困難である。
本発明は、上記のような課題に鑑みてなされたものであり、実験誤差などに起因する例外データが生じる場合においても、クラスタリングを適切に実施することができるデータ解析装置を提供することを目的とする。
本発明に係るデータ解析装置は、実験誤差データが記述している実験誤差の範囲に応じて、クラスタ境界を緩和するためのクラスタ範囲パラメータをあらかじめ定めておく。クラスタリングの過程において、いずれのクラスタにも属さない例外データについては、例外データからさらにクラスタ範囲パラメータによって定められる距離だけ離れた領域がいずれかのクラスタに含まれる場合は、例外データがそのクラスタに属するものと判定し、前記距離だけ離れてなおいずれのクラスタにも含まれない場合は、独立したクラスタを構成するものと判定する。
本発明に係るデータ解析装置によれば、単一細胞の解析結果を用いて細胞を適切に分類することができる。また、分類された細胞の種類数を精度よく決定することができる。
<従来のデータ解析手法>
以下では本発明の理解を促進するため、まず初めに従来のデータ解析手法およびその課題について、具体例を挙げて説明する。その後に、本発明に係るデータ解析装置の具体的な構成について説明する。
以下では本発明の理解を促進するため、まず初めに従来のデータ解析手法およびその課題について、具体例を挙げて説明する。その後に、本発明に係るデータ解析装置の具体的な構成について説明する。
単一細胞解析で得られた細胞データは、たとえば主因子分析法を用いて解析することができる。主因子分析法によって得られたデータは、目視によってグループを判定するために用いられることが多い。従来のクラスタリング方法の課題を具体的に見るため、以下では模擬サンプルデータを用いることにする。具体的には、4つの遺伝子、180個の細胞について単一細胞解析を実施したことを想定し、コンピュータ上で乱数を用いて細胞データを生成した。
図1は、模擬サンプルデータのクラスタリング結果を示す図である。図1(a)は模擬データ内の各クラスタの構成を示す。εは乱数であり、平均値0標準偏差σの分布に従うものとした。クラス2とクラスタ3の間の距離を制御するためのパラメータαを設定し、α=6×σに設定した。乱数εは4次元空間において等方的なガウス分布に従うと仮定しており、その片側標準偏差を上記σとした。
図1(b)は、図1(a)に示す模擬データについて主因子分析を実施し、上位2個の主因子について可視化したデータである。第1主因子に対応する軸をPC1とし、第2主因子に対応する軸をPC2とした。各データ集合に添えている数字は、図1(a)に記載している各クラスタ番号である。図1(b)に示すように、6つのクラスタを容易に目視確認することができる。しかし、主因子分析法はクラスタリングを実施するためのアルゴリズムではないため、目視確認以外の手段によって明確にクラスタリング結果を与えるものではない。また、クラスタリング結果の信頼度も定量的に与えられるものではない。
クラスタリング結果の信頼度は、複数の被験者から得た細胞データを比較する際に用いることができる。例えば、健常な被験者から採取した細胞データのクラスタリング結果が信頼できる場合、そのクラスタリング結果を基準として他の被験者の細胞データを検証することができる。しかしクラスタリング結果の信頼度が不明であれば、これを基準として用いることが適切であるのか否か判断することができない。したがって、細胞解析においてクラスタリング結果の信頼度を得ることは、非常に有用である。
図2は、階層的クラスタリング法について説明する図である。階層的クラスタリング法は、主因子分析法と同様にデータを分類する手法としてよく用いられる。ここでは、図1と同じデータを階層的クラスタリング法によってクラスタリングした結果を示す。各細胞のデータは4次元ユークリッド空間内のベクトルとして表現されているため、各細胞間の距離はユークリッド距離を用いて評価した。
階層的クラスタリング法においては、データ間の距離が最も近い2つのデータをペアにして、その距離に応じた高さのトーナメント図に現れる矩形を設定する。以下、ペアの平均値の位置に2つのデータを代表するデータがあるものと仮定して次のデータに進み、すべてのデータがペアとして結びつけられるまで同様の処理を繰り返す。トーナメント図における矩形の高さが高いほど、クラスタ間の距離が大きいことを示す。ある高さにおいてトーナメント図を水平方向に切断した時に得られる長い鉛直線との交点の数がクラスタ数に相当すると考えられる。
しかし、階層的クラスタリング法においては、最適なクラスタ数は求めることができない。図2に示す例においては、5個のクラスタあるいは6個のクラスタいずれも妥当なクラスタ数であるように思われる。すなわち、階層的クラスタリング法によって得られるクラスタ数は、主因子分析法において目視確認により得られるクラスタ数とは必ずしも一致しない。したがって、階層的クラスタリング法は、最適なクラスタ数およびその時のクラスタリングの信頼度を求めることは困難である。
クラスタリングの信頼度を評価するためには、サンプルデータを確率分布へフィッティングさせることが最も自然である。確率分布へフィッティングさせる手法としては、混合ガウシアンモデルと呼ばれる方法が最もよく知られている。混合ガウシアンモデルにおいては、各細胞データはガウス分布していると仮定し、各細胞データはいずれかのクラスタに所属していると見なす。次に、ガウシアン確率密度関数に対する対数尤度を計算し、最尤法によってクラスタ数、クラスタ位置(平均値)、分布(標準偏差)を決定する。
一般に、単純に対数尤度を用いてクラスタ数などを決定する場合においては、クラスタ数とデータ数が一致するまでクラスタ総数を増やせば良いフィッティングが得られる。したがって、単一細胞解析のように最適なクラスタ総数を決定したい場合においては不向きである。
図3は、赤池情報量基準を図1と同じ模擬サンプルデータに対して適用した結果を示す図である。上記のような混合ガウシアンモデルにおける問題を回避するため、各種情報量規準を混合ガウシアンモデルに対して適用することが考えられる。赤池情報量基準は、基本的な情報量基準としてよく知られている。
赤池情報量基準を最適化において用いる場合は、多くの情報量を与えて最適化を実施するほど多くのペナルティ値を課す。これを混合ガウシアンモデルに対して適用する場合、クラスタ数を増やすとその分だけペナルティを課すので、必ずしもクラスタ数を増やすほど良いフィッティングが得られるとは限らないことになる。したがって、評価値が極値を有するクラスタ数が最適であろうと推測することができる。
MatLab2008aが搭載しているEMアルゴリズムを用いて、赤池情報量基準を混合ガウシアンモデルに対して適用した最適化計算を実施したところ、図3に示すように明確な極値を得ることはできなかった。したがって、混合ガウシアンモデル、および赤池情報量基準をこれに適用した修正手法のいずれも、最適なクラスタ数を決定することは困難であることが分かった。
その他のクラスタリングアルゴリズムとしては、サポートベクトルマシン法、k−means法などがある。しかしこれらの手法は、事前にクラスタ数が分からないデータに対して適用しても、有効な結果を得ることは困難である。仮に、これら手法によって最適なクラスタ数が得られたとしても、そのクラスタリングの信頼度を数値によって評価することはやはり困難である。
事前に情報を与える必要がないデータ解析手法として、多くのデータマイニング手法が知られている。たとえば、非特許文献2において採用されている自己組織化マップ(Self Organizing Maps)がある。しかし、自己組織化マップを用いてクラスタリングを実施したとしても、クラスタリング結果の信頼度を得ることができない。
上記課題の他にも、単一細胞解析において得られるサンプルデータは実験誤差が無視できないという課題がある。実験誤差を多く含むサンプルデータは、誤差の少ないサンプルデータのクラスタリング結果から外れてしまうため、いずれのクラスタに属するのか、あるいはそもそも独立のクラスタとみなすべきであるのかを判断することが難しい。したがって、実験誤差を含むサンプルデータについて、有意な分解能でクラスタリングを実行することも、細胞解析/診断装置として重要であると考えられる。
<実施の形態1>
以上、従来のデータ解析手法およびその課題について説明した。以下では本発明の実施形態1に係るデータ解析装置について説明する。単一細胞ごとの生体分子に関する解析、特に遺伝子発現解析において得られるデータは、各細胞の生体分子の定量値を要素値とする行列を用いて表現することができる。個々の細胞に対する各遺伝子の発現量のサンプルデータは、遺伝子数をn、細胞数をmとすると、m行n列の行列で表示できる。以下ではこの形式によって記述されたサンプルデータを前提とする。
以上、従来のデータ解析手法およびその課題について説明した。以下では本発明の実施形態1に係るデータ解析装置について説明する。単一細胞ごとの生体分子に関する解析、特に遺伝子発現解析において得られるデータは、各細胞の生体分子の定量値を要素値とする行列を用いて表現することができる。個々の細胞に対する各遺伝子の発現量のサンプルデータは、遺伝子数をn、細胞数をmとすると、m行n列の行列で表示できる。以下ではこの形式によって記述されたサンプルデータを前提とする。
図4は、本実施形態1に係るデータ解析装置がサンプルデータをクラスタリングする処理の概略フローを示す図である。以下、図4に示す各ステップを説明する。各ステップの詳細については後述の図5〜図6を用いて改めて説明する。
(図4:ステップS401〜S402)
データ解析装置は、図1で説明した形式のサンプルデータを取得する(S401)。データ解析装置は、サンプルデータを取得したときの実験誤差に関するデータを取得する(S402)。これらのデータは、オペレータがデータ解析装置へ入力してもよい。実験誤差に関するデータとは、サンプルデータが実験誤差に起因してどの程度ばらついているかを数値的に示すデータである。たとえば、各遺伝子データの標準偏差ベクトル値を実験誤差データとすることができる。
データ解析装置は、図1で説明した形式のサンプルデータを取得する(S401)。データ解析装置は、サンプルデータを取得したときの実験誤差に関するデータを取得する(S402)。これらのデータは、オペレータがデータ解析装置へ入力してもよい。実験誤差に関するデータとは、サンプルデータが実験誤差に起因してどの程度ばらついているかを数値的に示すデータである。たとえば、各遺伝子データの標準偏差ベクトル値を実験誤差データとすることができる。
(図4:ステップS403)
データ解析装置は、実験誤差に起因していずれのクラスタにも属さなくなっている例外データをいずれかのクラスタに包含させるため、以後のクラスタリング処理においてクラスタ境界を緩和する。具体的には、クラスタリング空間において例外データから所定範囲離れた位置にいずれかのクラスタが存在する場合は、その例外データは当該クラスタに属するものとみなす。上記所定範囲のことを、本発明においてはCR(Clustering Resolution)値と呼ぶことにする。例外データは実験誤差によって生じると考えられるため、CR値は実験誤差データが記述している実験誤差以上の値として設定する。例えば実験誤差データを誤差の標準偏差σによって表す場合、CR値はσ〜4σ程度の値とすることができる。
データ解析装置は、実験誤差に起因していずれのクラスタにも属さなくなっている例外データをいずれかのクラスタに包含させるため、以後のクラスタリング処理においてクラスタ境界を緩和する。具体的には、クラスタリング空間において例外データから所定範囲離れた位置にいずれかのクラスタが存在する場合は、その例外データは当該クラスタに属するものとみなす。上記所定範囲のことを、本発明においてはCR(Clustering Resolution)値と呼ぶことにする。例外データは実験誤差によって生じると考えられるため、CR値は実験誤差データが記述している実験誤差以上の値として設定する。例えば実験誤差データを誤差の標準偏差σによって表す場合、CR値はσ〜4σ程度の値とすることができる。
(図4:ステップS404)
データ解析装置は、CR値を掃引しながら各CR値について以下のステップS405を実施する。CR値の掃引範囲としては、実験誤差データを誤差の標準偏差σによって表す場合、例えば上述のようにσ〜4σなどとすればよい。
データ解析装置は、CR値を掃引しながら各CR値について以下のステップS405を実施する。CR値の掃引範囲としては、実験誤差データを誤差の標準偏差σによって表す場合、例えば上述のようにσ〜4σなどとすればよい。
(図4:ステップS405)
データ解析装置は、クラスタ数kを2〜予想最大値まで仮設定して実際にクラスタリングを実施しながら、各クラスタリング結果についてクラスタ数の最適度を評価する。具体的には、サンプルデータが各クラスタに属する確率の対数尤度および属さない確率の対数尤度を用いて、現在のクラスタ数の尤もらしさを評価する。
データ解析装置は、クラスタ数kを2〜予想最大値まで仮設定して実際にクラスタリングを実施しながら、各クラスタリング結果についてクラスタ数の最適度を評価する。具体的には、サンプルデータが各クラスタに属する確率の対数尤度および属さない確率の対数尤度を用いて、現在のクラスタ数の尤もらしさを評価する。
(図4:ステップS406)
データ解析装置は、ステップS405において求めたクラスタ数の尤もらしさが極値を取るとき、そのクラスタ数が最適であるものとみなし、そのクラスタ数を最終的なクラスタ数として採用する。
データ解析装置は、ステップS405において求めたクラスタ数の尤もらしさが極値を取るとき、そのクラスタ数が最適であるものとみなし、そのクラスタ数を最終的なクラスタ数として採用する。
(図4:ステップS407)
データ解析装置は、ステップS406において決定したクラスタ数に基づくクラスタリング結果とともに、そのクラスタリング結果の信頼度を出力する。クラスタリング結果の信頼度としては、ステップS405において求めたクラスタ数の尤もらしさの値を用いることができる。
データ解析装置は、ステップS406において決定したクラスタ数に基づくクラスタリング結果とともに、そのクラスタリング結果の信頼度を出力する。クラスタリング結果の信頼度としては、ステップS405において求めたクラスタ数の尤もらしさの値を用いることができる。
図5は、ステップS405の詳細を説明する処理フロー図である。以下、図5の各ステップについて説明する。
(図5:ステップS501)
データ解析装置は、与えられた仮のクラスタ数kについて、サンプルデータを暫定的にクラスタリングする。本ステップにおけるクラスタリング手法は、例えば階層的クラスタリング法やk−means法など任意の手法でよい。
データ解析装置は、与えられた仮のクラスタ数kについて、サンプルデータを暫定的にクラスタリングする。本ステップにおけるクラスタリング手法は、例えば階層的クラスタリング法やk−means法など任意の手法でよい。
(図5:ステップS502)
データ解析装置は、サンプルデータをクラスタリングしたデータセットをk個分複製する。複製された各データセットは、それぞれ以下のステップにおいて各クラスタリング結果の尤もらしさを評価するために用いる。データ解析装置は、以下のステップにおいて用いるカウンタiを初期化する(i=1)。
データ解析装置は、サンプルデータをクラスタリングしたデータセットをk個分複製する。複製された各データセットは、それぞれ以下のステップにおいて各クラスタリング結果の尤もらしさを評価するために用いる。データ解析装置は、以下のステップにおいて用いるカウンタiを初期化する(i=1)。
(図5:ステップS503)
データ解析装置は、i番目(i=1〜k)のクラスタについて、i番目の複製データセットを用いて以下のステップを実施する。i番目の複製データセットに関してはi番目のクラスタのみを残し、その他データは全てi番目のクラスタ外に属するものと仮定し、i番目のクラスタ以外のクラスタについては削除する。すなわち、i番目のクラスタに属さないデータについては、全て単一のクラスタに属するものと仮定する。
データ解析装置は、i番目(i=1〜k)のクラスタについて、i番目の複製データセットを用いて以下のステップを実施する。i番目の複製データセットに関してはi番目のクラスタのみを残し、その他データは全てi番目のクラスタ外に属するものと仮定し、i番目のクラスタ以外のクラスタについては削除する。すなわち、i番目のクラスタに属さないデータについては、全て単一のクラスタに属するものと仮定する。
(図5:ステップS504)
データ解析装置は、i番目のクラスタが例外データによって構成されているものであるか否かを判定する。例外データの例については後述の図6で説明する。例外データでなければステップS505〜S507を実施し、例外データであればステップS508〜S509を実施する。
データ解析装置は、i番目のクラスタが例外データによって構成されているものであるか否かを判定する。例外データの例については後述の図6で説明する。例外データでなければステップS505〜S507を実施し、例外データであればステップS508〜S509を実施する。
(図5:ステップS504:例外データであるか否かの判定基準その1)
i番目のクラスタがその内部に十分なサンプルデータ数を保持している場合は、そのクラスタは例外データによるものではないと判定する。この場合は、クラスタに属するサンプルデータから計算される相関行列などを用いて決定したデータ構造は、信頼度が高いと認められるからである。サンプルデータ数が十分であるか否かの閾値thは、あらかじめ定めておく。例えばクラスタ内のサンプルデータ数が2個以下であれば例外データとみなす、などの判定基準が考えられる。
i番目のクラスタがその内部に十分なサンプルデータ数を保持している場合は、そのクラスタは例外データによるものではないと判定する。この場合は、クラスタに属するサンプルデータから計算される相関行列などを用いて決定したデータ構造は、信頼度が高いと認められるからである。サンプルデータ数が十分であるか否かの閾値thは、あらかじめ定めておく。例えばクラスタ内のサンプルデータ数が2個以下であれば例外データとみなす、などの判定基準が考えられる。
(図5:ステップS504:例外データであるか否かの判定基準その2)
閾値thは、例えばある範囲内でランダムに選択することもできる。あるいは、適当な確率分布を仮定して、その確率分布を前提としてランダムに選択することもできる。この場合は、確率分布の中央に位置するクラスタ数が最も選択され易くなる。確率分布のパラメータは任意に定めることもできるし、最適化計算によって望ましい確率分布を決定することもできる。
閾値thは、例えばある範囲内でランダムに選択することもできる。あるいは、適当な確率分布を仮定して、その確率分布を前提としてランダムに選択することもできる。この場合は、確率分布の中央に位置するクラスタ数が最も選択され易くなる。確率分布のパラメータは任意に定めることもできるし、最適化計算によって望ましい確率分布を決定することもできる。
(図5:ステップS505〜S507の総括)
データ解析装置は、i番目のクラスタ内に属するサンプルデータの分布を用いてそのサンプルデータの確率分布を決定する。データ解析装置は、各サンプルデータがi番目のクラスタに属する、あるいは属さないと判定した判断が正しいか否かについて、上記確率分布を用いてその妥当性を評価する。具体的な手法は以下の通りである。
データ解析装置は、i番目のクラスタ内に属するサンプルデータの分布を用いてそのサンプルデータの確率分布を決定する。データ解析装置は、各サンプルデータがi番目のクラスタに属する、あるいは属さないと判定した判断が正しいか否かについて、上記確率分布を用いてその妥当性を評価する。具体的な手法は以下の通りである。
(図5:ステップS505〜S507)
データ解析装置は、i番目のクラスタのクラスタ中心(クラスタ内に属するサンプルデータの平均)とクラスタ内のサンプルデータの標準偏差を算出し、サンプルデータを規格化する(S505)。データ解析装置は、i番目のクラスタ内のサンプルデータの相関行列の逆行列を計算する(S506)。データ解析装置は、各サンプルデータとi番目のクラスタ中心との間のマハラノビス距離を計算する(S507)。クラスタ中心からのマハラノビス距離を計算する理由については後述の図6〜図7で説明する。
データ解析装置は、i番目のクラスタのクラスタ中心(クラスタ内に属するサンプルデータの平均)とクラスタ内のサンプルデータの標準偏差を算出し、サンプルデータを規格化する(S505)。データ解析装置は、i番目のクラスタ内のサンプルデータの相関行列の逆行列を計算する(S506)。データ解析装置は、各サンプルデータとi番目のクラスタ中心との間のマハラノビス距離を計算する(S507)。クラスタ中心からのマハラノビス距離を計算する理由については後述の図6〜図7で説明する。
(図5:ステップS508〜S509)
データ解析装置は、i番目のクラスタのクラスタ中心を算出し、クラスタ中心とCR値を用いてサンプルデータを規格化する(S508)。データ解析装置は、各サンプルデータとi番目のクラスタ中心との間のユークリッド距離を計算する(S509)。クラスタ中心からのユークリッド距離を計算する理由については後述の図6〜図7で説明する。
データ解析装置は、i番目のクラスタのクラスタ中心を算出し、クラスタ中心とCR値を用いてサンプルデータを規格化する(S508)。データ解析装置は、各サンプルデータとi番目のクラスタ中心との間のユークリッド距離を計算する(S509)。クラスタ中心からのユークリッド距離を計算する理由については後述の図6〜図7で説明する。
(図5:ステップS510:ステップ1)
データ解析装置は、ステップS501においてi番目のクラスタに属さないと判定したサンプルデータについて、クラスタ中心からの距離が離れるほど当該サンプルデータがi番目のクラスタに属さない確率が高くなるような確率分布関数を用いて、当該サンプルデータがi番目のクラスタに属さない確率を算出する。同様にステップS501においてi番目のクラスタに属すると判定したサンプルデータについて、クラスタ中心からの距離が離れるほど当該サンプルデータがi番目のクラスタに属する確率が低くなるような確率分布関数を用いて、当該サンプルデータがi番目のクラスタに属する確率を算出する。例えば、前者については自由度nのx2乗分布の累積確率分布関数にしたがって確率値を計算し、後者については1から上記累積確率分布関数を引いた関数にしたがって確率値を計算する。この関数の例については後述の図6〜図7で例示する。
データ解析装置は、ステップS501においてi番目のクラスタに属さないと判定したサンプルデータについて、クラスタ中心からの距離が離れるほど当該サンプルデータがi番目のクラスタに属さない確率が高くなるような確率分布関数を用いて、当該サンプルデータがi番目のクラスタに属さない確率を算出する。同様にステップS501においてi番目のクラスタに属すると判定したサンプルデータについて、クラスタ中心からの距離が離れるほど当該サンプルデータがi番目のクラスタに属する確率が低くなるような確率分布関数を用いて、当該サンプルデータがi番目のクラスタに属する確率を算出する。例えば、前者については自由度nのx2乗分布の累積確率分布関数にしたがって確率値を計算し、後者については1から上記累積確率分布関数を引いた関数にしたがって確率値を計算する。この関数の例については後述の図6〜図7で例示する。
(図5:ステップS510:ステップ2)
データ解析装置は、上記関数にしたがって計算した確率値の対数尤度を計算することにより、各サンプルデータがi番目のクラスタに属することの尤もらしさ、または属さないことの尤もらしさを計算する。すべてのサンプルデータおよびすべてのクラスタについてこの対数尤度の和を求め、クラスタ数kで除算することにより、現在のクラスタ数kの尤もらしさを計算する。本ステップの次は、i+1番目のクラスタについてステップS503から同様の処理を実施する。
データ解析装置は、上記関数にしたがって計算した確率値の対数尤度を計算することにより、各サンプルデータがi番目のクラスタに属することの尤もらしさ、または属さないことの尤もらしさを計算する。すべてのサンプルデータおよびすべてのクラスタについてこの対数尤度の和を求め、クラスタ数kで除算することにより、現在のクラスタ数kの尤もらしさを計算する。本ステップの次は、i+1番目のクラスタについてステップS503から同様の処理を実施する。
(図5:ステップS510:補足)
対数尤度の評価に用いているのはクラスタリングが妥当である確率の評価値であるため、最適パラメータが得られた時の対数尤度の値から、クラスタリングの信頼度を確率値で出力することが可能である。
対数尤度の評価に用いているのはクラスタリングが妥当である確率の評価値であるため、最適パラメータが得られた時の対数尤度の値から、クラスタリングの信頼度を確率値で出力することが可能である。
(図5:ステップS511:最適化ループおよびクラスタ数掃引ループ)
データ解析装置は、例えばモンテカルロ法などのような最適化手法を用いて、ステップS510で計算した対数尤度が小さくなるようにクラスタリングを繰り返し実施することにより、最適なクラスタリング結果と最適なクラスタ数を決定する。最適化手法として例えばモンテカルロ法を用いる場合は、各クラスタに属するサンプルデータをランダムに入れ替えながらステップS503から同様の処理を実施する。最適化ループの終了条件は、例えばステップS510で計算した現在のクラスタ数kの尤もらしさが所定閾値に達した時点などとすればよい。最適化ループを完了した後はクラスタ数kをインクリメントしてステップS501に戻り、同様の処理を実施する。
データ解析装置は、例えばモンテカルロ法などのような最適化手法を用いて、ステップS510で計算した対数尤度が小さくなるようにクラスタリングを繰り返し実施することにより、最適なクラスタリング結果と最適なクラスタ数を決定する。最適化手法として例えばモンテカルロ法を用いる場合は、各クラスタに属するサンプルデータをランダムに入れ替えながらステップS503から同様の処理を実施する。最適化ループの終了条件は、例えばステップS510で計算した現在のクラスタ数kの尤もらしさが所定閾値に達した時点などとすればよい。最適化ループを完了した後はクラスタ数kをインクリメントしてステップS501に戻り、同様の処理を実施する。
(図5:ステップS511:CR値掃引ループ)
データ解析装置は、現在のCR値について最適化ループおよびクラスタ数掃引ループを終了した後は、CR値をインクリメントしてステップS501に戻り、同様の処理を実施する。CR値をインクリメントする幅は、想定しているCR値の最小値と最大値の差に応じて適宜定める。
データ解析装置は、現在のCR値について最適化ループおよびクラスタ数掃引ループを終了した後は、CR値をインクリメントしてステップS501に戻り、同様の処理を実施する。CR値をインクリメントする幅は、想定しているCR値の最小値と最大値の差に応じて適宜定める。
図6は、ステップS505〜S507の処理イメージを示す図である。ここでは図1(b)と同様に、2つの主因子に対応する軸を設定したクラスタリング空間を例示した。図6に示すクラスタリング例においては、クラスタ中心に比較的多くのサンプルデータが集まっているため、当該クラスタ内に含まれるデータは実験誤差の影響を受けた例外データではないと仮定する。
クラスタの形状は必ずしもクラスタリング空間において円形ではないため、図6(a)の左に示すデータ分布から右に示すデータ分布へと線形変換を実行し、変換後のクラスタ形状においてクラスタ中心から各サンプルデータまでの距離を求める。これは、各データのクラスタの中心からのマハラノビス距離を求めることに相当する。換言すると、クラスタがマハラノビス田口法における単位空間を形成していると仮定することに対応する。
ステップS505〜S507においては、i番目のクラスタに属さないと仮判定したサンプルデータ(例えば図6(a)の左上のデータ群)については、図6(b)に示すx2累積確率分布を用いて、そのサンプルデータがi番目のクラスタに属さない確率を計算する。i番目のクラスタに属すると仮判定したサンプルデータ(図6(a)の円内のデータ群)については、図6(b)に示す1−x2累積確率分布を用いて、そのサンプルデータがi番目のクラスタに属する確率を計算する。これにより、i番目のクラスタに属すると仮判定したサンプルデータについては、クラスタ中心に近いほど確率値が高くなり、i番目のクラスタに属さないと仮判定したサンプルデータについては、クラスタ中心から離れるほど確率値が高くなる。
図7は、ステップS508〜S509の処理イメージを示す図である。ここでは図1(b)と同様に、2つの主因子に対応する軸を設定したクラスタリング空間を例示した。図6に示すクラスタリング例においては、中央のクラスタに属するデータ数が少ないため(例外1においては2個、例外2においては1個)、例外データであると仮判定される。これはすなわち、i番目のクラスタに属するデータ数が少ないため、そのクラスタ内のデータ構造を十分に決定することができない場合に相当する。
ステップS508〜S509においては、CR値を当該クラスタのサイズであると仮定する。i番目のクラスタの周辺に存在する、同クラスタに属さないと判定されたデータ個数が所定個数未満である場合(例外1)、例外データは独立したクラスタを構成している可能性が高い。この場合、例外クラスタに属さないサンプルデータはその確率値が大きくなる。結果として、例外データが独立したクラスタであることの尤もらしさが高く評価されることになる。
一方、i番目のクラスタの周辺に、同クラスタに属さないと判定されたデータ個数が所定個数以上存在する場合(例外2)、例外データは本来他のクラスタに属していたが実験誤差によってクラスタから外れた可能性が高い。この場合、近傍のクラスタが例外クラスタに属する確率値が相応に高くなる。結果として、例外データが独立したクラスタであることの尤もらしさが低く評価されることになる。
図7に示すように、本発明におけるCR値は、実験誤差によってクラスタから外れたサンプルデータを元のクラスタに戻す機能を備える。そのためCR値は、実験誤差よりも大きな値とし(実験誤差より小さいと誤差を修正できないため)、他の医学的、生物学的な知見から得られる制限範囲内で設定することが望ましい。
図8は、図1で説明した模擬サンプルデータについてCR値を1、2、4(それぞれ4×σ、8×σ、16×σに対応)に掃引し、クラスタ数の尤もらしさを表す対数尤度を計算した結果を示す。クラスタ数k=6において対数尤度が最小値を示し、その最小値はCR値を変えても安定であることが観察される。すなわち、クラスタ数の対数尤度が極値となるクラスタ数が最適であると判定することができる。
<実施の形態1:まとめ>
以上のように、本実施形態1に係るデータ解析装置は、実験誤差に基づき定めた、クラスタ境界を緩和するクラスタ範囲パラメータ(CR値)を用いて、いずれのクラスタにも属さないと仮判定された例外データが他のクラスタに属するか否かを判定する。これにより、実験誤差に起因する例外データが生じる場合であっても、精度よくクラスタリングを実施することができる。
以上のように、本実施形態1に係るデータ解析装置は、実験誤差に基づき定めた、クラスタ境界を緩和するクラスタ範囲パラメータ(CR値)を用いて、いずれのクラスタにも属さないと仮判定された例外データが他のクラスタに属するか否かを判定する。これにより、実験誤差に起因する例外データが生じる場合であっても、精度よくクラスタリングを実施することができる。
また、本実施形態1に係るデータ解析装置は、クラスタ中心からのマハラノビス距離またはユークリッド距離を基準として、CR値を掃引しながらクラスタリング結果の尤もらしさを評価し、その値が極値を取るクラスタ数を最適とみなす。これにより、最適なクラスタ数があらかじめ分かっていない場合であっても、良好なクラスタリング結果を得ることができる。
また、本実施形態1に係るデータ解析装置は、クラスタリング結果の信頼度をクラスタリング結果と併せて出力する。これにより、特定種類に属する細胞数やその種類に属する遺伝子発現マーカによる診断精度を向上させることができる。すなわち、従来は単一細胞解析以外の方法においては複数種類の細胞に対して医学的評価を実施していたところ、本実施形態1に係るデータ解析装置のクラスタリング結果に基づき特定集団に対するバイオマーカを評価することによって、バイオマーカによる疾患の種類や状態の判定の精度が向上するものと期待される。
本実施形態1に係るデータ解析装置によって決定したクラスタは、それぞれサンプルデータから導かれた細胞の種類に対応する。最適なクラスタ数を決定できるデータ解析装置は、ひとつの細胞種と別の細胞種との間の境界を明示できる能力を有すると考えられる。それゆえ、細胞の解析や診断で得られたサンプルデータから細胞種毎の個数を出力することができる。また、個々の細胞種において、特定の遺伝子発現量などに代表されるバイオマーカの平均値などの統計量を出力するための母集団を決定することもできる。さらに、その決定には信頼度が付与されているため、ある信頼度以下のデータを採用しないなどの判定も容易である。
本実施形態1に係るデータ解析装置は、具体的には、個々の細胞の特性を1細胞ずつ解析するバイオ・医療分野の解析・診断装置において適用することができる。特に、血球を解析対象とする解析/診断装置、または尿中の細胞を対象とする解析/診断装置、または組織切片を対象とする解析/診断装置において適用することができる。以下の実施形態についても同様である。
<実施の形態2:装置構成>
本発明の実施形態2では、実施形態1で説明したデータ解析装置の具体的な適用例として、単一細胞中の遺伝子発現解析によるデータによって細胞を分類するデータ解析装置について説明する。
本発明の実施形態2では、実施形態1で説明したデータ解析装置の具体的な適用例として、単一細胞中の遺伝子発現解析によるデータによって細胞を分類するデータ解析装置について説明する。
本実施形態2においては、単一細胞中の遺伝子発現解析を実施するため、非特許文献1に記載されている方法に基づき、磁性ビーズ上に単一細胞中のcDNAライブラリを作製し、単一細胞中の0.5pgという微量なmRNAをqPCR装置によって定量する方法を用いた。
図9は、本実施形態2に係るデータ解析装置906の機能ブロック図である。データ解析装置906は、演算部904、データ入出力部905、サンプルデータ入力部908、実験誤差データ入力部909を備える。各データを測定する測定系についても図9内に併記した。
図9は、実験誤差データを取得する動作と、クラスタリング対象となる単一細胞解析とを平行して実施する場合における測定系の構成を示した。機能ブロック901は、クラスタリング対象となるサンプルデータを取得する機能部である。機能ブロック902は、実験誤差を取得する機能部である。機能ブロック903は、両者に共通する処理を実施する機能部である。正確な実験誤差を評価するためには、実験誤差を評価するためのデータとクラスタリング対象とするサンプルデータは共通である部分が多い方が望ましい。実験誤差に関する全部または一部のデータを事前に取得している場合は、そのデータをデータベース907に保存しておき、クラスタリングを実施するとき演算部904がこれを読み出すようにしてもよい。
機能ブロック901は、細胞を1つずつ採取し、個々の細胞ごとに反応容器に分注し、この細胞を分注した反応容器に、mRNAを抽出するためのポリTプローブ付きビーズとリシスバッファを導入する。
機能ブロック902は、多数の細胞を反応容器に導入した後、機能ブロック901と同様にmRNA抽出を実施し、その後mRNA溶液を希釈し、単一細胞相当量を採取して、ポリTプローブ付ビーズ溶液にこの溶液を分注する。
機能ブロック903は、リシスバッファを除去し、逆転写酵素を含む反応溶液を分注し、逆転写反応後に反応溶液を除去し、mRNAの分解酵素を加えて洗浄後、qPCR試薬を導入し、PCRのための温度サイクルを印加しながら蛍光測定を実施することによって定量を実施する。
個々のサンプル中の遺伝子発現量は、蛍光強度がある閾値を超えたサイクル数Ct値に基づき算出される。分子数が既知の複数のDNAに対してqPCR定量を実行することによって、Ct値は分子数に読み替えられる。実験誤差は、定量を実施するまでにおける単一細胞のサンプリング以外のすべてのプロセス内の誤差を含む。定量値はガウス分布に従っているほうが望ましいので、分子数の対数をとった値をサンプルデータとして演算部904に出力する。実験誤差データについても同様である。
サンプルデータ入力部908は、機能ブロック901、903を介してサンプルデータを取得する。実験誤差データ入力部909は、機能ブロック902、903を介して実験誤差データを取得する。演算部904は、これらのデータを用いて、実施形態1で説明した処理フローにしたがってクラスタリングを実施する。データ入出力部905は、これらデータの入出力を制御する。
本実施形態2においては、仮クラスタ内の要素数が3より小さい(2以下)場合に、ステップS508〜S509を実施することとした。生物学的あるいは医学的に特定の遺伝子の発現量(サンプルデータの値)が変動すると知られている場合、その変動幅をCR値として用いてもよい。このような、遺伝子とCR値の対応関係やサンプルデータのカテゴリとCR値の対応関係をデータベース907に保存しておき、データ解析装置906に入力された遺伝子やサンプルデータのカテゴリのなかでデータベース907に保存しておいたものと合致する該当するものがあれば、演算部906はそれをデータベース907から読み出す。あるいは、実験誤差データに代えてステップS504における判定閾値をデータベース907にあらかじめ格納しておき、データ解析装置904に入力されたサンプルデータや遺伝子情報に基づいて適当な値を用いるようにしてもよい。
本実施形態2では、ステップS504における判定閾値は、クラスタ内のサンプルデータ数に対応する固定した値としたが、判定閾値を複数設けておき、各閾値をそれぞれ用いてクラスタリングを実施し、もっとも対数尤度の低い(最も尤もらしい)閾値を選択するようにしてもよい。さらには、閾値をランダムに選択するか、あるいは閾値がある確率分布にしたがう前提の下でランダムに選択してもよい。このときの確率分布関数はあらかじめ複数種類をデータベース907に格納しておき、サンプルデータの内容に応じて適当な確率分布関数を選択すればよい。
<実施の形態2:サンプリング結果>
以下では、本実施形態2に係るデータ解析装置906によるサンプリング結果の例を説明する。サンプルとしてマウスの間葉系幹細胞(C3H10T1/2)を92個用いた。すなわち、細胞を取得したマウスの骨の分化誘導状態を知るためにこのような細胞を採取した。もちろん、他の目的のための他の生体(人を含む)から、別の細胞(たとえば、血中の免疫細胞やがんの組織切片など)を採取してもよい。サンプルデータとして用いる遺伝子発現データは、Bglap1、Col1a1、Pparg、Col2a1、Eef1gの5種類の遺伝子についての定量値とした。遺伝子の種類およびその数は1例であり、可能性のある遺伝子すべてを計測対象にしてもかまわない。
以下では、本実施形態2に係るデータ解析装置906によるサンプリング結果の例を説明する。サンプルとしてマウスの間葉系幹細胞(C3H10T1/2)を92個用いた。すなわち、細胞を取得したマウスの骨の分化誘導状態を知るためにこのような細胞を採取した。もちろん、他の目的のための他の生体(人を含む)から、別の細胞(たとえば、血中の免疫細胞やがんの組織切片など)を採取してもよい。サンプルデータとして用いる遺伝子発現データは、Bglap1、Col1a1、Pparg、Col2a1、Eef1gの5種類の遺伝子についての定量値とした。遺伝子の種類およびその数は1例であり、可能性のある遺伝子すべてを計測対象にしてもかまわない。
次に、実験誤差を評価するため、5×105個程度以上の多数の細胞から抽出したmRNAを単一細胞相当量である0.5pgずつ96個分サンプリングし、5×105個の細胞の遺伝子発現量の平均値を計測するためqPCR装置を用いて定量した。定量値のばらつきはcDNAライブラリを作製するときのハンドリング誤差や、qPCR定量時の定量誤差から構成される実験誤差の総量である。この誤差は遺伝子ごとに異なる。この誤差を5次元ベクトルとして評価した。
実験誤差は以下のような方法によって取得することが望ましい。まず、多数の細胞から抽出したmRNAサンプルを作成するとき、細胞数の異なる複数のサンプルを準備し、これらのサンプルから単一細胞相当量のmRNAを採取し、cDNAを構築後、各細胞数および各遺伝子に関するqPCR定量時の誤差を定量する。そして、細胞数無限大の場合の値を外挿によって推定する。実際には、細胞数の逆数を横軸に実験誤差(標準偏差)を縦軸にとって、y切片の値を実験誤差の推定値とする。
得られた実験誤差に基づき、許容できるCR値を定める。このCR値(ベクトル値)をσとする。本実施形態2においては、実験誤差と最小CR値が一致するものとする。ただし、特定の遺伝子の発現は、同じ細胞状態であっても時間的に変化しており、実験誤差によるデータのばらつき以外に生物学的にもデータの変動幅が分かっている場合がある。これらの変動をクラスタリングにおいて用いたくない場合は、σの値をサンプルデータに基づいて一部変更しても構わない。具体的には、たとえば、特定の遺伝子について10倍程度の変動がある場合、この遺伝子のみについて、実験誤差ではなくこの変動幅の数値をCR値として用いてもよい。ただし、実験誤差が生物学的あるいは医学的変動より十分小さい場合に限られる。
図10は、本実施形態2において、クラスタ数の尤もらしさを表す対数尤度を計算した結果を示す。実験誤差をσとして、CR値は1σから2.1σまでについて評価した。すべてのCR値について、クラスタ数k=15のとき、安定的に極値が得られた。
図11は、階層的クラスタリングの結果と本実施形態2に係るデータ解析装置906によるクラスタリング結果を比較する図である。15個のクラスタは階層的クラスタリング結果に対応するトーナメント図の下に四角の箱で表現している。15個のクラスタは、5個の有意なクラスタと10個の例外クラスタから構成されることが分かった。5個の有意なクラスタは階層的クラスタリングの結果と一致する。本実施形態2により、有意なクラスタと例外クラスタの境界が明確になったことが分かる。
図12は、各クラスタ(1〜5の番号で表示)のBglap1、PParg, Col2a1の遺伝子発現量を示す図である。縦軸は、単一細胞中の各遺伝子の分子数である。横軸はクラスタ番号を示している。黒いバーは分化誘導剤あり、灰色のバーは分化誘導剤なしに対応する細胞についての遺伝子発現量である。分化誘導剤あり・なしの区別は事前に分かっている細胞情報の例である。
図12に示すように、分化誘導剤あり・なしは、個々のクラスタのなかでは細胞を区別するための重要な情報でないことが分かる。一方、3つの遺伝子の発現プロファイル、すなわち、(Bglap1,Pparg,Col2a1)の発現量を+−で表現すると、クラスタ1は(−,+,+)、クラスタ2は(+,+,+)、クラスタ3は(−,−,+)、クラスタ4は(−,+,−)、クラスタ5は(+,−,−)という明確な特徴を持っていることがわかる。同時にクラスタごとの遺伝子発現量の変動量は、上記特徴を識別するには十分でないことがわかる。事前情報との対応は各クラスタに属する細胞数に反映されるので、その意味においても正確なクラスタリングが重要であることがわかる。
図13は、図12で説明した3つの遺伝子について、クラスタリングを実施せず、分化誘導ありの場合となしの場合について発現量を示した図である。図13に示すようにBglap1については事前情報に対応した変化が見られるが、PpargおよびCol2a1については事前情報に対応した発現量の変化がなく、これら2つの遺伝子に関しては情報が得られない。以上から、クラスタリングによって、より詳細な生体の遺伝子発現に関する情報が得られることが分かる。
以上のように、本実施形態2によれば、個々の細胞中の遺伝子発現量を計測し、クラスタリングによって、細胞を構成する生体の情報を得ることができる。すなわち本実施形態2に係るデータ解析装置906は、生体中にどのような種類(クラスタ)に属する細胞がどの程度の数だけ存在するかを推定することによって生体の状態を推定するための装置である。計測対象となっている生体の健康状態が変化すると、クラスタの種類とクラスタに属する細胞の数が変化する場合において、本実施形態2は有効である。
<実施の形態3>
図14は、本発明の実施形態3に係るデータ解析装置906の機能ブロック図である。本実施形態3では、遺伝子発現の定量方法として、大規模DNAシーケンサを用いる構成例を説明する。これにより、計測対象とする遺伝子数を大幅に増やすことができる。
図14は、本発明の実施形態3に係るデータ解析装置906の機能ブロック図である。本実施形態3では、遺伝子発現の定量方法として、大規模DNAシーケンサを用いる構成例を説明する。これにより、計測対象とする遺伝子数を大幅に増やすことができる。
機能ブロック901は、クラスタリング対象サンプルをプレート上のポリTプローブ付ビーズの入った反応容器に1つずつ分注し、反応容器内で細胞を破砕し、mRNAをビーズ表面でトラップすることによって抽出する。
機能ブロック902は、実験誤差測定のため、複数の遺伝子について既知量のmRNAを個々の反応容器に分注する。既知量RNAを導入した、反応容器に対応するシーケンシングデータからその標準変偏差を算出することによって、サンプル処理とシーケンシングに関する実験誤差を定量することができる。
機能ブロック903は、逆転写反応とmRNA分解を実施する。このとき、cDNAの末端に一括増幅用プライマを導入し、次にこのプライマを用いて一括増幅を実施する。さらに、断片化し、各反応槽ごとに配列のことなる細胞認識タグをもった増幅プライマ(プライマ配列はすべての容器で共通)を用いて一括増幅を実施し、シーケンシングライブラリを構築する。シーケンシングライブラリの末端には容器ごとすなわち細胞ごとに異なる配列のタグが挿入されているので、以下の処理においてはサンプルを混合することができる。混合したサンプルを大規模シーケンサでシーケンシングするために、エマルジョンPCRやブリッジPCRなどの個別増幅を用いてシーケンシングを実行する。シーケンシングは、蛍光計測を用いる装置、FETを用いる装置、ナノポアを用いる装置などいずれの装置を用いてもよい。断片化したサンプルから得られたシーケンシングデータを既知の遺伝子配列にマッピングすることによって、どの遺伝子のどの領域の配列が計測できたかを判定する。その後、データを遺伝子毎に集計し、遺伝子毎の発現量データを算出する。このときの算出アルゴリズムは当該専門家が一般的に用いるアルゴリズムを用いてよい。その結果、細胞毎・遺伝子毎の発現量データが得られる。
クラスタリング手順については実施形態1〜2と同様であるが、計測遺伝子数が数万程度と非常に大きいため、細胞をより詳細に分類することができる。
同様の処理を用いて、各細胞のゲノムデータを解析することもできる。データ解析装置906に入力されるデータとしては、ゲノムを複数の領域に分割し、各分割領域についてカウントした変異数である。この場合の計測目的は、例えば癌の組織切片からの細胞を計測することによって、癌の発生や転移のメカニズムを解明すること、あるいは分子標的薬を選択するための診断、などが考えられる。
ゲノムデータを解析することによって得られるデータについて補足する。ゲノム全体またはその一部を単一細胞ごとに配列解析し、(mRNAの配列解析におけるデータを用いることを想定)たとえば50k塩基ごとに領域を分けて、変異を計測する。計測対象は、たとえば一塩基置換、欠失、挿入、遺伝子コピー数異常、などである。入力データはそれぞれの変異数である。実験誤差は、変異がないサンプルを人工的に作製し、それを評価することによって、評価することができる。
ゲノムを直接配列決定するためには、mRNA抽出の替わりにRNA分解、DNA抽出後、断片化、ポリAテイリングを酵素試薬の添加によって実行する必要がある。その後はmRNAの処理と同様である。
<実施の形態4>
本発明の実施形態4では、細胞中の遺伝子の発現量(分子数)を定量することによって細胞を特徴付ける代わりに、免疫染色法などで得られた細胞サンプルを蛍光顕微鏡でイメージングして、蛍光強度と分子数の対応データまたは単分子蛍光カウンティングによって細胞中の蛋白質量を定量することにより細胞を分類する構成例を説明する。
本発明の実施形態4では、細胞中の遺伝子の発現量(分子数)を定量することによって細胞を特徴付ける代わりに、免疫染色法などで得られた細胞サンプルを蛍光顕微鏡でイメージングして、蛍光強度と分子数の対応データまたは単分子蛍光カウンティングによって細胞中の蛋白質量を定量することにより細胞を分類する構成例を説明する。
本実施形態4においては、遺伝子の種類が蛋白の種類に対応する。細胞ごとの蛋白量(分子数)を、サンプルデータとしてデータ解析装置906に入力する。免疫染色などのサンプル作製時は、蛍光計測時の誤差を評価して実験誤差としてデータ解析装置906に入力する。これにより、細胞のクラスタリングを実行することができる。
図15は、本実施形態4に係るデータ解析装置906の機能ブロック図である。機能ブロック901は、細胞サンプルとして組織切片を採取し、脱パラフィン処理後に免疫染色を実施する。このとき、免疫染色と異なる蛍光波長の既知量の蛍光色素を細胞内に導入すると同時に核染色と細胞膜染色を実施する。7色の多色の蛍光計測を実施することによって、ターゲット蛋白(3種類)と実験誤差計測用色素、核、細胞膜に対応する像を得る。ターゲット蛋白の種類を増やす場合には計測色を増やす必要がある。細胞膜と核を認識し、細胞膜中に核が1つある物体を細胞と認識する。個々の細胞と認識された領域のターゲット蛋白に相当する蛍光色の強度、または単分子蛍光時の分子数を定量値として記録する。同時に実験誤差用色素の定量も実行する。実験誤差は、認識した細胞の面積あたりの強度の標準偏差を求め、細胞の平均面積から算出する。
本発明は上記した実施形態に限定されるものではなく、様々な変形例が含まれる。上記実施形態は本発明を分かりやすく説明するために詳細に説明したものであり、必ずしも説明した全ての構成を備えるものに限定されるものではない。また、ある実施形態の構成の一部を他の実施形態の構成に置き換えることもできる。また、ある実施形態の構成に他の実施形態の構成を加えることもできる。また、各実施形態の構成の一部について、他の構成を追加・削除・置換することもできる。
上記各構成、機能、処理部、処理手段等は、それらの一部や全部を、例えば集積回路で設計する等によりハードウェアで実現してもよい。また、上記の各構成、機能等は、プロセッサがそれぞれの機能を実現するプログラムを解釈し、実行することによりソフトウェアで実現してもよい。各機能を実現するプログラム、テーブル、ファイル等の情報は、メモリ、ハードディスク、SSD(Solid State Drive)等の記録装置、ICカード、SDカード、DVD等の記録媒体に格納することができる。
906:データ解析装置、904:演算部、905:データ入出力部、908:サンプルデータ入力部、909:実験誤差データ入力部。
Claims (10)
- 複数のサンプルデータをクラスタリング解析する装置であって、
複数のサンプルデータを受け取るサンプルデータ入力部と、
前記サンプルデータの実験誤差についての情報を記述した実験誤差データを受け取る実験誤差データ入力部と、
前記複数のサンプルデータをクラスタリング空間内においてクラスタリングする演算部と、
前記クラスタリングの結果を出力する出力部と、
を備え、
前記演算部は、
前記実験誤差データが記述している前記実験誤差の範囲に応じて、前記クラスタリングを実施する際におけるクラスタ境界を緩和するためのクラスタ範囲パラメータをあらかじめ取得しておき、
仮設定したクラスタ総数に準じて前記複数のサンプルデータをクラスタリングし、
前記複数のサンプルデータのうちいずれのクラスタにも属さない例外データについては、
前記クラスタリング空間上で前記例外データからさらに前記クラスタ範囲パラメータによって定められる距離だけ離れた領域がいずれかのクラスタに含まれる場合は、前記例外データがそのクラスタに属するものと判定し、前記領域がいずれのクラスタにも含まれない場合は、前記例外データが独立したクラスタを構成するものと判定し、
前記演算部は、
各前記サンプルデータが前記クラスタリングの結果得られた各クラスタに属する確率の尤もらしさを表す第1対数尤度と、各前記サンプルデータが前記クラスタリングの結果得られた各クラスタに属さない確率の尤もらしさを表す第2対数尤度とをそれぞれ算出する処理を、前記第1対数尤度と前記第2対数尤度を用いて算出される前記クラスタリング結果の尤もらしさが所定閾値に達するまで繰り返すことにより、最適なクラスタ総数を求め、
得られた最適なクラスタ総数に準じて、前記複数のサンプルデータの最終的なクラスタリング結果を決定する
ことを特徴とするデータ解析装置。 - 前記演算部は、
前記仮設定したクラスタ総数に準じて前記複数のサンプルデータをクラスタリングする過程において、前記サンプルデータが仮設定したクラスタに属すると仮定して前記第1対数尤度と前記第2対数尤度を算出し、
前記仮設定したクラスタに属すると仮定した前記サンプルデータについては、前記クラスタリング空間上において、前記仮設定したクラスタの中心からの距離が離れるにしたがって、前記仮設定したクラスタに属する確率を低く評価し、
前記仮設定したクラスタに属さないと仮定した前記サンプルデータについては、前記クラスタリング空間上において、前記仮設定したクラスタの中心からの距離が離れるにしたがって、前記仮設定したクラスタに属さない確率を高く評価する
ことを特徴とする請求項1記載のデータ解析装置。 - 前記演算部は、
クラスタ内に属する前記サンプルデータの個数が所定個数以上であるか否かを基準として、前記サンプルデータが前記例外データであるか否かを判定する
ことを特徴とする請求項1記載のデータ解析装置。 - 前記演算部は、前記所定個数をランダムに選択する
ことを特徴とする請求項3記載のデータ解析装置。 - 前記演算部は、所定の確率分布にしたがって前記所定個数をランダムに選択する
ことを特徴とする請求項3記載のデータ解析装置。 - 前記演算部は、
前記クラスタ範囲パラメータを掃引し、各前記クラスタ範囲パラメータの値を採用した場合における前記クラスタリングによって得られたクラスタ総数を取得し、
前記第1対数尤度と前記第2対数尤度に基づき算出した前記クラスタリング結果の尤もらしさの値が極値となるクラスタ総数を、最適なクラスタ総数として採用する
ことを特徴とする請求項1記載のデータ解析装置。 - 前記演算部は、
前記クラスタリングを実施する過程において得られる情報を用いて、前記クラスタリングの結果の信頼度指標を算出し、
前記出力部は、
前記信頼度指標を前記クラスタリングの結果と併せて出力する
ことを特徴とする請求項1記載のデータ解析装置。 - 前記演算部は、
前記第1対数尤度と前記第2対数尤度に基づき算出した前記クラスタリング結果の尤もらしさの値を、前記クラスタリングの結果の信頼度指標として算出する
ことを特徴とする請求項7記載のデータ解析装置。 - 前記サンプルデータ入力部と前記実験誤差データ入力部は、細胞の解析結果に関するデータをそれぞれ前記サンプルデータおよび前記実験誤差データとして受け取り、
前記演算部は、前記クラスタリングによって前記細胞をグループ化する
ことを特徴とする請求項1記載のデータ解析装置。 - 複数のサンプルデータをクラスタリング解析する方法をコンピュータに実行させるデータ解析プログラムであって、前記コンピュータに、
複数のサンプルデータを受け取るサンプルデータ入力ステップ、
前記サンプルデータの実験誤差についての情報を記述した実験誤差データを受け取る実験誤差データ入力ステップ、
前記複数のサンプルデータをクラスタリング空間内においてクラスタリングする演算ステップ、
前記クラスタリングの結果を出力する出力ステップ、
を実行させ、
前記演算ステップにおいては、前記コンピュータに、
前記実験誤差データが記述している前記実験誤差の範囲に応じて、前記クラスタリングを実施する際におけるクラスタ境界を緩和するためのクラスタ範囲パラメータをあらかじめ取得するステップ、
仮設定したクラスタ総数に準じて前記複数のサンプルデータをクラスタリングするステップ、
前記複数のサンプルデータのうちいずれのクラスタにも属さない例外データについては、
前記クラスタリング空間上で前記例外データからさらに前記クラスタ範囲パラメータによって定められる距離だけ離れた領域がいずれかのクラスタに含まれる場合は、前記例外データがそのクラスタに属するものと判定し、前記領域がいずれのクラスタにも含まれない場合は、前記例外データが独立したクラスタを構成するものと判定するステップ、
を実行させ、
前記演算ステップにおいては、前記コンピュータにさらに、
各前記サンプルデータが前記クラスタリングの結果得られた各クラスタに属する確率の尤もらしさを表す第1対数尤度と、各前記サンプルデータが前記クラスタリングの結果得られた各クラスタに属さない確率の尤もらしさを表す第2対数尤度とをそれぞれ算出する処理を、前記第1対数尤度と前記第2対数尤度を用いて算出される前記クラスタリング結果の尤もらしさが所定閾値に達するまで繰り返すことにより、最適なクラスタ総数を求めるステップ、
得られた最適なクラスタ総数に準じて、前記複数のサンプルデータの最終的なクラスタリング結果を決定するステップ、
を実行させる
ことを特徴とするデータ解析プログラム。
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
PCT/JP2012/080003 WO2014080447A1 (ja) | 2012-11-20 | 2012-11-20 | データ解析装置、データ解析方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
JP6029683B2 true JP6029683B2 (ja) | 2016-11-24 |
JPWO2014080447A1 JPWO2014080447A1 (ja) | 2017-01-05 |
Family
ID=50775654
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
JP2014548349A Expired - Fee Related JP6029683B2 (ja) | 2012-11-20 | 2012-11-20 | データ解析装置、データ解析プログラム |
Country Status (3)
Country | Link |
---|---|
US (1) | US20150302042A1 (ja) |
JP (1) | JP6029683B2 (ja) |
WO (1) | WO2014080447A1 (ja) |
Families Citing this family (14)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN104903957A (zh) * | 2013-01-10 | 2015-09-09 | 富士通株式会社 | 控制方法、控制程序以及控制装置 |
JP6336932B2 (ja) * | 2015-03-03 | 2018-06-06 | 富士フイルム株式会社 | 細胞群検出装置および方法並びにプログラム |
JP7041614B6 (ja) * | 2015-08-17 | 2022-05-31 | コーニンクレッカ フィリップス エヌ ヴェ | 生体データにおけるパターン認識のマルチレベルアーキテクチャ |
EP3361416A1 (en) | 2017-02-09 | 2018-08-15 | Tata Consultancy Services Limited | Statistical modeling techniques based neural network models for generating intelligence reports |
US10846311B2 (en) * | 2017-11-01 | 2020-11-24 | Mad Street Den, Inc. | Method and system for efficient clustering of combined numeric and qualitative data records |
JP6655595B2 (ja) | 2017-12-21 | 2020-02-26 | 三菱日立パワーシステムズ株式会社 | 単位空間生成装置、プラント診断システム、単位空間生成方法、プラント診断方法、及びプログラム |
US20200012890A1 (en) | 2018-07-06 | 2020-01-09 | Capital One Services, Llc | Systems and methods for data stream simulation |
CN109325059A (zh) * | 2018-12-03 | 2019-02-12 | 枘熠集成电路(上海)有限公司 | 一种数据比较方法及装置 |
CN110287456A (zh) * | 2019-06-30 | 2019-09-27 | 张家港宏昌钢板有限公司 | 基于数据挖掘的大盘卷轧制表面缺陷分析方法 |
CN110490229A (zh) * | 2019-07-16 | 2019-11-22 | 昆明理工大学 | 一种基于spark和聚类算法的电能表检定误差诊断方法 |
WO2022009342A1 (ja) * | 2020-07-08 | 2022-01-13 | 富士通株式会社 | 情報処理プログラム、情報処理方法および情報処理装置 |
CN113418885A (zh) * | 2021-07-22 | 2021-09-21 | 合肥学院 | 一种用于紫外分光光度计实验数据的分析方法 |
CN113673582B (zh) * | 2021-07-30 | 2023-05-09 | 西南交通大学 | 基于系统聚类分析的铁路动态基准点多层级分群方法 |
CN116796214B (zh) * | 2023-06-07 | 2024-01-30 | 南京北极光生物科技有限公司 | 一种基于差分特征的数据聚类方法 |
Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JP2004272350A (ja) * | 2003-03-05 | 2004-09-30 | Nec Corp | クラスタリング装置、クラスタリング方法、クラスタリングプログラム |
JP2006163894A (ja) * | 2004-12-08 | 2006-06-22 | Hitachi Software Eng Co Ltd | クラスタリングシステム |
JP2011520183A (ja) * | 2008-04-25 | 2011-07-14 | コーニンクレッカ フィリップス エレクトロニクス エヌ ヴィ | サンプルデータの分類 |
-
2012
- 2012-11-20 US US14/443,118 patent/US20150302042A1/en not_active Abandoned
- 2012-11-20 WO PCT/JP2012/080003 patent/WO2014080447A1/ja active Application Filing
- 2012-11-20 JP JP2014548349A patent/JP6029683B2/ja not_active Expired - Fee Related
Patent Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JP2004272350A (ja) * | 2003-03-05 | 2004-09-30 | Nec Corp | クラスタリング装置、クラスタリング方法、クラスタリングプログラム |
JP2006163894A (ja) * | 2004-12-08 | 2006-06-22 | Hitachi Software Eng Co Ltd | クラスタリングシステム |
JP2011520183A (ja) * | 2008-04-25 | 2011-07-14 | コーニンクレッカ フィリップス エレクトロニクス エヌ ヴィ | サンプルデータの分類 |
Non-Patent Citations (1)
Title |
---|
JPN6016016161; Rajesh N. Dave: 'Characterization and detection of noise in clustering' Pattern Recognition Letters Vol.12, NO.11, 199111, p.657-664 * |
Also Published As
Publication number | Publication date |
---|---|
WO2014080447A1 (ja) | 2014-05-30 |
US20150302042A1 (en) | 2015-10-22 |
JPWO2014080447A1 (ja) | 2017-01-05 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
JP6029683B2 (ja) | データ解析装置、データ解析プログラム | |
JP6313757B2 (ja) | 統合デュアルアンサンブルおよび一般化シミュレーテッドアニーリング技法を用いてバイオマーカシグネチャを生成するためのシステムおよび方法 | |
JP2019531700A5 (ja) | ||
JP2022533003A (ja) | バイオ医薬品生産における細胞株選択のためのデータ駆動予測モデリング | |
US20050159896A1 (en) | Apparatus and method for analyzing data | |
US20200227168A1 (en) | Machine learning in functional cancer assays | |
CN109411015A (zh) | 基于循环肿瘤dna的肿瘤突变负荷检测装置及存储介质 | |
CN112634987B (zh) | 一种单样本肿瘤dna拷贝数变异检测的方法和装置 | |
EP2883179A2 (en) | Population classification of genetic data set using tree based spatial data structure | |
JP5854346B2 (ja) | トランスクリプトーム解析方法、疾病判定方法、コンピュータプログラム、記憶媒体、及び解析装置 | |
Raza et al. | A novel anticlustering filtering algorithm for the prediction of genes as a drug target | |
CN107463797B (zh) | 高通量测序的生物信息分析方法及装置、设备及存储介质 | |
KR101067352B1 (ko) | 생물학적 네트워크 분석을 이용한 마이크로어레이 실험 자료의 작용기작, 실험/처리 조건 특이적 네트워크 생성 및 실험/처리 조건 관계성 해석을 위한 알고리즘을 포함한 시스템 및 방법과 상기 방법을 수행하기 위한 프로그램을 갖는 기록매체 | |
JP5963198B2 (ja) | 動的ネットワークバイオマーカーの検出装置、検出方法及び検出プログラム | |
Padmanaban et al. | Between-tumor and within-tumor heterogeneity in invasive potential | |
CN101517579A (zh) | 蛋白质查找方法和设备 | |
WO2012149107A2 (en) | Stratifying patient populations through characterization of disease-driving signaling | |
Wang et al. | Learning dynamics by computational integration of single cell genomic and lineage information | |
KR20240046481A (ko) | 지문 분석을 이용하여 화합물을 생리학적 조건과 연관시키는 시스템 및 방법 | |
JP6517933B2 (ja) | 検査システム、検査装置、及び検査方法 | |
WO2022056478A2 (en) | Automated classification of immunophenotypes represented in flow cytometry data | |
JP2016201123A (ja) | 動的ネットワークバイオマーカーの検出装置、検出方法及び検出プログラム | |
CN110751983A (zh) | 一种筛选特征mRNA用于诊断早期肺癌的方法 | |
CN110475874A (zh) | 脱靶序列在dna分析中的应用 | |
CN116825367A (zh) | 一种组织硬度预测模型的建立、组织硬度预测方法及装置 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
TRDD | Decision of grant or rejection written | ||
A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 Effective date: 20161011 |
|
A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20161018 |
|
R151 | Written notification of patent or utility model registration |
Ref document number: 6029683 Country of ref document: JP Free format text: JAPANESE INTERMEDIATE CODE: R151 |
|
LAPS | Cancellation because of no payment of annual fees |