JP6079511B2 - ピーク抽出方法及びプログラム - Google Patents

Description

本発明は、クロマトグラムやスペクトルからピークを抽出するピーク抽出方法及びプログラムに関する。

液体クロマトグラフ（ＬＣ）やガスクロマトグラフ（ＧＣ）では、試料をカラムに導入して該試料に含まれる成分を時間的に分離し、それらを順に検出することによって、横軸を保持時間、縦軸を信号強度とするクロマトグラムを取得する。そして、クロマトグラムからピークを抽出し、ピークトップの位置に対応する保持時間から成分を同定し、また、ピーク強度やピーク面積から該成分を定量する。

特許文献１に記載されているように、ＧＣにおいて昇温分析を行ったり、ＬＣにおいてグラジエント分析を行ったりしてクロマトグラムを取得すると、ベースラインが時間的に変動する、いわゆるベースラインドリフトが生じることが多い。そこで、クロマトグラムを２次微分することによってベースラインドリフトの影響を排除した後で、２次微分波形からピークを抽出する方法が提案されている（非特許文献１）。

クロマトグラムを２次微分すると、直線的に増加あるいは減少するベースラインドリフトの影響を排除できる。仮に、ベースラインドリフトが直線的に増加あるいは減少するものでなくても、ベースラインドリフトは一定の傾向で緩やかに変動するため、２次微分を行うことでその影響を低減できる。

クロマトグラムを２次微分した後は、クロマトグラムを２次微分して得た波形において、クロマトグラムに重畳するノイズの大きさを考慮して予め設定した閾値以上の絶対値を持つ極値を抽出し、その位置をピークとして検出する。

ここでは、クロマトグラムを一例に挙げて説明したが、他の測定により得られる波形でも同様にベースラインのドリフトが生じる場合がある。例えば、特許文献２に記載のように、分光測定において、試料を配置しないで取得したリファレンスデータと、試料を配置して取得した測定データとの差分を計算して得たスペクトルでは、リファレンスデータ取得時と測定データ取得時との測定環境の違い（光学系の変化、温度変化等）に起因するベースラインのドリフトが生じる場合がある。従って、クロマトグラムの場合と同様に、スペクトルのピークを検出する際にも２次微分によりベースラインドリフトの影響を排除している。

特開２００４−２７１４２２号公報 特開２００１−３４３３２４号公報

Thermo (DIONEX) Software Postres, "DIONEX Taking the Pain Out of Chromatographic Peak Integration", [online], [平成25年8月12日検索],インターネット<URL:http://www.dionex.com/en-us/webdocs/77494-PO-HPLC-LPN2297-01-Chrome.pdf>

クロマトグラムやスペクトル（以下、クロマトグラム等とする。）を２次微分すると、広範囲にわたる緩やかな信号強度の変化を抑制しつつ、局所的な信号強度の変化を強調した波形が得られるため、ベースラインドリフトの影響を軽減してピークを抽出することができる。ところが、クロマトグラムに重畳するノイズも局所的な信号強度の変化であるため、２次微分後の波形では真のピークだけでなくノイズも強調されてしまう。特に、急峻な波形のノイズが重畳していると、２次微分後の波形ではノイズが大きく強調されてしまい、低Ｓ／Ｎ比のピークの検出が困難になってしまう。

本発明が解決しようとする課題は、クロマトグラフィや分光分析により取得した測定波形からピークを抽出する際に、ベースラインドリフトの影響を低減しつつ、ノイズの中から低Ｓ／Ｎ比のピークでも正しく抽出することのできる方法、及び該方法をコンピュータに実行させるためのプログラムを提供することである。

