Claude Codeでクロマトグラムのピークを検出し面積を検算する

クロマトグラムのピーク検出と面積計算は、SciPyのfind_peaks、peak_widths、trapezoidの3関数で組めます。Claude Codeにスクリプトを書かせ、実行して出力を読ませるところまでは数分です。難しいのは、出てきた保持時間と面積を信じてよいかの判定です。

この記事では、面積の真値が分かる合成クロマトグラムを作り、検出したピークの面積と突き合わせる型を示します。検出の閾値や積分の窓を変えると面積がどれだけ動くかも、実際に走らせた数値で並べます。

先に答えを持つ — 合成データで検算する理由

実測データのピーク面積には、正解がありません。スクリプトが間違っていても、もっともらしい数字が出るだけです。そこで最初に、真値を自分で決めたデータを作ります。

ガウス型のピークの面積は、高さA、標準偏差sとして A × s × √(2π) で求まります。これは手計算できる式なので、Claudeに計算させず、人間が電卓や別の経路で確かめられます。

今回の合成データは次の条件です。0〜10分を0.01分刻み(1000点)で並べ、3本のピークと直線ベースライン、ランダムノイズを足します。

ベースラインは 2.0 + 0.5 × t、ノイズは標準偏差0.3の正規乱数(乱数シード固定)です。検出と積分のコードがこの表に戻ってこられれば、少なくとも「真値がある条件では動く」ことは言えます。

検出から面積までのスクリプト

次のコードを、Claude Codeに書かせる対象とします。手順は、ピークを探す、各ピークの左右の端を決める、ベースラインを引く、台形則で積分する、の4段です。

import numpy as np
from scipy.signal import find_peaks, peak_widths
from scipy.integrate import trapezoid
 
rng = np.random.default_rng(42)
t = np.arange(0, 10.0, 0.01)  # 分
def g(t, A, mu, s):
    return A * np.exp(-0.5 * ((t - mu) / s) ** 2)
 
true = [(100.0, 2.0, 0.08), (60.0, 4.5, 0.10), (30.0, 7.0, 0.12)]
y = sum(g(t, *p) for p in true) + 2.0 + 0.5 * t + rng.normal(0, 0.3, t.size)
 
for A, mu, s in true:
    print("true", mu, round(A * s * np.sqrt(2 * np.pi), 3))
 
peaks, props = find_peaks(y, prominence=5, width=3)
w = peak_widths(y, peaks, rel_height=0.99)
for i, p in enumerate(peaks):
    l, r = int(np.floor(w[2][i])), int(np.ceil(w[3][i]))
    seg_t, seg_y = t[l:r + 1], y[l:r + 1]
    base = np.interp(seg_t, [seg_t[0], seg_t[-1]], [seg_y[0], seg_y[-1]])
    area = trapezoid(seg_y - base, seg_t)
    print(f"peak{i+1} rt={t[p]:.2f} h={y[p]:.2f} area={area:.3f} "
          f"window=({t[l]:.2f},{t[r]:.2f})")

SciPy 1.18.1、NumPy 2.5.4で実行すると、次の出力になりました。

python chrom.py
# true 2.0 20.053
# true 4.5 15.04
# true 7.0 9.024
# peak1 rt=2.00 h=103.10 area=20.033 window=(1.74,2.28)
# peak2 rt=4.50 h=64.21 area=15.125 window=(4.14,4.88)
# peak3 rt=7.00 h=35.83 area=9.048 window=(6.66,7.38)

保持時間は3本とも真値と一致し、面積の差は次のとおりです。

ノイズが乗ったデータで1%未満に収まっています。ピーク高さが103.10なのは、ベースライン分の約3(2.0 + 0.5×2.0)とノイズが上乗せされているためで、Aの100とずれていても不思議ではありません。

find_peaksの引数が何を絞るか

