JP2020000197A - 計数方法、濃度測定装置及び濃度測定システム - Google Patents

Abstract

【課題】検体となる溶液を分割した小区画の大きさにばらつきがあっても、比較的少ない数の小区画を用いて精度良く注目対象物の数及び濃度を推定できる計数方法を提供する。【解決手段】複数の第１の小区画のうち、第１の注目対象を内包する前記第１の小区画を第１の対象区画とし、該第１の対象区画に含まれる当該第１の注目対象の数を取得する計数方法であって、第１の注目対象の数に関する第１の候補値を、アルゴリズムを用いて算出する第１の算出工程と、第２の小区画の体積に関する情報を用いて、第２の注目対象を第２の小区画に分配する計算を行う分配計算工程と、第１の算出工程と共通のアルゴリズムを一部に用いて、第２の注目対象の数に関する第２の候補値を算出する第２の算出工程と、を有し、第２の候補値と第２の注目対象の数とを比較し、所定の収束条件を満たす第２の候補値を第１の注目対象の数として決定する決定工程と、を含む。【選択図】図１２

Description

本発明は、計数方法、濃度測定装置及び濃度測定システムに関する。

病原の特定や癌の治療方針決定等において、標的核酸濃度測定の利用が進みつつある。測定方法として、核酸数を増幅するＰＣＲ（Ｐｏｌｙｍｅｒａｓｅ Ｃｈａｉｎ Ｒｅａｃｔｉｏｎ）法と、標的核酸に結合する際に蛍光を発する蛍光プローブとを組み合わせた方法が、精力的に開発されている。しかしながら、医療応用に対しては、測定精度が不足しているのが現状である。

通常、蛍光色素濃度（つまり標的核酸濃度）は蛍光強度から測定するが、蛍光強度が微弱であり階調を正確に測定することはできない。ＰＣＲ法による核酸数の増幅率の再現性は高くないため、蛍光強度（核酸数）を増幅してからの測定によっても測定精度は向上しない。その中で、濃度が既知の核酸に対し、標的核酸と同一環境でＰＣＲ法を適用してリファレンスを作成する方法（リアルタイムＰＣＲ法）が広く利用されている。蛍光強度の階調を測定する手法の中では高い精度が得られるが、リファレンス作成等の処理が複雑で、リファレンスの増幅率との差異が推定誤差として残る問題がある。

この課題に対する解として、蛍光強度の階調の測定を必要としない、ｄｉｇｉｔａｌ ＰＣＲ（以下、ｄＰＣＲと表記する）法が知られている。ｄＰＣＲ法では、標的核酸の含まれる溶液を多数の小区画に分割し、区画毎にＰＣＲ法を行って、蛍光発光の有無を調べる。全区画数に対する蛍光を発しない区画数の比率から、標的核酸濃度を推定する。

非特許文献１に紹介されているように、ｄＰＣＲ法を基にした核酸濃度測定装置が多数提案されている。その中で、オイル内に生成した多数の小さい液滴を小区画として使用する方式は、ランニングコストの低減や小型化で有利であり、主流となっている。液滴を小区画として使用する方式では、ポアソン分布を用いて標的核酸濃度を推定するが、上記推定方法は、液滴の大きさが同じであることを前提としている。

液滴の大きさがばらついた際には、測定に利用可能な液滴数が減少し、標的核酸濃度の推定精度が低下する。特許文献１には、小区画の体積を測定する工程を有し、液滴の体積が一定でない場合に適用可能なｄＰＣＲ法が開示されている。

米国特許第９，３２２，０５５号明細書

Ｌ．Ｃａｏ ｅｔ ａｌ．，"Ａｄｖａｎｃｅｓ ｉｎ ｄｉｇｉｔａｌ ｐｏｌｙｍｅｒａｓｅ ｃｈａｉｎ ｒｅａｃｔｉｏｎ（ｄＰＣＲ） ａｎｄ ｉｔｓ ｅｍｅｒｇｉｎｇ ｂｉｏｍｅｄｉｃａｌ ａｐｐｌｉｃａｔｉｏｎｓ"，Ｂｉｏｓｅｎｓｏｒｓ ａｎｄ Ｂｉｏｅｌｅｃｔｒｏｎｉｃｓ，ｖｏｌ．９０，ｐａｇｅｓ４５９−４７４（２０１７）．

小区画として液滴を利用するｄＰＣＲ型濃度測定装置においては、液滴生成にかかる時間の短縮が望まれている。従来技術において、液滴数は２万個程度であり、今後高精度化、広ダイナミックレンジ化のためにさらなる増加が予想されている。一般に、液滴生成速度を上げると、液滴の大きさがばらつく傾向がある。ポアソン分布を利用した標的核酸濃度推定は、液滴の大きさが均一であることを前提としているが、特許文献１で開示される方法のように、各液滴体積を測定して、濃度推定値を補正することで対処可能と考えられている。

しかしながら、現在知られている濃度推定値の補正方法には、いくつかの課題がある。例えば、大きさが同程度である液滴をグループ化し、グループ毎にポアソン分布により推定した濃度を求め、その平均値を全体の濃度とする方法がある。ポアソン分布を使った場合、濃度の推定精度は液滴数に依存するが、全液滴数に比べてグループ毎の液滴数は少ないため、グループ毎の推定精度が低下する。精度の低い数値の平均を求めても、推定精度は向上しない。

特許文献１において挙げられているように、事後確率や尤度を最大化する濃度推定方法も考えられる。事後確率等を利用する濃度推定方法では、標的核酸を液滴に散布した場合に、蛍光を発する液滴数が、装置の測定値と一致する確率（標的核酸数を条件とした条件付き確率）を計算する。結果として確率が最大となる標的核酸数、つまり濃度を求める。この手法では、確率を計算する際に、一つの液滴内に複数の標的核酸が入る場合を考慮する必要がある。例えば、１０３個の標的核酸により、１００個の液滴が蛍光を発する確率を考える。実際の測定では、標的核酸数１０３個は少ない部類であるが、各液滴が標的核酸を２個までしか含まないとしても、液滴の組み合わせは、１６１７００（＝１００×９９×９８÷（２×３））通りまで増加する。多数の組み合わせに対し、確率を求めて加算する必要があるため、標的核酸数が多い場合には長時間の計算を必要とする。また、液滴の大きさの測定精度が低い場合には、各組み合わせの確率に含まれる微小誤差が蓄積されて、推定精度が低下するおそれがある。

本発明は、上記に鑑み、粒子を用いた注目対象物を含む、検体となる溶液を分割した小区画の大きさにばらつきがあっても、比較的少ない数の小区画を用いて精度良く注目対象物の数、つまり濃度を推定することができる計数方法を提供することを目的とする。

本実施形態に係る方法は、複数の第１の注目対象を含む検体を分割した複数の第１の小区画のうち、前記第１の注目対象を内包する前記第１の小区画を第１の対象区画とし、該第１の対象区画に含まれる当該第１の注目対象の数に関する情報を取得する計数方法であって、複数の第１の小区画の数と、第１の対象区画の数と、に基づいて、第１の注目対象の数に関する情報の第１の候補値を、アルゴリズムを用いて算出する第１の算出工程と、複数の第２の小区画の体積に関する情報を用いて、複数の第２の注目対象を複数の第２の小区画に分配する計算を行う分配計算工程と、複数の第２の小区画の数と、分配計算工程により得られた、第２の注目対象を内包する第２の対象区画の数と、に基づいて、第１の算出工程と共通のアルゴリズムを一部に用いて、第２の注目対象の数に関する情報の第２の候補値を算出する第２の算出工程と、第２の候補値と第２の注目対象の数とを比較し、比較の結果が所定の条件を満たした場合の第２の候補値を、第１の注目対象の数に関する情報として決定する決定工程と、を含むことを特徴とする。

本発明によれば、粒子を用いた注目対象物を含む、検体となる溶液を分割した小区画の大きさにばらつきがあっても、比較的少ない数の小区画を用いて精度良く注目対象物の数、つまり濃度を推定することができる。

粒子を含む画像データの特徴を説明するための図である。図１（ａ）は、粒子をオイルの入ったシャーレ内に配置し、下面から照明光で照射して、上方から撮像した際の像を示す図である。図１（ｂ）は、光軸に平行な断面における粒子表面と光線の例を示す図である。図１（ｃ）は、明暗状態を表わす二種類の画素値を粒子の像に割りあてて生成した二値化画像データを示す図である。図１（ｄ）は、図１（ｃ）に示した二値化画像データの一部の領域に対して領域分割工程を適用した例を示す図である。 本発明の実施形態に係る粒子測定方法の処理手順の一例を示すフローチャートである。 本発明の実施形態に係る粒子測定方法の粒子位置取得工程における処理手順の一例を示すフローチャートである。 形状評価工程における形状評価の一例を説明するための図である。 粒子間の白領域が検出できる場合の一例を説明するための図である。図５（ａ）は、隣接粒子抽出工程における処理手順の一例を示すフローチャートである。図５（ｂ）は、二値化画像における白領域と黒領域との関係の一例を説明するための図である。 隣接粒子グラフのデータを説明するための図である。図６（ａ）は、ノードテーブルの一例を示す図である。図６（ｂ）は、エッジテーブルの一例を示す図である。 粒子間の白領域が検出できない場合の一例を説明するための図である。図７（ａ）は、隣接粒子抽出工程における処理手順の一例を示すフローチャートである。図７（ｂ）は、二値化画像における白領域と黒領域との関係の一例を説明するための図である。 本発明の実施形態に係る粒子測定方法の算出工程における処理手順の一例を示すフローチャートである。 粒子の大きさの推定に用いる連立方程式の一例を説明するための図である。図９（ａ）は、球体が配置されている二次元面に対して、垂直な方向から観察したときの球体の一例を示す図である。図９（ｂ）は、球体が配置されている面に平行で、球の接点と球の中心を含む面に垂直な方向から観察したときの球体の一例を示す図である。図９（ｃ）は、粒子半径差が大きい場合の連立方程式の一例を示す図である。図９（ｄ）は、粒子半径差が無視できる場合の連立方程式の一例を示す図である。 本発明の実施形態に係る標的核酸計数方法の処理手順の一例を示すフローチャートである。 比較例としてポアソン分布法及び重み付き平均法の推定精度を検証したシミュレーションの結果を示す図である。 本発明の実施形態に係る標的核酸計数方法の計数工程における処理手順の一例を示すフローチャートである。 計数工程内の分配計算工程におけるモンテカルロシミュレーションの処理手順の一例を示すフローチャートである。 本発明の実施形態に係る標的核酸計数方法の計数工程の推定精度を検証した結果の一例を示す図である。 第一の実施形態に係る濃度測定システムの構成の一例を示す図である。 第一の実施形態に係る濃度測定システムにおいて、粒子を保持する検査プレートの構造の一例を示す図である。図１６（ａ）は、検査プレートの断面図である。図１６（ｂ）は、検査プレートの平面図である。 第一の実施形態に係る濃度測定システムにおいて、ＰＣＲ法と蛍光撮像を実施する撮像装置の構成の一例を示す図である。 第二の実施形態に係る濃度測定装置の構成の一例を示す図である。 第二の実施形態に係る濃度測定装置が備える画像処理ボードの構成の一例を示すブロック図である。 第三の実施形態に係るデジタルＥＬＩＳＡ装置の構成の一例を示す図である。

