(57)【要約】【課題】検出限界がある測定方法で得られた測定データから、正確に母集団の大きさを推定する母集団大きさ推定方法を提供する。【解決手段】  コンピュータに、検出限界がある測定方法で得られた測定データの総データ数の推定を実行させる。まず、大きい順／小さい順等に、測定データをソートする。次に、測定データから、前記総データ数の候補値の範囲を設定する。さらに、範囲内の全ての候補値について、標準化した前記測定データと、想定される測定データの分布との確率分布の散布図を描画する。そして、これらの散布図の直線性をピアソンのｒ値等の相関関係又は線形関係に係る関数により評価し、最適な候補値を母集団の大きさとして推定する。【選択図】図２
2014035653.tif 000002
【特許請求の範囲】
【請求項１】
コンピュータに、検出限界がある測定方法で得られた測定データの総データ数の推定を実行させる母集団大きさ推定方法において、
前記測定データをソートし、
前記測定データにより算出された範囲で、前記総データ数の大きさの候補値を複数用意し、
前記候補値毎に、前記測定データの母集団の分布に係る理論値のセットを用意し、
ソートされた前記測定データと前記理論値のセットとの一致を相関関係又は線形関係に係る関数により評価し、前記関数による評価が最適な前記候補値を母集団の大きさとして推定する
ことを特徴とする母集団大きさ推定方法。
【請求項２】
前記関数は、上位からの同順の前記測定データの値と前記理論値のセットの値との直線性又は距離を評価する
ことを特徴とする請求項１に記載の母集団大きさ推定方法。
【請求項３】
前記関数として、ピアソンの積率相関係数を用い、
最適な前記候補値が前記範囲の端であった場合、又は用意した候補値の精度が低い場合には、前記候補値の前記範囲を調整し、
十分な精度が得られたら評価が最適な前記候補値を確定する
ことを特徴とする請求項１に記載の母集団大きさ推定方法。
【請求項４】
請求項１乃至３のいずれか１項に記載の、評価が最適な前記候補値の前記理論値のセットと、ソートされた前記測定データとを用いて、
上位からの同順の前記測定データの値と前記理論値のセットの値とを対比させ、一次近似した切片と傾きを算出し、前記母集団の分布の特性を示す母数を算出する
ことを特徴とする母数推定方法。
【請求項５】
測定データとして、遺伝子発現量の定量化データを用い、
請求項１乃至３のいずれか１項に記載の母集団大きさ推定方法により推定された前記総データ数により前記遺伝子発現量を標準化する
ことを特徴とする遺伝子発現量解析方法。
【請求項６】
前記遺伝子発現量の定量化データは、トランスクリプトの発現量、ｃＤＮＡの読み取り回数、及び派生データである
ことを特徴とする請求項５に記載の遺伝子発現量解析方法。
【請求項７】
測定データとして、タンパク質量又は物質量の定量化データを用い、
請求項１乃至３のいずれか１項に記載の母集団大きさ推定方法により推定された前記総データ数により前記タンパク質量又は物質量を標準化する
ことを特徴とするプロテオーム又はメタボロームのデータ分析方法。
【請求項８】
測定データとして、ＲＮＡｓｅｑ及びＤＮＡマイクロアレイの遺伝子発現量の定量化データを用い、
請求項１乃至３のいずれか１項に記載の母集団大きさ推定方法により推定された前記総データ数によりそれぞれ標準化し、
前記ＤＮＡマイクロアレイの各スポットの感度を、標準化された前記ＲＮＡｓｅｑの遺伝子発現量で補正して求める
ことを特徴とするＤＮＡマイクロアレイの感度測定方法。
【請求項９】
コンピュータに、検出限界がある測定方法で得られた測定データの総データ数を推定する母集団大きさ推定方法を実行させるプログラムにおいて、
前記コンピュータに、
前記測定データをソートさせ、
前記測定データにより算出された範囲で、前記総データ数の大きさの候補値を複数用意させ、
前記候補値毎に、前記測定データの母集団の分布に係る理論値のセットを用意させ、
ソートされた前記測定データと前記理論値のセットとの一致を相関関係又は線形関係に係る関数により評価し、前記関数による評価が最適な前記候補値を母集団の大きさとして推定させる母集団大きさ推定方法を実行させる
ことを特徴とするプログラム。
【請求項１０】
コンピュータに、検出限界がある測定方法で得られた測定データの総データ数を推定する母集団大きさ推定方法を実行させるプログラムを記録したコンピュータ読み取り可能な記録媒体において、
前記コンピュータに、
前記測定データをソートさせ、
前記測定データにより算出された範囲で、前記総データ数の大きさの候補値を複数用意させ、
前記候補値毎に、前記測定データの母集団の分布に係る理論値のセットを用意させ、
ソートされた前記測定データと前記理論値のセットとの一致を相関関係又は線形関係に係る関数により評価し、前記関数による評価が最適な前記候補値を母集団の大きさとして推定させる母集団大きさ推定方法を実行させるためのプログラムを記録した
ことを特徴とする記録媒体。
【請求項１１】
検出限界がある測定方法で得られた測定データの総データ数の推定を実行させる母集団大きさ推定装置において、
前記測定データをソートする、測定データソート手段と、
前記測定データにより算出された範囲で、前記総データ数の大きさの候補値を複数用意する候補値用意手段と、
前記候補値毎に、前記測定データの母集団の分布に係る理論値のセットを用意する分布用意手段と、
ソートされた前記測定データと前記理論値のセットとの一致を相関関係又は線形関係に係る関数により評価し、前記関数による評価が最適な前記候補値を母集団の大きさとして推定する候補値評価手段とを備える
ことを特徴とする母集団大きさ推定装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、母集団大きさ推定方法、母数推定方法、遺伝子発現量解析方法、データ分析方法、感度測定方法、プログラム、記録媒体、及び母集団大きさ推定装置に係り、特に、検出感度に限界があるデータを用いた際の母集団大きさ推定方法、母数推定方法、遺伝子発現量解析方法、データ分析方法、感度測定方法、プログラム、記録媒体、及び母集団大きさ推定装置に関する。
【背景技術】
【０００２】
遺伝する生物情報はゲノムにＤＮＡの塩基配列の形で記録されている。ＤＮＡは細胞内でＲＮＡへと転写され、これらＲＮＡは細胞内に蓄積する。
それらＲＮＡは特定のタンパク質の情報を保持しているが、頻繁に翻訳されるＲＮＡは細胞内の濃度（発現量）が高く保たれる。つまり、これらＲＮＡの濃度は、遺伝情報がはじめて物理的な形に変わる段階である。細胞に含まれるＲＮＡの種類と発現量に係るコンテンツ（内容）が、トランスクリプトームと称されている。各ＲＮＡの発現量は、その細胞の機能や状態を知る上で重要な手掛かりとなる。
【０００３】
ＲＮＡの種類は数万に登るため、それらの種類を見分けた上で定量するのは容易ではない。しかし、１９９０年代にＤＮＡチップを使ったＤＮＡマイクロアレイ法が登場し、その網羅的な解析が可能になった。これはハイブリダイゼーションの画像等から、遺伝子ごとに感度が異なる相対値を測定値として得る方法である。
また、ＤＮＡチップを用いないトランスクリプトームの測定方法として、定量ＰＣＲを使う方法と、ＤＮＡシーケンサーを使う方法がある。
さらに、細胞内の代謝産物である、有機化合物やタンパク質や代謝産物等を、クロマトグラフィ等の方法で分離・定量しながら質量分析し、物質の同定を行う方法も実用化されている。これらの物質の同定方法の例として、ＬＣ−ＭＳ（Ｌｉｑｕｉｄ  Ｃｈｒｏｍａｔｏｇｒａｐｈ  Ｍａｓｓ  Ｓｐｅｃｔｏｍｅｔｅｒ）等が存在する。
【０００４】
ここで、複数のサンプルからなる相対値の測定データを解析する際に、最初に行うことは、当該測定データの標準化である。測定に用いられるサンプルは、先に何らかの方法でその量を揃えることが多いが、揃えられた状態にて取得される数値は、絶対値にはならない。これは、測定の原理上、サンプルごとに感度が異なるからである。この感度差のために、測定間においてデータを直接に比較することはできない。そのため、測定データを別途、標準化することが必要である。この標準化により、各サンプルが比較可能になる。
測定データの標準化を行うためには、データの母集団の分布をあらわす「母数」を推定する。ここでは、あるサンプルに内包されている被測定物質のすべてを、母集団であると考える。たとえば、ある特定の強度域から失われた「強め」に測定されたサンプルのデータからでも、その母集団の性質を表す母数が分かれば、そのサンプルがどの程度に強めに測定されたのかを推定することができる。そのため、母数を使うことで、測定間の比較が可能になる。
たとえば、正規分布モデルであるなら、平均で推定される分布の中心の母数μと、標準偏差で推定される分布の巾の母数σがある。これらの母数の推定は、その後のあらゆる解析の基盤となるために重要である。
【０００５】
図９の例を参照して、測定データの標準化と母集団の関係について、より詳しく説明する。
たとえば、図９（ａ）のように、本来の母集団が平均μ＝０，標準偏差σ＝１の正規分布に従うデータである場合、そのデータのヒストグラムは、図９（ａ）のようなものになる。
さて、従来から、感度限界等により、検出できたものだけがデータとして得られるような検出限界のある測定方法は多い。たとえば電気泳動、クロマトグラフィ、質量分析のように、データを経時的に記録するものは、必ず検出限界がある。また、次世代型シーケンサーを用いたＲＮＡ−ｓｅｑのように、データが基本的に「読まれる・読まれない」というベルヌーイ試行を用いた測定の場合も、読まれなかったデータは得られないために、検出限界が存在する。さらに、蛍光や電気的な検知方法を用いたハイブリダイゼーション等も、機器の感度限界に係る検出限界がある。
図９（ｂ）は検出限界があるデータの測定値のヒストグラムの例である。このように、検出限界より小さい分布の一部が特異的に失われることが多い。逆に、検出限界を超えて大きいシグナルが無視されたり、測定機の不全を招くことで、当該データが失われることもある。
【０００６】
従来から、このように検出限界がある測定方法を用いる場合、測定データの母集団の大きさが不明となり、このため分布特性や母数が求められないことが多かった。
ここで、「母集団の大きさ」は、被測定物の種類の総数を示す値である。従来は、データが母集団からのランダムサンプルだと仮定し、たとえば正規分布を仮定して母数を推定して、母手段の大きさを標準化してきた。
しかしながら、この図９（ｂ）のデータから、例えば、平均を算出すると０．３で、標準偏差は０．８程度となる。これは本来の母集団の平均０、標準偏差１から大きく外れている。こうした違いは、データが、ある特定の強度域から失われ、ランダムなサンプリングであるという前提が崩れるために引き起こされる。
【０００７】
また、このように分布の一部を特異的に失う場合、いわゆるロバストな計算方法を用いても、正しく母数を算出できない。これは、撹乱の原因が少数の外れ値ではないためである。たとえば、図９（ｂ）のデータで、よく使われるロバストなμの推定方法である中央値を計算すると０．２２となり、同じくロバストなσの推定方法であるＭＡＤ（ｍｅｄｉａｎ  ａｂｓｏｌｕｔｅ  ｄｅｖｉａｔｉｏｎ）は０．８４となる。
加えて、実際には、検出限界はサンプルや測定によって異なる。これは、たとえ測定機器が同じであっても、測定されている個々のサンプル毎に、そのシグナル総量が一定しないためである。
つまり、図９（ｃ）、（ｄ）、（ｅ）の例で灰色（塗りつぶし）で示すような、感度限界以下で検出できない消失領域の割合は、一定にはならない。このため、図９（ｃ）では、平均が０．１５で標準偏差が０．８８となり、図９（ｄ）では、平均が０．３０で標準偏差が０．８０となり、図９（ｅ）では平均が０．５１で標準偏差が０．７０となる。
【０００８】
この消失領域の割合のばらつきは、とくに母集団が対数正規分布である場合には顕著になる。たとえば細胞の熱力学的な特性のために、ｍＲＮＡは対数正規分布する。この性質から、多くの量的なオミックスデータは同様な分布をするものと予想される。
対数正規分布では、もっとも下位のデータがもっとも頻度よく現れる。そのため、僅かな検出感度の違いが、得られるデータの数を大きく変えることになる。当然、分布の母数の推定値は失われるデータ数によって大きく影響される。
【０００９】
ここで、従来のＤＮＡマイクロアレイデータの標準化方法として、特許文献１を参照すると、一つのＤＮＡマイクロアレイから得られたデータセットにおいて、各データを大きさ順に並べ替えてデータの順位を決定するステップと、あるデータの順位を、標準正規累積分布関数の逆関数に入力して、当該データの標準化後の予測値を出力するステップと、当該データを前記予測値で置き換えるステップと、を備えたＤＮＡマイクロアレイデータの処理方法が開示されている（以下、従来技術１とする。）。
従来技術１の処理方法によれば、従来のパラメトリック手法の標準化における直線性からの歪みと、パラメータ推定に必要な繰り返し演算の必要性とを回避する標準化方法を提供することができる。
【先行技術文献】
【特許文献】
【００１０】
【特許文献１】国際公開第２００８／０５６６９３号
【発明の概要】
【発明が解決しようとする課題】
【００１１】
従来の検出限界がある測定方法では、いずれも感度に限界があり、その限界が測定ごとに異なる。またトランスクリプトーム等に係る発現量データは、基本的に対数正規分布をする。このため、僅かな測定限界の違いに影響される測定データのデータ数（母集団の大きさ）の違いが大きくなる傾向がある。
しかしながら、従来技術１は、検出限界のためにデータを失う場合における、測定データの母集団の大きさの正確な推定は想定していなかった。
【００１２】
本発明は、このような状況に鑑みてなされたものであり、上述の課題を解消することを課題とする。
【課題を解決するための手段】
【００１３】
本発明の母集団大きさ推定方法は、コンピュータに、検出限界がある測定方法で得られた測定データの総データ数の推定を実行させる母集団大きさ推定方法において、前記測定データをソートし、前記測定データにより算出された範囲で、前記総データ数の大きさの候補値を複数用意し、前記候補値毎に、前記測定データの母集団の分布に係る理論値のセットを用意し、ソートされた前記測定データと前記理論値のセットとの一致を相関関係又は線形関係に係る関数により評価し、前記関数による評価が最適な前記候補値を母集団の大きさとして推定することを特徴とする。
本発明の母集団大きさ推定方法は、前記関数は、上位からの同順の前記測定データの値と前記理論値のセットの値との直線性又は距離を評価することを特徴とする。
本発明の母集団大きさ推定方法は、前記関数として、ピアソンの積率相関係数を用い、最適な前記候補値が前記範囲の端であった場合、又は用意した候補値の精度が低い場合には、前記候補値の前記範囲を調整し、十分な精度が得られたら評価が最適な前記候補値を確定することを特徴とする。
本発明の母数推定方法は、前記評価が最適な前記候補値の前記理論値のセットと、ソートされた前記測定データとを用いて、上位からの同順の前記測定データの値と前記理論値のセットの値とを対比させ、一次近似した切片と傾きを算出し、前記母集団の分布の特性を示す母数を算出することを特徴とする。
本発明の遺伝子発現量解析方法は、測定データとして、遺伝子発現量の定量化データを用い、前記母集団大きさ推定方法により推定された前記総データ数により前記遺伝子発現量を標準化することを特徴とする。
本発明の遺伝子発現量解析方法は、前記遺伝子発現量の定量化データは、トランスクリプトの発現量、ｃＤＮＡの読み取り回数、及び派生データであることを特徴とする。
本発明のプロテオーム又はメタボロームのデータ分析方法は、測定データとして、タンパク質量又は物質量の定量化データを用い、前記母集団大きさ推定方法により推定された前記総データ数により前記タンパク質量又は物質量を標準化することを特徴とする。
本発明のＤＮＡマイクロアレイの感度測定方法は、測定データとして、ＲＮＡｓｅｑ及びＤＮＡマイクロアレイの遺伝子発現量の定量化データを用い、前記母集団大きさ推定方法により推定された前記総データ数によりそれぞれ標準化し、前記ＤＮＡマイクロアレイの各スポットの感度を、標準化された前記ＲＮＡｓｅｑの遺伝子発現量で補正して求めることを特徴とする。
本発明のプログラムは、コンピュータに、検出限界がある測定方法で得られた測定データの総データ数を推定する母集団大きさ推定方法を実行させるプログラムにおいて、前記コンピュータに、前記測定データをソートさせ、前記測定データにより算出された範囲で、前記総データ数の大きさの候補値を複数用意させ、前記候補値毎に、前記測定データの母集団の分布に係る理論値のセットを用意させ、ソートされた前記測定データと前記理論値のセットとの一致を相関関係又は線形関係に係る関数により評価し、前記関数による評価が最適な前記候補値を母集団の大きさとして推定させる母集団大きさ推定方法を実行させることを特徴とする。
本発明の記録媒体は、コンピュータに、検出限界がある測定方法で得られた測定データの総データ数を推定する母集団大きさ推定方法を実行させるプログラムを記録したコンピュータ読み取り可能な記録媒体において、前記コンピュータに、前記測定データをソートさせ、前記測定データにより算出された範囲で、前記総データ数の大きさの候補値を複数用意させ、前記候補値毎に、前記測定データの母集団の分布に係る理論値のセットを用意させ、ソートされた前記測定データと前記理論値のセットとの一致を相関関係又は線形関係に係る関数により評価し、前記関数による評価が最適な前記候補値を母集団の大きさとして推定させる母集団大きさ推定方法を実行させるためのプログラムを記録したことを特徴とする。
本発明の母集団大きさ推定装置は、検出限界がある測定方法で得られた測定データの総データ数の推定を実行させる母集団大きさ推定装置において、前記測定データをソートする、測定データソート手段と、前記測定データにより算出された範囲で、前記総データ数の大きさの候補値を複数用意する候補値用意手段と、前記候補値毎に、前記測定データの母集団の分布に係る理論値のセットを用意する分布用意手段と、ソートされた前記測定データと前記理論値のセットとの一致を相関関係又は線形関係に係る関数により評価し、前記関数による評価が最適な前記候補値を母集団の大きさとして推定する候補値評価手段とを備えることを特徴とする。
【発明の効果】
【００１４】
本発明によれば、候補値の範囲を設定して想定される分布と散布図で比較することで、検出限界のある測定データの母集団の大きさを正確に推定して標準化に用いることができる母集団大きさ推定方法を提供することができる。
【図面の簡単な説明】
【００１５】
【図１】本発明の実施の形態に係る解析装置１０の制御構成を示すブロック図である。
【図２】本発明の実施の形態に係る母集団大きさ推定処理のフローチャートである。
【図３】本発明の実施の形態に係る測定データ２５１と散布図の描画の例を示すグラフである。
【図４】本発明の実施の形態に係る実施例１にて、母集団大きさ推定方法をＲＮＡｓｅｑに適用した結果を示すグラフである。
【図５】本発明の実施の形態に係る実施例２にて、母集団大きさ推定方法によるＲＮＡｓｅｑデータの標準化と、従来の標準化との比較を行った結果を示すグラフである。
【図６】本発明の実施の形態に係る実施例３にて、母集団大きさ推定方法を用いた際の、ＲＮＡｓｅｑデータとＤＮＡマイクロアレイデータとの一致を調べた結果を示すグラフである。
【図７】本発明の実施の形態に係る実施例４にて、母集団大きさ推定方法を用いて、ＲＮＡｓｅｑデータを用いてＤＮＡマイクロアレイの感度測定を行った結果を示すグラフである。
【図８】本発明の実施の形態に係る実施例５にて、実際にＲ言語のコードを使用して、作成した計算機実験用のヒストグラムである。
【図９】従来の検出限界のある測定データの標準化の手法を説明する概念図である。
【発明を実施するための形態】
【００１６】
＜実施の形態＞
〔解析装置１０の制御構成〕
まず、図１を参照して、本発明の実施の形態に係る解析装置１０（コンピュータ、母集団大きさ推定装置、解析手段）の制御構成について説明する。
解析装置１０は、例えばＰＣ／ＡＴ互換機や汎用機等である計算装置であって、Ｌｉｎｕｘ（登録商標）、Ｗｉｎｄｏｗｓ（登録商標）等のＯＳがインストールされている。
解析装置１０の主な構成要素としては、ＣＰＵ（Ｃｅｎｔｒａｌ  Ｐｒｏｃｅｓｓｉｎｇ  Ｕｎｉｔ）やＧＰＵ（Ｇｒａｐｈｉｃｓ  Ｐｒｏｃｅｓｓｉｎｇ  Ｕｎｉｔ）等の制御・演算装置である制御部１００と、ＲＡＭ（Ｒａｎｄｏｍ  Ａｃｃｅｓｓ  Ｍｅｍｏｒｙ）やＲＯＭ（Ｒｅａｄ  Ｏｎｌｙ  Ｍｅｍｏｒｙ）やＨＤＤ（Ｈａｒｄ  Ｄｉｓｋ  Ｄｒｉｖｅ）やフラッシュメモリやＳＳＤ（Ｓｏｌｉｄ  Ｓｔａｔｅ  Ｄｒａｉｖｅ）等の記憶装置／記録媒体である記憶部１１０と、キーボードやマウス等のポインティングデバイスやタッチパネル等やマイクロアレイ解析装置等の外部機器からのＩ／Ｏインターフェイス等を含む入力部１３０と、液晶ディスプレイや有機ＥＬディスプレイや印刷を行うプリンタ等である表示部１４０と、１０００Ｂａｓｅ−Ｔ等の規格のＬＡＮボードや無線ＬＡＮボードやシリアル／パラレル／ＵＳＢインターフェイス等である入出力部１５０とを備えている。
解析装置１０は、主に記憶部１１０に記憶された各種プログラムと、データベース等を含むデータとを用いて制御部１００が実行することで、本発明の実施の形態に係る母集団大きさ推定方法をハードウェア資源を用いて実現することができる。
【００１７】
記憶部１１０は、本発明の実施の形態に係る母集団大きさ推定方法を実現するためのプログラム及びデータが記憶されている。この記憶部１１０のプログラム及びデータは、制御部１００により、ハードウェア資源を用いて実行／処理することができる。
このプログラム及びデータは、データベース２００、測定データソート部２１０（測定データソート手段）、候補値範囲設定部２２０（候補値用意手段）、分布作成部２２５（分布用意手段）、候補値描画部２３０（候補値描画手段）、及び候補値評価部２４０（候補値評価手段、母数算出手段）とを含んで構成される。
また、このプログラム及びデータは、別途記録媒体に記録され、解析装置１０にインストール可能に構成されていてもよい。また、インターネット等からダウンロードしてインストールすることも可能である。
【００１８】
データベース２００は、ＳＱＬ等のデータベースや各種データを記憶する部位である。データベース２００には、主に、測定データ２５１、分布データ２５２（理論値のセット）、設定データ２５３、及び元データ２５４を記憶している。
【００１９】
測定データ２５１は、元データ２５４の値を、測定データソート部２１０により、大きい順にソート等して作成されたデータである。
分布データ２５２は、分布作成部２２５で作成された、元データ２５４や測定データ２５１と同じと想定される分布の理論値のセットに係るデータである。なお、分布データ２５２には、想定される分布の乱数値を、候補値分記憶しておいてもよい。
設定データ２５３は、元データ２５４から測定データ２５１を作成するための設定、測定データ２５１の想定される分布の種類、候補値の範囲を設定するためのパラメータ、候補値の範囲のパラメータ、候補値の範囲、候補値のセット等を記憶するデータである。
元データ２５４は、検出限界がある測定方法を用いる実験機器から取得した、元のデータである。この元データ２５４は、例えば、トランスクリプトームに係る遺伝子発現量の定量化データや、タンパク質や低分子有機化合物等の質量計測のデータ等を用いることができる。遺伝子発現量の定量化データとしては、例えば、解析する遺伝子のＲＮＡや翻訳されたタンパク質の量等に対応する、ＲＮＡｓｅｑのリード数、ＤＮＡマイクロアレイのシグナル強度、ＲＴ−ＰＣＲの発現量の値、リアルタイムＲＴ−ＰＣＲの蛍光値（蛍光強度比）、ＳＡＧＥ法（Ｓｅｒｉａｌ  Ａｎａｌｙｓｉｓ  ｏｆ  Ｇｅｎｅ  Ｅｘｐｒｅｓｓｉｏｎ）のタグの計測数、及びこれらの派生データ等を用いることができる。また、検出限界がある質量測定のデータとしては、ＬＣ−ＭＳのデータのような各種質量計測のデータも用いることができる。さらに、計測値が大きさの定まる数値のセットとして得られる、各種の生物学的、非生物学的な値のデータを元データ２５４として用いることが可能である。
【００２０】
測定データソート部２１０は、実験機器から取得した元データ２５４を、大きさによりソートする等の初期化を行って測定データ２５１として記憶部１１０に記憶する。
また、測定データソート部２１０は、測定データ２５１から、候補値の範囲を決定するため、最大値、最小値、平均、メディアン、分散等の統計的な値を算出する。
【００２１】
候補値範囲設定部２２０は、測定データ２５１の本来の完全なデータ数として想定される理論値（以下、「候補値」という。）の範囲を設定する部位である。
また、候補値範囲設定部２２０は、候補値評価部２４０の評価に従って新たな候補値の範囲を設定する。この候補値の範囲の設定については後述する。
【００２２】
分布作成部２２５は、設定データ２５３の想定される分布の種類の設定に従って、設定された候補値に対応する分布の値（理論値）のセットを、分布データ２５２として作成する。ここでは、分布作成部２２５は、所定のｚスコアにおける分布の理論値のセットを算出し、分布データ２５２に記憶することができる。
たとえば、分布作成部２２５は、測定データ２５１が正規分布と想定される場合、標準正規分布を用意して理論値のセットを算出し、分布データ２５２に記憶する。つまり、分布作成部２２５は、例えば、正規分布ならμ＝０，σ＝１の分布を用意する。
また、例えば、測定データ２５１が対数正規分布と想定される場合、分布作成部２２５は、作成した正規分布を対数正規分布に変換して、対数正規分布の値のセットを理論値のセットとして分布データ２５２に記憶する。
なお、分布作成部２２５は、乱数等を用いて、候補値に対応した正規分布の値のセットを直接、作成することもできる。この場合、分布作成部２２５は、作成した分布の値のセットをソートして、分布データ２５２として記憶することもできる。
また、分布作成部２２５は、想定される分布の種類に係る分布関数や確率密度関数自体を分布データ２５２に記憶し、各候補値における所定のｚスコアにおける分布の理論値をその都度算出してもよい。
【００２３】
候補値描画部２３０は、分布データ２５２と測定データ２５１との散布図を作成する部位である。
候補値描画部２３０は、例えば、用意した分布データ２５２の分布の値のセットと、測定データ２５１を各候補値で標準化した値のセットとを、それぞれ大きい順に対応させ、散布図を作成、描画する。
【００２４】
候補値描画部２３０は、本実施形態においては、データが正規分布である場合に、いわゆるｚスコアを算出し、推定された母数を用いた標準化を行う。すなわち、下記の式（１）により、ｚスコアを算出し標準化する：
ｚ＝（測定値−ｍ）／ｓ  ……  式（１）
ただし、上述の式（１）において、ｍは母平均μの推定値、ｓは母標準偏差の推定値である。
【００２５】
候補値評価部２４０は、分布データ２５２と測定データ２５１との関係性を評価する部位である。
候補値評価部２４０は、ソートされた測定データ２５１と、設定された候補値に係る分布データ２５２との一致を、「直線性」または「順位の理論値−ソートデータ間の距離」等で評価する。
候補値評価部２４０は、例えば、直線性の評価関数として、ピアソンの積率相関係数ｒ値を算出して評価することができる。
また、候補値評価部２４０は、ソートデータ間の距離として、分布データ２５２と測定データ２５１の各順位の差の絶対値の総和や平均、二乗和の総和や平均等のマハラノビス汎距離を用いることができる。
候補値評価部２４０は、もっとも評価の高い候補値を、母集団の大きさの推定値として確定する。
また、候補値評価部２４０は、母集団の大きさの推定値を使って、母集団の分布の特性を示す母数を算出することもできる。
【００２６】
〔母集団大きさ推定処理〕
次に、図２〜図３を参照して、本発明の実施の形態に係る母集団大きさ推定方法を実行する母集団大きさ推定処理について説明する。
ここで、本実施形態の母集団大きさ推定処理では、被測定物の種類の数を「母集団の大きさ」として推定する。この「母集団の大きさ」として、例えば、電気泳動の一種であるＳＤＳ−ＰＡＧＥの実験データの場合、疎水性が大きいタンパクや、電荷の強いタンパクはそもそも電気泳動されず、検出できないため、母集団に入れないことが好適である。すなわち、この例では、ＳＤＳ−ＰＡＧＥで検出することが原理的に可能なタンパクの全てが母集団となる。
つまり、被測定物として多くの種類があり、そのうちいくらかを所定の実験手法方法で測定可能である場合、その測定可能な種類の数が母集団の大きさとなる。
【００２７】
本実施形態の母集団大きさ推定処理では、全体の総データ数ｎを未知数とおいて、この総データ数ｎに様々な候補値を仮に代入しながら、最もふさわしい候補値を探していく。この候補値のふさわしさの確認のためには、理論値とデータの一致を相関関係又は線形関係に係る関数により評価する。この際、相関関係又は線形関係としては、「直線性」又は「順位の理論値−ソートデータ間の距離」等で評価する。このソートデータ間の距離としては、マハラノビス距離や「ノルム」等を用いることができる。
ここで、理論値と評価された候補値とをＱＱプロット（Ｑｕａｎｔｉｌｅ−Ｑｕａｎｔｉｌｅ  ｐｌｏｔ）の変形である散布図等で描画した場合に、直線に近づくほど評価が高く、総データ数ｎに対して最もふさわしい値となる。すなわち、望ましい「母集団の大きさ」の推定結果となる。
【００２８】
つまり、本実施形態では、想定される分布の種類が分かっており、検出限界がある測定方法で得られた測定データ２５１の母集団の大きさを推定する。これにより、測定データ２５１の標準化を行う。このために、測定データ２５１の本来の完全なデータ数として想定される総データ数ｎを未知とし、その値を解析的に解く処理を行う。
より具体的に説明すると、まず解析装置１０は、測定データソート部２１０を用いて、元データ２５４をソートし、測定データ２５１を作成する。
次に、解析装置１０は、候補値範囲設定部２２０を用いて、総データ数ｎの候補値の範囲を決定する。この上で、解析装置１０は、候補値範囲設定部２２０を用いて、範囲内から複数の候補値を用意する。
そして、解析装置１０は、候補値描画部２３０により、想定される分布である分布データ２５２と、範囲内の各候補値で標準化した測定データ２５１とを、それぞれ大きい順に対応させ、ＱＱプロットの変形である散布図を描画する。
次に、解析装置１０は、候補値評価部２４０の評価関数により、この総データ数ｎに係る理論値と測定データ２５１データとの一致を評価する。なお、散布図を描画せずに、直接データ間で一致を評価することもできる。
解析装置１０は、候補値範囲設定部２２０により、範囲内の端が最適な候補値と算出されたり、用意した候補値の精度が低かったりした場合、候補値の範囲を調整しつつ、計算を繰り返す。そして、解析装置１０は、十分な精度が得られたら候補値を確定する。この確定した候補値が、測定データ２５１の本来の完全な総データ数ｎの期待値となる。
以下で、図２のフローチャートを参照して、各母集団大きさ推定処理について、ステップ毎に詳しく説明する。
【００２９】
（ステップＳ１０１）
制御部１００は、測定データソート部２１０を用いて、データ読み込みソート処理を行う。
この処理においては、まず、制御部１００は、検出限界がある測定方法で得られた元データ２５４を取得する。この元データ２５４は、他の記録媒体から、記憶部１１０に記憶されていてもよい。また、元データ２５４は、入出力部１５０を介して、ネットワーク上の他のサーバや実験機器等（図示せず）から取得してもよい。
制御部１００は、取得した元データ２５４の発現量等、標準化したい値のセットを、大きさにより、小さい順又は大きい順になるよう並べ替え（ソート）する。
そして、制御部１００は、ソートされた値のセットを、測定データ２５１として記憶部１１０に記憶する。
また、制御部１００は、測定データ２５１の最大値、最小値、平均、分散（標準偏差）等の統計的な値を求め、記憶部１１０に記憶する。
また、図３（ａ）を参照すると、制御部１００は、測定データ２５１を標準正規分布とみて標準化した場合の、ｚスコアを横軸、値の頻度を縦軸としたヒストグラムを描画し、又は度数分布表を作成して記憶部１１０に記憶することもできる。図３（ａ）は、図９（ｂ）と同様の正規分布のデータを描画した例を示している。
【００３０】
データの並べ替えに関しては、例えば、統計用プログラミング言語である「Ｒ言語」を用いた場合、下記のようなコードを用いることができる：
sorted_data <- sort(data, na.last=T, decreasing=T) 
ここでは、配列ｄａｔａの並び替えを行い、配列ｓｏｒｔｅｄ＿ｄａｔａに代入している。
【００３１】
（ステップＳ１０２）
次に、制御部１００は、候補値範囲設定部２２０を用いて、候補値範囲初期設定処理を行う。
この処理では、制御部１００は、測定データ２５１を用いて、測定データ２５１の総データ数ｎの推定値となる「候補値」の範囲の初期値を決定する。
制御部１００は、この候補値の範囲として、まず、設定データ２５３から、想定される分布の種類を取得する。この上で、測定データ２５１の統計的な値を用いて、想定される分布からあり得る範囲を初期値として設定し、設定データ２５３に記憶する。
この候補値の範囲としては、例えば、測定データ２５１を標準化した場合には、ｚスコアで−４から−１までの分布になるような値の積分値を算出することができる。また、データの種類毎に所定値を用意したり、以前に計測した同種のデータから、候補値の範囲の初期値を用意してもよい。
【００３２】
（ステップＳ１０３）
制御部１００は、候補値範囲設定部２２０を用いて、候補値用意処理を行う。
ここでは、制御部１００は、設定された候補値の範囲から、候補値を複数個用意して記憶部１１０の設定データ２５３に記憶する。
制御部１００は、例えば、設定された候補値の範囲を所定値で等分して、この所定値＋１の値分の候補値を用意できる。つまり、上述の例の場合、この所定値が１００の場合、制御部１００は、０〜６５０の範囲を１００等分して、１０１個の候補値を用意することができる。
なお、測定データ２５１がリード数のように整数の場合は、総データ数ｎを整数と仮定し、整数の候補値を用意することもできる。
【００３３】
この候補値用意処理について、Ｒ言語のコードとして、例えば、９０００〜１１０００の範囲を、２０刻みで１０１個の候補値を用意する場合、下記のようなコードを用いることができる：
n_candidates <- seq(from=9000, to=11000, by=20)
ここで、n_candidatesは、候補値の配列である。
また、後述する直線性評価処理において、散布図の直線性を、例えばピアソンのｒ値を評価関数の値として算出し、判定する。このため、各候補値に係る評価関数の値を記憶するためのベクターｒｓも、下記のＲ言語のコードのように用意する：
rs <- 1:101 + NA
【００３４】
（ステップＳ１０４）
制御部１００は、分布作成部２２５を用いて、分布用意処理を行う。
この処理において、まず、制御部１００は、用意された候補値から、ひとつの候補値を選択する。
次に、制御部１００は、設定データ２５３を参照し、測定データ２５１として想定される分布の理論値のセットを用意する。
制御部１００は、正規分布なら、例えば、Ｎ（０，１）、すなわち平均０で分散σ2＝１の標準正規分布を用意し、分布データ２５２として、記憶部１１０に記憶する。
また、制御部１００は、異なる分布モデルを使う場合には、それぞれの標準状態のものを使う。制御部１００は、例えば、対数正規分布なら、γ＝０、μ＝０、σ＝１の対数正規分布を用意し、分布データ２５２として、記憶部１１０に記憶する。
【００３５】
（ステップＳ１０５）
制御部１００は、候補値描画部２３０を用いて、散布図描画処理を行う。
具体的には、まず、制御部１００は、測定データ２５１と、用意した分布である分布データ２５２とを、それぞれ大きい順に対応させて散布図を描く。
ここで、従来のＱＱプロットでは、それぞれ同じ分位点にあるデータを対応させた図である。これに対して、制御部１００は、同じ順位にある数を対応させた図を、本実施形態の散布図として描画する。
このため、本実施形態の散布図は、最も低い領域にあっては、理論値又はソートしたデータが、対応させる相手を持たないという状況が起こる。この場合、制御部１００は、対応がつくペアだけを用いて描画する。
【００３６】
図３（ｂ）〜（ｄ）を参照して詳細に説明すると、制御部１００は、ある確率ｐを与えたときに、２つの確率点（ｑｕａｎｔｉｌｅ）となるｑ１とｑ２とを、例えば、それぞれＸ軸（横軸）、Ｙ軸（縦軸）にとってプロットした散布図（確率プロット）を作成する。ここでは、分布データ２５２をｑ１として描画する。また、制御部１００は、測定データ２５１を、選択された候補値により標準化し、ｑ２として描画する。これは、ＱＱプロット（「Ｇｎａｎａｄｅｓｉｋａｎ，  Ｒ．；  Ｗｉｌｋ，  Ｍ．Ｂ．（１９６８），  "Ｐｒｏｂａｂｉｌｉｔｙ  ｐｌｏｔｔｉｎｇ  ｍｅｔｈｏｄｓ  ｆｏｒ  ｔｈｅ  ａｎａｌｙｓｉｓ  ｏｆ  ｄａｔａ"，  Ｂｉｏｍｅｔｒｉｋａ  ５５  （１）：  １〜１７」を参照）の変形となる。つまり、この散布図では、Ｘ軸の座標に用意した分布の理論値である分布データ２５２を、Ｙ軸の座標にソートしたデータである測定データ２５１をとって比較している。この際、分布データ２５２を大きい順に、また測定データ２５１を大きい順に対応させて、その関係を調べている。
【００３７】
たとえば、この散布図を表示するＲ言語のコードとしては、下記のようなコードを用いることができる：
ideal <- -qnorm(ppoints(n_candidates[i]))
s_limit <- min(length(data[!is.na(data)]), n_candidates[i])
plot(sorted_data[1:s_limit], ideal[1:s_limit])
【００３８】
図３を再び参照すると、図３（ｂ）は、図３（ａ）のデータにおいて、候補値が７０００の際の散布図の例を示している。
図３（ｃ）は、図３（ｂ）と同様に、候補値が１００００の際の散布図の例を示している。
図３（ｄ）は、図３（ｂ）（ｃ）と同様に、候補値が１４０００の際の散布図の例を示している。
【００３９】
なお、制御部１００は、分布作成部２２５を用いて、上述の分布用意処理において、選択された候補値に対応する分布の値のセットを（疑似）乱数列から直接作成し、これを分布データ２５２として記憶することもできる。この際、制御部１００は、正規分布を、ボックス＝ミューラー法（Ｂｏｘ−Ｍｕｌｌｅｒ  ｔｒａｎｓｆｏｒｍ）や、中心極限定理を用いた方法等を用いて作成可能である。また、制御部１００は、対数正規分布も、正規分布から作成することができる。
また、制御部１００は、候補値描画部２３０を用いて、このような分布をもつ乱数列である分布データ２５２と、測定データ２５１とを対応させた散布図を、直接描画するような処理も可能である。
また、制御部１００は、この散布図を描画せずに、直接下記の候補値評価処理を行うこともできる。
【００４０】
（ステップＳ１０６）
制御部１００は、候補値評価部２４０を用いて、候補値評価処理を行う。
この処理において、制御部１００は、分布データ２５２と測定データ２５１との一致について、相関関係又は線形関係に係る関数で評価する。この相関関係又は線形関係に係る関数としては、「直線性」又は「順位の理論値−ソートデータ間の距離」等に係る所定の評価関数を用いることができる。
直線性により評価する場合、この評価関数として、制御部１００は、例えば、ピアソンのｒ値を算出して用いることができる。
また、制御部１００は、このピアソンのｒ値の算出の際、正規分布の場合には、μとσを、それぞれ直線のｙ切片と傾きとして推定することもできる。
【００４１】
図３を再び参照すると、図３（ｂ）の散布図では、ピアソンのｒ値は、０．９８４４３６３であった。また、標準化した測定データ２５１のμとσの推定値は（０．４６，０．７１）であった。
図３（ｃ）では、ピアソンのｒ値は０．９９９８９９６であり、標準化に用いた測定データ２５１のμとσの推定値は（−０．０１，０．９９）であった。
図３（ｄ）では、ピアソンのｒ値は０．９９８１４１６であり、標準化に用いた測定データ２５１のμとσの推定値は（−０．４３，１．１４）であった。
このように、図３（ａ）のデータは、図３（ｃ）のように候補値が１００００の際に、ピアソンのｒ値が大きく、直線性が高くなり、測定データ２５１の総データ数ｎの推定値として良好となることが分かる。
【００４２】
なお、制御部１００は、分布データ２５２と測定データ２５１とを「順位の理論値−ソートデータ間の距離」により評価する場合、この「距離」として、各順位の値の差の絶対値の総和や平均、二乗和の総和や平均等のマハラノビス汎距離等を算出して用いることができる。
【００４３】
（ステップＳ１０７）
制御部１００は、候補値評価部２４０を用いて、範囲内で用意されたすべての候補値を評価したか否かを判定する。
Ｙｅｓ、すなわち全ての候補値を評価し終わった場合、制御部１００は、処理をステップＳ１０８に進める。
Ｎｏ、すなわち範囲内で候補値がある場合、制御部１００は、処理をステップＳ１０４に戻して、分布の用意、散布図の描画、直線性の評価を続ける。この際、例えば、候補値の範囲内で、大きい値、小さい値について交互に評価を行い、徐々にもっとも確からしい候補値を求めるように範囲を狭め、計算を高速に行っても良い。また、ニュートン法を用いてもよい。
【００４４】
このような直線性の評価のＲ言語のコードとして、例えば、下記のようなコードを用いることができる：
for (i in 1:101){
n <- n_candidates[i]
ideal <- -qnorm(ppoints(n))  # 理想値
s_limit <- min(length(data[!is.na(data)]), n)
rs[i]<- cor(sorted_data[1:s_limit], ideal[1:s_limit]) } # Peasronのr値
このコードでは、上述のコードで用意したベクターｒｓについて、単純に全ての候補値を評価して代入する例を示している。
この上で、下記のようなコードを用いて、候補値の中から最適なものを選抜する。ここでは、制御部１００は、最も大きい評価値を与えた候補を選択する。
selected <- n_candidates[which(rs==max(rs))] 
変数ｓｅｌｅｃｔｅｄは、ベクターｒｓ中の評価が最尤な候補値となる。
【００４５】
図３（ｅ）を参照すると、図３（ａ）とは別の測定データ２５１について、０〜６５０を候補値の範囲として、散布図を描いた際のピアソンのｒ値の変化の例を示している。この図３（ｅ）のＸ軸は各候補値、Ｙ軸はピアソンのｒ値である。
このように、制御部１００は、候補値の値を変化させながら、関係の直線性を、上述したようにピアソンのｒ値のような評価関数でモニターして、もっとも確からしい候補値を探すことができる。この関係は、単純に上に凸か下に凸かであると考えられる。
このため、もっとも確からしい候補値を求める作業に、ニュートン法を適用することができ、計算資源を節約することが可能である。
【００４６】
なお、もしこの図が上に凸でなければ、候補値が不適であると考えられるため、制御部１００は、後述する候補値範囲ずらし処理を行う。または、制御部１００は、候補をもっと狭い範囲から選択する候補値範囲狭め処理を行う。
【００４７】
また、ピアソンのｒ値ではない、別の評価関数を用いて、上述の散布図を評価することもできる。
図３（ｆ）（ｇ）を参照すると、例えば、上述のＲ言語のコードのうち、ベクターｒｓに代入する箇所を下記のように変更することもできる：
rs[i]<- -mean((sorted_data[1:s_limit]-ideal[1:s_limit])^2) # 別の評価関数
図３（ｆ）のＸ軸は各候補値、Ｙ軸はこの別の評価関数の値である。このように、評価関数が異なると値が変わるため、データの特性に合わせて評価関数を選択することができる。
図３（ｇ）は、図３（ｃ）と同様に、この別の評価関数の値が最大値となった候補値を代入した変数ｓｅｌｅｃｔｅｄについて、下記コードで散布図を描画した例である：
ideal <- -qnorm(ppoints(selected))
s_limit <- min(length(data[!is.na(data)]), selected)
plot(sorted_data[1:s_limit], ideal[1:s_limit])
【００４８】
（ステップＳ１０８）
制御部１００は、候補値評価部２４０を用いて、用意した候補値の外側（端）により最適な候補値があるか否かを判定する。
ここでは、制御部１００は、上述した評価関数のグラフを、極値がほぼなく、単調増加／減少の場合、もっと最適な候補値がある、すなわちＹｅｓと判定する。つまり、制御部１００は、用意した候補値の端で、評価関数の値が最高値の場合、Ｙｅｓと判定する。それ以外の場合は、Ｎｏと判定する。
Ｙｅｓの場合、制御部１００は、処理をステップＳ１０９に進め、候補値の範囲を変更する。
Ｎｏの場合、制御部１００は、処理をステップＳ１１０に進める。
【００４９】
（ステップＳ１０９）
制御部１００は、候補値の範囲を変更する必要がある場合、候補値範囲設定部２２０を用いて、候補値範囲ずらし処理を行う。
ここでは、用意した候補値のいずれかの端が、評価関数の値が最高値であった場合、候補値の範囲をそちら側に調整してずらす。つまり、例えば、用意した候補値の範囲が、５００〜１０００であり、候補値１０００の評価関数の値が最大値であった場合、候補値の範囲を１００１〜１５０１に再設定するような処理を行う。
制御部１００は、その後処理をステップＳ１０３に戻して、ずらされた候補値の範囲で再度計算を進める。
なお、候補値の再設定の際、候補値の範囲を同一の幅でずらすのではなく、評価関数のグラフの増加率／減少率を基に、範囲の幅を設定することもできる。
【００５０】
（ステップＳ１１０）
制御部１００は、候補値の範囲内にもっとも確からしい候補値があった場合、候補値評価部２４０を用いて、より精度が必要か否かを判定する。
制御部１００は、この精度については、例えば、用意された候補値の範囲の幅と、設定データ２５３の候補値の範囲のパラメータを参照し、より狭い幅で検索する必要がある場合にはＹｅｓと判定する。それ以外の場合、たとえば候補値が整数で用意され、評価した連続する整数内に最適な候補値があった場合には、Ｎｏと判定する。
Ｙｅｓの場合、制御部１００は、処理をステップＳ１１１に進める。
Ｎｏの場合、制御部１００は、処理をステップＳ１１２に進める。
【００５１】
（ステップＳ１１１）
制御部１００は、候補値範囲設定部２２０を用いて、候補値範囲狭め処理を行う。
制御部１００は、候補値の範囲を狭めて再度設定する。ここでは、制御部１００は、上述の直線性の評価関数の値が最高値であった候補値の前後の値を、新たな候補値の範囲と設定することができる。
制御部１００は、その後処理をステップＳ１０３に戻して、狭められた候補値の範囲で再度計算を進める。
【００５２】
（ステップＳ１１２）
制御部１００は、充分な精度で最適な候補値が得られた場合、候補値評価部２４０を用いて、候補値確定処理を行う。
ここでは、制御部１００は、十分な精度が得られたら候補値を確定する。
これが、測定データ２５１の本来の完全な総データ数ｎの期待値、すなわち母集団の大きさと推定される。
【００５３】
（ステップＳ１１３）
次に、制御部１００は、候補値評価部２４０を用いて、確定した候補値を基に母数の推定を行い、この母数により測定データ２５１を標準化して記憶する。
この際に、制御部１００は、失われたデータ数を勘案しながら、例えば、直線近似を用いて算出することが可能である。
この際、制御部１００は、ソートデータと理論値を上位からの同順の値でもって対比させ、一次近似して切片と傾きを算出する。制御部１００は、算出された切片と傾きを用いて、母数を推定する。制御部１００は、この一次近似において、最小二乗法や、ロバストな別の方法等を用いることができる。
【００５４】
このような母数の推定のためのパラメータの設定に関するＲ言語のコードとして、例えば、下記のようなコードを用いることができる：
u_limit <- selected-length(data[!is.na(data)])
params <- coef ( lm(sorted_data[u_limit: s_limit] ~ ideal[u_limit: s_limit]))
この上で、μとσの推定値は、それぞれ：
(mu <- params[1]) # μ
(sigma <- params[2]) # σ
のようなコードで算出できる。すなわち、正規分布の場合、この例のようにμが切片、σが傾きとなる。
制御部１００は、これらの算出された母数を用いて、測定データ２５１を標準化して記憶する。
【００５５】
（ステップＳ１１４）
制御部１００は、候補値評価部２４０を用いて、分析処理を行う。
次に、必要な場合、制御部１００は、実際の測定データと、他に標準化された測定データ２５１との比較を行う。たとえば、制御部１００は、同じ発現量データであるＲＮＡｓｅｑとＤＮＡマイクロアレイのデータを比較することができる。これにより、ＤＮＡマイクロアレイの各スポットの感度を、標準化されたＲＮＡｓｅｑの遺伝子発現量で補正して求めることが可能となる。
また、制御部１００は、主成分分析等のデータ解析を行うこともできる。
以上により母集団大きさ推定処理を終了する。
【００５６】
以上のように構成することで、以下のような効果を得ることができる。
上述したように、検出限界のためにデータを失う測定方法のデータでは、母数の正確な推定は困難であった。
これに対して、本発明の実施の形態に係る母集団大きさ推定方法は、検出限界がある測定方法の測定データ２５１の分布モデルと、想定される分布との関係から母集団の大きさである総データ数ｎを推定し、μやσのような分布の母数も求めることもできる。このため、標準化した値の精度と確度が高まる。特に、対数正規分布をする測定データ２５１においては、総データ数ｎは、母数を正確に求める際に重要である。本発明はこれにより、検出限界がある測定方法のデータについて、標準化を正確に行うことができる。
すなわち、本実施形態の母集団大きさ推定方法により母集団の大きさを推定して標準化することで、ＲＮＡｓｅｑ、メタボローム解析、ｑＰＣＲを用いるマイクロアレイ等のトランスクリプトーム解析、次世代シーケンサーのデータ分析等を、従来より正確に行うことができる。
また、測定では検出限界が変化することがあるが、本実施形態に係る母集団大きさ推定方法は、こうした場合でも、標準化した値を比較することができる。よって、例えば、有意差が得られる測定項目を増やすことができる。
また、本発明の実施の形態に係る母集団大きさ推定方法で、疾病等に係るトランスクリプトームのデータを標準化して、主成分分析等に用いることで、より正確な分析を行うことができる。
【００５７】
また、従来、ＤＮＡマイクロアレイのそれぞれのスポットの感度は、測定方法がなかった。そのため、異なる設計のＤＮＡマイクロアレイのデータは、たとえ同じ遺伝子を測っていたとしても、直接に比較することが困難であった。これは、遺伝子のどの配列を測定するかによって、プローブに用いられる塩基配列が違うため、そのプローブの感度が異なるからであった。
本実施形態の母集団大きさ推定方法は、こういったＤＮＡマイクロアレイのデータに関しても、別の測定方法と対応させて標準化することで、各スポットの感度を推定することができる。このように、各スポットの感度を求めることで、他の臓器や患者等のＤＮＡマイクロアレイのデータを補正することが可能となる。これによって、マイクロアレイのデータを一種の絶対値として用いることができる。このため、異なる設計のチップのデータを、直接に比較できるようになる。
【００５８】
またＤＮＡマイクロアレイは、チップの設計にないＲＮＡは検出できない。このため、特に近年その重要性がクローズアップされてきた、タンパクをコードしないＲＮＡのような、新規／未知の対象を測定することはできないという問題があった。
これに対して、本実施形態の母集団大きさ推定方法は、ＲＮＡｓｅｑのような新規／未知の対象を測定できるトランスクリプトーム等の解析にも適用できる。また、よく知られているＤＮＡマイクロアレイのデータと比較することで、当該新規／未知の対象の発現量等を確実に測定できる。また、ＤＮＡマイクロアレイの測定結果の感度差に関わらず、ＲＮＡｓｅｑと結果を直接比較することができる。
【００５９】
なお、本実施形態の母集団大きさ推定方法は、生物学的データ以外にも適用することができる。
たとえば、食品、土壌、水、空気中に含まれる物質を、クロマトグラフィ等の物質の分離方法と質量分析等の検出方法と組み合わせて同定・定量する場合等に応用することが可能である。その場合、これらの物質の分布形式に関する情報が必要である。たとえば、食品なら、生物由来なので対数正規分布を用いることができる。
【実施例１】
【００６０】
次に図面に基づき本発明を実施例によりさらに説明するが、以下の具体例は本発明を限定するものではない。
【００６１】
〔ＲＮＡｓｅｑのデータの標準化への適用〕
ＲＮＡｓｅｑは、ＲＮＡを逆転写してＤＮＡにして、そのＤＮＡ断片の配列を決定するトランスクリプトームの測定方法である。ＲＮＡｓｅｑでは、より多量に存在したＲＮＡの相補配列がより高い頻度で観測されるため、この頻度から各ＲＮＡの濃度を推定する。ＲＮＡｓｅｑは、スプライシング等が異なる未知のＲＮＡでも測定可能である。またデータが観測される／されないというデジタルな測定方法であり、画像のような仲介物もないため、より正確な測定結果が期待されている。
しかしながら、ＲＮＡｓｅｑは、測定項目が大きい多変量データであり、解析は容易ではない。また観測数を数えるベルヌーイ試行に由来する、二項分布のノイズを持つ。当初は標準化が不要ではないかとさえ言われたが、総リード数の違いをあわせるためにも、データの比較には当然、標準化が必要である。
このため、本実施例において、ＲＮＡｓｅｑデータについて、本実施形態の母集団大きさ推定方法により母集団の大きさを推定して、標準化を行った。
【００６２】
（使用データ）
文献「Ｊｏｈｎ  Ｃ．  Ｍａｒｉｏｎｉ  ｅｔ．  ａｌ．，  ＲＮＡ−ｓｅｑ：  Ａｎ  ａｓｓｅｓｓｍｅｎｔ  ｏｆ  ｔｅｃｈｎｉｃａｌ  ｒｅｐｒｏｄｕｃｉｂｉｌｉｔｙ  ａｎｄ  ｃｏｍｐａｒｉｓｏｎ  ｗｉｔｈ  ｇｅｎｅ  ｅｘｐｒｅｓｓｉｏｎ  ａｒｒａｙｓ，  Ｇｅｎｏｍｅ  Ｒｅｓ．  ２００８，  ８（９）：  １５０９−１５１７．」に報告されたデータを元データ２５４として使用した。
これは単一のヒト検体からの腎臓および肝臓のサンプルをイルミナ社製の次世代シーケンサーによるＲＮＡｓｅｑ、及びアフィメトリクス社製のＤＮＡマイクロアレイで測定したデータである。
【００６３】
（結果）
肝臓（Ｒ１Ｌ４Ｌｉｖｅｒ）と腎臓（Ｒ１Ｌ１Ｋｉｄｎｅｙ）のサンプルについて、初期範囲として、１２０００〜２３０００を指定し、元データ２５４の分布として対数正規分布を用い、本実施形態の母集団大きさ推定方法で各候補値の散布図を描画し、ピアソンのｒ値を測定した。対数の底には１０を用いた。
結果として、腎臓のサンプルでは、最適な候補値は、１２６２０であった。つまり、総データ数ｎは、１２６２０と推定される。この腎臓のサンプルの母数は、平均μ＝１．７７、標準偏差σ＝０．５１２であった。
また、肝臓のサンプルでは、最適な候補値は１３０６０であった。つまり、総データ数ｎは、１３０６０と推定される。この肝臓のサンプルの母数は、平均μ＝１．５１５、標準偏差σ＝０．６１０であった。
図４を参照して、この測定データ２５１と分布データ２５２の分布モデルの関係とを、下記で説明する。
【００６４】
図４（ａ）は、肝臓のサンプルについて描画した、上述のＱＱプロットの変形である散布図である。図４（ａ）のＸ軸は用意した正規分布のｚスコア、Ｙ軸は肝臓のサンプルから作成した測定データ２５１を標準化したもののｚスコアである。また、丸印で示される各点は、算出した確率点を示す。
図４（ｂ）は、この肝臓のサンプルを標準化した際のヒストグラムである。Ｘ軸は用意した正規分布のｚスコア、Ｙ軸のバーは各サンプルの頻度を示す。また、点線は、最適な候補値の分布から期待される頻度を示す。
この図４（ａ）（ｂ）の結果のように、ｚスコアにして、下が０．５付近までのデータが、算出した最適な候補値により求められた母数のモデルと一致することがわかる。これは全体の約４０％に相当する。ｚスコアで０．５以下のデータは、測定限界以下であり、サンプル由来のノイズ、及び二項分布に従う測定機器／実験方法由来のノイズに影響されると考えられる。
【００６５】
図４（ｃ）は、図４（ａ）（ｂ）と同じ肝臓のサンプルを散布図ではなく、大きさ順に直接プロットした図である。この図でも、Ｘ軸は用意した正規分布のｚスコア、Ｙ軸は測定データ２５１を標準化したもののｚスコアである。この図では、分布データ２５２に合わせて算出した最適な候補値分のデータを作成して大きさ順に並べ、測定データ２５１とペアで描画して比較する。
この結果は上述の図４（ａ）（ｂ）と同様に、０．５以下のデータはむしろノイズによる影響が大きく、データに再現性が期待できない。
このため、このサンプルにおいて、大きさがｚスコアで０．５以下のデータを用いないで発現量を比較することも考えられる。
【００６６】
図４（ｄ）（ｅ）は、上述の図４（ａ）（ｂ）と同様に、腎臓のサンプルについて、本実施の形態の母集団大きさ推定方法で標準化した測定データ２５１と正規分布との散布図による比較、及びヒストグラムを示す。
この結果として、ｚスコアにして−１程度のところで、算出した最適な候補値により求められた母数のモデルからの乖離が観測される。これ以下のデータは、サンプル由来のノイズ、又は二項分布のノイズの影響を受けていることが示唆される。
【実施例２】
【００６７】
〔ＲＮＡｓｅｑデータの従来の標準化法との比較〕
次に、図５を参照して、上述の実施例１のＭａｒｉｏｎｉらの文献と同じＲＮＡｓｅｑのデータを用いて、よく用いられている従来の標準化法（以下、「総和一定」と称する。）と、本実施形態の母集団大きさ推定方法とで、発現量を比較した。
従来の「総和一定」の標準化法は、サンプル（リード）数の総数で全体を除して、測定した値の強度を合わせる標準化法である。
【００６８】
（遺伝子の選択）
７回繰り返し測定をした肝臓及び腎臓のＲＮＡｓｅｑのデータを、それぞれの遺伝子ごとにｔ検定し、肝臓又は腎臓のいずれかで両側１％の有意水準で発現量（リード数）の差がある、すなわちポジティヴ（陽性）となる遺伝子を選んだ。
これら遺伝子に着目し、肝臓及び腎臓で対になる発現量の組み合わせを得た。つまり、統計検定で１％有意となった遺伝子において、肝臓と腎臓の発現量の対数比（ログレシオ）を求めて比較した。
【００６９】
（結果）
図５（ａ）は、従来の、総数で全体を除した総和一定の標準化法により標準化した各遺伝子の発現量を用い、２組のデータから腎臓と肝臓の差を求めて、それぞれを比較したプロットである。Ｘ軸とＹ軸はそれぞれ、肝臓と腎臓の遺伝子発現の差を、ログレシオを用いて示している。各点は、各遺伝子の腎臓及び肝臓の発現量の差を、Ｘ座標に一方のデータ組み合わせ，Ｙ座標に他方の組み合わせとしてプロットしたものである。図中の「ｒ」は、ピアソンの相関係数ｒ値を示す。
図５（ａ）のように、総和一定の標準化だと、発現変化の再現性がさほど高くないことがわかる。特に、原点（Ｘ座標，Ｙ座標）＝（０，０）付近でｙ＝−ｘ方向にばらつく、即ち逆相関を示す遺伝子が散見される。これは、特に発現量の差が少ない遺伝子群において、検定が偽陽性を生じやすいという傾向を示している。この現象は、測定にノイズ成分が多いときに起きる。
【００７０】
図５（ｂ）は、本実施形態の母集団大きさ推定方法で発現量を標準化した測定データ２５１の腎臓と肝臓での発現量の違いを、（ａ）と同じ２組の測定の組み合わせにおいて比較したプロットである。Ｘ軸とＹ軸はそれぞれ、肝臓と腎臓の遺伝子発現の差をログレシオを用いて示している。各点は、各遺伝子の腎臓及び肝臓の発現量を図５（ａ）と同様にプロットしたもので、図中の「ｒ」はピアソンの相関係数ｒ値である。
図５（ｂ）では、原点を境に、各プロットが直線ｙ＝ｘ上で２つに分離されている。つまり、逆相関をする遺伝子はほとんど検出されず、検定が偽陽性を生じにくいことを示している。これは、測定結果のノイズ成分が、シグナルに比べて小さいときにおきる現象である。測定の再現性は高くなっており、ｒ値も従来の０．８９７から０．９４９程度に増えている。
また、またｙ＝ｘで大きな値をとる遺伝子の数が、本実施形態の母集団大きさ推定方法のプロットのほうが明らかに多いことがわかる。これは、ノイズ成分が減少したために検定で有意差が検出される遺伝子が増加したからだと考えられる。すなわち、結果的に相関係数を高めているだけでなく、データのもつ価値を高めることができる。
【実施例３】
【００７１】
〔ＲＮＡｓｅｑデータとＤＮＡマイクロアレイデータとの一致〕
次に、上述の実施例１、２と同じＭａｒｉｏｎｉらのＲＮＡｓｅｑとＤＮＡマイクロアレイのデータについて、同じ遺伝子の発現量を測定したＲＮＡｓｅｑと、ＤＮＡマイクロアレイのデータとが、一致するか比較した。具体的には、本実施形態の母集団大きさ推定方法が影響するか否かかを調べた。
実施例３では、上述の実施例２と同じ、肝臓又は腎臓のいずれかのＲＮＡｓｅｑの発現量のデータにおいて、両側１％の有意水準でｔ検定をかけ、発現量（リード数）の差があった遺伝子の発現量を元データ２５４として用いた。また、これに加えて、ＤＮＡマイクロアレイでも、１％有意水準で発現量の差があったものを選択して元データ２５４として用いた。
【００７２】
図６（ａ）は、元々のＭａｒｉｏｎｉらのオリジナルのデータをプロットしたコントロールである。Ｙ軸がＲＮＡｓｅｑでの発現量の対数比であり、Ｘ軸がアフィメトリクス社製のＤＮＡマイクロアレイの発現量データである。これは、ＤＮＡマイクロアレイの発現データを、従来のＲＭＡ（Ｒｏｂｕｓｔ  Ｍｕｌｔｉ−Ａｒｒａｙ  Ａｖｅｒａｇｅ）法で標準化処理して算出した値である。ＲＮＡｓｅｑで２５０リード以上であった遺伝子を濃い色（濃い塗りつぶし）の点で示す。それより少ないデータをものを薄い色（薄い塗りつぶし）の点で示す。２５０リード以上の遺伝子は、上位６．４％のデータに限られている。ピアソンのｒ値は、０．８７程度であった。
【００７３】
図６（ｂ）は、本実施の形態の母集団大きさ推定方法で標準化を行った後、同じ遺伝子について図６（ａ）と同様にプロットを行ったものである。特に、薄い色（薄い塗りつぶし）の点で示す、２５０リード未満の比較的低いリード数において、上下へのバラツキが少なくなっている。具体的には、ピアソンのｒ値が０．９１４程度に向上している。
すなわち、少ない発現量の遺伝子においても、ＲＮＡｓｅｑ及びＤＮＡマイクロアレイの発現量を再現性高く、誤差を少なくして分析できる。また、ＲＮＡｓｅｑ及びＤＮＡマイクロアレイの発現量を、比較することも可能になる。
【実施例４】
【００７４】
〔ＤＮＡマイクロアレイの感度測定〕
次に、上述の実施例３と同じデータを用いて、ＲＮＡｓｅｑとＤＮＡマイクロアレイとの比較から、ＤＮＡマイクロアレイの各スポットの感度の測定を行った。
従来のＤＮＡマイクロアレイの発現量データの解析では、通常、コントロールの組織の測定量と、調べたい組織の発現量を測定した際の測定量との「倍率」を算出して、遺伝子発現の変化量を求めている。このため、ＤＮＡマイクロアレイの各スポットに形成されたＤＮＡ配列のＴｍ値（ｍｅｌｔｉｎｇ  ｔｅｍｐｅｒａｔｕｒｅ）や収率等により生じる感度差を相殺することができている。しかしながら、より一般性の高い定量値を得るために、各スポットの感度を測定したいという要望があった。
これに対して、本実施例では、同一遺伝子のＲＮＡｓｅｑのデータを用いて、ＤＮＡマイクロアレイの各スポット自体の感度を測定した。
【００７５】
（感度差の推定）
上述の実施例３の肝臓のデータを、ＲＮＡｓｅｑとＤＮＡマイクロアレイを本実施形態の母集団大きさ推定方法で標準化し、各ＤＮＡマイクロアレイの各遺伝子に対応するスポットの感度を推定する。
この各スポットの感度を用いて、各ＤＮＡマイクロアレイの各スポットの感度差を補正する。具体的には、各スポットについて、肝臓のＲＮＡｓｅｑとＤＮＡマイクロアレイでで得られた感度差を推定し、この感度差の値を腎臓のデータに掛けて標準化する。
【００７６】
図７（ａ）（ｂ）は、それぞれＹ軸にマイクロアレイの測定値を標準化した値、Ｘ軸にＲＮＡ−ｓｅｑを本実施形態の母集団大きさ推定方法で標準化した値をとって比較したものである。Ｘ軸、Ｙ軸ともに、感度の倍率ではなく、対数値のｚスコアを示している。また、図７（ａ）は肝臓のデータ、図７（ｂ）は、腎臓のデータをそれぞれ示している。
結果として、正の相関が見られるものの、大きくばらついている。すなわち、実施例３では、倍率を観察して感度が相殺されるため、結果が一致した。しかしながら、上述したように、マイクロアレイにはスポットごとに異なる感度があるため、絶対的な感度で比較すると、必ずしも結果が一致しないことが分かる。
【００７７】
図７（ｃ）は、図７（ａ）の結果を基にして、「ＲＮＡ−ｓｅｑの標準化データ−マイクロアレイの標準化データ」として感度を推定し、この感度を図７（ｂ）のマイクロアレイの結果に和することで補正したものである。すなわち、図７（ｃ）は、図７（ａ）の結果を基にして、肝臓のデータから推定された各スポット（遺伝子）の感度差を使って補正したプロットである。図７（ａ）（ｂ）と同様に、Ｘ軸、Ｙ軸とも絶対値（ログスケール）を示している。
結果として、ＲＮＡ−ｓｅｑとの相関が、明らかに図７（ｂ）よりも向上している。このように、ＤＮＡマイクロアレイの感度測定により、結果の散らばりを抑えて、各遺伝子毎に、発現量の測定を正確に行うことができる。また、感度推定のためのデータ数を増やせば、さらに精度を高めることができる。
【実施例５】
【００７８】
〔Ｒ言語での実装例〕
次に、上述の実施の形態に係る母数の推定を実際にＲ言語で実装して、計算機実験を行った実施例を以下で示す。
まず、以下のコードで、正規分布する乱数を一万個を用意して、配列ｒａｎｄに記憶した。
また、配列ｒａｎｄから、ｚスコアで−０．５より小さいデータを消去して作成した配列ｄａｔａを用意して、測定限界をシミュレートした：
rand <- rnorm(10000) # 
data <- rand ; data[ data< -0.5] <-NA  
【００７９】
この用意した乱数ｒａｎｄと配列ｄａｔａとについてヒストグラムの作成を行った：
hr <- hist(rand, breaks=-50:50/10) #ヒストグラムの作成
hd <- hist(data, breaks=-50:50/10)
【００８０】
この上で、以下のコードで、ヒストグラムを重ね書きして描画した：
plot(hr, col="gray50", xlab="", ylab="", main="")
par(new=T)
plot(hd, col="white", xlab="", ylab="", main="")
【００８１】
図８は、描画したヒストグラムの例である。図８の影がついている部分が、図９と同様に、検出されなかったデータに相当する。
次に、得られた配列ｒａｎｄの母数を、以下のコードで推定した：
# 平均
> mean(rand, na.rm=T)
[1] -0.0009032517
# 中央値（ロバストなμの推定）
> median(rand, na.rm=T)
[1] 0.0079797
# 標準偏差（ロバストなμの推定）
> sd(rand, na.rm=T)
[1] 0.9975465
# 標準偏差（ロバストなＭＡＤの推定）
> mad(rand, na.rm=T)
[1] 0.9996847
【００８２】
これらの配列ｒａｎｄから求めた母数の値は、どれもよく作成した正規分布の設定を反映していた。
次に、従来方法である下記のコードで、得られた配列ｄａｔａの母数を推定した：
# 平均
> mean(data, na.rm=T)
[1] 0.5033331
# 中央値（ロバストなμの推定）
> median(data, na.rm=T)
[1] 0.3849361
# 標準偏差（ロバストなμの推定）
> sd(data, na.rm=T)
[1] 0.6977276
# 標準偏差（ロバストなＭＡＤの推定）
> mad(data, na.rm=T)
[1] 0.7157046
【００８３】
このように、配列ｄａｔａは、従来方法で計算すると、ロバストな計算方法でも、本来の正規分布の設定から離れている。
これに対して、上述の実施形態のＲ言語のコードにより、失われたデータ数を勘案しつつ、直線近似を用いて母数を推定した：
(mu <- params[1]) # μ
-0.0003264361 
(sigma <- params[2]) # σ
0.9854621 
【００８４】
結果として、それぞれ、設定値に近い母数の推定に成功した。
【００８５】
なお、上記実施の形態の構成及び動作は例であって、本発明の趣旨を逸脱しない範囲で適宜変更して実行することができることは言うまでもない。
【産業上の利用可能性】
【００８６】
本発明の母集団大きさ推定方法は、次世代シーケンサーやＤＮＡマイクロアレイ等を用いたトランスクリプトームの解析に幅広く利用でき、産業上に利用することができる。
【符号の説明】
【００８７】
１０  解析装置
１００  制御部
１１０  記憶部
１３０  入力部
１４０  表示部
１５０  入出力部
２００  データベース
２１０  測定データソート部
２２０  候補値範囲設定部
２２５  分布作成部
２３０  候補値描画部
２４０  候補値評価部
２５１  測定データ
２５２  分布データ
２５３  設定データ
２５４  元データ
【図１】
2014035653.tif 000003
【図２】
2014035653.tif 000004
【図３】
2014035653.tif 000005
【図４】
2014035653.tif 000006
【図５】
2014035653.tif 000007
【図６】
2014035653.tif 000008
【図７】
2014035653.tif 000009
【図８】
2014035653.tif 000010
【図９】
2014035653.tif 000011