find_peaksは、両隣より大きいサンプルをすべて極大として拾い、指定した条件で絞り込みます。ノイズだらけの信号では、極大の数が桁違いに増えます。上のデータで条件を変えて、拾われた本数を数えました。

高さだけで絞ると、ベースラインが2〜7のあたりにあるため、ノイズの極大が大量に残ります。突出度(prominence)は、ピークが周囲の最低線からどれだけ突き出しているかを表す量です。ベースラインが傾いていても効くので、クロマトグラムでは高さよりこちらを主役にします。

使った条件の評価順にも意味があります。ドキュメントによれば、条件はplateau_size、height、threshold、distance、prominence、widthの順に評価されます。distanceは小さいピークから先に落とすので、近接する2本のうち小さいほうが消えることがあります。

返り値のpropsには、prominenceやwidthを指定したときだけprominencesやwidthsなどのキーが入ります。値を見たいだけで絞り込みたくないときは、(None, None)を渡します。

wlenは長いクロマトグラムの処理と値の両方に効く

wlenは、prominenceを計算する範囲をピークの前後の窓(サンプル数)に限る引数です。ピークは窓の真ん中に置かれ、長さは次の奇数に切り上げられます。ドキュメントは、データが長いときや極大が多いときの計算時間を減らす用途を挙げています。

窓がピークの裾より狭いと、窓の中に本当の谷底が入らず、prominenceが小さく出ます。上のデータでprominence=5のままwlenだけを変えると、3本のprominenceは次のようになりました。

ピーク1は標準偏差0.08分で、0.01分刻みなら約8サンプルです。窓が21サンプルだと片側0.1分しか見えず、突出度が半分近くに縮みます。検出本数は3本のままなので、閾値が高いと気づきにくい変化です。窓はピークの裾の幅より十分広く、目安としては最も太いピークの半値幅の数倍を取ります。

処理時間も変わります。乱数で作った20万点の信号(極大約2万9千個)でprominence=0.5を指定すると、wlenなしで0.033秒、wlen=2001で0.007秒でした。手元の1回の計測値なので傾向の目安です。

面積の窓をどこで切るか — rel_heightの落とし穴

面積の結果をいちばん大きく動かすのは、積分の窓です。peak_widthsのrel_heightは、幅を測る高さをprominenceに対する割合で指定します。1.0なら最も低い等高線、0.5なら突出度の半分の高さです。

ガウス型では、半値幅の区間だけでピーク面積の約76%が入ります。それでもrel_height=0.5で面積が3分の1ほどに落ちるのは、窓の両端が半値の高さにあるからです。上のコードは窓の両端を結んだ直線をベースラインとして引くので、半値の高さより下の部分が丸ごと差し引かれます。同じデータでrel_heightだけを変えて面積を出すと、次のようになりました。

0.5のままだと、ピーク1は真値の約33%しか出ません。コードがエラーなく動き、保持時間も正しいので、見た目では気づけません。表の最下行の真値と並べて初めて分かる不一致です。

ただし1.0に近づけるほど、隣のピークやノイズの裾まで窓に入ります。0.99はこの合成データで通った値であって、実試料の最適値ではありません。実データでは、装置ソフトが採用している積分の区切り方に窓を揃えて比べます。

実データのCSVを読み込む

合成データを実データに置き換えるには、tとyを作る部分だけを差し替えます。検出から面積までのコードは変えません。列が「時間(分)、検出器の信号」の2列のCSVを想定した例です。

data = np.genfromtxt("run1.csv", delimiter=",", skip_header=1)
t, y = data[:, 0], data[:, 1]
ok = np.isfinite(y)          # NaN と inf を除く
t, y = t[ok], y[ok]
print(len(t), np.median(np.diff(t)))   # 点数とサンプリング間隔(分)

最後の1行で、点数とサンプリング間隔を確かめます。合成データの0.01分刻みと違う場合は、width=3のようにサンプル数で指定した引数の意味が変わります。時間軸が秒なら、面積の単位も秒になります。装置ソフトの出力と比べるときは、単位を揃えてから比べます。