上記課題を解決するために成された本発明は、クロマトグラフィや分光分析により取得した測定波形から真のピークを抽出するピーク抽出方法であって、
a) 該測定波形を２次微分して２次微分波形を取得する２次微分波形取得ステップと、
b) 前記２次微分波形の極大値及び／又は極小値に基づいて前記測定波形上のピークを抽出して仮ピークとする仮ピーク抽出ステップと、
c) 前記真のピークに対して予め設定されたモデルピーク関数の理論式に基づいて前記仮ピークのピーク幅を決定するピーク幅決定ステップと、
d) 前記モデルピーク関数の理論式に基づいて、前記ピーク幅の両端に対応する前記測定波形上又は前記２次微分波形上の２点を用いて前記仮ピークの高さの理論値を計算するピーク高さ理論値計算ステップと、
e) 前記測定波形上のノイズのばらつきの指標値を前記２次微分波形から計算するノイズ指標値計算ステップと、
f) 前記仮ピーク高さの理論値と前記指標値の比であるＳ／Ｎ比を計算し、該Ｓ／Ｎ比の値が予め設定された値以上である仮ピークを真のピークとして抽出するピーク抽出ステップと、
を有することを特徴とする。

本発明に係るピーク抽出方法により抽出されるピークは１つに限定されない。複数の真のピークを抽出する場合には、仮ピーク抽出ステップにおいて複数の仮ピークを抽出し、各仮ピークについて、ピーク幅の決定、ピーク高さの理論値の計算等を行う。

モデルピーク関数は、典型的にはガウス関数であるが、ローレンツ関数等としてもよい。また、ピークトップを挟んで対称な関数に限らず、測定系の状態を考慮し、例えばクロマトグラムのピークがテーリングあるいはリーディングしている場合には、Exponentially Modified Gaussian(EMG)や、BI-Gaussian等の非対称関数を用いてもよい。

前記ピーク高さ理論値計算ステップでは、例えば、後述するように、前記２点と仮ピークのピークトップの高さの実測値を用いて仮ピークの高さの理論値を計算することができる。この場合には、前記仮ピーク抽出ステップにおいて、仮ピークのピークトップの位置を特定し、該ピークトップの高さを決定しておく。
また、前記２次微分波形上でピーク幅の範囲内でモデルピーク関数を２次微分した関数をフィッティングすることによりモデルピーク関数の係数を設定し、それを元にモデルピークの高さの理論値を計算するようにしてもよい。

モデルピーク関数がガウス関数である場合のピーク抽出手順の一例を説明する。
まず、測定波形を２次微分して得た２次微分波形の極大値／極小値に基づいて仮ピークのピークトップの位置と高さH（実測値）を特定する。続いて、そのピークトップの両側に位置する変曲点の位置をピーク幅の両端に対応する２点とする。そして、それら２点の高さとピークトップの高さの差(H-h)を計算する。ガウス関数の理論式によれば、変曲点の高さはピーク高さのexp(-0.5)倍であるため、ピーク高さの理論値H₀を以下の式(1)により計算することができる。
H₀=(H-h)/(1-exp(-0.5)) …(1)

前記ノイズのばらつきの指標値は、典型的には２次微分波形の標準偏差であるが、必ずしもこれに限らず、２次微分波形を平均化して算出した値などとしてもよい。
ノイズの標準偏差を計算する方法としては、従来用いられている種々の方法を用いることができる。例えば、２次微分波形のMedian Absolute Deviation (MAD)の1.4826倍を標準偏差とする方法や、２次微分波形の絶対値のMedianの1.4826倍をノイズの標準偏差とする方法等を用いることができる。

本発明に係るピーク抽出方法では、まず、測定波形の２次微分波形の極大値及び／又は極小値に基づいて仮ピークを抽出する。続いて、モデルピーク関数の理論式に基づいて各仮ピークのピーク高さの理論値を計算する。そして、ピーク高さの理論値に基づいてＳ／Ｎ比を計算し、その結果から真のピークを抽出する。これにより、ブロードで低いピークについても、そのピーク高さを正確に計算して、真のピークとして抽出することができる。従って、本発明に係るピーク抽出方法を用いると、ベースラインのドリフトの影響を低減しつつ、ノイズの中から、低Ｓ／Ｎ比のピークでも正しく抽出することができる。

また、上述の一例のように、仮ピークのピーク高さの理論値を計算する際に、ピーク幅に対応する２点の高さと、仮ピークのピーク高さの実測値を用いると、モデルピーク関数に従う形状のピークに対しては正確なピーク高さが計算され、これにより、真のピークを強調して低Ｓ／Ｎ比のピークをより正確に抽出することができる。

