peak deconvolution の覚え書き②～Savitzky-Golay filter ごっこ～
はじめに
本稿では HPLC のクロマトデータ処理として有名な Savitzky-Golay 法について勉強したことをまとめます．本稿は高校の数学定期試験で零点を連発していた人間が書いています．そのつもりでお読みください．間違っているところを見つけたら優しく教えてくださると幸いです．
前回の続き
前回 5 点を通る三次関数の算出法を求め，Excel で計算するところまで実施した (前回)．今回はそれを踏まえてクロマトグラムの平滑化に汎用される Savitzky-Golay 法 (SG 法) を小規模ながら Excel で実施し，クロマトグラムの微分でピーク分離ができるか確認しようと思う．
クロマトグラムを微分してみる
なぜ微分するか
今，ある試料を HPLC 分析したところ Fig. 1 A のようなクロマトグラムが得られたとする．なおこれはガウス関数で近似した二つのピーク (Fig. 1 B) を足し合わせて作った仮想的なクロマトグラムである (クロマトをガウス関数で近似してよいかはまたどこかで)．
さて今回は「二つのピークが重なっている」という「正解」のある問題だが，HPLC 分析を実施して Fig. 1 A のような (いわゆる「肩のある」) ピークが出た時，即座に 2 化合物であると断ずるのは難しい．イマドキの市販の HPLC は操作用 PC にデータ解析ソフトが搭載されており，自動でこうしたピークを分割してくれるものが多いが，これはクロマトグラムを微分することによりピーク数を検出している．その原理についてはおおよそ下記の通りである．
まず理想的なクロマトグラフィーでは化合物のピークをガウス関数で近似できることはすでに述べた．Fig. 2 に示す通り，ガウス関数は二階微分でピークを上下反転させたような形になる．これは Fig. 1 A のようなガウス関数を足し合わせた場合も同様であり，重なったピークでも二階微分により二峰性がより明瞭になることが知られている．
たとえば Fig. 1 のクロマトの場合，二階微分すると FIg. 3 下段のようになる．生データだと「ショルダーピークがありそう」くらいに思えるが，二階微分だと明確に二峰性であると見てとれる．ピークトップの位置もおおむね元のピークと一致し，
微分の弱点
さてクロマトデータの解析で優れた効果を示す微分であるが，弱点もある．それはノイズに弱いということである．例として Fig. 1 のデータ (Fig. 4 灰色，理論値と記載) に平均 0, 分散 300 のガウス白色ノイズを付加したデータを作成し，これを実測値ということにする (Fig. 4 赤色)．この時，微分する前のデータは理論値とほぼ変わらないように見えるが，二階微分した後のデータはかなり激しく波打っていることがわかる．Excel でやっている都合上かなりデータ数が少ないのもあるが，これだけを見て二峰性のピークと断ずる勇気のある人はあまりいないだろう．そこで一旦ノイズの影響を低減する，言い換えるとギザギザのデータを滑らかにするような平滑化法が必要になってくる．
クロマトグラムを平滑化してみる
中央移動平均
COVID 19 に際して陽性者数の推移がニュースで取り沙汰されていたのは記憶に新しいが，そこで週移動平均を使用するニュースが多かったはずである．あれは「土日に検査数が減ることで陽性者が過小に見積もられ，月曜にその反動が来る」という周期性を相殺するために，当日 + 過去 6 日の一週間で平均をとり，大まかな流れをとらえていた．これを単純移動平均という．
クロマトデータのように，すでに全工程のデータがそろっている場合は対象のデータ点を中心に前後 2n+1 点のデータの平均をとり，その平均値を代入する方法で平滑化する方法があり，これを中央移動平均 (CMA: Centered Moving Average) という．
ノイズがランダムに生じるのであれば，複数点を取って割り算すれば相殺されてほぼゼロになるはずであり，また COVID 19 陽性者数のように周期性が明白な場合もこれを相殺することができる．したがって，CMA である程度データは滑らかになると予想できる．この時に平均値をとるデータ幅 (2n+1 点) を窓と呼ぶらしい．今回は窓を 5 点で取っている (つまり対象の点と前後 2 点)．実際に Fig. 4 と同一のデータで中央移動平均をとり，二階微分した結果が Fig. 5 である．
確かに微分する前のデータも二階微分した後のデータも滑らかでギザギザしていないことが見てとれる．………見てとれるのだが，しかしこれは完全にシングルピークであろう．
CMA はもっとも単純な平滑化法でありわかりやすくもあるが，このような弱点もあるという好例である．すなわち，クロマトピークのような大きく変動する条件において周辺の値に引きずられることでピークが「鈍って」しまうのである．Fig. 6 に理想的なピーク (理論値) に分散 0.5 平均値 0 のガウス白色ノイズを付加したデータ (実測値) と窓 5 点の中央移動平均をとった結果 (CMA) を示しているが，CMA で微分値が滑らかになったかわりにピークが幅広かつピークトップが小さくなっていることが見てとれる．これを解決するいろいろな手法の1つが重みつき平均をとる方法であり，代表的なものが Savitzky-Golay 法である．
Savitzky-Golay 法
重み付き平均をとる代表的な方法として Savitzky-Golay 法がある．詳細は上記 wikipedia リンクなどに譲るが，A. Savitzky と M. J. E. Golay が考案した平滑化法で，アメリカ化学会の Analytical Chemistry 誌で最も多く引用された論文として有名な手法である．
SG 法では 2n+1 点の単純な平均をとるのではなく，窓の 2n+1 点にフィッティングする多項式を最小二乗法により求め，その値を対象の点に代入することで平滑化を行う，ある種の重み付き移動平均，らしい．今回は先の CMA と同様に窓を 5 点とし前回にならって三次関数をフィッティングすることとした．
Fig. 6 と同じデータを上記条件で SG 法により平滑化した結果が Fig. 7 である．図に示す通り，ピークの「鈍り」なく理論値に近い挙動を再現できていることが見てとれる．少なくとも「理想的な」ピークから作成したダミーデータでは CMA より優秀そうである．
同様にして Fig. 4 と同一のデータで Savitzky-Golay 平滑化を行い，二階微分した結果が FIg. 8 である．
CMA と同様に二階微分した後のデータは処理前よりギザギザしていないことが見てとれる．かなり怪しいが，微分前のデータも若干肩があるように見え，二階微分した結果も二峰性に見えなくもない．少なくとも CMA よりは優秀そうである．
おわりに
なんか微妙な感じになってしまったが，おそらくデータが少なすぎるのが一因と思われる．微分なのでデータの取得幅が大きいと $${\cfrac{\Delta y}{\Delta x}}$$ の分母も分子も大きくなって二峰性とかの鋭敏な動きが読めなくなっている気がする．前回のを見た方はご存知と思うが，あの方法だとすべての点で四行四列の行列計算をして多項式近似を行う必要がある．そういうわけで取れるデータ点は 20 個が (やる気的に) 限界だった．半手動 SG filter とかアホなことせんでも R や C++ でまともなパッケージがありそうなので実務的にはそちらを使えばそれでよいはず．どこかで余力があれば数値計算プログラムもキチンと勉強してみたい．
この記事が気に入ったらサポートをしてみませんか？