本発明の実施の形態に係る粒子測定方法及び標的核酸計数方法の概要について説明する。ただし、本発明はこれらに限定されるものではない。

［粒子測定方法］
以下、図面を参照しながら、本発明の実施形態に係る粒子測定方法について説明する。本実施形態に係る粒子測定方法は、粒子を用いた注目対象物を計数することを含む。注目対象物としては、例えば、特定の塩基配列を持つ標的核酸、標的タンパク質、標的細胞、細菌、酵母、及び、微生物などが挙げられるが、これらに限定されるものではない。これらの注目対象物の中でも、特に、標的核酸及び標的タンパク質が好ましい。

＜粒子測定方法の測定対象＞
まず、本実施形態に係る粒子測定方法の測定対象について説明する。本実施形態に係る粒子測定方法により測定する粒子とは、例えば、注目対象物を含む検体を小区画に分割した微小な粒子であり、具体的には、標的核酸を含む溶液を分割して生成された液滴やエマルジョン、あるいはジェル等である。測定する粒子の大きさとは、球形に近い形状を持つ粒子、あるいは、一つの粒子を内部に含む球（以下、内包球と呼ぶ）の体積、断面積、直径・半径を表わす。計算の説明の都合上、上記球形粒子あるいは上記内包球を粒子と呼ぶ場合もある。本実施形態に係る粒子測定方法は、このような粒子を複数含む検体の画像データから、当該画像に含まれる粒子の個々の大きさを測定することができる。なお、上記粒子（小区画）のうち、注目対象物を内包する粒子を内包粒子（対象区画）と呼ぶ場合がある。

＜粒子測定方法で扱う画像データ＞
本実施形態に係る粒子測定方法に使用する画像データの見方について、図１を用いて説明する。図１は、粒子を含む画像データの特徴を説明するための図である。図１（ａ）は、粒子をオイルの入ったシャーレ内に配置し、下面から照明光で照射して、上方から撮像した際の像を示す図である。図１（ｂ）は、光軸に平行な断面における粒子表面と光線の例を示す図である。図１（ｃ）は、明暗状態を表わす二種類の画素値を粒子の像に割りあてて生成した二値化画像データを示す図である。図１（ｄ）は、図１（ｃ）に示した二値化画像データの一部の領域に対して領域分割工程を適用した例を示す図である。

以後、図１（ａ）に示されるような画像データを、粒子像データと呼ぶ。粒子像データから、各粒子の中央部、周辺部で濃淡が生じていることがわかる。濃淡の発生原因であるが、図１（ｂ）に示すように、粒子の中央部を通過する光線αが屈折せずにカメラの入射瞳に到達するのに対し、周辺部を通過する光線βは屈折や全反射により入射瞳に到達しないために発生する。

図１（ａ）に示される粒子周辺部の暗画素領域と、粒子外にあるオイルに相当する明画素領域との境界は、粒子と同じ直径を持つ円とは一致しない。例えば、撮影レンズの被写界深度が図１（ｂ）のΔｚで与えられる場合、カメラのイメージセンサ上では、Δｚの範囲にある断面（例えばｚ_０、ｚ_１）の光分布が加算されて受光される。断面ｚ_０には粒子直径に近い大きさの円形の光分布が形成されるが、断面ｚ_１には径が小さい光分布が形成される。これらの光分布が加算された場合、粒子直径よりも小さい円が形成される。

ただし、断面ｚ_０、ｚ_１は完全拡散面ではないので、上記説明は厳密なものではない。正確には、光線毎に物体面の高さｚが異なり、物体面ｚと光線の交点ｘに相当するイメージセンサ上位置の画素に光線が到達して光強度が加算される。隣接する粒子がある場合は、粒子表面の多重反射等も加わって複雑な分布が形成される。このように、撮影レンズの被写界深度が大きく粒子が隣接している状況では、粒子像データから輪郭情報を正確に得ることは困難である。

本実施形態に係る粒子測定方法では、粒子像データとは別に、図１（ｃ）に示すように、明暗状態を表わす二種類の画素値を粒子の像に割りあてて生成した二値化画像データを使用する。画像データに用いる二つの画素値は任意の値を選択できるが、明状態を白色、暗状態を黒色で表わし、各々の画素で構成される領域を白領域、黒領域と呼ぶこととする。図１（ｃ）に示した画像データ内では、粒子の中心部及び粒子間領域が白色の画素で表された白領域であり、粒子の周辺部及び粒子間領域が黒色の画素で表された黒領域である。

図１（ａ）に示した粒子像データは、透過照明を用いて撮像したが、例えば色素を混ぜた粒子を照明し、粒子からの散乱光を撮影することも可能である。この場合、粒子中央部、周辺部の濃淡は無くなって粒子全体が明るくなり、粒子間のオイルに相当する画素は暗くなる。濃淡分布と、粒子・オイルとの対応関係は変化するものの、粒子の輪郭の検出が困難であるのは透過照明の場合と同様である。二値化画像についても、明暗状態と粒子・オイルの対応関係は変化するが、後述する粒子位置取得工程の実施には影響しないため、本発明を適用できる。

＜粒子測定方法全体の流れ＞
本発明の実施形態に係る粒子測定方法０１０の手順を、図２を用いて説明する。図２は、本実施形態に係る粒子測定方法０１０の処理手順の一例を示すフローチャートである。

粒子位置取得工程Ｓ１００は、粒子像データから、各粒子の位置データを取得する工程である。ここで、各粒子の位置データとは、粒子像データ内の座標系、あるいは、物体面上・像面上の座標系における各粒子の中心に相当する座標値である。

隣接粒子抽出工程Ｓ１０１は、複数の粒子の中から、隣接している二つの粒子を抽出する工程である。隣接粒子抽出工程Ｓ１０１の詳細については後述する。

算出工程Ｓ１０２は、隣接した二つの粒子の位置データから距離を計算し、距離を基にした連立方程式を解いて粒子の大きさを算出する工程である。

＜粒子測定方法における粒子位置取得工程＞
粒子位置取得工程Ｓ１００の処理手順の一例を、図３を用いて説明する。図３は、本発明の実施形態に係る粒子測定方法０１０の粒子位置取得工程Ｓ１００における処理手順の一例を示すフローチャートである。本手順は、図１（ａ）に示した透過照明による粒子像データを入力とするが、後述する形状評価工程Ｓ１１３の評価値を変えることにより、色素からの散乱光を受光して生成したような、他の種類の粒子像データにも適用可能である。

二値化工程Ｓ１１０は、粒子像データから二値化画像データを生成する工程である。ここで、粒子の中心部及び粒子間領域を表わす白色の画素と、粒子の周辺部及び粒子間領域を表わす黒色の画素とに分類される。
領域分割工程Ｓ１１１は、二値化画像データの白領域をグループ化する工程である。例えば、図１（ｃ）に示した二値化画像データの一部の領域に対して領域分割工程を適用した場合、図１（ｄ）に示すように、白領域のグループが１８個あり、各々に０〜１７の番号が割り当てられる。

重心算出工程Ｓ１１２は、各々の白領域グループに対し、以下の数式１で与えられる重心（ｘ_ｗ，ｙ_ｗ）を算出する工程である。
ただし、ｘ_ｉ、ｙ_ｉ（ｉ＝０，…，ｎ−１）は、一つの白領域グループを構成する画素の座標値である。このようにして、粒子領域の画素の重心位置を、粒子位置として求めることができる。

形状評価工程Ｓ１１３では、各々の白領域グループに対し、形状評価値を求める。形状評価値の例として、扁平率について図４を用いて説明する。図４は、形状評価工程Ｓ１１３における形状評価の一例を説明するための図である。扁平率は、白領域の扁平度を表わす値であり、０に近ければ領域は円に近く、値が大きくなるほど偏りが大きいことを示す。図４に示すように、扁平率の計算では、対象の白領域グループの重心ｗから、複数の方向に沿って進み、一番端にある画素までの距離を求める。この例では上方向に沿った時の端の画素がuであり、左方向の端がvであり、uｗ間、vｗ間の距離を求める。得られた距離の中で最大の値をｄ_Ｌ、最小の値をｄ_Ｓとしたとき、扁平率Ｌは、以下の数式２で与えられる。
ただし、図１（ｄ）に示す白領域グループ０のように大きく湾曲した白領域の場合、重心が白領域の外に存在する場合がある。白領域グループ０の重心は白領域グループ１の近辺に存在するため、右方向に沿って調べても白領域グループ０の構成画素が存在しない。このような場合は、逆方向に沿って重心に近い方の端画素を求め、距離に対し、負の符号を与える。

形状評価値の種類は複数であっても良く、例えば白領域グループを構成する画素数や、重心から構成画素までの距離の平均値や標準偏差を用いることができる。

抽出工程Ｓ１１４では、形状評価値を基にして、粒子に相当する白領域グループを抽出する。粒子の白領域グループは必ずしも真円にはならず、扁平度から直接抽出することは難しいが、形状評価値を用いて粒子以外の白領域を除外し、消去法で残すことは容易である。例えば、図１（ｄ）に示す白領域グループ０の扁平率は１に近い値であり、粒子以外の領域として除外できる。また、白領域グループ２についても、属する画素の数は、粒子の白領域より明らかに少ないため、除外できる。このようにして、粒子の中心部を表わす部分領域である粒子領域を抽出することができる。以上から、粒子に相当する白領域グループの番号と、位置データである重心位置を得ることができる。

