WO2020054292A1

WO2020054292A1 - スペクトル校正装置及びスペクトル校正方法

Info

Publication number: WO2020054292A1
Application number: PCT/JP2019/031581
Authority: WO
Inventors: 洋平川口; ジョージチャルキディス; 安田　知弘
Original assignee: 株式会社日立ハイテクノロジーズ
Priority date: 2018-09-10
Filing date: 2019-08-09
Publication date: 2020-03-19
Also published as: US20210310951A1; JP2020041876A; DE112019003466T5; DE112019003466B4; JP7022670B2; GB202102022D0; GB2590833B; US11692937B2; CN112513619B; GB2590833A; CN112513619A

Abstract

蛍光色素間でピークが同時に出現する場合であってもスペクトル校正を高精度に行えるスペクトル校正装置及びスペクトル校正方法を提供する。色変換処理に用いられる色変換行列を算出するスペクトル校正装置であって、経時的に検出された蛍光の分光信号を取得する分光信号取得部と、前記分光信号に基づいて、蛍光色素の蛍光のピークが同時に出現する頻度に依存するパラメータの値毎に色変換行列の候補を算出する候補算出部と、各候補に対して算出される評価値に基づいて色変換行列を選択する選択部を備えることを特徴とする。

Description

スペクトル校正装置及びスペクトル校正方法

　本発明はキャピラリ電気泳動型シーケンサのスペクトル校正装置及びスペクトル校正方法に関する。

　キャピラリ電気泳動型シーケンサは、キャピラリと呼ばれる細管の中で核酸試料を電気泳動させながら、励起光の照射により生じる蛍光を分光して経時的に検出し、検出された信号を用いて核酸を分析する装置である。具体的には、核酸を構成する塩基配列を決定するシーケンシング処理やＤＮＡ断片の鎖長を比較するフラグメント解析処理が行われる。シーケンシング処理やフラグメント解析処理の過程では色変換処理が実施される。色変換処理とは、分光された信号の時系列データを、核酸試料を標識する蛍光色素の蛍光の時系列データに変換する処理であり、各蛍光色素のスペクトル応答特性に対応する色変換行列の擬似逆行列を検出された信号の時系列データに乗じて行われる。

　色変換行列を取得する処理はスペクトル校正と呼ばれ、特許文献１にはスペクトル識別可能な複数の分子化学種から生成される電気泳動図を用いたスペクトル校正方法が開示されている。すなわち、電気泳動図から生成されるデータマトリクスの部分をランク解析することによりスペクトル応答特性を有する領域を識別し、識別された領域から判定される純成分スペクトル応答をクラスタに分け、各クラスタの代表的なスペクトル応答と分子化学種を相関付ける。

国際公開ＷＯ０１／０５５６９９号公報

　しかしながら特許文献１では、電気泳動図の時間領域において、蛍光色素を表す各ピークが十分に分離されていると仮定しており、蛍光色素間でピークが同時に出現する場合にスペクトル校正の精度が低下する。キャピラリ電気泳動型シーケンサでは装置の小型化にともないキャピラリが短くなるので、キャピラリ内で塩基が分離されにくくなり、検出される信号のピークの重なりが増している。

　そこで本発明は、蛍光色素間でピークが同時に出現する場合であってもスペクトル校正を高精度に行えるスペクトル校正装置及びスペクトル校正方法を提供することを目的とする。

　上記目的を達成するために本発明は、色変換処理に用いられる色変換行列を算出するスペクトル校正装置であって、経時的に検出された蛍光の分光信号を取得する分光信号取得部と、前記分光信号に基づいて、蛍光色素の蛍光のピークが同時に出現する頻度に依存するパラメータの値毎に色変換行列の候補を算出する候補算出部と、各候補に対して算出される評価値に基づいて色変換行列を選択する選択部を備えることを特徴とする。

　また本発明は、色変換処理に用いられる色変換行列を算出するスペクトル校正方法であって、経時的に検出された蛍光の分光信号を取得する取得ステップと、前記分光信号に基づいて、蛍光色素の蛍光のピークが同時に出現する頻度に依存するパラメータの値毎に色変換行列の候補を算出する候補算出ステップと、各候補に対して算出される評価値に基づいて色変換行列を選択する選択ステップを備えることを特徴とする。

　本発明によれば、蛍光色素間でピークが同時に出現する場合であってもスペクトル校正を高精度に行えるスペクトル校正装置及びスペクトル校正方法を提供することが可能になる。

スペクトル校正装置のハードウェア構成図である。検出された分光信号の例と、分光信号と色変換行列、蛍光の時間強度の関係を説明する図である。実施例１の機能ブロック図である。実施例１の処理の流れの一例を示す図である。実施例１の初期値の算出処理の流れの一例を示す図である。実施例１の色変換行列の候補の算出処理の流れの一例を示す図である。実施例１の候補の選択処理の流れの一例を示す図である。実施例２の候補の選択処理の流れの一例を示す図である。実施例３の候補の選択処理の流れの一例を示す図である。移動度補正処理を説明する図である。プルアップやプルダウンを含む色変換信号に移動度補正処理を施したときの例を説明する図である。