上記課題を解決するために成された本発明の別の態様は、クロマトグラフィや分光分析により取得した測定波形から真のピークを抽出する処理をコンピュータに実行させるピーク抽出プログラムであって、前記真のピークに対して設定されるモデルピーク関数の理論式を含むモデルピーク関数情報及び後記Ｓ／Ｎ比に関する閾値の情報が保存された記憶部にアクセス可能なコンピュータを、
a) 該測定波形を２次微分して２次微分波形を取得する２次微分波形取得部、
b) 前記２次微分波形の極大値及び／又は極小値に基づいて前記測定波形上のピークを抽出して仮ピークに設定する仮ピーク抽出部、
c) 前記モデルピーク関数の理論式に基づいて前記仮ピークのピーク幅を決定するピーク幅決定部、
d) 前記モデルピーク関数の理論式に基づいて、前記ピーク幅の両端に対応する前記測定波形上又は前記２次微分波形上の２点を用いて前記仮ピークの高さの理論値を計算するピーク高さ理論値計算部、
e) 前記測定波形上のノイズのばらつきの指標値を前記２次微分波形から計算するノイズ指標値計算部、
f) 前記仮ピーク高さの理論値と前記指標値の比であるＳ／Ｎ比を計算し、該Ｓ／Ｎ比の値が前記閾値以上である仮ピークを真のピークとして抽出するピーク抽出部、
として動作させることを特徴とする。

本発明に係るピーク抽出方法では、測定波形の２次微分波形の極大値及び／又は極小値に基づいて仮ピークを抽出し、そのピーク高さの理論値を計算して正確なピーク高さを得る。そして、該理論値に基づいてＳ／Ｎ比を計算して真のピークを抽出する。従って、ベースラインのドリフトの影響を低減しつつ、ノイズの中から、低Ｓ／Ｎ比のピークでも正しく抽出することができる。本発明に係るピーク抽出プログラムにおいても同様の効果を得ることができる。

本発明に係るピーク抽出方法の一実施例を実行するピーク抽出装置の概略構成図。 本実施例におけるピーク抽出方法のフローチャート。 本実施例における、(a)仮ピークを２次微分した波形の例、及び(b)クロマトグラム上の仮ピーク波形の例。

本発明に係るピーク抽出方法及びピーク抽出プログラムの実施例について、以下、図面を参照して説明する。図１は本実施例のピーク抽出方法を実行するピーク抽出装置１の概略構成図、図２は本実施例のピーク抽出方法のフローチャート、図３は本実施例における、(a)仮ピークを２次微分した波形の例、及び(b)クロマトグラム上の仮ピーク波形の例である。

ピーク抽出装置の実体は汎用のコンピュータ（ＰＣ）であり、中央演算装置（ＣＰＵ）１１、メモリ１２、ハードディスクドライブ（ＨＤＤ）やソリッドステートドライブ（ＳＳＤ）等から構成される大容量記憶装置（記憶部）１５、外部記憶媒体読取部１６、等を備えている。また、ピーク抽出装置１には、液晶ディスプレイ等から構成され、各種の情報を表示する表示部１３と、マウスやキーボードから構成され、ユーザが各種の命令を入力する入力部１４とが接続されている。

記憶部１５にはピーク抽出プログラム１７が保存されており、該プログラム１７の起動により、ＣＰＵ１１が、後述する２次微分波形取得部１７１、仮ピーク抽出部１７２、ピーク幅決定部１７３、ピーク高さ理論値計算部１７４、ノイズ指標値計算部１７５、及びピーク抽出部１７６として動作する。本実施例では、記憶部１５にピーク抽出プログラム１７を保存して実行するが、ピーク抽出プログラムを記憶したＤＶＤ−ＲＯＭ等の記憶媒体２３を、記録媒体読取部１６に挿入して該プログラムを実行するようにしても良い。

また、記憶部１５には、後述するモデルピーク関数情報１８及び閾値情報１９が保存されている。本実施例では、クロマトグラムの真のピークがガウス関数で表されるピーク形状を有していると想定されており、記憶部１５にはモデルピーク関数情報１８として、ガウス関数の理論式、ガウス関数に従う形状のモデルピークにおけるピーク高さとピーク幅の両端位置の高さの比に関する情報が保存されている。なお、複数のモデルピーク関数に関するモデルピーク関数情報１８を保存しておき、入力部１４を通じた操作によってユーザに使用するモデルピーク関数を選択させるようにしてもよい。