Claude Codeへの頼み方 — 実行させて、出力を読ませる

Claudeの暗算や目視の推定を、そのまま結果にしません。数値はコードの実行結果だけを根拠にします。この方針はCLAUDE.mdに書いておくと、セッションをまたいで効きます。公式のメモリのページも、指示は「確認できるくらい具体的に」書くよう勧めています。

## クロマトグラム解析のルール
- ピーク面積・保持時間は、必ず python で実行した出力から読む。暗算で補完しない
- ピーク検出は scipy.signal.find_peaks の prominence を使う。height だけで絞らない
- 面積は peak_widths(rel_height=0.99) の窓で台形則(scipy.integrate.trapezoid)で積分する
- 検算用に tests/synthetic.py の合成データ(真値あり)を先に流し、差が1%を超えたら止まって報告する
- NaN が含まれる列は find_peaks に渡す前に除去か補間をする

依頼は、例えば次のような文面になります。

chrom.py を実行して、合成データの3ピークの保持時間と面積を表にして。
真値(20.053 / 15.040 / 9.024)との差を%で出して、1%を超える行があれば
rel_height を変えた場合の面積も並べて。数値は実行結果から写すこと。

この型の要点は、Claudeに結論を言わせるのではなく、検算の表を作らせて人間が読むことです。真値との差が表に出ていれば、設定が妥当かどうかの判断材料が手元に残ります。

CLAUDE.mdの末尾の1行は、公式のドキュメントにある警告に由来します。find_peaksはNaNを含むデータで予期しない結果を返すことがあり、事前に取り除くか置き換えるよう書かれています。

実データに移すときの切り分け

合成データで1%以内に収まっても、実データで同じ精度になるとは限りません。差が出たときは、原因を次の順に切り分けます。

台形則と単純な総和の差も、最初に確認しておく項目です。合成データの最初の100点で、y.sum() × 0.01は2.232、trapezoid(y, t)は2.211でした。差は約1%で、端点の扱いの違いから出ます。面積の定義は1つに決めて、手計算側もそれに合わせます。

重なったピークは別問題

隣り合うピークが重なると、find_peaksの挙動も面積の意味も変わります。ピーク高さ100・80を保持時間4.00分と4.25分に置いた合成データ(sはどちらも0.10)で試しました。

検出されたのは4.01分と4.23分の2本です。ただし2本目のprominenceは5.2しかなく、prominence=5の閾値すれすれで拾われています。閾値を少し上げれば、2本目は消えます。

面積を、2本の真値(25.07と20.05)と並べます。peak_widthsの窓で、本編と同じ台形則を2本それぞれに当てた値です。

0.99では、ピーク1の窓(3.66〜4.61分)が2本をまとめて覆い、合計に近い45.21が出ます。ピーク2は突出度が5.2しかなく、窓が4.15〜4.28分と狭いまま、面積は0.48に留まります。2本の合計は真値の45.12に近くても、1本ずつの値は分離されていません。

ピークごとの面積が必要なら、ガウス関数などの当てはめや、谷の位置で垂直に切る積分を別途検討する局面です。find_peaksとpeak_widthsの窓だけでは、重なった2本の面積は分けられません。

まとめ

find_peaksは検出、peak_widthsは窓、trapezoidは面積と役割が分かれています。検算で効いたのは、真値のある合成データを先に流すことと、rel_heightの既定値0.5を面積の窓に流用しないことの2点です。

Claude Codeには、スクリプトを書かせて実行させ、真値との差を表で返させます。数値を写すのは実行結果だけにする、という規約をCLAUDE.mdに入れておけば、次のセッションでも同じ検算が再現できます。

別分野の同じ型は、CalculiXの片持ち梁の検算やRDKitによるSMILESの検算にあります。MATLABの解析コードを移すなら、数値結果の一致を確かめる手順が使えます。測定値の合成はGUMの不確かさバジェットにつながります。