　以下、添付図面に従って本発明に係るスペクトル校正装置及びスペクトル校正方法の好ましい実施例について説明する。なお、以下の説明及び添付図面において、同一の機能構成を有する構成要素については、同一の符号を付することにより重複説明を省略することにする。

　図１はスペクトル校正装置１のハードウェア構成を示す図である。スペクトル校正装置１は、ＣＰＵ(Central Processing Unit)２、主メモリ３、記憶装置４、表示メモリ５、表示装置６、コントローラ７、入力装置８、ネットワークアダプタ１０がシステムバス１１によって信号送受可能に接続されて構成される。スペクトル校正装置１は、ネットワーク１２を介してキャピラリ電気泳動型シーケンサ１３や分光信号データベース１４と信号送受可能に接続される。ここで、「信号送受可能に」とは、電気的、光学的に有線、無線を問わずに、相互にあるいは一方から他方へ信号送受可能な状態を示す。

　ＣＰＵ２は、各構成要素の動作を制御する装置である。ＣＰＵ２は、記憶装置４に格納されるプログラムやプログラム実行に必要なデータを主メモリ３にロードして実行する。記憶装置４は、ＣＰＵ２が実行するプログラムやプログラム実行に必要なデータを格納する装置であり、具体的にはハードディスクやＳＳＤ（Solid State Drive）等の記録装置や、ＩＣカード、ＳＤカード、ＤＶＤ等の記録媒体に読み書きする装置である。各種データはＬＡＮ(Local Area Network)等のネットワーク１２を介して送受信される。主メモリ３は、ＣＰＵ２が実行するプログラムや演算処理の途中経過を記憶するものである。

　表示メモリ５は、液晶ディスプレイ等の表示装置６に表示するための表示データを一時格納するものである。入力装置８は、操作者がスペクトル校正装置１に対して操作指示を行う操作デバイスであり、具体的にはキーボードやマウス等である。マウスはトラックパッドやトラックボールなどの他のポインティングデバイスであっても良い。コントローラ７は、キーボードやマウスの状態を検出して、検出した情報等をＣＰＵ２へ出力するものである。また表示装置６がタッチパネルである場合には、タッチパネルが入力装置８としても機能する。ネットワークアダプタ１０は、スペクトル校正装置１をＬＡＮ、電話回線、インターネット等のネットワーク１２に接続するためのものである。

　キャピラリ電気泳動型シーケンサ１３は、キャピラリの中で核酸試料を電気泳動させながら、励起光の照射により生じる蛍光を分光して経時的に検出し、検出された信号である分光信号を用いて核酸を分析する装置である。核酸の分析には、分光信号に色変換行列の擬似逆行列を乗じることにより、核酸試料を標識する蛍光色素の蛍光の時系列データに変換する色変換処理が実施される。分光信号データベース１４は、キャピラリ電気泳動型シーケンサ１３によって取得された分光信号を保管するデータベースシステムである。

　図２を用いて、分光信号と色変換行列、蛍光の時間強度の関係について説明する。図２（ａ）は検出された分光信号の一例であり、縦軸が信号強度、横軸が時間であり、分光された波長に対応する検出チャンネル毎の波形が示されている。

　図２（ｂ）は、分光信号Ｆと色変換行列Ｓ、蛍光の時間強度Ｃの関係を示す図である。キャピラリの中を電気泳動により移動する核酸試料には、四種の塩基Ａ（アデニン）、Ｇ（グアニン）、Ｃ（シトシン）、Ｔ（チミン）が含まれ、各塩基の移動順を知ることにより核酸試料中の塩基配列を把握できる。各塩基には異なる蛍光色素が対応するので、各蛍光色素の蛍光の時間強度Ｃから各塩基の移動順を知ることにより塩基配列が決定できる。

　キャピラリ電気泳動型シーケンサ１３によって検出される分光信号Ｆは、蛍光の時間強度Ｃを蛍光色素毎に分けずにまとめて検出したものを分光した結果であるので、色変換行列Ｓを用いるとＦ＝Ｓ×Ｃとなる。つまり分光信号Ｆから蛍光の時間強度Ｃを得るには、Ｓの擬似逆行列であるＳ^＋を両辺に乗じてＣ＝Ｓ^＋×Ｆを算出すれば良いので、色変換行列Ｓを取得する処理であるスペクトル校正が必要である。なお、Ｆは波長チャンネル番号ｎの行と時間番号ｔの列の行列、Ｓは波長チャンネル番号ｎの行と蛍光色素番号ｋの列の行列、Ｃは蛍光色素番号ｋの行と時間番号ｔの列の行列である。