＜粒子測定方法における隣接粒子抽出工程＞
隣接粒子抽出工程Ｓ１０１の処理手順の一例を、図５を用いて説明する。図５は、粒子間に白領域が存在する場合の一例を説明するための図である。図５（ａ）は、隣接粒子抽出工程Ｓ１０１における処理手順の一例を示すフローチャートである。図５（ｂ）は、二値化画像における白領域と黒領域との関係の一例を説明するための図である。被写界深度が浅い撮像レンズと透過照明を用いて粒子像データを撮像した場合には、図５（ｂ）に示すような二値化画像が得られる。

粒子対選択工程Ｓ１２０は、粒子位置取得工程Ｓ１００により得られた粒子相当の白領域グループ（以後、粒子白領域と呼ぶ）の中から、距離が近い（位置データの差が小さい）二つの粒子白領域を選択する工程である。図５（ｂ）では、粒子白領域ｗ_０、ｗ_２が選択されたとしている。

ラインサーチ工程Ｓ１２１では、二値化画像内で、選択された二つの粒子白領域の中心位置（位置データに相当する座標）を結ぶ直線上の画素が属する白領域グループを調べる。二つの粒子白領域以外の白領域グループに属する画素がない場合は、二つの粒子は隣接していると判断する。図５（ｂ）では、粒子白領域ｗ_０、ｗ_２の中心位置がγ、δであり、γとδを結ぶ直線上の白色画素は、白領域グループｗ_１にも属するため、隣接していないと判断する。

ラインサーチ工程Ｓ１２１で隣接すると判断された白領域グループは、図６に示すようなグラフ構造を表現するテーブルに登録する。図６は、隣接粒子グラフのデータを説明するための図である。図６（ａ）はノードテーブル１３０の一例を示す図である。図６（ｂ）はエッジテーブル１３５の一例を示す図である。

グラフ構造とは、多数の要素（以後、ノードと呼ぶ）と、二つの要素間の接続（以後、エッジと呼ぶ）を表わすデータ構造であり、粒子の隣接関係を表現できる。図６（ａ）に示すノードテーブル１３０は、ノードである粒子を表現するテーブルであり、図６（ｂ）に示すエッジテーブル１３５はエッジである隣接関係を表わすテーブルである。以下、この２つのテーブルを合わせてグラフ構造テーブルと呼ぶ場合がある。

図６（ａ）及び図６（ｂ）では、一例として、各テーブルの要素に、図１（ｄ）に示した白領域グループの番号と隣接関係を当てはめている。白領域グループの登録では、まず、ノードテーブル１３０あるいはエッジテーブル１３５に、すでに行として記録されていないことを確認する。記録されていない場合には、記録済み行の次の行に白領域グループ情報を記録する。このようにして、隣接粒子抽出工程Ｓ１０１により隣接粒子のグラフ構造のデータが得られる。

粒子対の有無判定Ｓ１２２では、グラフ構造テーブルに登録されていない粒子白領域の有無を判定し、存在する場合は粒子対選択工程Ｓ１２０に戻り、存在しない場合はグラフ構造テーブルを隣接粒子抽出工程Ｓ１０１の出力とする。

隣接粒子抽出工程Ｓ１０１の手順の異なる例である隣接粒子抽出工程Ｓ１０３について、図７を用いて説明する。図７は、粒子間に白領域が存在しない場合、すなわち粒子間の白領域が検出できない場合の一例を説明するための図である。図７（ａ）は、隣接粒子抽出工程Ｓ１０３における処理手順の一例を示すフローチャートである。図７（ｂ）は、二値化画像における白領域と黒領域との関係の一例を説明するための図である。被写界深度が深い撮像レンズを用いて粒子像データを撮像した場合や、色素の散乱光を撮像した場合には、図７（ｂ）に示すような二値化画像が得られる。

図７（ａ）に示すフローチャートの粒子対選択工程Ｓ１２３及び粒子対の有無判定Ｓ１２６は、図５（ａ）に示す粒子間に白領域が存在する場合の粒子対選択工程Ｓ１２０及び粒子対の有無判定Ｓ１２２と同じである。白領域径算出工程Ｓ１２４では、二値化画像内で、選択された二つの粒子白領域の中心位置（位置データに相当する座標）を結ぶ直線上において、中心位置から離れた画素までの最大距離（以後、白領域径と呼ぶ）を算出する。図７（ｂ）では、粒子白領域ｗ_３、ｗ_４の中心位置がｓ、ｔであり、白領域径は、ｄ_３、ｄ_４である。

粒子距離算出工程Ｓ１２５では、二つの粒子白領域の中心位置の間の距離を算出し、白領域径算出工程Ｓ１２４で得られた白領域径の合計との比較を行う。白領域の大きさ次第ではあるが、中心位置の距離が白領域径の合計の２倍以下である場合は、隣接していると判断する。例えば、二つの粒子白領域の中心位置を結ぶ線上にある白領域が二つの粒子白領域のみに属し、かつ、二つの中心位置の間の距離との差が、中心位置の間の距離の半分より小さい場合に隣接していると判断する。なお、図７（ｂ）では、中心間距離Ｄが白領域径ｄ_３、ｄ_４の合計の２倍以上離れており、隣接していないと判断できる。粒子対の有無判定Ｓ１２６では、二つの粒子位置の間の距離の差が、粒子位置の間の距離の半分より小さい場合に隣接していると判断し、二つの粒子を抽出する。粒子距離算出工程Ｓ１２５で隣接すると判断された白領域グループは、図５（ａ）に示すラインサーチ工程Ｓ１２１と同様に、図６に示すグラフ構造を表現するテーブルに登録する。

＜粒子測定方法における算出工程＞
算出工程１０２の処理手順の一例を、図８を用いて説明する。図８は、本発明の実施形態に係る粒子測定方法０１０の算出工程Ｓ１０２における処理手順の一例を示すフローチャートである。

部分集合削除工程Ｓ１４０では、隣接粒子抽出工程Ｓ１０１の出力データであるノードテーブル１３０とエッジテーブル１３５から、粒子半径の算出に適さない部分グラフ構造を削除する。粒子半径の算出は、後述する連立方程式により行うが、連立方程式が優決定になる場合に半径を求めることができる。粒子の数が粒子対の数より多い場合に、部分集合に含まれる粒子に係る係数を連立方程式から除くことで、粒子半径の算出に適さない部分グラフ構造を削除することができる。つまり、互いに接続されているノードで構成される部分集合の内、エッジ数がノード数より少ない部分集合は劣決定になり、半径が求まらないため、削除対象となる。

削除対象の部分集合の抽出について、図６に戻って説明する。図６に示したノードテーブル１３０からノード番号の列のノード番号を一つ選択する。このノードと接続関係にあるノードが存在する場合には、エッジテーブル１３５の小ノード番号及び大ノード番号に選択されたノード番号が含む行が存在する。これらの行から接続関係にあるノード番号を取得できる。接続関係にあるノード番号に対して、エッジテーブル１３５の行を検索する処理を繰り返すことで、接続関係にあるノードとエッジを抽出できる。得られたノードの数とエッジの数とを比較し、エッジの数がノードの数より小さい場合には、削除対象と判定する。削除対象となったノードとエッジについては、関連する行をノードテーブル１３０とエッジテーブル１３５から除去し、ノード番号及びエッジ番号が連番となるように、番号を更新する。

解算出工程Ｓ１４１では、ノードテーブル１３０とエッジテーブル１３５の情報を基に、連立方程式を構成し、連立方程式を解くことで粒子半径を求める。

連立方程式の例について、図９を用いて説明する。図９は、粒子の大きさの推定に用いる連立方程式の一例を説明するための図である。図９（ａ）は、球体が配置されている二次元面に対して、垂直な方向から観察したときの球体の一例を示す図である。図９（ｂ）は、球体が配置されている面に平行で、球の接点と球の中心を含む面に垂直な方向から観察したときの球体の一例を示す図である。図９（ｃ）は、粒子半径差が大きい場合の連立方程式１５０の一例を示す図である。図９（ｄ）は、粒子半径差が無視できる場合の連立方程式１５１の一例を示す図である。

図９（ａ）及び図９（ｂ）に示したこれらの球体は、球形の粒子、あるいは、内包球に相当する。図９（ｂ）に示す球の半径ｒ_０，ｒ_１と中心間距離ｄ_０１の間の関係は、ピタゴラスの定理から、以下の数式３のように導出できる。

この方程式を複数の球体毎に作成することで、図９（ｃ）に示す連立方程式１５０を作成する。仮に球体毎の半径の差が小さく、球体毎の高さの差（図９（ｂ））の影響を無視できる場合は、中心間距離が半径の和に等しいと仮定できる。この場合は、図９（ｄ）に示される連立方程式１５１を用いることも可能である。

次に、図６に示したノードテーブル１３０とエッジテーブル１３５から、図９（ｃ）に示す連立方程式１５０を構成して、解を得るまでの手順の一例を説明する。連立方程式１５０を解くために必要とされるパラメータは、連立方程式１５０の左辺にあるベクトル（以下、距離ベクトルと呼ぶ）と右辺の行列（以下、隣接行列と呼ぶ）であり、最初に作成する。

距離ベクトルは、エッジの数と等しい長さを持つベクトルである。エッジテーブル１３５を上の行から順番にたどりながらエッジに対応する二つのノード番号を求め、二つのノード番号に対応する粒子位置１３３をノードテーブル１３０から取得する。粒子位置１３３の間の距離を基に距離ベクトルの要素を算出できる。

隣接行列は、列と行の要素数が、各々ノード数、エッジ数に対応している。隣接行列の作成では、距離ベクトルと同様に、エッジテーブル１３５を上の行から順番にたどりながらエッジに対応する二つのノード番号を求める。その中で、エッジ番号に対応する行において、ノード番号に対応する列の要素を１とし、他の要素を０とすることで、隣接行列を作成することができる。隣接行列は、大部分の要素が０であるスパース行列であるので、圧縮行格納方式（ＣＲＳ；Ｃｏｍｐｒｅｓｓｅｄ Ｒｏｗ Ｓｔｏｒａｇｅ）のように値が１の要素のインデクスを記録する方式を採用することも可能である。

距離ベクトルと隣接行列が求まった後には、連立方程式のソルバーを適用することで、連立方程式１５０の右辺のベクトルを求めることができる。球体の数が１万を超える場合、隣接行列のサイズが大きくなり、ＬＵ分解を用いたような一般的なソルバーは適用できない。このような場合は、ＬＳＱＲ法のような大規模スパース行列に適用可能なソルバーを使用する。得られた右辺のベクトルの要素を指数項とした指数関数を計算することで、各球体の半径を得ることができる。

