用 alignDE 做 HPLC 層析圖對位、峰積分與檢量線

1. 為什麼需要 peak alignment

本文使用的是食品中苯甲酸與己二烯酸的反相 HPLC 資料——這是最常見的防腐劑檢驗項目之一， 兩者可在同一次沖提中以 UV 230 nm 同時定量^[15]。 本批資料共 5 個標準品（0.25 / 10 / 25 / 50 / 100 ppm）與 6 個樣品， 由 Shimadzu LCsolution 匯出為 ASCII 檔。

同一支管柱、同一套方法，不同時間跑出來的層析圖，滯留時間 (retention time, RT) 不會完全一樣。 造成漂移的原因包括管柱溫度、流動相組成微變、幫浦壓力、管柱老化等。

這批資料就是活生生的例子：sample-1 ~ sample-3 是上午跑的，整條層析圖比下午的標準品左移約 0.08 ~ 0.10 分鐘。下圖左為標準品、右為樣品，虛線是標準品的平均峰位置—— 可以清楚看到樣品的兩個峰都往左偏了。

Peak alignment 就是把每張圖的峰「推」回共同的時間軸上，讓後續可以用同一個積分視窗 處理所有樣品，也讓不同批次的資料可以直接比較。

2. alignDE 方法總覽

alignDE 出自 Zhang, Chen & Liang (2011)^{[1, 2]}，核心想法是 「用小波找出峰在哪裡，再用差分演化搜尋每個峰該平移多少」。 和常見的 COW (correlation optimised warping)^[3] 相比，它最大的優點是 只平移峰、只拉伸峰與峰之間的基線區段，因此峰形、峰高、峰面積都不會被扭曲—— 這對定量分析非常關鍵^[4]。

STEP 1連續小波轉換與波峰偵測

層析圖的峰近似高斯形。Mexican Hat 小波（高斯函數的二階導數）長得就像一個峰， 所以拿它去和訊號做連續小波轉換 (CWT)，在「峰的位置、且尺度 s 接近峰寬」的地方 係數會出現極大值。

下圖上半是參考層析圖（100 ppm 標準品），下半是 CWT 係數的絕對值 |W(t,s)| （顏色越紅越大），黑線是脊線 (ridge line)——把不同尺度下的局部極大值串起來的軌跡。 一條夠長、且訊噪比夠高的脊線，就代表一個真正的波峰。

偵測結果

值得注意的是，sample-3 與 sample-5 只偵測到 1 個峰——小波演算法認為它們的 SA 位置 沒有夠強的訊號。這和後面定量算出的「低於偵測極限」結論完全一致，等於是兩個獨立方法互相驗證。

STEP 2峰寬估計 = 積分邊界

知道峰頂在哪裡還不夠，要積分就必須知道峰從哪裡開始、到哪裡結束。 alignDE 改用 Haar 小波再做一次轉換：Haar 小波是一個階梯函數，它的轉換係數在 訊號斜率變化最劇烈的地方（也就是峰的起點與終點）會出現局部極值。 函式 widthEstimationCWT() 就是在峰頂兩側搜尋這些位置。

共同積分視窗

本流程的做法是：只從參考層析圖推導一組積分視窗，然後套用到所有對位後的圖。 由於 BA 與 SA 兩峰之間留有空隙，程式再依層析慣例把共用邊界移到谷底 (valley-to-valley)：

STEP 3基線校正

UV 偵測器的基線會隨梯度沖提、管柱溫度而漂移。若不扣掉，積分面積會多算一塊梯形， 而且不同樣品多算的量還不一樣。

alignDE 用 Whittaker 平滑器^{[6, 7]}（懲罰最小平方法）擬合基線： 在「貼近資料」與「基線要夠平滑」之間取平衡， 並且只讓非峰區的點參與擬合（峰區的權重設得很低），避免基線被峰拉高。

替代方案：airPLS（不需事先偵測波峰）

上述做法必須先知道峰在哪裡（步驟 1–2 的結果）才能設定權重。 同一團隊後續提出的 airPLS（自適應迭代重新加權懲罰最小平方法）^[8] 則不需要任何先驗資訊：先用全部點擬合一條基線，把高於基線的點視為波峰、權重歸零， 把低於基線的點依偏離量給予指數權重，反覆迭代直到殘餘負值夠小為止。