　色変換行列Ｓを取得するために、分光信号Ｆに対して公知の独立成分分析を行う方法はあるものの、得られる色変換行列Ｓでは精度が不十分である。ところで、各塩基には異なる蛍光色素が標識されるため、各蛍光色素の蛍光のピークが同時に出現する頻度は低い。そこで本実施例では、蛍光の時間強度Ｃにおいて各蛍光色素の蛍光のピークが同時に出現する頻度が低くなる色変換行列Ｓの候補を複数算出し、各候補に対して算出される評価値に基づいて複数の候補の中から、色変換行列Ｓを選択する。

　図３を用いて本実施例の要部について説明する。なおこれらの要部は、専用のハードウェアで構成されても良いし、ＣＰＵ２上で動作するソフトウェアで構成されても良い。以降の説明では本実施例の要部がソフトウェアで構成された場合について説明する。本実施例は、分光信号取得部３１と初期値算出部３２と候補算出部３３と選択部３４を備える。
以下、各部について説明する。

　分光信号取得部３１はスペクトル校正に用いられる分光信号３０を取得する。分光信号３０は記憶装置４に記憶される信号であっても良いし、ネットワークアダプタ１０を介してキャピラリ電気泳動型シーケンサ１３や分光信号データベース１４から取得される信号であっても良い。

　初期値算出部３２は色変換行列Ｓと蛍光の時間強度Ｃの初期値を算出する。初期値の算出には分光信号３０が用いられ、初期値の算出に先立ち、いくつかの前処理が行われても良い。前処理として、例えばスパイクの除去やベースラインの除去、波長チャンネル間の正規化等が行われることにより、以降の処理での精度を向上できる。色変換行列Ｓの初期値の算出には、公知の独立成分分析が用いられる。蛍光の時間強度Ｃの初期値は、色変換行列Ｓの初期値と分光信号Ｆを用いて算出される。

　候補算出部３３は色変換行列Ｓの候補を算出する。色変換行列Ｓの候補の算出には、蛍光の時間強度Ｃにおいて、蛍光色素の蛍光のピークが同時に出現する頻度である同時出現頻度に依存するパラメータλが用いられる。パラメータλは、同時出現頻度に対して負の相関関係にある。すなわち、パラメータλが小さいときには同時出現頻度は多くなり、パラメータλが大きくなると同時出現頻度が少なくなる。色変換行列Ｓの候補はパラメータλの値毎に算出される。

　選択部３４はパラメータλの値毎に算出された候補の中から、各候補に対して算出される評価値に基づいて色変換行列を選択する。各候補に対して算出される評価値は、例えば色変換行列の各候補を用いて算出される蛍光の時間強度Ｃの要素に基づいて算出される。蛍光の時間強度Ｃの要素が負の値となることは物理的に発生しえないので、そのような色変換行列の候補は低い評価値となるコスト関数が用いられる。

　図４を用いて、以上の各部を備えるスペクトル校正装置１が実行する処理の流れの一例を説明する。以下、各ステップについて説明する。

　（Ｓ４０１）
　分光信号取得部３１が分光信号３０を取得する。分光信号３０は記憶装置４から取得されても良いし、キャピラリ電気泳動型シーケンサ１３や分光信号データベース１４から取得されても良い。なお、分光信号３０は波長チャンネル番号ｎの行と時間番号ｔの列からなる行列であり、各要素をｘ（ｔ，ｎ）とする。また、１≦ｔ≦Ｔであり、Ｔは検出時刻の点数を表す整数であって例えばＴ＝１０００である。さらに、１≦ｎ≦Ｎであり、Ｎは波長チャンネル数を表す整数であって、例えばＮ＝２０である。

　（Ｓ４０２）
　初期値算出部３２が分光信号３０から色変換行列Ｓと蛍光の時間強度Ｃの初期値を算出する。なお色変換行列Ｓは波長チャンネル番号ｎの行と蛍光色素番号ｋの列からなる行列であり、各要素をｓ（ｎ，ｋ）とする。また、１≦ｋ≦Ｋであり、Ｋは蛍光色素の個数を表す整数であって例えばＫ＝４である。さらに蛍光の時間強度Ｃは蛍光色素番号ｋの行と時間番号ｔの列からなる行列である。

　色変換行列Ｓと蛍光の時間強度Ｃの初期値を算出する処理の流れの例を図５に示し、図５の各ステップについて以降で説明する。

　（Ｓ５０１）
　初期値算出部３２が分光信号３０からスパイク信号を除去する。スパイク信号の除去は、例えば以下の手順で実行される。
（１）まず［数１］により、ｘ（ｔ，ｎ）の時間方向の差分信号ｖ（ｔ，ｎ）を計算する。