本実施例のピーク抽出装置１によりクロマトグラムから真のピークを抽出する手順を説明する。クロマトグラムは、図示しないクロマトグラフ装置において取得され、記憶部１５に保存されている。

はじめに、２次微分波形取得部１７１は、記憶部１５からクロマトグラムを読み出し、その波形に対して、平滑化微分フィルタの一種であるSavitzky-Goleyフィルタを用いた２次微分を行うことにより２次微分波形を取得する（ステップＳ１）。そして、クロマトグラムと２次微分波形を表示部１３の画面上に表示する。２次微分の方法としては、その他、ガウシアン２次微分フィルタ等の方法を用いることもできる。いずれの場合においてもフィルタ係数の二乗和が１となるような正規化を行っておく。

次に、仮ピーク抽出部１７２は、２次微分波形の極大値と極小値を検出してピークの位置を特定する。そして、クロマトグラム上で特定した位置に存在する複数のピークをそれぞれ仮ピークに設定する（ステップＳ２）。また、クロマトグラム上での該複数の仮ピークのそれぞれのピーク高さH、即ち各仮ピークの高さの実測値を決定する（ステップＳ３）。ステップＳ２及びＳ３は、全ての仮ピークに対して行う。

続いて、ピーク幅決定部１７３は、複数の仮ピークのうちの１つについて、２次微分波形における０クロス点（２次微分波形が横軸と交差する点）を検出し、ピークトップを挟んで両側に位置する０クロス点をピーク幅の両端に対応する位置に決定する（図３(a)参照）。そして、２次微分波形における２つの０クロス点に対応する、クロマトグラム上の２点をピーク幅の両端に対応する位置に設定する。即ち、２次微分波形に基づき、クロマトグラム上でピーク幅を決定する（ステップＳ４）。

なお、本実施例では、０クロス点に対応する位置からピーク幅（±１σ）を決定したが、ピーク幅の両端として設定する点はモデルピーク関数の性質に応じて適宜に変更する。モデルピーク関数がガウス関数である場合には、０クロス点のほか、２次微分波形においてピークトップの位置に対応する極小値（あるいは極大値）の両端に位置する極大値（あるいは極小値）の位置に基づいてピーク幅（±√３σ）に対応する位置を決定してもよい。

クロマトグラムのピーク幅が決定すると、ピーク高さ理論値計算部１７４は、クロマトグラム上における該ピーク幅の両端に対応する位置の高さhを検出する。そして、上述のピーク高さの実測値Hとの差H-hを計算する（ステップＳ５）。ガウス関数に従うモデルピークでは、ピーク幅の両端に対応する位置の高さは、モデルピークのピーク高さのexp(-0.5)倍となる。そこで、次式(1)により、仮ピークのピーク高さの理論値H₀を計算する（ステップＳ６）。
H₀=(H-h)/(1-exp(-0.5)) …(1)

次に、ノイズ指標値計算部１７５が、測定波形上のノイズのばらつきの指標値であるノイズ指標値を計算する（ステップＳ７）。本実施例ではノイズのばらつきの指標値を２次微分波形の標準偏差により計算する。ノイズのばらつきの指標値は、これに限らず、２次微分波形を平均化して算出した値など、適宜に求めることができる。

本実施例では、２次微分波形のMedian Absolute Deviation (MAD)を計算し、その値の1.4826倍を標準偏差、即ちノイズのばらつきの指標値とする。ノイズのばらつきの指標値の計算についても、その他、２次微分波形の絶対値のMedianの1.4826倍をノイズの標準偏差とする方法など、適宜の方法を用いることができる。

本実施例では、ノイズ指標値計算部１７５によるノイズのばらつきの指標値の計算（ステップＳ７）を、ピーク高さ理論値計算部１７４による仮ピークのピーク高さの理論値の計算（ステップＳ６）の後に行っている。しかし、ノイズのばらつきの指標値の計算は、仮ピークに対する計算等を行うステップＳ２〜Ｓ５と独立に実行することができる。従って、ステップＳ７はステップＳ１において２次微分が取得された後であって、後述のステップＳ８において仮ピークのＳ／Ｎ比の値を計算するまでの任意の時点（ステップＳ２〜Ｓ６と並行して行うことも含む。）で実行することができる。