它可以直接建構在 alignDE 內建的 WhittakerSmooth() 之上，不需額外套件：

airPLS <- function(x, lambda = 100, porder = 1, itermax = 20, tol = 0.001) {
  m <- length(x); w <- rep(1, m); z <- x
  for (i in seq_len(itermax)) {
    z <- WhittakerSmooth(x, w, lambda, porder)
    d <- x - z; neg <- d < 0
    dssn <- sum(abs(d[neg]))
    if (dssn < tol * sum(abs(x))) break        # 收斂
    w[!neg] <- 0                                # 高於基線 → 視為波峰
    w[neg]  <- exp(i * abs(d[neg]) / dssn)      # 低於基線 → 依偏離量加權
    w[1] <- exp(i * max(abs(d[neg])) / dssn); w[m] <- w[1]
  }
  z
}

本流程實際跑了這個比較：同一組積分視窗、同一套加權檢量線，只換基線演算法， 並掃描 λ = 10 / 100 / 1000 / 10000。

結論一：對可定量的結果，兩種方法一致

8 個高於報告下限的結果中，CWT 與 airPLS (λ = 10–100) 的差異 除 sample-1 的 BA 為 4.0% 外，其餘全部在 1.8% 以內。 檢量線品質 airPLS 甚至略優（BA 的 R² 由 0.99707 提升到 0.99773， 回收率最大偏差由 8.90% 降到 8.01%）。換句話說，這批資料用不用 airPLS 都可以。

結論二：λ 不能設太大，否則會製造假訊號

灰色列為低於報告下限 0.25 ppm 的結果（判定皆為 N.D.，差異不影響結論）。「低 λ 群一致性」= CWT / airPLS_10 / airPLS_100 三者的全距相對中位數。

結論三：實務建議

本流程維持 CWT 法為主（它能利用步驟 1–2 已經算出的峰位置，在本資料這種 「微量峰坐在大峰拖尾上」的情境有實質優勢），並把 airPLS 當作獨立的交叉驗證： 兩法一致的結果可以放心報告；兩法分歧的結果（都落在報告下限附近）就該人工覆核。 若要改用 airPLS 當主方法，在腳本頂端把 BASELINE_METHOD 改成 "airpls" 即可。

STEP 4差分演化對位

到這裡我們已經知道每張圖有哪些峰、峰的起訖在哪。最後一步是問： 每個峰要平移幾個取樣點，才能和參考圖最像？

alignDE 把這個問題交給差分演化 (Differential Evolution, DE)^[9] ——一種族群式的全域最佳化演算法。 目標函數是「平移後的訊號與參考訊號的線性相關係數」，要最大化：

DE 的好處是不需要目標函數可微、也不容易卡在局部最佳解，代價是要跑很多次評估 （本例族群 60、150 世代）。alignDE 3.0.0 把最佳化交給 DEoptim 套件的 C 實作^[10]， 比舊版內建的純 R 版本快 1–2 個數量級，因此每張圖只需約 0.2–0.4 秒。

對位效果

用與參考圖的相關係數當指標，改善最大的是 sample-3.txt（0.356 → 0.822）。 套用的位移量也和肉眼觀察到的漂移吻合：sample-1 ~ sample-3 都被往右推了約 0.08 ~ 0.10 分鐘。

STEP 5峰積分

對位完成後，用梯形法在共同視窗內對基線校正後的訊號積分：

驗證一：和 LCsolution 原始峰表比對

本流程完全獨立於儀器軟體重新算了一次面積。五個標準品的十個峰， 和 LCsolution 報告的面積差異最大只有 3.02%—— 可以確認積分邏輯是正確的。

驗證二：對位真的不會改變峰面積嗎？

表格中的「未對位偏差」欄位，是拿沒有對位的訊號套用同一個視窗積分，再和對位後的結果比較。 要看這個欄位，必須只挑實際被平移過的層析圖才有意義（位移量為 0 的圖偏差必然是 0，不能當證據）。 真正被移動的幾張圖中，25ppm（位移 −0.025 min）偏差為 +0.02% / −0.45%， sample-1（位移 +0.083 min）為 +1.13% / −0.44%—— 都在積分誤差的量級內，證實 alignDE 確實保留了峰面積。