このように、本発明の粒子測定方法０１０は、粒子の輪郭の測定が困難な画像データを基にして、粒子の大きさを得ることができる。

［標的核酸計数方法］
以下、図面を参照しながら、本発明の実施形態に係る標的核酸計数方法について説明する。

＜標的核酸計数方法で扱う画像データと測定対象＞
まず、本実施形態に係る標的核酸計数方法において扱う画像データについて説明する。粒子に励起光を照射した際に発生する蛍光の中で撮像した、標的核酸に関連する分光成分の光強度分布データ（以下、蛍光撮像データと呼ぶ）を入力データとして扱う。蛍光撮像データ内では、標的核酸を一つ以上含む粒子の画素値は、標的核酸を含まない粒子の画素値よりも大きくなるため、閾値処理等で両者を区別することができる。以下の説明では、蛍光撮像データ内で画素値が大きい粒子をポジ粒子、画素値が小さい粒子をネガ粒子と呼ぶ。

本実施形態に係る標的核酸計数方法の測定対象は、粒子生成の基になった溶液中に含まれる（基準体積あたりの）標的核酸数である。ＰＣＲ法適用後に増幅された標的核酸数や、蛍光発光しているポジ粒子数とは異なることに注意されたい。

＜標的核酸計数方法全体の流れ＞
本発明の実施形態に係る標的核酸計数方法０１１の処理手順の一例を、図１０を用いて説明する。図１０は、本実施形態に係る標的核酸計数方法０１１の処理手順の一例を示すフローチャートである。

ポジ粒子数算出工程Ｓ１６０は、蛍光撮像データからポジ粒子を抽出し、ポジ粒子数を算出する工程である。ポジ粒子の抽出では、粒子測定方法０１０により得られた粒子位置にある画素値、あるいは粒子位置近辺の領域の平均画素値を閾値と比較することで判定する。閾値は、事前の実験において、蛍光撮像装置を使用して得た蛍光撮像データの画素値の統計量から算出した値である。

計数工程Ｓ１６１は、ポジ粒子数算出工程Ｓ１６０により得られたポジ粒子数と、粒子測定方法０１０により測定された粒子数・大きさと、粒子の大きさに関する基準値とを基に、標的核酸の濃度（単位体積当たりの標的核酸数）を推定する工程である。

上記基準値は、事前の実験により得られた統計量であり、粒子の大きさの平均値、あるいは測定可能な粒子の総体積等が用いられる。粒子測定方法０１０では、各粒子の周囲を囲い、互いに接する球体の大きさが得られるが、粒子を隙間無しに接するよう配置することは困難である。粒子形状も完全な球ではないことから、粒子測定方法０１０のみでは、粒子の大きさの精度が過度に低下するおそれがある。以上の理由で、上記基準値による補正が必要とされる場合がある。

＜比較例によるシミュレーション＞
ここで、本発明の実施形態に係る標的核酸計数方法０１１における計数工程Ｓ１６１の詳細について説明する前に、比較例として、一般に用いられる二つの標的核酸濃度推定手法の課題について説明する。

＜ポアソン分布法（比較例１〜３）＞
通常のｄＰＣＲ法では、以下の数式４に示すポアソン分布を用いて標的核酸の濃度ｆの推定を行う（以下、ポアソン分布法と呼ぶ）。
ここで、Ｎ_Ａは粒子数、Ｎ_Ｐはポジ粒子数、Ｖ_ａは粒子の体積である。ポアソン分布法は、全ての粒子の体積が一定の場合に正常に機能し、粒子の大きさがばらついている場合には誤差を伴うことが知られている。

ポアソン分布法の推定誤差の影響を示すために行ったシミュレーションについて説明する。このシミュレーションは、後述する二つの標的核酸濃度推定方法（比較例である重み付き平均法及び実施例である計数工程Ｓ１６１）において、共通のパラメータを用いて実行し、各々の手法の性能比較を可能にする。シミュレーションでは、粒子の大きさに従った確率で、事前に決めた濃度の標的核酸を粒子に分配するモンテカルロシミュレーション（以下、モンテカルロ法と呼ぶ場合もある）を実行し、粒子数と、シミュレーションにより得たポジ粒子数から濃度を推定する。モンテカルロシミュレーションの結果であるポジ粒子数は試行毎にばらついているため、複数回実行することが好ましい。例えば、分配と推定を５０回試行し、濃度の平均推定値と標準偏差を算出することができる。ただし、試行回数については特に限定されるものではない。濃度の平均推定値と、シミュレーションで事前に与えた濃度真値の差が０に近く、濃度の標準偏差が０に近いほど、標的核酸濃度推定方法の性能が高いと言える。そこで、両者の評価値として、平均値誤差割合Ｔ（下記数式５）と、９５％信頼区間割合Ｓ（下記数式６）とを定義し、シミュレーション結果の比較に用いることができる。
ここで、ｇ_ｅｓｔは濃度の平均推定値、ｇ_ｏｒｇはシミュレーションで与えた濃度の真値である。
ただし、σは濃度の標準偏差である。

シミュレーションにおいて共通に用いたパラメータについて述べる。粒子数は１５００個に設定し、粒子のばらつきについては、粒子の半径の標準偏差を半径の２０％に設定する。上記分配と推定の実行回数は、５０回とする。粒子数に対するポジ粒子数の割合（以下、ポジ率と呼ぶ）を変化させて精度検証のシミュレーションを行った。図１１にシミュレーションの結果を示す。

図１１は、比較例としてポアソン分布法及び重み付き平均法の推定精度を検証したシミュレーションの結果を示す図である。それぞれの濃度推定法について、濃度目標値に対応するポジ率（％）の値を１％、１０％、５０％と変化させた場合における、濃度の目標値（ｃｏｐｉｅｓ／ｕｌ）、平均推定値（ｃｏｐｉｅｓ／ｕｌ）、標準偏差（ｃｏｐｉｅｓ／ｕｌ）、平均値誤差割合（％）、及び、９５％信頼区間の割合（％）の各数値が得られた。なお、重み付き平均法の詳細については後述する。

図１１に示すように、ポアソン分布法では、ポジ率が１％及び１０％の場合、真の標的核酸濃度である目標値に対する平均値誤差割合がそれぞれ１％未満及び約２％であり適用可能な範囲であるが、ポジ率が５０％の場合には、約１０％となり適用可能な範囲を超えていた。

なお、ｄＰＣＲ法を用いた従来技術による測定誤差は、９５％信頼区間で真値の±１０％以下に抑えられている。しかし、ポアソン分布法を用いた場合、測定誤差を含まない値である平均値誤差割合のみで、すでに１０％の誤差が発生することから、大きさがばらついている粒子に対しては、計数工程Ｓ１６１の代わりにポアソン分布法を適用することは困難である。

＜重み付き平均法（比較例４〜６）＞
次に、重み付き平均法について説明する。重み付き平均法は、一般に用いられる標的核酸濃度推定手法であり、大きさが同程度である粒子をグループ化し、グループ毎にポアソン分布法を適用する手法である。重み付き平均法の濃度推定値は、以下の数式７で与えられる。

ここで、Ｎ_Ａ ^（ｉ）、Ｎ_Ｐ ^（ｉ）、Ｖ_ａ ^（ｉ）（ｉ＝０，１，…，Ｍ−１）は、Ｍ個存在するグループの内、番号ｉのグループにおける粒子数、ポジ粒子数、粒子の平均体積である。ｗ^（ｉ）は、各グループにより得られる濃度に積算される重みであり、番号ｉのグループにおける粒子体積合計値を、全粒子の体積合計値で割った値である。

重み付き平均法の問題は、グループ数が多くなると、各グループに属する粒子数が減ることである。粒子数が減ると、数式７におけるグループ毎に得られるポジ率（Ｎ_Ｐ ^（ｉ）／Ｎ_Ａ ^（ｉ））の値が連続的では無くなり、量子化誤差の影響で推定精度が低下する。グループ毎に求められた精度の低い標的核酸推定値が合算され、平均濃度の推定精度も低下することになる。

重み付き平均法の推定誤差の程度を示すために行ったシミュレーション結果について、図１１を用いて説明する。図１１に示すように、重み付き平均法（比較例４〜６）では、ポジ率にかかわらず、平均値誤差割合は４％以下に抑えられていたが、９５％信頼区間割合はいずれも８％以上の値であり、ポジ率１％においては４２．９％と非常に大きな数値となった。重み付き平均法では、例えば、粒子数を増やすことで９５％信頼区間割合における誤差の影響を抑制することはできるが、粒子数１５００個と少ない場合には精度の低下が生じるため、計数工程Ｓ１６１の代わりに重み付き平均法を適用することは困難である。

＜標的核酸計数方法における計数工程＞
以下に、本発明の実施形態に係る標的核酸計数方法０１１における計数工程Ｓ１６１について詳細に説明する。

計数工程Ｓ１６１にて実行される濃度推定法は、ポアソン分布法と重み付き平均法に発生する精度低下の問題を解決するよう、新規開発した手法である。計数工程Ｓ１６１の処理手順の一例について、図１２を用いて説明する。図１２は、本発明の実施形態に係る標的核酸計数方法０１１の計数工程Ｓ１６１における処理手順の一例を示すフローチャートである。

第１候補算出工程Ｓ１８０は、ポジ粒子数、粒子数、粒子の大きさ（体積）、及び、上記大きさの基準値（体積基準値）から、数式４を用いて濃度の第１候補値を算出する工程である。数式４で必要とされる粒子の体積Ｖ_ａには、粒子の大きさ、あるいは基準値を用いて算出した平均粒子体積を使用する。

分配計算工程Ｓ１８１は、濃度の第２候補値を基にして算出された個数の標的核酸を、各粒子に分配するモンテカルロシミュレーションを実行する工程である。モンテカルロシミュレーションにより得られた標的核酸の分布から、ポジ粒子数を計算する。濃度の第２候補値が決まっていない場合は、任意の正の整数を初期値として採用することができる。真の濃度と離れた数値を初期値として採用すると収束が遅くなるため、第１候補値やポジ粒子数算出工程Ｓ１６０により得たポジ粒子数を初期値とするのが好ましい。モンテカルロシミュレーション方法の詳細については後述する。