仮ピークのピーク高さの理論値と、ノイズのばらつきの指標値の計算が完了すると、ピーク抽出部１７６は、該理論値と該指標値の比を計算して仮ピークのＳ／Ｎ比を求める（ステップＳ８）。そして、記憶部１５に保存されている閾値と該Ｓ／Ｎ比の値を比較し、Ｓ／Ｎ比の値が閾値以上である仮ピークを真のピークとして抽出する（ステップＳ９）。

ピーク抽出部１７６が、仮ピークのＳ／Ｎ比の値に基づくピークの抽出を完了すると、ステップＳ２において設定した全ての仮ピークに対する処理が完了したか否かを判定し、未だ処理が行われていない仮ピークが存在する場合には（ステップＳ１０でＮＯ）、ステップＳ４に戻って、未処理の仮ピークに対する処理を実行する。ただし、２回目以降の処理である場合には、ノイズのばらつきの指標値の計算（ステップＳ７）をスキップする。そして、全ての仮ピークに対する処理が完了すると（ステップＳ１０でＹＥＳ）、表示部１３の画面上で、クロマトグラムに真のピークを重畳表示させるとともに、真のピークのプロファイル（溶出時間、強度、ピーク面積等）を表示して動作を完了する。

以上、説明したように、本実施例のピーク抽出方法では、まず、測定波形の２次微分波形の極大値及び／又は極小値に基づいて仮ピークを抽出する。続いて、モデルピーク関数の理論式に基づいて各仮ピークのピーク高さの理論値を計算する。そして、ピーク高さの理論値に基づいてＳ／Ｎ比を計算し、その結果から真のピークを抽出する。これにより、ブロードで低いピークを含む全ての真のピークについて、そのピーク高さを正確に計算することができる。そのため、低Ｓ／Ｎ比のピークについても、真のピークとして正しく抽出することができる。

上記実施例は一例であって、本発明の趣旨に沿って適宜に変更することができる。
上記実施例では、ピーク抽出装置と、クロマトグラフ装置を別体としたが、ピーク抽出プログラムをクロマトグラフ装置にインストールして実行することもできる。
また、入力部１４を通じたユーザの操作に基づいて閾値情報１９を設定あるいは変更するようにしてもよい。
上記実施例では、複数の仮ピークから真のピークを抽出する場合を例に挙げて説明したが、その他、例えば、測定波形上にある１つのピークが真のピークであるか否かを判定する場合にも本実施例の方法を用いることができる。この場合には、仮ピークの設定（ステップＳ２）や全仮ピークに対する処理が完了した否かの判定（ステップＳ１０）を行う必要はない。
上記実施例では、仮ピークのピーク幅に対応する位置の高さhと、仮ピークのピーク高さHの実測値の差から、仮ピークのピーク高さの理論値を計算したが、クロマトグラムあるいはその２次微分波形上でピーク幅に対応する２点を決定し、その範囲内でモデルピーク関数をフィッティングすることによってピーク高さの理論値を求めるようにしてもよい。具体的には、クロマトグラム上で仮ピークのピーク幅に対応する２点を用い、該２点を通るようにモデルピーク関数の係数を設定したり、あるいは、ピーク幅の両端に対応する２次微分波形上の２点を用い、２次微分波形上で該２点を通るようにモデルピーク関数の２次微分関数の係数を設定することによって、ピーク高さの理論値を求めることができる。

上記実施例では、クロマトグラムから真のピークを抽出する場合を説明したが、分光測定により取得したスペクトル等の波形に対しても上記同様の方法及びプログラムを用いてピークを抽出することができる。

１…ピーク抽出装置
１１…ＣＰＵ
１２…メモリ
１３…表示部
１４…入力部
１５…記憶部
１６…外部記憶媒体読取部
１７…ピーク抽出プログラム
１７１…２次微分波形取得部
１７２…仮ピーク抽出部
１７３…ピーク幅決定部
１７４…ピーク高さ理論値計算部
１７５…ノイズ指標値計算部
１７６…ピーク抽出部
１８…モデルピーク関数情報
１９…閾値情報