偏差特別大的三個（sample-2 的 BA、sample-3 與 sample-5 的 SA）面積都小於 2 萬 µV·s， 屬於雜訊等級的訊號，本來就沒有定量意義。

STEP 6檢量線

把五個標準品的面積對濃度做線性迴歸。這裡有一個容易踩的陷阱： 本組標準品濃度跨 400 倍（0.25 → 100 ppm）。 一般未加權的最小平方法會被最高濃度點主導，低濃度端嚴重失真。

為什麼是 1/x²，而不是 1/x？權重不該憑感覺挑。判準是「儀器響應的標準差 σ 隨濃度 x 怎麼變」^[12]：σ 若不隨濃度改變就不必加權；σ² 與 x 成正比時用 1/x；σ 與 x 成正比（即相對標準差固定）時用 1/x²。層析峰面積多屬最後一種—— 濃度愈高、絕對誤差愈大而相對誤差大致穩定，所以 1/x² 是合理起點。 實務上的做法是三種都跑一次、比較低濃度端的回收率再定案， 本流程的 results/weighting_comparison.csv 就是這個比較的輸出。 若要判斷資料是否真的具異質變異，可用殘差分析加統計檢定確認^[13]。

本文的 LOD 與 LOQ 依 ICH Q2(R2)^[14] 的定義計算： LOD = 3.3σ/S、LOQ = 10σ/S，其中 S 為檢量線斜率、σ 為響應的標準差 （加權迴歸下取最低標準點處的殘差標準差）。

線性檢核：回算回收率

檢查檢量線好不好，比 R² 更靈敏的指標是把標準品的面積代回檢量線、看能不能還原成原本的濃度。 R² 對高濃度點的誤差不敏感（0.99 以上看起來都很漂亮），但回收率會誠實地暴露問題。

STEP 7樣品定量

把樣品的積分面積代入檢量線反推濃度：

10. 參數速查與重現方式

本流程使用的參數

重現步驟

# 1. 安裝套件（僅需一次）
install.packages(c("DEoptim", "Matrix", "remotes"))
remotes::install_github("Tai-ShengYeh/alignDE")

# 2. 執行分析：產生 results_align/ 的 10 張圖與 6 個 CSV
Rscript hplc_alignDE.R

# 3. 產生這份教學網頁
python build_teaching_html.py

檔案清單

11. 參考文獻

依在本文中首次出現的順序編號；內文的上標數字可點選跳至對應條目。 所有 DOI 均經查證。

A. 峰對位方法

1. Zhang, Z.-M., Chen, S., & Liang, Y.-Z. (2011). Peak alignment using wavelet pattern matching and differential evolution. Talanta, 83(4), 1108–1117. doi:10.1016/j.talanta.2010.08.008alignDE 的原始論文，本流程步驟 1–4 的方法依據。
2. Zhang, Z.-M., Chen, S., & Liang, Y.-Z. alignDE: Peak Alignment Using Wavelet Pattern Matching and Differential Evolution. R 套件版本 3.0.1. https://github.com/Tai-ShengYeh/alignDE本流程實際呼叫的軟體實作。
3. Nielsen, N.-P. V., Carstensen, J. M., & Smedsgaard, J. (1998). Aligning of single and multiple wavelength chromatographic profiles for chemometric data analysis using correlation optimised warping. Journal of Chromatography A, 805(1–2), 17–35.COW（相關性最佳化扭曲）的原始論文，是 alignDE 的主要比較對象。
4. Tomasi, G., van den Berg, F., & Andersson, C. (2004). Correlation optimized warping and dynamic time warping as preprocessing methods for chromatographic data. Journal of Chemometrics, 18(5), 231–241. doi:10.1002/cem.859COW 與 DTW 的系統性比較，說明為何扭曲類方法會改變峰面積。

B. 小波峰偵測