第２候補算出工程Ｓ１８２は、分配計算工程Ｓ１８１により得られたポジ粒子数を用いて、第２候補値の更新値を算出する工程である。上記更新値は、以下の数式８で与えられる。
ここで、ｇ_ｉは分配計算工程Ｓ１８１で使用した第２候補値、ｇ_ｉ＋１は第２候補値の更新値である。ｆ、Ｎ_Ａ、Ｖ_ａは数式４と同じであるが、Ｎ_Ｐ’は分配計算工程Ｓ１８１で得られたポジ粒子数である。ｆ_０は第１候補値である。

収束判定Ｓ１８３は、分配計算工程Ｓ１８１で使用した第２候補値が所定の収束条件を満たすかを判定する。つまり、第２候補値と、第２候補算出工程Ｓ１８２により得られた更新値の差の絶対値が、許容値未満であるかを判定する。モンテカルロシミュレーションにより、上記更新値にばらつきが生じるため、必要以上に小さい許容値を採用することは困難である。通常は、第１候補値の０．１％程度の値を採用することが好ましい。

収束判定Ｓ１８３により許容値以上と判定された場合は、更新値を第２候補値として、分配計算工程Ｓ１８１、第２候補算出工程Ｓ１８２、及び、収束判定Ｓ１８３を繰り返す。収束判定Ｓ１８３により許容値未満と判定された場合、あるいは、繰り返し回数が一定回数以上の場合には、第２候補値を濃度推定値とする。計数工程Ｓ１６１の収束性と精度については後述するが、標的核酸濃度や粒子数に影響されずに、ポアソン分布法及び重み付き平均法よりも高い精度が得られる。

＜モンテカルロシミュレーション＞
計数工程Ｓ１６１内の分配計算工程Ｓ１８１において行うモンテカルロシミュレーションの処理手順について、図１３を用いて説明する。図１３は、計数工程Ｓ１６１内の分配計算工程Ｓ１８１におけるモンテカルロシミュレーションの処理手順の一例を示すフローチャートである。

配列初期化工程Ｓ１８５は、粒子に分配される核酸の個数を保持する整数配列（以下、核酸保持配列と呼ぶ）を準備し、配列要素の値を０に初期化する工程である。以下では、各粒子に０〜（粒子数−１）の番号（以下、粒子番号と呼ぶ）を割り当て、核酸保持配列の要素番号を粒子番号に対応づける方式、について説明する。この方式の核酸保持配列の数値は、配列要素と同じ粒子番号の粒子に含まれる核酸の数を表わす。

核酸分配工程Ｓ１８６は、粒子の大きさに関連した確率に従って、粒子番号（核酸保持配列の要素番号）を選択（確率選択）し、核酸保持配列要素（数値配列要素）の数値に１を加算（要素変更）する工程である。粒子番号の選択では、以下の数式９で与えられる確率関数Ｐ（ｑ）に関する方程式を解くことで、０〜１の値を持つ一様乱数ｓに対応する粒子番号ｑを求める。
ここで、Ｔ_ｉは粒子番号ｉの粒子の大きさから求めた比率で、球形粒子（又は、内包球）の体積を、体積合計相当値で割った値（体積比）である。体積合計相当値は、全球形粒子（又は、全内包球）の体積の合計に相当する値である。各球形粒子（又は、内包球）の体積の測定精度が低いときには、合計値が大きくずれる可能性がある。そのような場合には、上記基準値を使って体積合計相当値を求めてもよい。上記基準値が、統計的に求めた粒子の平均体積の場合は、基準値に粒子数を掛けた値を体積合計相当値として採用できる。基準値が粒子生成前の溶液量等の場合は、基準値から粒子生成時の損失を差し引いた値を、体積合計相当値として採用することができる。

確率関数Ｐ（ｑ）は離散値を返すため、連続値である乱数ｓを使って方程式を解くことはできず、近似解を求める必要がある。例えば、以下の数式１０で与えられる粒子番号ｑ_ａを近似解とすることができる。

近似精度は低下するが、ルックアップテーブル（以下、ＬＵＴと略記する）を作成して高速に近似解を求めることも可能である。例えば、十分に大きな要素数Ｎを持つ数値配列を使ってＬＵＴを作成する場合、以下の数式１１で与えられる要素番号ｐ_ｓから要素番号ｐ_ｌまでの要素の値を粒子番号ｑに設定する。
ここで、ｆｌｏｏｒ（ｘ）は、浮動小数点数ｘの小数点以下の切り捨てを行う関数である。０からＮ−１の間の要素番号を、一様乱数を使って生成することで、ＬＵＴの要素を粒子番号とすることができる。

核酸数判定Ｓ１８７は、核酸分配工程Ｓ１８６の実行回数と、第２候補値との比較を行う。核酸分配工程Ｓ１８６の実行回数が第２候補値と等しくなるまで、核酸分配工程Ｓ１８６を繰り返し実行する（対象散布）。

ポジ粒子数算出工程Ｓ１８８は、核酸保持配列から、ポジ粒子数を算出する（要素数算出）工程である。
試行回数判定Ｓ１８９は、配列初期化工程Ｓ１８５からポジ粒子数算出工程Ｓ１８８までの実行回数と、あらかじめ決められた試行回数との比較を行う。核酸分配工程Ｓ１８６において乱数を使用していることにより、試行毎に異なるポジ粒子数が得られる。
ポジ粒子数平均工程Ｓ１９０は、試行毎に得られたポジ粒子数の平均値を求める工程である。十分な回数の試行を実施することにより、ばらつきの少ないポジ粒子数平均値を得ることができる。

以上の説明では、粒子を小区画とし、注目対象物である標的核酸を内包する粒子を内包粒子（対象区画）として、核酸保持配列の要素番号を粒子番号に対応づける方式について述べたが、他の実装方式を用いてもよい。例えば、粒子よりも十分に小さい微小区画を考え、核酸保持配列の要素番号を微小区画の番号に対応づける方式も採用できる。この場合、複数の微小区画の集合が一つの粒子に対応することになる。核酸分配工程Ｓ１８６は、粒子の大きさに関する確率ではなく、一様確率に従って微小区画の番号を選択し、核酸保持配列の要素の値に１を加算する工程となる。粒子に対応する微小区画の要素値の合計が、粒子に含まれる核酸数となる。

＜計数工程の収束性と効果＞
計数工程Ｓ１６１で行う反復処理の収束性は、バナッハの不動点定理（縮小写像の定理）を用いて説明できる。バナッハの不動点定理は、以下の数式１２を満たす関数ｇ（ｘ）には、不動点ｘ_ｆ＝ｇ（ｘ_ｆ）が存在し、反復処理ｘ_ｉ＋１＝ｇ（ｘ_ｉ）（ｉ＝０，１，２，…）を実行することで不動点ｘ_ｆが得られることを示している。
ａ，ｂは、反復処理を実行する区間にある任意の二点である。

分配計算工程Ｓ１８１の入力値である第２候補値をｘ、分配計算工程Ｓ１８１と第２候補算出工程Ｓ１８２を連続して行う関数をｇとすると、以下に示す数式１３のように式を展開することができる。
Ｈ（ｘ）は分配計算工程Ｓ１８１により、第２候補値をｘからポジ粒子数ｙを算出する関数であり、ｈ（ｙ）は第２候補算出工程Ｓ１８２により、ポジ粒子数ｙから第２候補値の更新値を算出する関数である。Ｈ^−１（ｙ）は、分配計算工程Ｓ１８１であるＨ（ｘ）の仮想的な逆関数である。計数工程Ｓ１６１のアルゴリズムは、ｈ（ｙ）がＨ^−１（ｙ）になるように、つまりポジ粒子数ｙが測定で得られたポジ粒子数と一致するように、第２候補値ｘを調整すると考えてよい。ｆ、Ｎ_Ａ、Ｖ_ａは数式4と同じで、ｆ_０は第１候補値である。

一方、理想的な濃度推定式Ｈ^−１（ｙ）は、ポアソン分布法の数式４から、ポジ粒子数に依存した誤差を減算した以下の数式１４に示す関数で近似することができる。
ここで、Ｄ（ｙ）は、ポアソン分布法により発生する誤差を表わした関数である。Ｒ（ｙ）は、関数Ｈ（ｙ）中で実行されるモンテカルロシミュレーションで発生するポジ粒子数のばらつきの影響（濃度推定値の変動）を表わす。図１１に示した表において、Ｄ（ｙ）は平均推定値と目標値の差に相当し、Ｒ（ｙ）は標準偏差に相当する。

数式１４を数式１３に代入すると、以下の数式１５に示すように、ｆ（Ｎ_Ａ，ｙ，Ｖ_ａ）が相殺され、誤差と定数のみの式に変換できる。
Ｈ（ｘ）とＤ（ｙ）の合成関数をＤ_ｈ（ｘ）、Ｒ（ｘ）とＤ（ｙ）の合成関数をＲ_ｈ（ｘ）と表記した。

数式１５が、バナッハの不動点定理の条件である数式１２を満たすことを示す。数式１２に数式１５を代入すると、以下の数式１６に示す関係式となる。
２点ａ，ｂの値に差がある場合、数式１２の右辺では、大きい方の値（仮に、ａが大きいとする）が支配的になる。図１１に示した表のポアソン分布法に関する行（すなわち比較例１〜３）において、Ｄ_ｈ（ｘ）は平均推定値と目標値の差に相当し、Ｒ_ｈ（ｘ）は標準偏差に相当する。そのため、数式１２の左辺では、Ｄ_ｈ（ａ）が支配的になる。図１１の平均推定値の誤差割合を見ると、Ｄ_ｈ（ａ）の大きさはａの１０．６％以下の値とされており、数式１５は数式１２の条件を満たす。

一方、２点ａ，ｂの値がほぼ同じである場合、Ｄ_ｈ（ｘ）はｘの微小変化に対して急激に変化しないため、Ｄ_ｈ（ａ）−Ｄ_ｈ（ｂ）はほぼ０となる。その結果、モンテカルロシミュレーションによる値のばらつきであるＲ_ｈ（ａ）−Ｒ_ｈ（ｂ）を小さくすることで、数式１５は数式１２の条件を満たすことができる。モンテカルロシミュレーションの値のばらつきは、図１３を用いて説明した試行回数判定Ｓ１８９の試行回数を増やすことで抑制することができる。つまり、計数工程Ｓ１６１で行う反復処理の収束が進むにつれて、試行回数判定Ｓ１８９の試行回数を増やすことで精度を向上できる。