Claims (6)

1. クロマトグラフィや分光分析により取得した測定波形から真のピークを抽出するピーク抽出方法であって、
   a) 該測定波形を２次微分して２次微分波形を取得する２次微分波形取得ステップと、
   b) 前記２次微分波形の極大値及び／又は極小値に基づいて前記測定波形上のピークを抽出して仮ピークとする仮ピーク抽出ステップと、
   c) 前記真のピークに対して予め設定されたモデルピーク関数の理論式に基づいて前記仮ピークのピーク幅を決定するピーク幅決定ステップと、
   d) 前記モデルピーク関数の理論式に基づいて、前記ピーク幅の両端に対応する前記測定波形上又は前記２次微分波形上の２点を用いて前記仮ピークの高さの理論値を計算するピーク高さ理論値計算ステップと、
   e) 前記測定波形上のノイズのばらつきの指標値を前記２次微分波形から計算するノイズ指標値計算ステップと、
   f) 前記仮ピーク高さの理論値と前記指標値の比であるＳ／Ｎ比を計算し、該Ｓ／Ｎ比の値が予め設定された値以上である仮ピークを真のピークとして抽出するピーク抽出ステップと、
   を有することを特徴とするピーク抽出方法。
2. 前記ピーク幅決定ステップにおいて、前記２次微分波形において極小値又は極大値が現れる位置を挟んで両側に位置する２つの０クロス点、あるいは、前記２次微分波形において極小値又は極大値が現れる位置を挟んで両側に存在する極大値又は極小値に対応する２点に基づいて前記仮ピークのピーク幅を決定することを特徴とする、請求項１に記載のピーク抽出方法。
3. 前記ピーク高さ理論値計算ステップにおいて、前記測定波形上で、前記仮ピークのピーク幅に対応する位置の高さと、前記仮ピークのピーク高さの実測値の差を用いて前記仮ピークのピーク高さの理論値を計算することを特徴とする請求項１又は２に記載のピーク抽出方法。
4. 前記ピーク高さ理論値計算ステップにおいて、前記２次微分波形上で、前記仮ピークのピーク幅に対応する範囲内でモデルピーク関数をフィッティングすることによって前記仮ピークのピーク高さの理論値を計算することを特徴とする請求項１又は２に記載のピーク抽出方法。
5. 前記ノイズ指標値計算ステップにおいて、前記２次微分波形のMedian Absolute Deviation又は前記２次微分波形の絶対値のMedianを計算することにより、ノイズのばらつきの指標値を計算することを特徴とする請求項１から３のいずれかに記載のピーク抽出方法。
6. クロマトグラフィや分光分析により取得した測定波形から真のピークを抽出する処理をコンピュータに実行させるピーク抽出プログラムであって、前記真のピークに対して設定されるモデルピーク関数の理論式を含むモデルピーク関数情報及び後記Ｓ／Ｎ比に関する閾値の情報が保存された記憶部にアクセス可能なコンピュータを、
   a) 該測定波形を２次微分して２次微分波形を取得する２次微分波形取得部、
   b) 前記２次微分波形の極大値及び／又は極小値に基づいて前記測定波形上のピークを抽出して仮ピークに設定する仮ピーク抽出部、
   c) 前記モデルピーク関数の理論式に基づいて前記仮ピークのピーク幅を決定するピーク幅決定部、
   d) 前記モデルピーク関数の理論式に基づいて、前記ピーク幅の両端に対応する前記測定波形上又は前記２次微分波形上の２点を用いて前記仮ピークの高さの理論値を計算するピーク高さ理論値計算部、
   e) 前記測定波形上のノイズのばらつきの指標値を前記２次微分波形から計算するノイズ指標値計算部、
   f) 前記仮ピーク高さの理論値と前記指標値の比であるＳ／Ｎ比を計算し、該Ｓ／Ｎ比の値が前記閾値以上である仮ピークを真のピークとして抽出するピーク抽出部、
   として動作させることを特徴とするピーク抽出プログラム。

Images