（２）Δ＝Δ＿ｍａｘとする。なお、Δ＿ｍａｘは時間番号の最大値、Δ＿ｍｉｎは時間番号の最小値、Δ＿ｔｏｌは時間番号の許容値、θ１は分光信号ｘ（ｔ，ｎ）の閾値、θ２は差分信号ｖ（ｔ，ｎ）の閾値である。
（３）［数２］を満たすｐを探索し、そのようなｐがあれば、ｔ ∈ [ｐ－Δ＿ｔｏｌ，ｐ＋Δ＋Δ＿ｔｏｌ]の範囲のｘ(ｔ，ｎ)を０に置き換える。

（４）Δから１を減じ、Δ ＜ Δ＿ｍｉｎであれば終了し、そうでなければ（３）に戻る。

　以上の手順により、スパイク信号が除去された分光信号ｘ（ｔ，ｎ）が取得される。分光信号にスパイク信号が含まれる場合、以降の処理の精度を低下させるため、スパイク信号は除去されることが好ましい。

　（Ｓ５０２）
　初期値算出部３２が、Ｓ５０１の算出結果ｘ（ｔ，ｎ）からベースラインを除去する。ベースラインの除去は、例えば各時刻ｔの近傍区間の最小値をｘ（ｔ，ｎ）から差し引くなどの公知の方法が用いられる。ベースラインが除去された信号をｚ（ｔ，ｎ）とする。

　（Ｓ５０３）
　初期値算出部３２が、Ｓ５０２の算出結果ｚ（ｔ，ｎ）に対して、波長チャンネル間強度の正規化処理をする。正規化処理には［数３］と［数４］が用いられ、正規化信号ｆ（ｔ，ｎ）が取得される。

波長チャンネル間の強度に大きく差異がある場合、強度が大きい波長チャンネルを過剰評価し、強度が小さい波長チャンネルを過小評価することになり、以降の処理の精度を低下させる。そこで、波長チャンネル間の強度を正規化しておくことが好ましい。

　なおＳ５０１～Ｓ５０３は、以降の処理のための前処理であって、必ずしも実行されなくても良い。ただし、以降の処理の算出精度を向上させるには実行されることが好ましい。

　（Ｓ５０４）
　初期値算出部３２が、独立成分分析により色変換行列Ｓの初期値を算出する。独立成分分析の入力には、Ｓ５０３の算出結果ｆ（ｔ，ｎ）を用いることが好ましい。また独立成分分析には、Fast ICA、Auxiliary-function-based ICA、Natural-gradient-based ICAなど公知の方法が用いられる。

　なお、独立成分分析に先立ち、正規化信号ｆ（ｔ，ｎ）から時間方向の差分信号ｇ（ｔ，ｎ）＝ｆ（ｔ，ｎ）－ｆ（ｔ－１，ｎ）を算出し、差分信号ｇ（ｔ，ｎ）を独立成分分析に用いても良い。独立成分分析は入力信号が零平均の分布に従うことを仮定している。差分信号ｇ（ｔ，ｎ）は零平均の分布に従うので、独立成分分析を有効に機能させることができる。同じ理由により、正規化信号ｆ（ｔ，ｎ）から正規化信号ｆ（ｔ，ｎ）の平均値を差し引いて得られる信号を独立成分分析に用いても良い。

　また独立成分分析により得られた色変換行列Ｓの各要素ａ（ｎ，ｋ）に対して非負化処理を施して得られた各要素ｂ（ｎ，ｋ）を色変換行列Ｓの初期値としても良い。非負化処理には例えば［数５］が用いられる。

また、他の非負化処理として、ａ（ｎ，ｋ）の最小値を各ａ（ｎ，ｋ）から差し引く処理や、ａ（ｎ，ｋ）の中の負の値を零に置換する処理を施して、色変換行列Ｓの初期値を算出しても良い。算出された色変換行列Ｓ（ｎ，ｋ）の初期値は、記憶装置４等に保管される。

　（Ｓ５０５）
　初期値算出部３２が、Ｓ５０４の算出結果を用いて蛍光の時間強度Ｃの初期値を算出する。Ｃの初期値の算出には［数６］が用いられる。

なお、Ｓ^＋はＳの擬似逆行列であり、Ｆは正規化信号ｆ（ｔ，ｎ）を要素とする行列である。算出された蛍光の時間強度Ｃ（ｋ，ｔ）の初期値は、記憶装置４等に保管される。

　図４の説明に戻る。

　（Ｓ４０３）
　候補算出部３３が色変換行列の候補を算出する。色変換行列の候補は、蛍光の時間強度Ｃにおいて蛍光色素の蛍光のピークが同時に出現する頻度である同時出現頻度に依存するパラメータλの値毎に算出される。またパラメータλは例えば１．０１、１．０２、・・・、１．１０に設定される。

