--- title: "HPLC 分離最佳化：從滯留資料到最佳條件" subtitle: "R 互動教學（Quarto 版）" lang: zh-Hant format: html: toc: true toc-depth: 3 number-sections: true code-fold: false code-tools: true theme: cosmo fig-width: 8 fig-height: 5 execute: warning: false message: false knitr: opts_chunk: dev: ragg_png dpi: 130 --- ## 開場 這份教學會帶你走過一次完整的 HPLC（高效液相層析，High-Performance Liquid Chromatography）分離最佳化流程：從一張「滯留時間」實驗表格出發，建立數學模型、 模擬層析圖、計算「解析度」（resolution），最後掃過整個實驗範圍找出最佳的移動相 比例。讀完並跑過所有程式碼之後，你會能夠： - 說明為什麼分析化學家用「滯留因子 k」而不是直接用滯留時間 $t_R$ 來建模； - 用 `nest()` + `map()` 對多個化合物同時擬合線性模型，並用 `broom::tidy()` 整理結果； - 畫出模擬層析圖，理解波峰的高度、寬度怎麼隨滯留時間改變； - 計算相鄰兩個波峰之間的解析度 $R_s$，並找出「最難分開的一對」； - 掃描一個實驗參數（此處是有機溶劑比例 f），畫出解析度與最後滯留時間的關係， 在時間預算內選出最佳條件。 本教學重現的方法出自： > Zisi, Ch., Pappa-Louisi, A., & Nikitas, P. (2020). Separation optimization > in HPLC analysis implemented in R programming language. *Journal of > Chromatography A*, 1617, 460823. **版權說明**：上述論文本文與其補充材料（`Data.xlsx`、`RChromOptim.RData` 等） 皆受著作權保護，**並未**隨本教學檔案附上或散布。本文件僅嵌入一小段具明確 出處標註的事實性數據摘錄（滯留時間表），目的僅在於讓讀者能夠自行重現計算 結果；若你需要完整資料集或原始程式碼，請直接向原期刊或原作者取得。 ```{r} #| label: setup library(tidyverse) # 讓圖上的中文正確顯示：自動挑一個系統上真的存在的中文字型。 # 沒有這一步，Windows 的預設繪圖裝置會把中文畫成空白方框（俗稱「豆腐」）。 cjk_font <- local({ cands <- c("Microsoft JhengHei", "Noto Sans TC", "PingFang TC", "Heiti TC", "Noto Sans CJK TC", "SimSun") avail <- if (requireNamespace("systemfonts", quietly = TRUE)) systemfonts::system_fonts()$family else character(0) hit <- cands[cands %in% avail] if (length(hit)) hit[1] else "" # 找不到就用預設值 }) theme_set(theme_minimal(base_size = 12, base_family = cjk_font)) ``` ## 資料 {#sec-data} 下表是九個化合物（以 A1–A9 標示）在四種有機溶劑比例 f（acetonitrile 的體積 分率）下量測到的滯留時間 $t_R$（單位：分鐘）。層析管柱是 Kinetex 2.6 µm XB-C18，150 x 4.6 mm；移動相為 acetonitrile / 水（pH 5.7）；已知的無滯留 時間（死時間）$t_0 = 1.4$ 分鐘。 資料出處：Zisi et al. (2020) 補充材料 `Data.xlsx`，工作表 `i-ret.fit` （僅摘錄用於教學重現，非完整資料集）。 ```{r} #| label: data t0 <- 1.4 dat <- tribble( ~f, ~A1, ~A2, ~A3, ~A4, ~A5, ~A6, ~A7, ~A8, ~A9, 0.40, 3.735, 4.671, 6.646, 7.662, 8.498, 10.760, 11.090, 15.990, 19.000, 0.45, 3.000, 3.650, 5.200, 5.600, 5.900, 6.850, 8.000, 10.380, 12.690, 0.50, 2.537, 3.010, 4.208, 4.349, 4.430, 4.842, 6.051, 7.239, 8.921, 0.60, 2.013, 2.265, 3.000, 2.940, 2.910, 3.061, 3.866, 4.198, 5.085 ) dat ``` 每一列是一次實驗（固定 f），每一欄（A1–A9）是該次實驗中，某個化合物從進樣 到偵測器出現訊號的時間。f 越大，代表移動相中「洗提能力較強」的有機溶劑比例 越高，通常會讓化合物更快沖出管柱（$t_R$ 變小）——你可以在表格中看到， f 從 0.40 增加到 0.60 時，每一欄的數字都在下降。 ## 滯留因子 k {#sec-k} 在層析學裡，我們很少直接對滯留時間 $t_R$ 建模，而是先換算成**滯留因子** （retention factor）： $$k = \frac{t_R - t_0}{t_0}$$ 為什麼要多做這一步轉換？因為 $t_R$ 本身包含了「化合物在管柱裡待多久」與 「化合物根本沒有跟固定相作用、純粹被流動相帶著走要花多久」（也就是 $t_0$） 兩部分。$t_0$ 會隨管柱長度、內徑、流速等儀器設定而改變，跟化合物本身的 化學性質無關。扣掉 $t_0$ 並除以 $t_0$ 之後得到的 k，去除了管柱長度與流速的 影響，只反映化合物與固定相之間真正的作用強弱，因此不同實驗室、不同管柱量到 的 k 值才具有可比較性、也才適合拿來建立化學模型。 ```{r} #| label: add-k long <- dat |> pivot_longer(-f, names_to = "solute", values_to = "tR") |> mutate( k = (tR - t0) / t0, lnk = log(k) ) long ``` ## 模型 1：ln k 對 f 的線性關係 {#sec-model1} Zisi et al. (2020) 使用的第一個模型非常簡單：對每個化合物，$\ln k$ 與有機 溶劑比例 f 呈線性關係： $$\ln k = c_0 - c_1 f$$ $c_0$ 可以理解為「f 趨近於 0 時的 $\ln k$」（即化合物在幾乎純水的移動相下 理論上有多強的滯留），$c_1$ 則描述滯留因子對 f 的敏感程度——$c_1$ 越大， 代表這個化合物的滯留隨溶劑比例改變得越劇烈。 我們用 `nest(.by = solute)` 把資料依化合物分組、each 組收成一個小表格， 再用 `map()` 對每組各自跑一次線性迴歸 `lm(lnk ~ f)`，最後用 `broom::tidy()` 把每個模型的係數整理成規整的表格。 ```{r} #| label: fit-model1 fits <- long |> nest(.by = solute) |> mutate( model = map(data, \(d) lm(lnk ~ f, data = d)), tidy = map(model, broom::tidy) ) params <- fits |> select(solute, tidy) |> unnest(tidy) |> select(solute, term, estimate) |> pivot_wider(names_from = term, values_from = estimate) |> rename(c0 = `(Intercept)`, slope = f) |> mutate(c1 = -slope) |> select(solute, c0, c1) params ``` 自己動手檢查一下：`c0` 欄位的範圍應該落在約 **3.1404 到 5.6030** 之間， `c1` 欄位的範圍應該落在約 **5.9113 到 8.5366** 之間。這兩個範圍剛好對應 原論文補充材料 `i-optim` 工作表中發表的擬合參數——代表我們這裡用 R 從頭 重新擬合出來的模型，跟原作者用他們自己的軟體算出來的結果是一致的。這是一個 你可以自己驗證正確性的檢查點，不需要相信任何人說的話，跑一次程式碼就能 看到答案。 ```{r} #| label: check-ranges range(params$c0) range(params$c1) ``` ## 診斷圖：ln k 對 f {#sec-diagnostic} 把實驗點（散點）跟擬合出來的直線畫在一起，可以直觀檢查模型合不合理——如果 每個化合物的點都乖乖落在自己的直線附近，代表線性模型是合理的近似。 ```{r} #| label: fig-diagnostic #| fig-cap: "ln k 對 f 的實驗點與模型 1 擬合直線" f_grid <- tibble(f = seq(0.30, 0.60, by = 0.01)) pred_lines <- params |> cross_join(f_grid) |> mutate(lnk = c0 - c1 * f) ggplot() + geom_line( data = pred_lines, aes(x = f, y = lnk, color = solute), linewidth = 0.8 ) + geom_point( data = long, aes(x = f, y = lnk, color = solute), size = 2 ) + labs( x = "f（acetonitrile 體積分率）", y = "ln k", color = "化合物", title = "模型 1：ln k = c0 - c1 f" ) + theme_minimal() ``` ## 模擬層析圖 {#sec-sim} 有了每個化合物在任意 f 下的滯留時間，我們就可以模擬出一整張層析圖：每個 化合物在圖上是一個波峰（peak），中心位置在它的 $t_R$，高度和寬度則跟它的 $t_R$ 有關（滯留越久的化合物，波峰通常會更寬、更矮，這是層析學裡常見的 「波峰擴散」現象）。 我們用以下簡化的波峰模型： - 波峰高度：$h = \max(0.014,\ h_0 + h_1 t_R)$（設一個最低高度，避免高度變成 負值或太小以致無法辨識） - 波峰標準差：$s = s_0 + s_1 t_R$ - 波峰底寬：$w = 4s / \sqrt{2}$ - 波峰形狀：高斯函數 $y = h \exp\!\left(-\dfrac{(t - t_R)^2}{s^2}\right)$ ```{r} #| label: shape-constants SHAPE <- list(h0 = 0.138, h1 = -0.0038, s0 = 0.013, s1 = 0.0105) predict_tR <- function(f, c0, c1) { lnk <- c0 - c1 * f k <- exp(lnk) t0 * (1 + k) } peak_height <- function(tR) { pmax(0.014, SHAPE$h0 + SHAPE$h1 * tR) } peak_sigma <- function(tR) { SHAPE$s0 + SHAPE$s1 * tR } peak_width <- function(tR) { 4 * peak_sigma(tR) / sqrt(2) } ``` 現在選一個 f 值，算出每個化合物的滯留時間，再畫出模擬層析圖。 ```{r} #| label: fig-chromatogram #| fig-cap: "模擬層析圖" f_choice <- 0.45 # 練習：改這個數字（例如 0.40、0.50、0.55），再重新 render 看層析圖怎麼變 solutes_at_f <- params |> mutate( tR = predict_tR(f_choice, c0, c1), h = peak_height(tR), s = peak_sigma(tR) ) t_axis <- tibble(t = seq(0, max(solutes_at_f$tR) * 1.2, length.out = 2000)) curves <- solutes_at_f |> select(solute, tR, h, s) |> cross_join(t_axis) |> mutate(y = h * exp(-((t - tR) / s)^2)) ggplot(curves, aes(x = t, y = y, color = solute)) + geom_line(linewidth = 0.7) + labs( x = "時間（分鐘）", y = "訊號強度（模擬）", color = "化合物", title = paste0("模擬層析圖（f = ", f_choice, "）") ) + theme_minimal() ``` ## 解析度 Rs {#sec-rs} 有了波峰之後,我們最關心的是:**相鄰兩個波峰有沒有分開?** 這件事用 「解析度」（resolution）$R_s$ 來量化。對相鄰的兩個波峰（依滯留時間排序後 彼此相鄰的一對）： $$R_s = \frac{2 (t_{R2} - t_{R1})}{w_1 + w_2}$$ 其中 $w_1, w_2$ 是兩個波峰各自的底寬。$R_s$ 越大代表兩個波峰分得越開； 業界慣用的「基線分離」（baseline separation）門檻是 $R_s \geq 1.5$—— 低於這個值，兩個波峰的訊號會有肉眼可見的重疊。 要注意的是：一張層析圖裡通常有很多個波峰，我們真正在乎的不是「平均解析度」 或「隨便哪一對的解析度」，而是**所有相鄰對之中最小的那個 $R_s$**——因為 只要有一對分不開，這個分離條件就不合格。這一對就是「最難分的一對」。 ```{r} #| label: compute-rs compute_min_rs <- function(f, params) { sol <- params |> mutate(tR = predict_tR(f, c0, c1)) |> arrange(tR) w <- peak_width(sol$tR) n <- nrow(sol) tibble( f = f, pair = paste0(sol$solute[-n], "/", sol$solute[-1]), Rs = 2 * (sol$tR[-1] - sol$tR[-n]) / (w[-1] + w[-n]) ) |> mutate(tR_max = max(sol$tR)) } rs_table <- compute_min_rs(f_choice, params) rs_table ``` 在上面這張表裡，找出 `Rs` 最小的那一列，就是目前這個 f 值下「最難分的 一對」化合物。 ## 掃描最佳化 {#sec-scan} 只看單一個 f 值不夠——我們想知道在哪一個 f 之下，最難分的一對也能達到 $R_s \geq 1.5$，同時整張層析圖跑完的時間（最後一個波峰的滯留時間 $t_{R,\max}$）還在可接受的時間預算內。做法是把 f 從 0.30 掃到 0.60， 每隔 0.005 算一次「最小 $R_s$」與「$t_{R,\max}$」。 ```{r} #| label: scan-f t_max <- 20 # 時間預算（分鐘）：層析圖跑完不能超過這個時間 f_seq <- seq(0.30, 0.60, by = 0.005) scan <- map(f_seq, \(f) compute_min_rs(f, params)) |> list_rbind() |> slice_min(Rs, by = f, n = 1) scan ``` 把最小解析度與最後滯留時間都對 f 畫圖，並加上 $R_s = 1.5$ 的參考虛線： ```{r} #| label: fig-scan #| fig-cap: "最小解析度與最後滯留時間隨 f 的變化" p1 <- ggplot(scan, aes(x = f, y = Rs)) + geom_line(linewidth = 0.8, color = "steelblue") + geom_hline(yintercept = 1.5, linetype = "dashed", color = "firebrick") + labs(x = "f", y = "最小解析度 Rs（最難分的一對）") + theme_minimal() p2 <- ggplot(scan, aes(x = f, y = tR_max)) + geom_line(linewidth = 0.8, color = "darkgreen") + labs(x = "f", y = "最後滯留時間 tR_max（分鐘）") + theme_minimal() library(patchwork) p1 / p2 ``` 留意曲線在 **f ≈ 0.35 附近有一個明顯的凹陷**：這裡最小 $R_s$ 掉到只有約 **0.235**，代表在這個 f 值下有兩個化合物幾乎完全重疊（共流出， co-elution）。這正是為什麼我們要對整個範圍逐點掃描，而不是只憑經驗挑 三、五個 f 值去試——如果剛好沒試到 0.35 附近，你可能完全不會發現這個 危險區域，誤以為分離狀況一路平順。 最後，在時間預算 `t_max` 之內，找出讓最小解析度最大的 f： ```{r} #| label: find-optimum optimum <- scan |> filter(Rs >= 1.5, tR_max <= t_max) |> slice_max(Rs, n = 1) optimum ``` 在預設的 `t_max <- 20` 之下，最佳條件是 **f = 0.410**，此時最小解析度 **$R_s$ = 2.4535**、最後滯留時間 **tR_max = 17.13** 分鐘，最難分的一對是 **A4/A5**。你可以在 @sec-scan 的程式碼與圖中確認：f 再往下（更接近 0.35 的凹陷區）雖然某幾個 f 值理論上 $R_s$ 更大，但對應的 $t_{R,\max}$ 超過了 20 分鐘的時間預算，因此被排除；f = 0.410 是「合乎時間限制」的條件中， 解析度最好的一個。 ## 練習 {#sec-exercises} 以下三個練習，請直接修改對應的空白程式碼區塊。修改完後重新 render 整份 文件，觀察結果如何改變。 **練習 (a)**：把時間預算 `t_max` 從 20 改成 12，重新跑一次最佳化， 新的最佳 f 是多少？此時的最小解析度和最後滯留時間各是多少？ *提示：你只需要複製 @sec-scan 裡「找出最佳條件」那段程式碼，把 `t_max <- 20` 改成 `t_max <- 12`，其餘不用動。* ```{r} #| label: exercise-a # 你的程式碼： ``` **練習 (b)**：從 `params` 裡把化合物 A6 移除，重新掃描 f，看看最佳的 f 是否改變。為什麼移除一個化合物可能會（或不會）影響最佳條件？ *提示：用 `params |> filter(solute != "A6")` 得到新的參數表，再把它傳入 `compute_min_rs()` 重新掃描一次。* ```{r} #| label: exercise-b # 你的程式碼： ``` **練習 (c)**：把 `SHAPE` 裡的 `s1`（波峰隨滯留時間變寬的速率）改大或改小 （例如從 0.0105 改成 0.02 或 0.005），重新計算 @sec-scan 的掃描結果， 描述 $R_s$ 曲線發生了什麼變化。 *提示：`s1` 越大，波峰越寬，同樣的滯留時間差距換算出來的 $R_s$ 就會越小； 留意這對「最佳 f」的位置有沒有影響，還是只影響 $R_s$ 的數值大小。* ```{r} #| label: exercise-c # 你的程式碼： ``` ## 延伸閱讀 {#sec-further} > Zisi, Ch., Pappa-Louisi, A., & Nikitas, P. (2020). Separation optimization > in HPLC analysis implemented in R programming language. *Journal of > Chromatography A*, 1617, 460823. 本教學檔案旁邊還提供了完整的腳本版本，涵蓋更完整的功能（例如更多模型、 更完整的最佳化流程）：`rchromoptim_modern.R`（R 版）、 `rchromoptim_modern.py`（Python 版），以及對應的 Python notebook 版本。 若你想深入研究原始方法的完整實作細節，可以參考這些檔案。