以上から、分配計算工程Ｓ１８１と第２候補算出工程Ｓ１８２とを連続して行う関数ｇ（ｘ）に反復処理を適用することで、不動点ｘ_ｆ＝ｇ（ｘ_ｆ）を求めることが可能である。得られる不動点ｘ_ｆは、モンテカルロシミュレーションにより得られたポジ粒子数が、測定時のポジ粒子数と一致するときの濃度推定値となる。

実際に計数工程Ｓ１６１の推定精度を検証したシミュレーション結果について、図１４を用いて説明する。図１４は、本発明の実施形態に係る標的核酸計数方法０１１の計数工程Ｓ１６１の推定精度を検証した結果の一例を示す図である。シミュレーション内容は、図１１に示したポアソン分布法と重み付き平均法のシミュレーションと同じである。実施例１〜３として、濃度目標値に対応するポジ率（％）の値を１％、１０％、５０％と変化させた場合における、濃度の目標値（ｃｏｐｉｅｓ／ｕｌ）、平均推定値（ｃｏｐｉｅｓ／ｕｌ）、標準偏差（ｃｏｐｉｅｓ／ｕｌ）、平均値誤差割合（％）、及び、９５％信頼区間の割合（％）の各数値が得られた。

図１４に示すように、平均値誤差割合は、ポジ率によらず１％以下に抑えられ、９５％信頼区間割合も７％以下に抑えられている。このように、本発明の標的核酸計数方法０１１を用いることで、粒子数が少ない場合でも、高い精度で標的核酸濃度を得ることができる。

［濃度測定システム及び装置］
以下に、本発明の実施形態に係る粒子測定方法又は核酸濃度測定方法を実装したソフトウェアを備える、濃度測定システム及び装置の実施形態について説明する。ただし、本発明はこれらに限定されるものではない。本発明の実施形態に係る濃度測定システム及び装置は、粒子を複数含む検体の画像データから、当該画像に含まれる粒子の個々の大きさを測定することができる。

［第一の実施形態］
本発明に係る第一の実施の形態である濃度測定システムについて、図１５を用いて説明する。図１５は、第一の実施形態に係る濃度測定システム２００の構成の一例を示す図である。図１５に示すように、本実施形態に係る濃度測定システム２００は、撮像装置２０１、画像処理装置２０２、タブレットＰＣ２０３、有線ＬＡＮ２０４、及び、無線ＬＡＮルーター２０６より構成されている。

撮像装置２０１は、例えば、イーサネット（登録商標）ケーブル等のネットワークケーブル２０５を介して、有線ＬＡＮ２０４に接続されている。撮像装置２０１はプレートホルダー２０８を備えており、プレートホルダー２０８には、複数の検査プレート２０９を設置することができる。撮像装置２０１の詳細な構成については後述する。

画像処理装置２０２は、ネットワークケーブル２１０を介して有線ＬＡＮ２０４に接続されている。また、画像処理装置２０２は、ネットワークケーブル２１１を介して無線ＬＡＮルーター２０６に接続されている。画像処理装置２０２は、無線ＬＡＮルーター２０６を介して無線通信２０７によりタブレットＰＣ２０３と通信することができる。

タブレットＰＣ２０３は、無線ＬＡＮルーター２０６を介して無線通信２０７により画像処理装置２０２と通信することができる。タブレットＰＣ２０３は、撮像装置の状態を表示する表示画面２１２を備える。表示画面２１２には、選択用リストボックス２１３、情報表示・修正部２１４、エディットボックス２１５、指示ボタン２１６を表示することができる。

＜濃度測定システムにおける処理の流れ＞
続いて、濃度測定システム２００で行われる処理の流れについて説明する。濃度測定システム２００の利用者により、まず、撮像装置２０１のプレートホルダー２０８に、検査プレート２０９が設置される。検査プレート２０９の詳細は後述する。撮像装置２０１は、設置後の検査プレート２０９の情報（プレートホルダー２０８内のスロット番号と、検査プレート２０９のバーコード情報）を画像処理装置２０２に送信する。その後、一定時間おきに画像処理装置２０２に作業指示の確認を行う待機モードに移行する。画像処理装置２０２は、Ｗｅｂサーバー機能を兼ねており、利用可能な撮像装置の情報（検査プレート２０９のバーコード情報から復元した識別情報、測定方法等）を利用者に公開している。

次に、濃度測定システム２００の利用者は、タブレットＰＣ２０３を利用して、画像処理装置２０２が利用者に公開しているＷｅｂサイトにログインし、撮像装置の状態を表示する表示画面２１２を開くことができる。利用者の操作によって、撮像装置２０１が選択用リストボックス２１３により選択された後、タブレットＰＣ２０３は、情報表示・修正部２１４に、検査プレート２０９の識別情報、測定方法、測定時間の予測値を表示する。ただし、検査プレート２０９にバーコードが添付されていない場合は、識別情報と測定方法のデフォルト値が表示されるので、操作により画面上で修正される。操作により、測定結果の識別情報（測定日時やラベル等）が入力用のエディットボックス２１５に入力された後、撮像開始・終了を指示する指示ボタン２１６がクリックされ、タブレットＰＣ２０３は、入力・修正された情報と撮像指示を画像処理装置２０２に送信する。

撮像装置２０１は、画像処理装置２０２から撮像指示を確認後、待機モードから撮像モードに移行する。撮像装置２０１の構成と、撮像モード時の動作については後述する。撮像装置２０１は、撮像を行いながら、検査プレート２０９の粒子像データと蛍光撮像データを画像処理装置２０２に送信し、全ての画像データの撮像・送信が完了後、待機モードに移行する。

画像処理装置２０２は、撮像装置２０１から送られてきた粒子像データに対して、図２に示した粒子測定方法０１０を実装したソフトウェアを実行し、粒子毎の位置と大きさを算出する。続いて、画像処理装置２０２は、撮像装置２０１から送られてきた蛍光撮像データと、粒子の位置・大きさに対して、図１０に示した標的核酸計数方法０１１を実装したソフトウェアを実行して濃度を算出する。ただし、上記ソフトウェアでは、粒子の大きさの基準値は使用していない。画像処理装置２０２は、算出された濃度値をＷｅｂサイトにて利用者に公開する。

濃度測定システム２００の利用者は、タブレットＰＣ２０３を利用して、画像処理装置２０２が利用者に公開しているＷｅｂサイトにログインし、情報表示・修正部２１４に新たに表示された濃度測定値を確認することができる。

＜検査プレート＞
検査プレート２０９の構造について、図１６を用いて説明する。図１６は、本実施形態に係る濃度測定システム２００において、粒子を保持する検査プレート２０９の構造の一例を示す図である。図１６（ａ）は検査プレート２０９の断面図である。図１６（ｂ）は検査プレート２０９の平面図である。検査プレート２０９は、標的核酸と試薬（プライマ、蛍光標識プローブ等のＰＣＲ法の実行に必要となる化学物質を含む薬剤）を混ぜた溶液２２０を、オイル２２１内で複数の粒子２２２として保持する器具である。検査プレート２０９には、注目対象を含む溶液を小区画に分割して生成した粒子を保持する検体を設置することができる。

次に、粒子２２２の生成方法について説明する。図１６（ａ）に示すように、まず、注射器等により溶液注入部２２３にオイル２２１を注入し、オイル吸出部２２５から空気を吸い出すことにより、オイル２２１を充填する。溶液２２０は、注射器等により加圧されて溶液注入部２２３に注入され、粒子生成部２２４が備える多孔質膜により多数の粒子２２２が生成される。本実施の形態では、多孔質膜を粒子生成に用いるが、マイクロ流路チップ等を用いてもよい。粒子２２２は、オイル２２１より比重が大きいため、下面のガラス基板２２６の上に粒子同士が重ならないように配置される。

粒子２２２は、ガラス基板２２６の変形等により容易に破壊されるため、検査プレート２０９の周囲は金属治具２２７で補強することが好ましい。ロボットにより検査プレート２０９を搬送する際には、金属治具２２７の窪み２２８にロボットハンドのツメを接触させて、ガラス基板２２６に力が加わらないようにすることができる。標的核酸（ＤＮＡ）の混入を防ぐため検査プレート２０９の再利用はできないが、金属治具２２７は、留め具２２９を外して分離することで再利用可能である。

検査プレート２０９には、溶液２２０の識別子や測定方法を符号化したバーコード２３０を添付することができる。バーコード２３０は、図１５に示した撮像装置２０１により撮像され、画像処理装置２０２で情報が復元される。バーコード２３０の情報は人間により識別できないことから、人間により識別可能なラベル２３１を添付することも可能である。

＜撮像装置の構成と撮像処理＞
図１５に示した撮像装置２０１の構成について、図１７を用いて説明する。図１７は、本実施形態に係る濃度測定システム２００において、ＰＣＲ法と蛍光撮像を実施する撮像装置２０１の構成の一例を示す図である。図１７に示すように、撮像装置２０１は、プレートホルダー２０８、ネットワークインターフェース２４０、コントローラー２４１、サーマルサイクラー２４２、カメラ２４３、撮像レンズ２４４、カメラ側フィルタホイール２４５、ｘｙステージ２４６、照明側フィルタホイール２４７、照明光学系２４８、ＬＥＤ光源２４９、単軸ロボット（ｘ軸方向駆動）２５０、単軸ロボット（ｚ軸方向駆動）２５１、単軸ロボット（ｙ軸方向駆動）２５２、及び、電源２５３により構成されている。電源２５３は、電源ケーブル２５４を介して外部の供給源に接続されている。電源２５３以外の構成要素は、それぞれネットワークケーブル２５５等を介してネットワークに接続されている。

続いて、検査プレート２０９がプレートホルダー２０８に設置された直後の動作について説明する。コントローラー２４１は、プレートホルダー２０８内のセンサーにより、検査プレート２０９が挿入されているスロットの位置を確認する。続いて、コントローラー２４１は、単軸ロボット２５０〜２５２を駆動してスロットに挿入されている検査プレート２０９を、ｘｙステージ２４６上に搬送し、バーコード２３０の撮影を行う。コントローラー２４１は、全ての検査プレート２０９のバーコード２３０を撮像後、検査プレート２０９が設置されているスロットの番号とバーコード画像を画像処理装置２０２に送信し、待機モードに移行する。