　色変換行列の候補を算出する処理の流れの例を図６に示し、図６の各ステップについて以降で説明する。

　（Ｓ６０１）
　候補算出部３３が、蛍光色素の蛍光のピークが同時に出現する頻度に依存するパラメータλを設定する。例えばλ＝１．０１が設定される。

　（Ｓ６０２）
　候補算出部３３が、色変換行列Ｓ（ｎ，ｋ）と蛍光の時間強度Ｃ（ｋ，ｔ）の初期値を記憶装置４等から読み出す。

　（Ｓ６０３）
　候補算出部３３が、色変換行列Ｓ（ｎ，ｋ）を波長チャネル番号ｎに関して正規化する。正規化には例えば［数７］が用いられ、正規化の結果としてＰ（ｎ｜ｋ）が算出される。

　また、［数７］の代わりに［数８］が用いられてもよい。

［数８］が用いられる場合、色変換行列Ｓ（ｎ，ｋ）の各列に相当する各蛍光色素の波長応答特性の互いの類似度をξにより調整することができる。

　（Ｓ６０４）
　候補算出部３３が、パラメータλを含む蛍光の時間強度Ｃ（ｋ，ｔ）を蛍光色素番号ｋに関して正規化する。正規化には例えば［数９］が用いられ、正規化の結果としてＰ（ｋ｜ｔ）が算出される。

１．０より大きいパラメータλが与えられ、［数９］が用いられることにより、パラメータλの大きさに応じて同時出現頻度が抑制されて、次ステップにおいて色変換行列Ｓ（ｎ，ｋ）と蛍光の時間強度Ｃ（ｋ，ｔ）が更新される。

　また、［数９］の代わりに［数１０］や［数１１］が用いられてもよい。

［数１０］や［数１１］が用いられる場合、強度が小さいピークが少ない蛍光の時間強度Ｃ（ｋ，ｔ）に更新されやすいので、雑音が大きい場合に精度の高い結果が得られる。

　（Ｓ６０５）
　候補算出部３３が、色変換行列Ｓ（ｎ，ｋ）と蛍光の時間強度Ｃ（ｋ，ｔ）を更新する。色変換行列Ｓ（ｎ，ｋ）の更新には［数１２］と［数１３］が用いられ、蛍光の時間強度Ｃ（ｋ，ｔ）の更新には［数１２］と［数１４］が用いられる。

　（Ｓ６０６）
　候補算出部３３が、規定回数に達したか否かを判定し、達した場合はＳ６０７に処理を進め、達してない場合はＳ６０３に処理を戻す。

　（Ｓ６０７）
　候補算出部３３が、Ｓ６０１で設定されたパラメータλに対する色変換行列Ｒ＿λ（ｎ，ｋ）＝σ（ｎ）Ｓ（ｎ，ｋ）を算出する。波長チャンネル間強度の正規化処理に用いられたσ（ｎ）をＳ（ｎ，ｋ）に乗じることにより、正規化処理前の信号ｚ（ｔ，ｎ）に対応する色変換行列Ｒ＿λ（ｎ，ｋ）が得られる。

　（Ｓ６０８）
　候補算出部３３が、パラメータλを更新するか否かを判定し、更新する場合はＳ６０１に処理を戻してパラメータλを再設定、例えばλ＝１．０２とし、更新しない場合はＳ４０３の処理が終了となる。パラメータλは、例えば１．０１、１．０２、・・・、１．１０の順に設定される。

　なお、Ｓ４０３において、候補算出部３３がベイズ推定を用いて色変換行列の候補を算出しても良い。すなわち［数１５］で定義される事後確率を最大化する色変換行列Ｓ（ｎ，ｋ）と蛍光の時間強度Ｃ（ｋ，ｔ）を、GibbsサンプリングやMetropolis-Hastingsサンプリングなどの公知のマルコフ連鎖モンテカルロを用いることにより求める。

ここで、αｋとβｋはガンマ分布のパラメータである。ベイズ推定を用いた場合、初期値依存性の問題を解決できるため、精度の高い色変換行列を得ることができる。

　図４の説明に戻る。

　（Ｓ４０４）
　選択部３４が、Ｓ４０３で算出された色変換行列の複数の候補Ｒ＿λ（ｎ，ｋ）の中から、各候補に対して算出される評価値に基づいて色変換行列を選択する。候補の中から色変換行列を選択する処理の流れの例を図７に示し、図７の各ステップについて以降で説明する。

　（Ｓ７０１）
　選択部３４がパラメータλ毎に色変換信号ｙ＿λ（ｔ，ｋ）を算出する。色変換信号ｙ＿λ（ｔ，ｋ）は、色変換行列の候補Ｒ＿λ（ｎ，ｋ）と正規化処理前の信号ｚ（ｔ，ｎ）を用いて、［数１６］により算出される。