5. Du, P., Kibbe, W. A., & Lin, S. M. (2006). Improved peak detection in mass spectrum by incorporating continuous wavelet transform-based pattern matching. Bioinformatics, 22(17), 2059–2065. doi:10.1093/bioinformatics/btl355CWT + 脊線 + 訊噪比的峰偵測框架（MassSpecWavelet），alignDE 的步驟 1 沿用此法。

C. 基線校正

6. Whittaker, E. T. (1922). On a new method of graduation. Proceedings of the Edinburgh Mathematical Society, 41, 63–75.懲罰最小平方平滑法的原始構想。
7. Eilers, P. H. C. (2003). A perfect smoother. Analytical Chemistry, 75(14), 3631–3636. doi:10.1021/ac034173t把 Whittaker 平滑器整理成現代化學計量學可直接使用的形式，即步驟 3 的基線擬合。
8. Zhang, Z.-M., Chen, S., & Liang, Y.-Z. (2010). Baseline correction using adaptive iteratively reweighted penalized least squares. Analyst, 135(5), 1138–1146. doi:10.1039/b922045cairPLS，同一團隊的基線校正改良版；不需事先偵測波峰，適合基線漂移嚴重的資料。

D. 最佳化演算法

9. Storn, R., & Price, K. (1997). Differential evolution — a simple and efficient heuristic for global optimization over continuous spaces. Journal of Global Optimization, 11(4), 341–359. doi:10.1023/A:1008202821328差分演化的原始論文，步驟 4 搜尋位移量所用的演算法。
10. Mullen, K., Ardia, D., Gil, D. L., Windover, D., & Cline, J. (2011). DEoptim: An R package for global optimization by differential evolution. Journal of Statistical Software, 40(6), 1–26. doi:10.18637/jss.v040.i06alignDE 3.0.0 改以此套件的 C 實作執行 DE，比舊版內建的純 R 版本快 1–2 個數量級。

E. 檢量線與方法確效

11. Almeida, A. M., Castel-Branco, M. M., & Falcão, A. C. (2002). Linear regression for calibration lines revisited: weighting schemes for bioanalytical methods. Journal of Chromatography B, 774(2), 215–222. doi:10.1016/S1570-0232(02)00244-1檢量線加權方案的實務依據；說明為何寬濃度範圍必須用 1/x 或 1/x² 加權。
12. Gu, H., Liu, G., Wang, J., Aubry, A.-F., & Arnold, M. E. (2014). Selecting the correct weighting factors for linear and quadratic calibration curves with least-squares regression algorithm in bioanalytical LC–MS/MS assays and impacts of using incorrect weighting factors on curve stability, data quality, and assay performance. Analytical Chemistry, 86(18), 8959–8966. doi:10.1021/ac5018265權重因子的選用準則出處：由響應標準差 σ 與濃度 x 的關係決定——σ 為常數用不加權、σ² 正比於 x 用 1/x、σ 正比於 x 用 1/x²。
13. Majer, D., & Finšgar, M. (2021). Single-drop analysis of epinephrine and uric acid on a screen-printed carbon electrode. Biosensors, 11(8), 285. doi:10.3390/bios11080285完整示範以殘差分析檢定異質變異、確認不適用 OLS 後改採加權迴歸，並直接比較兩者。
14. International Council for Harmonisation (2023). ICH Q2(R2) Validation of Analytical Procedures. ICH Harmonised Guideline, 採納於 2023 年 11 月 1 日（取代 Q2(R1)）. https://database.ich.org/sites/default/files/ICH_Q2(R2)_Guideline_2023_1130.pdfLOD = 3.3σ/S、LOQ = 10σ/S 的定義出處，以及線性、範圍、準確度的確效要求。

F. 苯甲酸／己二烯酸分析

15. Saad, B., Bari, M. F., Saleh, M. I., Ahmad, K., & Talib, M. K. M. (2005). Simultaneous determination of preservatives (benzoic acid, sorbic acid, methylparaben and propylparaben) in foodstuffs using high-performance liquid chromatography. Journal of Chromatography A, 1073(1–2), 393–397.食品中苯甲酸與己二烯酸同時定量的 RP-HPLC 方法，可對照本資料的分析條件。

本文引用方式