次に、撮像モードに移行後に行われる動作について説明する。コントローラー２４１は、タブレットＰＣ２０３により入力・修正された情報を、画像処理装置２０２から取得する。続いて、コントローラー２４１は、画像処理装置２０２から取得した情報と、すでに取得している検査プレート２０９の情報を基に、作業シーケンスを生成する。作業シーケンスは、単軸ロボットやカメラへ発行する命令と、命令の発行タイミングが記録されたデータである。命令の発行タイミングは、撮像装置２０１の部品毎の実行時間の差や、検査プレート２０９毎の測定方法の違いにより発生する待ち時間が最小になるように決定される。撮像・駆動命令だけでなく、ネットワークインターフェース２４０による画像データの送信命令も、作業シーケンスに含まれる。

理解しやすい例として、検査プレート２０９が一つで、蛍光撮像が一回行われるときの単純な作業シーケンスについて説明する。単軸ロボット２５０〜２５２は、検査プレート２０９をプレートホルダー２０８からサーマルサイクラー２４２に搬送する。単軸ロボット２５２の先には、フォークリフトのフォークのようなツメ（図示せず）が設けられている。単軸ロボット２５２は、金属治具２２７に設けられている窪み２２８に単軸ロボット２５２のツメを挿入し、上方に少し上げてから検査プレート２０９の抜き差しを行うことで、検査プレート２０９に振動を与えることなく搬送できる。次に、コントローラー２４１は、サーマルサイクラー２４２に対し、ペルチェ素子の加熱温度と冷却温度の制御命令を送信してＰＣＲ法を実行する。加熱温度と冷却温度は、出荷前に作成したテーブルの値を用いる。このように、ＰＣＲ法により粒子内の標的核酸の数を増幅させ、蛍光発光の強度を増大させることにより、標的核酸を含む粒子の検出及び測定が容易になる。

撮像装置２０１は、ＰＣＲ法の実行後、単軸ロボット２５２で検査プレート２０９をｘｙステージ２４６上に搬送させ、粒子像データの撮影を開始する。コントローラー２４１は、カメラ側フィルタホイール２４５に対し、干渉フィルタが挿入されていないスロットを選択するよう制御命令を送る。さらに、コントローラー２４１は、照明側フィルタホイール２４７に対し、ＮＤフィルタを選択するように制御命令を送信する。このとき、励起光は光強度が大きすぎて直接観察には適さないことから、ＮＤフィルタで減光することが好ましい。コントローラー２４１は、カメラ側フィルタホイール２４５及び照明側フィルタホイール２４７の設定後、ＬＥＤ光源２４９にシャッターを開放するための命令と、カメラ２４３に露光時間調整と撮像のための命令とを送る。コントローラー２４１は、撮影が完了してデータが保存されたタイミングで、カメラ２４３のメモリから粒子像データを取得し、ネットワークインターフェース２４０を通じて、画像処理装置２０２にデータを送信する。

次に、撮像装置２０１は、蛍光撮像データの撮像を行う。コントローラー２４１は、カメラ側フィルタホイール２４５に対し、励起光除去用干渉フィルタのスロットを選択するよう制御命令を送る。また、コントローラー２４１は、照明側フィルタホイール２４７に対し、励起光透過用干渉フィルタのスロットを選択するよう制御命令を送る。カメラ側フィルタホイール２４５及び照明側フィルタホイール２４７の設定後は、露光時間の数値が異なるものの、粒子像データの場合と同様である。

コントローラー２４１が、作業シーケンスに従い、命令の送信、及び画像データの送信を終えた後、撮像装置２０１は待機モードに移行する。

撮像装置２０１はサーマルサイクラー２４２を内部に備えているが、サーマルサイクラーは外部の別の装置であってもよい。例えば、外部のサーマルサイクラー装置（図示せず）により核酸数を増やした検査プレート２０９を、プレートホルダー２０８に挿入する場合は、撮像装置２０１は、撮像のみの構成でよい。

このように、本発明の第一の実施の形態は、被写界深度の影響等で、撮像による粒子体積の取得が困難でも、粒子のばらつきを補正できる標的核酸濃度測定システムである。つまり、本実施形態は、粒子輪郭を使わず、粒子の中心位置と隣接関係を用いて粒子の大きさを推定する手法と、モンテカルロ法により粒子補正精度を向上させた濃度推定手法により、推定精度を向上させた標的核酸濃度測定システムである。

［第二の実施形態］
本発明の第二の実施の形態である濃度測定装置について、図１８を用いて説明する。図１８は、本実施形態に係る濃度測定装置３００の構成の一例を示す図である。図１８に示すように、本実施形態に係る濃度測定装置３００は、第一の実施の形態である濃度測定システム２００の撮像装置２０１に、画像処理装置２０２と無線ＬＡＮルーター２０６の機能を追加した構成を備えている。機能的には両者はほぼ同じである。以下では、濃度測定システム２００と異なる点についてのみ説明する。

濃度測定装置３００は、撮像装置２０１が備えるネットワークインターフェース２４０の代わりに、無線ＬＡＮインターフェース３０２を備えている。画像処理装置２０２で担っていたＷｅｂサーバー機能を実行するために、コントローラー３０３は、例えば、ＣＰＵとメモリを備えた汎用ＰＣのような構成とすることができる。画像処理装置２０２で担っていた、粒子測定方法０１０と標的核酸計数方法０１１の実行は、画像処理ボード３０１で実行することができる。

次に、画像処理ボード３０１の構成について、図１９を用いて説明する。図１９は、本実施形態に係る濃度測定装置３００が備える画像処理ボード３０１の構成の一例を示すブロック図である。図１９に示すように、画像処理ボード３０１は、通信インターフェース３１０、ユニットコントローラー３１１、データバス３１２、メモリ３１３、及び、各種計算ユニットにより構成されている。画像処理ボード３０１は、計算ユニットとして、二値化部３１４、領域分割部３１５、重心算出部３１６、形状評価部３１７、抽出部３１８、隣接粒子抽出部３１９、部分集合削除部３２３、解算出部３２４、ポジ粒子数算出部３２５、及び、計数処理部３２６を備えている。これらの計算ユニットは、粒子測定方法０１０と標的核酸計数方法０１１とが備える複数の工程を実行する。

続いて、本実施形態において、図２に示した粒子測定方法０１０の実行に関する計算ユニットについて説明する。二値化部３１４、領域分割部３１５、重心算出部３１６、形状評価部３１７、及び、抽出部３１８は、粒子位置取得工程Ｓ１００を構成する５つの工程（図３を参照）を実行する計算ユニットである。隣接粒子抽出部３１９は、隣接粒子抽出工程Ｓ１０１を実行する計算ユニットである。部分集合削除部３２３、及び、解算出部３２４は、算出工程Ｓ１０２を構成する二つの工程（図８を参照）を実行する計算ユニットである。

本実施形態において、図１０に示した標的核酸計数方法０１１の実行に関する計算ユニットについて説明する。ポジ粒子数算出部３２５は、ポジ粒子数算出工程Ｓ１６０を実行する計算ユニットである。計数処理部３２６は、計数工程Ｓ１６１を実行する工程である。

次に、画像処理ボード３０１の構成要素の間において行うデータ交換について説明する。コントローラー３０３から送られてきた計算開始指示は、通信インターフェース３１０を介して、ユニットコントローラー３１１に送られる。同時に、計算対象となる画像データ等は、コントローラー３０３からデータバス３１２を介してメモリ３１３に送られ、メモリ３１３内で保持される。ユニットコントローラー３１１は、計算ユニットの実行開始指示・終了確認を行う。さらに、ユニットコントローラー３１１は、計算ユニット間での混信防止のための同期信号情報（データバス３１２のクロック）の割当処理（調停処理）も行う。具体的には、ユニットコントローラー３１１は、計算を開始する計算ユニットとメモリ３１３の組に対し、コントローラー・ユニット間通信路３２９を介して計算開始指示と、割り当てられた同期信号情報を送信する。

計算ユニットは、メモリ３１３内に保持された画像データや計算時に生成した数値配列中のデータに対して、読み取り・修正を行いながら計算を実行する。計算ユニットは、計算終了時に、終了確認情報をユニットコントローラー３１１に送信する。ユニットコントローラー３１１は、終了確認情報を受け取り後に、次の計算工程を実行する計算ユニットに実行開始指示を送ることを繰り返す。ユニットコントローラー３１１は、粒子測定方法０１０と標的核酸計数方法０１１の全処理を実行後、得られた濃度情報をコントローラー３０３に送信する。

濃度測定装置３００で行う処理内容については、濃度測定システム２００と異なる点が一つある。濃度測定装置３００の画像処理ボード３０１は、画像処理装置２０２に比べ処理能力が低い場合がある。大型の画像データを扱えないことから、粒子測定方法０１０の精度向上が困難である。そこで、計数処理部３２６では、過度な精度低下を防ぐために、統計的に測定した粒子の大きさの平均値を基準値として使用することができる。

このように、本発明の第二の実施の形態は、被写界深度の影響等で、撮像による粒子体積の取得が困難でも、粒子のばらつきを補正できる標的核酸濃度測定装置である。また、本実施形態は、画像処理機能を装置内に内蔵しており、ネットワーク接続が困難な場所でも使用できる標的核酸濃度測定装置である。

［第三の実施形態］
本発明の第三の実施の形態に係るデジタルＥＬＩＳＡ装置について説明する。本実施形態に係るデジタルＥＬＩＳＡ装置は、デジタルＥＬＩＳＡ法を用いた濃度測定装置である。デジタルＥＬＩＳＡ装置の構成について説明する前に、デジタルＥＬＩＳＡ法について述べる。

デジタルＥＬＩＳＡ法は、溶液中の標的タンパク質の濃度を測定する手法である。測定対象の溶液に、標的タンパク質に付着した際に蛍光物質を生成する抗体を混ぜて、タンパク質と抗体を反応させる。抗体には酵素が付けてあり、酵素反応によって生じた生成物が発色又は発光する。その光を測定することにより、抗体に結合したタンパク質を検出することができる。また、抗体さえあれば、どのようなバイオマーカーであっても検出することができる。次に、溶液を液滴に分割した後、蛍光像を撮像して、蛍光を発する液滴数からタンパク質濃度を算出することができる。デジタルＥＬＩＳＡ法では、ＰＣＲ法とは異なり、標的タンパク質の数の増幅手段は無いため、反応後の液滴を直接撮像することが好ましい。

本発明の第三の実施の形態について、図１を用いて説明する。図２０は、第三の実施形態に係るデジタルＥＬＩＳＡ装置４００の構成の一例を示す図である。図２０に示すように、デジタルＥＬＩＳＡ装置４００は、第二の実施の形態である濃度測定装置３００から、サーマルサイクラー２４２を除いた構成を持つ装置である。以下では、濃度測定装置３００と異なる点のみ説明する。