なお、Ｙ_λは色変換信号ｙ＿λ（ｔ，ｋ）を要素とする行列であり、Ｒ_λ ^＋は色変換行列の候補Ｒ＿λ（ｎ，ｋ）の擬似逆行列であり、Ｚは正規化処理前の信号ｚ（ｔ，ｎ）を要素とする行列である。

　（Ｓ７０２）
　選択部３４が、Ｓ７０１で算出された色変換信号ｙ＿λ（ｔ，ｋ）の中のプルダウンを抽出する。ある蛍光色素の色変換信号ｙ＿λ（ｔ，ｋ）がある時刻に十分に大きいピークを有するときに、他の蛍光色素の色変換信号ｙ＿λ（ｔ，ｋ）が同時刻に有する負の値のピークのことを「プルダウン」と呼ぶ。プルダウンは物理的には発生しえず、色変換信号ｙ＿λ（ｔ，ｋ）の算出に用いられた色変換行列の候補Ｒ＿λ（ｎ，ｋ）に含まれる誤差に起因するものである。そこで、本ステップでは色変換信号ｙ＿λ（ｔ，ｋ）に含まれるプルダウンを抽出し、抽出されたプルダウンに基づいて次ステップ以降で色変換行列の各候補Ｒ＿λ（ｎ，ｋ）を評価する。なお、色変換行列の候補Ｒ＿λ（ｎ，ｋ）に誤差が含まれる場合、プルダウンととともに、正の値のピークであるプルアップも生じるが、本実施例では扱わない。

　色変換信号ｙ＿λ（ｔ，ｋ）の中のプルダウンは、各蛍光色素の色変換信号ｙ＿λ（ｔ，ｋ）が負の値であって、ｙ＿λ（ｔ，ｋ）の時間方向の差分値ｙ＿λ（ｔ，ｋ）－ｙ＿λ（ｔ－１，ｋ）が負から正に変わる時刻を探索することにより抽出される。

　（Ｓ７０３）
　選択部３４が、Ｓ７０２で抽出されたプルダウンに基づきコスト関数ｃ＿λを算出する。コスト関数ｃ＿λの算出には例えば［数１７］が用いられる。

ここで、ｍは抽出されたプルダウンの番号であり、ｄｍは各プルダウンの強度の絶対値、ｅｍはプルダウンと同時刻における他の蛍光色素の最大強度、Ｄは予め定められたプルダウンの許容値である。［数１７］によれば、最大強度の総和による除算演算とすることにより、零除算による演算エラーを回避できる。

　なお、プルダウンに基づくコスト関数ｃ＿λは［数１７］に限定されず、プルダウンの面積の総和を評価可能な関数であっても良い。

　（Ｓ７０４）
　選択部３４が、Ｓ７０３で算出されたコスト関数ｃ＿λにより候補の中から色変換行列Ｒ＿λ（ｎ，ｋ）を選択する。コスト関数ｃ＿λが［数１７］の場合には、コスト関数ｃ＿λが最小となる色変換行列Ｒ＿λ（ｎ，ｋ）が選択される。

　以上説明した処理の流れにより、蛍光色素間でピークが同時に出現する場合であっても適切な色変換行列Ｒ＿λ（ｎ，ｋ）が求められる。すなわち、分光信号の時間分解能が低い場合や異なる蛍光色素間で蛍光のピークが重複する頻度が高い場合であってもスペクトル校正を高精度に行える。特にキャピラリが短い小型のキャピラリ電気泳動型シーケンサや測定時間が短い場合に本実施例は有効である。

　実施例１では色変換信号ｙ＿λ（ｔ，ｋ）の中のプルダウンに基づくコスト関数ｃ＿λを用いて色変換行列Ｒ＿λ（ｎ，ｋ）を選択することについて説明した。本実施例では色変換信号ｙ＿λ（ｔ，ｋ）に逆畳み込み処理を施して得られる逆畳み込み信号に基づくコスト関数ｃ＿λを用いることについて説明する。なお、本実施例の全体構成は実施例１と同じであるので、説明を省略する。

　図８を用いて本実施例の候補の選択処理について説明する。なお図８の処理の流れは、図４のＳ４０３において実行される。

　（Ｓ８０１）
　選択部３４が、Ｓ７０１と同様に、パラメータλ毎に色変換信号ｙ＿λ（ｔ，ｋ）を算出する。

　（Ｓ８０２）
　選択部３４が、Ｓ８０１で算出された色変換信号ｙ＿λ（ｔ，ｋ）に対して逆畳み込み処理をする。逆畳み込み処理には、公知のブラインドデコンボリューション処理等が用いられる。逆畳み込み処理により、パラメータλ毎に逆畳み込み信号が算出される。

　（Ｓ８０３）
　選択部３４が、Ｓ８０２で算出された逆畳み込み信号に基づいて、コスト関数ｃ＿λを算出する。

　（Ｓ８０４）
　選択部３４が、Ｓ８０３で算出されたコスト関数ｃ＿λにより、候補の中から色変換行列Ｒ＿λ（ｎ，ｋ）を選択する。

　以上説明した処理の流れにより、蛍光色素間でピークが同時に出現する場合であっても適切な色変換行列Ｒ＿λ（ｎ，ｋ）が求められる。なお、逆畳み込み信号では色変換信号に比べて分離度が高くなるため、より容易にプルダウンを検知することができる。

　実施例１では色変換信号の中のプルダウンに基づくコスト関数を用いて色変換行列を選択することについて説明した。本実施例では色変換信号に移動度補正処理を施して得られる移動度補正処理信号に基づくコスト関数を用いることについて説明する。なお、本実施例の全体構成は実施例１と同じであるので、説明を省略する。

　図９を用いて本実施例の候補の選択処理について説明する。なお図９の処理の流れは、図４のＳ４０３において実行される。

　（Ｓ９０１）
　選択部３４が、Ｓ７０１と同様に、パラメータλ毎に色変換信号ｙ＿λ（ｔ，ｋ）を算出する。

　（Ｓ９０２）
　選択部３４が、Ｓ９０１で算出された色変換信号ｙ＿λ（ｔ，ｋ）に対して移動度補正処理をする。移動度補正処理とは、蛍光色素の蛍光のピークが等間隔となるように、各蛍光色素の蛍光スペクトルを時間方向に移動させる処理である。四つの蛍光色素Ａ、Ｇ、Ｃ、Ｔの蛍光スペクトルに対する移動度補正処理の一例を図１０に示す。図１０（ａ）は移動度補正処理前の蛍光スペクトルの例であり、蛍光のピークが不等間隔である。図１０（ｂ）は移動度補正処理後の蛍光スペクトルの例であり、蛍光色素Ａと蛍光色素Ｃを時間の正方向に移動させ、蛍光色素Ｇをそのままとし、蛍光色素Ｔを時間の負方向に移動させることにより、蛍光のピークが等間隔になる。

　ところで色変換行列の候補Ｒ＿λ（ｎ，ｋ）に誤差が含まれる場合、色変換信号ｙ＿λ（ｔ，ｋ）には図１１（ａ）に示すようなプルアップ１１０１やプルダウン１１０２が含まれる。プルアップ１１０１やプルダウン１１０２が含まれる色変換信号ｙ＿λ（ｔ，ｋ）に対して移動度補正処理を施すと、図１１（ｂ）に示すようにピーク間隔が乱れる。本実施例ではピーク間隔の乱れを利用するために、本ステップにおいて移動度補正処理をする。移動度補正処理には公知の方法が用いられ、移動度補正処理よりパラメータλ毎に移動度補正信号ｈ＿λ（ｔ，ｋ）が算出される。

　（Ｓ９０３）
　選択部３４が、Ｓ９０２で算出された移動度補正信号ｈ＿λ（ｔ，ｋ）に基づくコスト関数ｃ＿λを算出する。コスト関数ｃ＿λの算出手順の一例について説明する。
（１）ｈ＿λ（ｔ，ｋ）をｋ方向に積算して時系列信号ｕ＿λ（ｔ）を得る。すなわち移動度補正信号ｈ＿λ（ｔ，ｋ）を蛍光色素で区別せずに扱う。
（２）ｕ＿λ（ｔ）に対し離散フーリエ変換を行い、Ｕ＿λ（ω）を得る。ただし、ωは離散化された周波数インデックスである。
（３）パワースペクトルＶ＿λ（ω）＝｜Ｕ＿λ（ω）｜^２を計算する。
（４）規定の範囲[ω_{ｂｅｇｉｎ}， ω_ｅｎｄ]内のωに対して、Ｖ＿λ（ω）の最大値Ｗ＿λを求める。Ｗ＿λは、乱れのないときのピーク間隔に相当する周波数成分の大きさを表す。
（５）規定の範囲[ω_{ｂｅｇｉｎ}， ω_ｅｎｄ]において、Ｖ＿λ（ω）のω方向にわたる積算値Ｅ＿λ＝Σ（Ｖ＿λ（ω））を計算する。Ｅ＿λは移動度補正信号ｈ＿λ（ｔ，ｋ）のピーク間隔の全周波数成分を含むパワーであり、ピーク間隔が乱れるほど大きくなる。
（６）コスト関数としてｃ＿λ＝Ｅ＿λ／Ｗ＿λを計算する。