検査プレート４０１には、注目対象を含む溶液を小区画に分割して生成した粒子を保持する検体を設置することができる。検査プレート４０１は、検査プレート２０９と同じ構造であるが、標的タンパク質と抗体を含む溶液を液滴（つまり粒子）の形で保持する。抗体には酵素が付けてあり、酵素反応によって生じた生成物が発色又は発光する。その光を測定することにより、注目対象物である標的タンパク質を含む粒子を検出することができる。

コントローラー４０２は、コントローラー３０３と異なり、サーマルサイクラー２４２関連の制御・搬送を除いた作業シーケンスを生成して実行する。デジタルＥＬＩＳＡ装置４００は、標的タンパク質の数の増幅手段を持たず、蛍光強度が弱いことから、高感度カメラ４０３と顕微鏡対物レンズ４０４を、カメラ２４３と撮像レンズ２４４の代わりに搭載することが好ましい。

このように、本発明の第三の実施の形態は、被写界深度の影響等で、撮像による粒子体積の取得が困難な状況でも、標的タンパク質の濃度の測定が可能な標的タンパク質濃度測定装置である。

以上、本発明の好ましい実施形態について詳述したが、本発明は係る特定の実施形態に限定されるものではなく、特許請求の範囲に記載された本発明の要旨の範囲内において、種々の変形・変更が可能である。

本発明の計数方法、濃度測定装置、及び濃度測定システムは、粒子を用いた注目対象物の計数及び濃度測定に広く利用可能であり、特に、医療分野における病理検査などの測定に有用である。

０１０…粒子測定方法
０１１…標的核酸計数方法

Ｓ１００…粒子位置取得工程
Ｓ１０１，Ｓ１０３…隣接粒子抽出工程
Ｓ１０２…算出工程

Ｓ１１０…二値化工程
Ｓ１１１…領域分割工程
Ｓ１１２…重心算出工程
Ｓ１１３…形状評価工程
Ｓ１１４…抽出工程

Ｓ１２０…粒子対選択工程
Ｓ１２１…ラインサーチ工程
Ｓ１２２…粒子対の有無判定

Ｓ１２３…粒子対選択工程
Ｓ１２４…白領域径算出工程
Ｓ１２５…粒子距離算出工程
Ｓ１２６…粒子対の有無判定

１３０…ノードテーブル
１３５…エッジテーブル

Ｓ１４０…部分集合削除工程
Ｓ１４１…解算出工程

１５０…粒子半径差が大きい場合の連立方程式
１５１…粒子半径差が無視できる場合の連立方程式

Ｓ１６０…ポジ粒子数算出工程
Ｓ１６１…計数工程

Ｓ１８０…第一候補算出工程
Ｓ１８１…分配計算工程
Ｓ１８２…第二候補算出工程
Ｓ１８３…収束判定

Ｓ１８５…配列初期化工程
Ｓ１８６…核酸分配工程
Ｓ１８７…核酸数判定
Ｓ１８８…ポジ粒子数算出工程
Ｓ１８９…試行回数判定
Ｓ１９０…ポジ粒子数平均工程

２００…濃度測定システム
２０１…撮像装置
２０２…画像処理装置
２０３…タブレットＰＣ
２０４…有線ＬＡＮ
２０５，２１０，２１１，２５５…ネットワークケーブル
２０６…無線ＬＡＮルーター
２０７…無線通信
２０８…プレートホルダー
２０９…検査プレート

２１２…表示画面
２１３…選択用リストボックス
２１４…情報表示・修正部
２１５…エディットボックス
２１６…指示ボタン

２２０…溶液
２２１…オイル
２２２…粒子
２２３…溶液注入部
２２４…粒子生成部
２２５…オイル吸出部
２２６…ガラス基板
２２７…金属治具
２２８…窪み
２２９…留め具
２３０…バーコード
２３１…ラベル

２４０…ネットワークインターフェース
２４１…コントローラー
２４２…サーマルサイクラー
２４３…カメラ
２４４…撮像レンズ
２４５…カメラ側フィルタホイール
２４６…ｘｙステージ
２４７…照明側フィルタホイール
２４８…照明光学系
２４９…ＬＥＤ光源
２５０…単軸ロボット（ｘ軸方向駆動）
２５１…単軸ロボット（ｚ軸方向駆動）
２５２…単軸ロボット（ｙ軸方向駆動）
２５３…電源
２５４…電源ケーブル

３００…濃度測定装置
３０１…画像処理ボード
３０２…無線ＬＡＮインターフェース
３０３…コントローラー

３１０…通信インターフェース
３１１…ユニットコントローラー
３１２…データバス
３１３…メモリ
３１４…二値化部
３１５…領域分割部
３１６…重心算出部
３１７…形状評価部
３１８…抽出部
３１９…隣接粒子抽出部
３２３…部分集合削除部
３２４…解算出部
３２５…ポジ粒子数算出部
３２６…計数処理部
３２９…コントローラー・ユニット間通信路

４００…デジタルＥＬＩＳＡ装置
４０１…検査プレート
４０２…コントローラー
４０３…高感度カメラ
４０４…顕微鏡対物レンズ

Claims (14)

1. 複数の第１の注目対象を含む検体を分割した複数の第１の小区画のうち、前記第１の注目対象を内包する前記第１の小区画を第１の対象区画とし、該第１の対象区画に含まれる当該第１の注目対象の数に関する情報を取得する計数方法であって、
   前記第１の小区画の数と、前記第１の対象区画の数と、に基づいて、前記第１の注目対象の数に関する情報の第１の候補値を、アルゴリズムを用いて算出する第１の算出工程と、
   複数の第２の小区画の体積に関する情報を用いて、複数の第２の注目対象を前記複数の第２の小区画に分配する分配計算を行う分配計算工程と、
   前記複数の第２の小区画の数と、前記分配計算工程により得られた、前記第２の注目対象を内包する第２の対象区画の数と、に基づいて、前記アルゴリズムと共通のアルゴリズムを一部に用いて、前記第２の注目対象の数に関する情報の第２の候補値を算出する第２の算出工程と、
   前記第２の候補値と前記第２の注目対象の数とを比較し、前記比較の結果が所定の条件を満たした場合の前記第２の候補値を、前記第１の注目対象の数に関する情報として決定する決定工程と、
   を含むことを特徴とする計数方法。
2. 前記分配計算工程は、
   前記第２の小区画の数以上の個数の要素を持つ数値配列に対し、前記要素の値を初期化する初期化工程と、
   前記第２の小区画の体積比に基づいた確率で、前記第２の小区画の一つを選択する確率選択工程と、
   選択された前記第２の小区画に対応した数値配列要素の値を変更する要素変更工程と、
   前記確率選択工程と前記要素変更工程を、前記第２の注目対象の数と等しい回数となるまで繰り返し実行する対象散布工程と、
   前記初期化された要素の初期値と異なる値を持つ前記数値配列要素の数を算出する要素数算出工程と、
   を含むことを特徴とする請求項１に記載の計数方法。
3. 前記分配計算工程は、
   前記対象散布工程と前記要素数算出工程とを複数回実行し、
   前記初期値と異なる値を持つ要素数の平均値を求める平均工程を、
   を更に含むことを特徴とする請求項２に記載の計数方法。
4. 前記比較の結果が前記所定の条件を満たさない場合には前記所定の条件を満たすまで、前記第２の候補値を前記複数の第２の注目対象の数として、前記分配計算工程および前記第２の算出工程を反復する反復工程を、
   更に含むことを特徴とする請求項１から３のいずれか一項に記載の計数方法。
5. 前記反復工程において、
   第２の注目対象の数と、第２の候補値との違いに基づいて、前記平均工程における前記複数回実行の回数、又は、前記平均工程の実行の有無、を切り換えることを特徴とする請求項４に記載の計数方法。
6. 前記分配計算工程で分配する第２の注目対象の数の初期値は、第１の候補値であることを特徴とする請求項１から５のいずれか一項に記載の計数方法。
7. 前記アルゴリズムは、ポアソン分布を使用したアルゴリズムであることを特徴とする請求項１から６にいずれか一項に記載の計数方法。
8. 画像データから、第１の粒子の数と、前記第１の粒子の体積比、及び、第１の注目対象を内包する第１の内包粒子の数を算出する粒子数算出手段と、
   前記第１の粒子の数と、前記第１の内包粒子の数とから、前記第１の注目対象の数に関する、第１の候補値を求める第１の算出手段と、
   第２の粒子の体積に関する情報を用いて、第２の注目対象を複数の第２の粒子に分配する分配計算を行う分配計算手段と、
   前記第２の粒子の数と、前記分配計算手段により得られた、前記第２の注目対象を内包する第２の内包粒子の濃度に関する、第２の候補値を求める第２の算出手段と、
   前記第２の候補値と前記第２の注目対象の数とを比較し、前記比較の結果が所定の条件を満たした場合の前記第２の候補値を、前記第１の注目対象の濃度に関する情報として決定する決定手段と、
   を備えることを特徴とする濃度測定装置。
9. 前記分配計算手段は、前記第２の粒子の体積比と、所定の体積基準値とに関する情報を用いて、前記第２の注目対象を複数の前記第２の粒子に分配する計算を行うことを特徴とする請求項８に記載の濃度測定装置。
10. 前記比較の結果が前記所定の条件を満たさない場合には前記所定の条件を満たすまで、前記第２の候補値を前記第２の注目対象の数として、前記分配計算及び前記第２の候補値を求める算出を繰り返し実行することを特徴とする請求項８又は９に記載の濃度測定装置。
11. 前記第１の注目対象を含む溶液を小区画に分割して生成した微小な粒子を前記第１の粒子として保持する検体と、
    前記検体を基にして前記画像データを撮像する撮像装置と、
    請求項８から１０のいずれか一項に記載の濃度測定装置を備え、かつ、前記画像データから、前記第１の注目対象の濃度を算出する画像処理装置と、
    を備えることを特徴とする濃度測定システム。
12. 前記撮像装置は、前記第１の粒子内の前記第１の注目対象の数を増幅させる機能を備えることを特徴とする請求項１１に記載の濃度測定システム。
13. 前記第１の注目対象は、標的タンパク質であることを特徴とする請求項１１に記載の濃度測定システム。
14. 前記第１の注目対象は、特定の塩基配列を持つ標的核酸であることを特徴とする請求項１１又は１２に記載の濃度測定システム。