　（Ｓ９０４）
　選択部３４が、Ｓ９０３で算出されたコスト関数ｃ＿λにより候補の中から色変換行列Ｒ＿λ（ｎ，ｋ）を選択する。本実施例のコスト関数ｃ＿λは、乱れのない時のピーク間隔に相当する周波数成分の大きさで、全周波数成分を含むパワーを除した値であるので、ピーク間隔が乱れるほど大きくなる。よって、コスト関数ｃ＿λが最小となる色変換行列Ｒ＿λ（ｎ，ｋ）が候補の中から選択される。

　以上説明した処理の流れにより、蛍光色素間でピークが同時に出現する場合であっても適切な色変換行列Ｒ＿λ（ｎ，ｋ）が求められる。なお、本実施例のコスト関数ｃ＿λではプルダウンやプルアップの検出が不要であるので、ノイズが多くプルダウンやプルアップの検出が困難な場合であっても適切な色変換行列Ｒ＿λ（ｎ，ｋ）を選択することができる。　
　以上、本発明の複数の実施例について説明した。本発明はこれらの実施例に限定されるものではなく、さまざまな変形例が含まれる。例えば、上記した実施例は本発明をわかりやすく説明するために詳細に説明したのであり、説明したすべての構成を備えるものに限定されるものではない。また、ある実施例の構成の一部を他の実施例の構成に置き換えることが可能である。さらに、ある実施例の構成に他の実施例の構成を加えることも可能である。また、各実施例の構成の一部について、他の構成の追加・削除・置換をすることが可能である。

１：スペクトル校正装置、２：ＣＰＵ、３：主メモリ、４：記憶装置、５：表示メモリ、６：表示装置、７：コントローラ、８：入力装置、１０：ネットワークアダプタ、１１：システムバス、１２：ネットワーク、１３：キャピラリ電気泳動型シーケンサ、１４：分光信号データベース、３０：分光信号、３１：分光信号取得部、３２：初期値算出部、３３：候補算出部、３４：選択部、３５：色変換行列、１１０１：プルアップ、１１０２：プルダウン

Claims

　色変換処理に用いられる色変換行列を算出するスペクトル校正装置であって、
　経時的に検出された蛍光の分光信号を取得する分光信号取得部と、
　前記分光信号に基づいて、蛍光色素の蛍光のピークが同時に出現する頻度に依存するパラメータの値毎に色変換行列の候補を算出する候補算出部と、
　各候補に対して算出される評価値に基づいて色変換行列を選択する選択部を備えることを特徴とするスペクトル校正装置。
　請求項１に記載のスペクトル校正装置であって、
　前記選択部は、前記パラメータの値毎に算出された色変換行列の各候補を用いて算出される色変換信号の中からプルダウンを抽出し、前記プルダウンに基づくコスト関数を用いて前記評価値を算出することを特徴とするスペクトル校正装置。
　請求項１に記載のスペクトル校正装置であって、
　前記選択部は、前記パラメータの値毎に算出された色変換行列の各候補を用いて算出される色変換信号に逆畳み込み処理を実行して逆畳み込み信号を取得し、前記逆畳み込み信号に基づくコスト関数を用いて前記評価値を算出することを特徴とするスペクトル校正装置。
　請求項１に記載のスペクトル校正装置であって、
　前記選択部は、前記パラメータの値毎に算出された色変換行列の各候補を用いて算出される色変換信号に移動度補正処理を実行して移動度補正信号を取得し、前記移動度補正信号に基づくコスト関数を用いて前記評価値を算出することを特徴とするスペクトル校正装置。
　色変換処理に用いられる色変換行列を算出するスペクトル校正方法であって、
　経時的に検出された蛍光の分光信号を取得する取得ステップと、
　前記分光信号に基づいて、蛍光色素の蛍光のピークが同時に出現する頻度に依存するパラメータの値毎に色変換行列の候補を算出する候補算出ステップと、
　各候補に対して算出される評価値に基づいて色変換行列を選択する選択ステップを備えることを特徴とするスペクトル校正方法。
　請求項５に記載のスペクトル校正方法であって、
　前記選択ステップは、前記パラメータの値毎に算出された色変換行列の各候補を用いて算出される色変換信号の中からプルダウンを抽出し、前記プルダウンに基づくコスト関数を用いて前記評価値を算出することを特徴とするスペクトル校正方法。
　請求項５に記載のスペクトル校正方法であって、
　前記選択ステップは、前記パラメータの値毎に算出された色変換行列の各候補を用いて算出される色変換信号に逆畳み込み処理を実行して逆畳み込み信号を取得し、前記逆畳み込み信号に基づくコスト関数を用いて前記評価値を算出することを特徴とするスペクトル校正方法。
　請求項５に記載のスペクトル校正方法であって、
　前記選択ステップは、前記パラメータの値毎に算出された色変換行列の各候補を用いて算出される色変換信号に移動度補正処理を実行して移動度補正信号を取得し、前記移動度補正信号に基づくコスト関数を用いて前記評価値を算出することを特徴とするスペクトル校正方法。