alignDE 使用教學

層析圖峰對齊：小波峰形匹配 + 差分演化。alignDE 3.0.0 ｜入門篇與進階篇｜本文中所有程式碼都可以直接複製執行。

入門篇

對應檔案 01_basic.R，整份約 15 秒跑完。從「為什麼需要對齊」開始，把五個步驟一個函式一個函式拆開看，最後重現論文的模擬與真實資料結果。

library(alignDE)
dir.create("figs", showWarnings = FALSE)

讓圖上的中文正常顯示
需要兩件事：(1) UTF-8 的 locale，(2) 一個有中文字的字型。
Windows 的 R 通常兩者都已經滿足，直接 png(...) 即可；
Linux / macOS 建議明確指定字型，並確認 locale 是 UTF-8（不是 POSIX/C）。

cjk_font <- if (.Platform$OS.type == "windows") "Microsoft JhengHei" else
            if (Sys.info()[["sysname"]] == "Darwin") "PingFang TC" else "Noto Sans CJK TC"
pngw <- function(f, w = 1700, h = 1000) {
  png(file.path("figs", f), w, h, res = 150,
      type = if (capabilities("cairo")) "cairo" else "Xlib",
      family = cjk_font)
}

0. 為什麼需要 peak alignment

同一種化合物在不同次分析裡的滯留時間會漂移，原因包括儀器溫度、
流動相配比、管柱老化等等。漂移雖然只有幾秒，但對「把整張層析圖
當成變數向量」的化學計量方法（PCA、PLS、PCR）是致命的：
這些方法假設「同一個欄位 = 同一個化合物」，漂移一破壞這個假設，
模型學到的就變成位移而不是組成差異。

套件內建兩組資料，先看第一組（模擬層析圖）：
  p1 = 參考層析圖 (reference, R)
  p2 = 待對齊層析圖 (to be aligned, C)，除了第二個峰以外都往右位移 50 點

data(simulate)

pngw("b01_problem.png", 1700, 700)
par(mar = c(4.2, 4.2, 3, 1))
plot(p1, type = "l", lwd = 1.5, xlab = "Sample intervals", ylab = "mAU",
     main = "對齊前：同樣的峰落在不同位置")
lines(p2, lty = 3, col = "#D9541E", lwd = 1.5)
legend("topleft", c("p1 參考", "p2 待對齊"), lty = c(1, 3),
       col = c("black", "#D9541E"), bty = "n")
dev.off()

cat("對齊前的相關係數 =", round(similarity(p1, p2), 4), "\n")

0.2401  —— 幾乎不相關，但這兩張圖其實是同一組峰

1. 安裝

alignDE 3.0.0 需要 R >= 3.5.0，並相依 DEoptim 與 Matrix：

  install.packages(c("DEoptim", "Matrix"))
  install.packages("path/to/alignDE_3.0.0.tar.gz", repos = NULL, type = "source")

檢查：
  packageVersion("alignDE")

2. alignDE 的五個步驟

(1) 長度均一化      不同取樣點數的層析圖用線性內插對齊長度
  (2) 峰匹配          連續小波轉換 (CWT) 找峰位置與峰寬
  (3) 基線校正        懲罰最小平方法 (Whittaker smoother)
  (4) 峰叢集          間距太近的峰合併成一個叢集，一起對齊
  (5) 峰對齊          差分演化 (DE) 找讓相關係數最大的位移量

對應的函式：
  cwt() -> getLocalMaximumCWT() -> getRidge() -> identifyMajorPeaks()
        -> widthEstimationCWT() -> baselineCorrectionCWT()
        -> peakClustering() -> alignDE()

3. 第二步：峰匹配，一個函式一個函式看

3.1　連續小波轉換

scales 是「尺度」，可以想成一系列寬度不同的尺，用來量峰有多寬。
峰寬大約等於 scale 的兩倍，所以要涵蓋最寬的峰。

scales <- seq(1, 56, 1)
wCoefs <- cwt(p1, scales = scales, wavelet = "mexh")   # mexh = 墨西哥帽小波
dim(wCoefs)          # 900 (資料點) x 56 (尺度)

3.2　每個尺度上找局部極大

localMax <- getLocalMaximumCWT(wCoefs)

3.3　把跨尺度的極大點串成脊線 (ridge)

真正的峰在所有尺度上都會出現，串起來就是一條從細尺度延伸到粗尺度的線；
雜訊只在最細的尺度出現，脊線很短。這就是用「脊線長度」濾雜訊的原理。

ridgeList <- getRidge(localMax, gapTh = 3, skip = 2)
length(ridgeList)    # 找到幾條脊線

3.4　依 SNR 與脊線長度挑出主峰

majorPeakInfo <- identifyMajorPeaks(p1, ridgeList, wCoefs,
                                    SNR.Th      = 3,   # 訊噪比門檻
                                    ridgeLength = 5)   # 脊線至少要跨過的尺度
majorPeakInfo$peakIndex

3.5　用 Haar 小波估計每個峰的起點與終點

peakWidth <- widthEstimationCWT(p1, majorPeakInfo)
data.frame(peak  = peakWidth$peakIndex,
           start = peakWidth$peakIndexLower,
           end   = peakWidth$peakIndexUpper)

pngw("b02_peak_detection.png", 1700, 1200)
par(mfrow = c(3, 1), mar = c(4.2, 4.2, 2.5, 1))
plot(p1, type = "l", xlab = "Sample intervals", ylab = "mAU",
     main = "(a) 原始層析圖與偵測到的峰")
points(majorPeakInfo$peakIndex, p1[majorPeakInfo$peakIndex], col = "red", cex = 1.3)
image(seq_len(nrow(wCoefs)), scales, wCoefs, col = terrain.colors(256),
      xlab = "index", ylab = "CWT scale", main = "(b) 2D 小波係數")
plotRidgeList(ridgeList, main = "(c) 脊線")
dev.off()

pngw("b03_peak_width.png", 1700, 700)
par(mar = c(4.2, 4.2, 3, 1))
plot(p1, type = "l", xlab = "Sample intervals", ylab = "mAU",
     main = "峰位置（圓圈）與起訖點（星號）")
points(peakWidth$peakIndex, p1[peakWidth$peakIndex], col = "red", cex = 1.4)
LR <- c(peakWidth$peakIndexLower, peakWidth$peakIndexUpper)
points(LR, p1[LR], pch = 8, col = "#2F6FEB", cex = 1.1)
dev.off()

4. 第三步：基線校正

基線會讓相關係數失真，所以對齊前要先扣掉。
baselineCorrectionCWT() 把峰的區間當成「遺失值」，用 Whittaker smoother
在峰底下內插出基線。
  lambda    越大基線越平滑（10~1000）
  threshold 峰形門檻，判斷峰是不是坐在斜的基線上（約 0.3）

backgr <- baselineCorrectionCWT(p1, peakWidth, threshold = 0.3, lambda = 100)
p1c <- p1 - backgr

pngw("b04_baseline.png", 1700, 800)
par(mfrow = c(2, 1), mar = c(4, 4.2, 2.5, 1))
plot(p1, type = "l", xlab = "", ylab = "mAU", main = "(a) 原始 + 擬合的基線")
lines(backgr, col = "#D9541E", lwd = 2)
plot(p1c, type = "l", xlab = "Sample intervals", ylab = "mAU",
     main = "(b) 基線校正後")
abline(h = 0, col = "#BBBBBB")
dev.off()

5. 把整套流程包成一個函式，兩張圖都跑一次

prepare <- function(x, scales = seq(1, 56, 1),
                    SNR.Th = 3, ridgeLength = 5,
                    lambda = 100, threshold = 0.3) {
  wCoefs    <- cwt(x, scales = scales, wavelet = "mexh")
  ridgeList <- getRidge(getLocalMaximumCWT(wCoefs), gapTh = 3, skip = 2)
  info      <- identifyMajorPeaks(x, ridgeList, wCoefs,
                                  SNR.Th = SNR.Th, ridgeLength = ridgeLength)
  peakWidth <- widthEstimationCWT(x, info)
  background <- baselineCorrectionCWT(x, peakWidth,
                                      threshold = threshold, lambda = lambda)
  list(corrected = x - background, peakWidth = peakWidth)
}

R <- prepare(p1)     # 參考
C <- prepare(p2)     # 待對齊

cat("基線校正後的相關係數 =", round(similarity(R$corrected, C$corrected), 4), "\n")

0.1188

6. 第四步：峰叢集

有肩峰、重疊峰的區域沒辦法個別對齊，間距小於 n 的峰要合併成一個叢集，
整個叢集一起位移。n 通常設 3~5。

cat("叢集前", length(C$peakWidth$peakIndex), "個峰 -> ")
pwC <- peakClustering(C$peakWidth, n = 5)
cat("叢集後", length(pwC$peakIndex), "個\n")

pngw("b05_clustering.png", 1700, 800)
par(mfrow = c(2, 1), mar = c(4, 4.2, 2.5, 1))
for (k in 1:2) {
  pw <- if (k == 1) C$peakWidth else pwC
  plot(C$corrected, type = "l", xlab = if (k == 2) "Sample intervals" else "",
       ylab = "mAU", main = sprintf("(%s) 峰叢集%s", letters[k],
                                    if (k == 1) "前" else "後"))
  plotwidthEstimation(C$corrected, pw)
}
dev.off()

7. 第五步：對齊

x         待對齊的（基線校正後）層析圖
  peakWidth 它自己的峰清單
  target    參考層析圖
  slack     每個峰最多可以移動幾點 —— 看圖估最大位移，抓寬一點
  n         幾個峰一組同時最佳化（預設 4）
  control   丟給 DEoptim 的參數
            NP      族群大小，慣例設 2 * slack
            itermax 世代數，150 以上比較保險
            trace   FALSE 就不印過程

set.seed(1)                      # DE 有隨機性，要可重現就設種子
aligned <- alignDE(C$corrected, pwC, R$corrected,
                   slack = 100, n = 4,
                   control = list(NP = 200, itermax = 150, trace = FALSE))

n <- min(length(aligned), length(R$corrected))
cat("對齊後的相關係數 =", round(similarity(R$corrected[1:n], aligned[1:n]), 4), "\n")

0.7819   （論文報告 0.7821）

pngw("b06_result.png", 1700, 1000)
par(mfrow = c(2, 1), mar = c(4, 4.2, 2.5, 1))
plot(R$corrected, type = "l", lwd = 1.5, xlab = "", ylab = "mAU",
     main = sprintf("(a) 對齊前  r = %.4f", similarity(R$corrected, C$corrected)))
lines(C$corrected, lty = 3, col = "#D9541E", lwd = 1.5)
legend("topright", c("參考", "待對齊"), lty = c(1, 3),
       col = c("black", "#D9541E"), bty = "n")
plot(R$corrected, type = "l", lwd = 1.5, xlab = "Sample intervals", ylab = "mAU",
     main = sprintf("(b) 對齊後  r = %.4f", similarity(R$corrected[1:n], aligned[1:n])))
lines(aligned, lty = 3, col = "#2F6FEB", lwd = 1.5)
legend("topright", c("參考", "已對齊"), lty = c(1, 3),
       col = c("black", "#2F6FEB"), bty = "n")
dev.off()

8. 換真實資料試一次

IBT11415 / IBT15670 是 HPLC-DAD 202 nm 的真菌萃取物層析圖。

data(IBT)

Rr  <- prepare(IBT11415)
Cr  <- prepare(IBT15670)
pwr <- peakClustering(Cr$peakWidth, n = 5)

set.seed(1)
alignedR <- alignDE(Cr$corrected, pwr, Rr$corrected, slack = 100, n = 4,
                    control = list(NP = 200, itermax = 150, trace = FALSE))
nr <- min(length(alignedR), length(Rr$corrected))

cat(sprintf("IBT 相關係數：對齊前 %.4f -> 對齊後 %.4f\n",
            similarity(Rr$corrected, Cr$corrected),
            similarity(Rr$corrected[1:nr], alignedR[1:nr])))

0.7867 -> 0.9185   （論文報告 0.7862 -> 0.9189）

pngw("b07_real.png", 1700, 1000)
par(mfrow = c(2, 1), mar = c(4, 4.2, 2.5, 1))
plot(Rr$corrected, type = "l", xlab = "", ylab = "mAU", main = "(a) IBT 對齊前")
lines(Cr$corrected, lty = 3, col = "#D9541E")
plot(Rr$corrected, type = "l", xlab = "Sample intervals", ylab = "mAU",
     main = "(b) IBT 對齊後")
lines(alignedR, lty = 3, col = "#2F6FEB")
dev.off()

9. 參數速查

函式                     參數           建議值      說明
 cwt()                    scales         1:56        要涵蓋最寬的峰（峰寬約 2x scale）
 identifyMajorPeaks()     SNR.Th         0~3         越小抓到的峰越多
                          ridgeLength    5~10        脊線長度門檻，濾雜訊
                          minNoiseLevel  1e-3        相對 max(wCoefs) 的雜訊下限
 baselineCorrectionCWT()  lambda         10~1000     越大基線越平滑
                          threshold      0.3         峰形門檻
 peakClustering()         n              3~5         間距門檻，合併相鄰峰
 alignDE()                slack          看資料      最大位移點數
                          n              4           一組幾個峰
                          control$NP     2 * slack   DE 族群大小
                          control$itermax >= 150     DE 世代數

下一步：02_advanced.R 講參數怎麼調、多樣本批次對齊、
        過載峰的處理，以及怎麼接到 PCA / PLS-DA。

cat("\n完成。圖在 figs/ 資料夾。\n")

展開：01_basic.R 的完整執行結果

1 
對齊前的相關係數 = 0.2401 
[1] 900  56
[1] 67
1_210 1_390 1_480 1_500 1_515 1_711 
  210   390   480   500   515   711 
      peak start end
1_210  210   190 226
1_390  390   378 403
1_480  480   470 490
1_500  500   486 510
1_515  515   510 524
1_711  711   684 738
          1 
          1 
          1 
基線校正後的相關係數 = 0.1188 
叢集前 6 個峰 -> 叢集後 4 個
          1 
對齊後的相關係數 = 0.7819 
          1 
IBT 相關係數：對齊前 0.7867 -> 對齊後 0.9185
          1 

完成。圖在 figs/ 資料夾。

進階篇

對應檔案 02_advanced.R，整份約 3 分鐘跑完。參數怎麼調、多樣本批次對齊、對齊品質怎麼驗證、兩個最常見的失敗模式，以及怎麼接到 PCA / PLS-DA。

library(alignDE)
library(parallel)
dir.create("figs", showWarnings = FALSE)

cjk_font <- if (.Platform$OS.type == "windows") "Microsoft JhengHei" else
            if (Sys.info()[["sysname"]] == "Darwin") "PingFang TC" else "Noto Sans CJK TC"
pngw <- function(f, w = 1700, h = 1000) {
  png(file.path("figs", f), w, h, res = 150,
      type = if (capabilities("cairo")) "cairo" else "Xlib", family = cjk_font)
}

data(simulate)
data(IBT)

prepare <- function(x, scales = seq(1, 56, 1), SNR.Th = 3, ridgeLength = 5,
                    minNoiseLevel = 0.001, lambda = 100, threshold = 0.3) {
  wCoefs    <- cwt(x, scales = scales, wavelet = "mexh")
  ridgeList <- getRidge(getLocalMaximumCWT(wCoefs), gapTh = 3, skip = 2)
  info      <- identifyMajorPeaks(x, ridgeList, wCoefs, SNR.Th = SNR.Th,
                                  ridgeLength = ridgeLength,
                                  minNoiseLevel = minNoiseLevel)
  peakWidth <- widthEstimationCWT(x, info)
  list(corrected = x - baselineCorrectionCWT(x, peakWidth, threshold = threshold,
                                             lambda = lambda),
       peakWidth = peakWidth, info = info, wCoefs = wCoefs)
}

1. 參數怎麼調

1.1　scales：尺度範圍要涵蓋最寬的峰

經驗法則：峰的半高寬大約是 2 x scale。所以最大的 scale 至少要有
「最寬的峰的一半」那麼大。設太小，寬峰的脊線會在尺度上端被截斷，
峰寬會低估；設太大只是浪費計算時間。

for (smax in c(8, 16, 32, 56)) {
  r <- prepare(p1, scales = seq(1, smax, 1))
  wd <- r$peakWidth$peakIndexUpper - r$peakWidth$peakIndexLower
  cat(sprintf("scales 1:%-3d -> %d 個峰，峰寬中位數 %3.0f 點\n",
              smax, length(r$peakWidth$peakIndex), median(wd)))
}

上限也不能無限加大：小波本身的長度會超過訊號長度。
這裡的訊號只有 900 點，scale 開到 64 就會報錯：

try(cwt(p1, scales = 64), silent = FALSE)

Error : scale 64 is too large!

遇到這個錯誤代表訊號太短或尺度設太大，兩者取其一調整。

1.2　SNR.Th 與 ridgeLength：抓多少峰

SNR.Th 低 -> 峰多（含雜訊峰）；ridgeLength 高 -> 只留跨越很多尺度的穩定峰。
對「對齊」而言，寧可少而準：錨定峰只要平均散布在整條層析圖上就夠了，
峰與峰之間的區域是用線性內插伸縮的。

grid <- expand.grid(SNR.Th = c(1, 3, 5, 10), ridgeLength = c(5, 10, 20))
grid$n_peak <- mapply(function(s, l) length(prepare(p1, SNR.Th = s,
                                                    ridgeLength = l)$peakWidth$peakIndex),
                      grid$SNR.Th, grid$ridgeLength)
print(grid)

1.3　minNoiseLevel：動態範圍很大時的關鍵參數

minNoiseLevel 預設 0.001，意思是「雜訊下限 = 0.001 x max(小波係數)」。
如果層析圖裡有一根巨大的峰，這個下限會連帶把所有小峰都壓掉。
真實樣品常常需要調到 1e-4 甚至 1e-5。

for (mn in c(1e-3, 1e-4, 1e-5)) {
  r <- prepare(IBT11415, minNoiseLevel = mn, SNR.Th = 3)
  cat(sprintf("minNoiseLevel = %.0e -> %3d 個峰\n", mn, length(r$peakWidth$peakIndex)))
}

1.4　lambda：基線的平滑度

pngw("a01_lambda.png", 1700, 1000)
par(mfrow = c(2, 2), mar = c(4, 4.2, 2.5, 1))
pw <- prepare(p1)$peakWidth
for (lam in c(1, 100, 1e4, 1e6)) {
  bg <- baselineCorrectionCWT(p1, pw, threshold = 0.3, lambda = lam)
  plot(p1, type = "l", xlab = "Sample intervals", ylab = "mAU",
       main = sprintf("lambda = %g", lam))
  lines(bg, col = "#D9541E", lwd = 2)
}
dev.off()

lambda 太小基線會跟著峰跑（把峰吃掉），太大則貼不上真實的基線起伏。
10~1000 之間通常沒問題。

1.5　peakClustering 的 n：什麼時候該合併

n 是「兩個峰之間允許的最大間距（點數）」。有肩峰、重疊峰時要調大。

pw2 <- prepare(p2)$peakWidth
for (g in c(1, 3, 5, 10, 30)) {
  cat(sprintf("gap = %2d -> %d 個峰叢集\n", g, length(peakClustering(pw2, n = g)$peakIndex)))
}

1.6　slack：最大位移量

slack 太小 -> 搜尋範圍不夠，對不齊；太大 -> 搜尋空間變大，需要更多世代，
而且可能把 A 峰對到 B 峰去。先用互相關估一下實際位移再決定。
用互相關估位移時，不要只算「整條圖一次」—— 不同的峰漂移量本來就不一樣，
一個全域位移量既不準也沒意義。正確做法是分段估，看漂移的範圍有多大。

lag_est <- function(x, y, window, lag.max = 200) {
  cc <- ccf(x[window], y[window], lag.max = lag.max, plot = FALSE)
  -cc$lag[which.max(cc$acf)]        # 負號：轉成「x 要移動多少才會對上 y」
}
Rtmp <- prepare(p1)$corrected; Ctmp <- prepare(p2)$corrected
wins <- list("150-320" = 150:320, "440-560" = 440:560, "600-800" = 600:800)
for (nm in names(wins))
  cat(sprintf("視窗 %s 的位移估計 = %4d 點\n", nm, lag_est(Ctmp, Rtmp, wins[[nm]])))

各段估出來的值不一樣（-41 / -70 ...），這正是需要「逐峰對齊」而不是
「整條圖平移」的理由。slack 要蓋得住其中最大的那個，再乘 1.5~2 倍。

R <- prepare(p1); C <- prepare(p2)
pwC <- peakClustering(C$peakWidth, n = 5)
slack_test <- sapply(c(20, 50, 100, 200), function(sl) {
  set.seed(1)
  a <- alignDE(C$corrected, pwC, R$corrected, slack = sl, n = 4,
               control = list(NP = 2 * sl, itermax = 150, trace = FALSE))
  n <- min(length(a), length(R$corrected))
  similarity(R$corrected[1:n], a[1:n])
})
print(data.frame(slack = c(20, 50, 100, 200), r = round(slack_test, 4)))

slack = 20 明顯不夠（實際位移 50 點），100 以上就穩定了。

2. 差分演化的隨機性與收斂

DE 是隨機演算法，不設種子每次結果都會略有不同。要可重現一定要 set.seed()。

reps <- sapply(1:8, function(s) {
  set.seed(s)
  a <- alignDE(C$corrected, pwC, R$corrected, slack = 100, n = 4,
               control = list(NP = 200, itermax = 150, trace = FALSE))
  n <- min(length(a), length(R$corrected))
  similarity(R$corrected[1:n], a[1:n])
})
cat(sprintf("8 次獨立執行：平均 %.4f，標準差 %.5f，範圍 %.4f-%.4f\n",
            mean(reps), sd(reps), min(reps), max(reps)))

這個例子夠簡單，8 次都收斂到同一個解。難一點的資料（峰多、slack 大、
itermax 小）就會出現差異，所以「有沒有設種子」和「itermax 夠不夠」
都要自己驗證，不要假設。

itermax 要設多少？跑收斂測試：

iters <- seq(20, 200, by = 20)
conv <- sapply(iters, function(it) {
  set.seed(1)
  a <- alignDE(C$corrected, pwC, R$corrected, slack = 100, n = 4,
               control = list(NP = 200, itermax = it, trace = FALSE))
  n <- min(length(a), length(R$corrected))
  similarity(R$corrected[1:n], a[1:n])
})

pngw("a02_convergence.png", 1500, 700)
par(mar = c(4.2, 4.2, 3, 1))
plot(iters, conv, type = "b", pch = 19, col = "#2F6FEB",
     xlab = "itermax（世代數）", ylab = "相關係數",
     main = "收斂測試：itermax 要設多少")
dev.off()

收斂之後再加世代數只是浪費時間。NP 建議 2 x slack，itermax 150 起跳。

3. 多樣本批次對齊

alignDE 一次只對齊「一張圖到一張參考圖」。多樣本時要先選參考。

選參考的三種常見做法：
  (a) medoid：與其他樣本相關係數中位數最高的那張（推薦，穩健）
  (b) 品質最好的那張（峰形最漂亮、基線最平）
  (c) QC 樣本或標準品

不要用「所有樣本的平均」當參考 —— 平均圖的峰已經被漂移抹平變寬了。

pick_reference <- function(X) {
  C <- cor(t(X))
  med <- apply(C, 1, function(v) median(v[-which.max(v)]))
  rownames(X)[which.max(med)]
}

先用「已知位移量」的人工資料驗證整套流程做得對不對。

shift_signal <- function(x, k) {
  n <- length(x); x[pmin(pmax(seq_len(n) - k, 1), n)]
}
true_shift <- c(base = 0, s30 = 30, sm25 = -25, s55 = 55, sm40 = -40, s18 = 18)
X1 <- t(sapply(true_shift, function(k) shift_signal(p1, k)))
rownames(X1) <- names(true_shift)

批次對齊的通用函式

align_all <- function(X, refName, slack = 100, n = 4, cores = 1, ...) {
  L    <- ncol(X)
  refc <- prepare(X[refName, ], ...)$corrected
  others <- setdiff(rownames(X), refName)
  res <- mclapply(others, function(s) {
    pr <- prepare(X[s, ], ...)
    pw <- peakClustering(pr$peakWidth, n = 5)
    set.seed(match(s, rownames(X)))
    a <- alignDE(pr$corrected, pw, refc, slack = slack, n = n,
                 control = list(NP = 2 * slack, itermax = 150, trace = FALSE))
    # 對齊後長度可能差幾點，補齊或截斷回原長度
    if (length(a) > L) a <- a[seq_len(L)]
    if (length(a) < L) a <- c(a, rep(a[length(a)], L - length(a)))
    list(sample = s, aligned = a, before = cor(refc, pr$corrected), after = cor(refc, a))
  }, mc.cores = cores)
  Xa <- X; Xa[refName, ] <- refc
  for (r in res) Xa[r$sample, ] <- r$aligned
  list(Xa = Xa, ref = refName,
       log = data.frame(sample = sapply(res, `[[`, "sample"),
                        cor_before = round(sapply(res, `[[`, "before"), 4),
                        cor_after  = round(sapply(res, `[[`, "after"),  4)))
}

ref1 <- pick_reference(X1)
cat("選到的參考層析圖：", ref1, "\n")
b1 <- align_all(X1, ref1, slack = 100, cores = max(1, detectCores() - 1))
print(b1$log)

已知位移的情況下，對齊後的相關係數應該非常接近 1。

4. 對齊品質怎麼驗證

不要只看相關係數 —— 它會被最大的那根峰主導。至少再看兩個指標：
  (a) 峰頂散布：挑幾個明確的峰，看頂點位置在樣本之間差多少
  (b) 目視疊圖：放大幾個區段直接看

apex_index <- function(M, centre, halfwin = 60) {
  apply(M, 1, function(v) {
    j <- max(centre - halfwin, 1):min(centre + halfwin, length(v))
    j[which.max(v[j])]
  })
}

markers <- prepare(p1)$peakWidth$peakIndex
qual <- t(sapply(markers, function(m) {
  c(range_before = diff(range(apex_index(X1, m))),
    range_after  = diff(range(apex_index(b1$Xa, m))))
}))
rownames(qual) <- paste0("peak@", markers)
cat("\n峰頂位置在樣本之間的差距（點數）：\n"); print(qual)

pngw("a03_batch_quality.png", 1700, 900)
par(mfrow = c(2, 1), mar = c(4, 4.2, 2.5, 1))
i <- 380:560
matplot(i, t(X1[, i]), type = "l", lty = 1, xlab = "", ylab = "mAU",
        main = "(a) 批次對齊前")
matplot(i, t(b1$Xa[, i]), type = "l", lty = 1, xlab = "Sample intervals",
        ylab = "mAU", main = "(b) 批次對齊後")
dev.off()

5. 疑難排解

5.1　過載峰讓 CWT 峰寬估計失效

症狀：基線校正後主峰被削掉一大截，或峰的起訖點落在上升沿的半山腰。
原因：widthEstimationCWT() 用 Haar 小波係數的局部極小值找峰界，
      過載（管柱超載）的峰很寬又極不對稱，這個估計會失準；
      上升沿被當成基線，接著就被 baselineCorrectionCWT() 扣掉。

判斷方法：算「峰界高度 / 峰高」。理想值應該很小（< 0.2）。

edge_ratio <- function(x, pw) {
  h <- sapply(seq_along(pw$peakIndex),
              function(k) max(x[pw$peakIndexLower[k]:pw$peakIndexUpper[k]]))
  pmax(x[pw$peakIndexLower], x[pw$peakIndexUpper]) / h
}
cat("\n模擬資料的峰界高度比中位數 =",
    round(median(edge_ratio(p1, prepare(p1)$peakWidth)), 3), "（正常）\n")

解法：不要用 CWT 估峰寬，改用門檻法直接從訊號取峰區間。
alignDE() 只需要 peakIndex / peakIndexLower / peakIndexUpper 三個欄位，
所以可以自己做一份 peakWidth 餵給它。

peakRegions <- function(x, block = 750, thresh_mult = 5, gap = 120, minw = 20,
                        lambda = 1e5) {
  n <- length(x)
  # (1) 粗基線：分段取最小值 -> 線性內插 -> 重度平滑
  nb  <- ceiling(n / block)
  idx <- rep(seq_len(nb), each = block)[seq_len(n)]
  bmin <- tapply(x, idx, min)
  bpos <- tapply(seq_len(n), idx, function(i) i[which.min(x[i])])
  base <- approx(as.numeric(bpos), as.numeric(bmin), xout = seq_len(n), rule = 2)$y
  base <- WhittakerSmooth(base, rep(1, n), lambda = lambda, differences = 2)

  resid <- x - base
  # (2) 雜訊水準用安靜區段的 MAD 估
  quiet <- resid[resid < quantile(resid, 0.75)]
  thr <- thresh_mult * mad(quiet)

  # (3) 超過門檻的連續區段，間距小於 gap 的合併
  r <- rle(resid > thr)
  ends <- cumsum(r$lengths); starts <- ends - r$lengths + 1L
  keep <- which(r$values); lo <- starts[keep]; hi <- ends[keep]
  LO <- HI <- numeric(0); ml <- lo[1]; mh <- hi[1]
  for (k in seq_along(lo)[-1]) {
    if (lo[k] - mh <= gap) mh <- hi[k]
    else { LO <- c(LO, ml); HI <- c(HI, mh); ml <- lo[k]; mh <- hi[k] }
  }
  LO <- c(LO, ml); HI <- c(HI, mh)
  ok <- (HI - LO) >= minw; LO <- LO[ok]; HI <- HI[ok]

  # (4) 邊界往外走到峰腳
  foot <- thr / 3
  for (k in seq_along(LO)) {
    i <- LO[k]; while (i > 1 && resid[i - 1] > foot &&
                       (k == 1 || i - 1 > HI[k - 1] + 1)) i <- i - 1
    LO[k] <- i
    j <- HI[k]; while (j < n && resid[j + 1] > foot &&
                       (k == length(HI) || j + 1 < LO[min(k + 1, length(LO))] - 1)) j <- j + 1
    HI[k] <- j
  }
  keep <- c(TRUE, LO[-1] > HI[-length(HI)] + 1); LO <- LO[keep]; HI <- HI[keep]

  apex <- sapply(seq_along(LO), function(k) LO[k] + which.max(x[LO[k]:HI[k]]) - 1)
  rng <- lapply(seq_along(apex), function(k) LO[k]:HI[k])
  names(rng) <- as.character(apex)
  c(list(peakIndex = apex, peakIndexLower = LO, peakIndexUpper = HI), rng)
}

pwCWT <- prepare(IBT11415, minNoiseLevel = 1e-4)$peakWidth
pwTHR <- peakRegions(IBT11415, thresh_mult = 5)
cat(sprintf("CWT 峰寬估計：%2d 個峰，峰界高度比中位數 %.3f\n",
            length(pwCWT$peakIndex), median(edge_ratio(IBT11415, pwCWT))))
cat(sprintf("門檻法      ：%2d 個峰，峰界高度比中位數 %.3f\n",
            length(pwTHR$peakIndex), median(edge_ratio(IBT11415, pwTHR))))

pngw("a04_peak_regions.png", 1700, 900)
par(mfrow = c(2, 1), mar = c(4, 4.2, 2.5, 1))
for (k in 1:2) {
  pw <- if (k == 1) pwCWT else pwTHR
  plot(IBT11415, type = "l", xlim = c(1, 1200), xlab = "Sample intervals", ylab = "mAU",
       main = sprintf("(%s) %s", letters[k],
                      if (k == 1) "widthEstimationCWT 的峰界" else "門檻法的峰界"))
  abline(v = pw$peakIndexLower, col = "#00808055")
  abline(v = pw$peakIndexUpper, col = "#80000055")
}
dev.off()

注意：改用門檻法取峰界時，通常就「不要」再跑 baselineCorrectionCWT()，
直接對原始訊號做對齊即可 —— 這樣比較的結果才單純只反映對齊的效果。

5.2　邊界峰：第一個峰太靠近訊號開頭

alignDE 的位移不是把峰整段搬走，而是「峰本身原封不動複製，靠伸縮峰與峰
之間的基線把它推到新位置」。所以第一個峰前面必須有足夠長的基線，
否則它沒辦法往右移；同理最後一個峰後面也要留空間。

經驗規則：lower[1] > slack 且 length(x) - upper[last] > slack。

IBT11415 的最大峰就在第 57 點，屬於典型的邊界峰。把它人工右移再對齊回來
就會失敗（往左移則沒問題，因為縮短前導基線永遠做得到）：

refI <- prepare(IBT11415, minNoiseLevel = 1e-4)
LI   <- length(refI$corrected)
fitlen <- function(a, L) if (length(a) > L) a[seq_len(L)] else
                         c(a, rep(a[length(a)], L - length(a)))

edge_demo <- function(k, pad = 0) {
  y  <- shift_signal(IBT11415, k)
  py <- prepare(y, minNoiseLevel = 1e-4)
  pw <- peakClustering(py$peakWidth, n = 5)
  tgt <- refI$corrected; sig <- py$corrected
  if (pad > 0) {                       # 前後各補一段平坦基線
    tgt <- c(rep(tgt[1], pad), tgt, rep(tgt[LI], pad))
    sig <- c(rep(sig[1], pad), sig, rep(sig[LI], pad))
    lo <- pw$peakIndexLower + pad; hi <- pw$peakIndexUpper + pad
    ap <- pw$peakIndex + pad
    rng <- lapply(seq_along(ap), function(j) lo[j]:hi[j]); names(rng) <- as.character(ap)
    pw <- c(list(peakIndex = ap, peakIndexLower = lo, peakIndexUpper = hi), rng)
  }
  set.seed(1)
  a <- alignDE(sig, pw, tgt, slack = 100, n = 4,
               control = list(NP = 200, itermax = 150, trace = FALSE))
  a <- fitlen(a, length(tgt))
  if (pad > 0) a <- a[(pad + 1):(pad + LI)]
  c(before = cor(refI$corrected, py$corrected), after = cor(refI$corrected, a))
}

edge <- t(sapply(c(-33, -22, 18, 40), function(k)
  c(k = k, no_pad = edge_demo(k)[["after"]], pad150 = edge_demo(k, pad = 150)[["after"]])))
cat("\n邊界峰：人工位移 k 點之後再對齊回來的相關係數\n")
print(round(edge, 4))

k 為負（峰被往左移，要往右推回去）時沒補基線就對不回來；
前後各補 150 點平坦基線之後就明顯改善。
實務上：確認方法的 dead time 有錄進去，或像上面一樣自己補。

6. 接到 PCA / PLS-DA

最重要的一課：對齊有沒有效，取決於後面用什麼 scaling。

  autoscaling（每個滯留時間點除以自己的標準差）
      -> 幾萬個基線雜訊點被放大到跟峰一樣重要，峰的資訊被稀釋，
         對齊的改善會被埋掉。
  mean-centering（只扣平均，不除標準差）
      -> 變異數大的峰主導模型，對齊的改善才看得出來。

建議流程：面積歸一化 -> 對齊 -> 中心化 -> PCA / PLS-DA

驗證一定要做：
  - LOOCV 正確率在小樣本下解析度很差（n = 14 時一個樣本就是 7%）
  - 加報 Q^2（交叉驗證的 R^2），它是連續的、比較敏感
  - 加做標籤置換檢定（permutation test），確認模型不是在記憶雜訊

pca_simple <- function(X, scale. = FALSE) {
  Xc <- scale(X, center = TRUE, scale = scale.)
  Xc[is.na(Xc)] <- 0
  s <- svd(Xc)
  list(scores = s$u %*% diag(s$d), var_exp = s$d^2 / sum(s$d^2) * 100)
}

pls1 <- function(X, y, ncomp = 2, scale. = FALSE) {
  n <- nrow(X); p <- ncol(X)
  xm <- colMeans(X); ym <- mean(y)
  xs <- if (scale.) sqrt(pmax(colSums(X^2) - n * xm^2, 0) / (n - 1)) else rep(1, p)
  xs[xs < .Machine$double.eps] <- 1
  E <- sweep(sweep(X, 2, xm), 2, xs, "/"); f <- y - ym
  W <- P <- matrix(0, p, ncomp); Tm <- matrix(0, n, ncomp); q <- numeric(ncomp)
  for (a in seq_len(ncomp)) {
    w <- crossprod(E, f); w <- w / sqrt(sum(w^2))
    t <- E %*% w; tt <- sum(t^2)
    pa <- crossprod(E, t) / tt; qa <- sum(f * t) / tt
    E <- E - t %*% t(pa); f <- f - qa * t
    W[, a] <- w; P[, a] <- pa; Tm[, a] <- t; q[a] <- qa
  }
  beta <- (W %*% solve(crossprod(P, W)) %*% q) / xs
  list(beta = beta, intercept = ym - sum(xm * beta), scores = Tm, W = W, q = q)
}

q2_loocv <- function(X, y, ncomp = 2, scale. = FALSE) {
  pr <- sapply(seq_len(nrow(X)), function(i) {
    f <- pls1(X[-i, , drop = FALSE], y[-i], ncomp, scale.)
    as.vector(X[i, , drop = FALSE] %*% f$beta) + f$intercept
  })
  1 - sum((y - pr)^2) / sum((y - mean(y))^2)
}

perm_test <- function(X, y, ncomp = 2, scale. = FALSE, nperm = 199, seed = 1) {
  set.seed(seed)
  obs <- q2_loocv(X, y, ncomp, scale.)
  null <- replicate(nperm, q2_loocv(X, sample(y), ncomp, scale.))
  list(Q2 = obs, p = (sum(null >= obs) + 1) / (nperm + 1))
}

兩組樣本（p1 家族 vs p2 家族），各自加上不同的位移

X2 <- rbind(a1 = p1, a2 = shift_signal(p1,  35), a3 = shift_signal(p1, -28),
            b1 = p2, b2 = shift_signal(p2,  42), b3 = shift_signal(p2, -19))
grp <- c(0, 0, 0, 1, 1, 1)
b2 <- align_all(X2, pick_reference(X2), slack = 100,
                cores = max(1, detectCores() - 1))

norm_area <- function(M) M / rowSums(M) * mean(rowSums(M))
for (nm in c("對齊前", "對齊後")) {
  M <- norm_area(if (nm == "對齊前") X2 else b2$Xa)
  p <- pca_simple(M)
  pt <- perm_test(M, grp, ncomp = 2)
  cat(sprintf("%s：PC1 %5.2f%%  PC2 %5.2f%%  Q2 = %6.3f  (置換檢定 p = %.3f)\n",
              nm, p$var_exp[1], p$var_exp[2], pt$Q2, pt$p))
}

7. 效能

alignDE 3.0.0 相對 2.0.1 的加速（3256 點 x 56 尺度）：
  cwt()                   1.9x   FFT 只算一次，重複用於所有尺度
  localMaximum()          7.6x   max.col() 取代 apply()
  widthEstimationCWT()    6.9x   Haar 局部極小圖依尺度快取
  WhittakerSmooth()       3.0x   帶狀稀疏矩陣 + 非置換 Cholesky
  alignDE()              >20x    改用 CRAN DEoptim 的 C 實作

大型資料的建議：
  - 用 mclapply()（Linux/macOS）或 parLapply()（Windows）平行處理樣本
  - 峰太多會讓 DEoptim 呼叫次數暴增：峰數 / n 就是要跑幾次最佳化，
    控制在 10~20 組通常就夠了
  - 先在一個樣本上把參數調好，再批次跑

if (requireNamespace("alignDE", quietly = TRUE)) {
  t1 <- system.time(cwt(IBT11415, scales = seq(1, 56, 1)))[["elapsed"]]
  t2 <- system.time(WhittakerSmooth(IBT11415, rep(1, length(IBT11415)), 100))[["elapsed"]]
  cat(sprintf("\n本機速度：cwt %.3f s，WhittakerSmooth %.3f s\n", t1, t2))
}

cat("\n完成。圖在 figs/ 資料夾。\n")

展開：02_advanced.R 的完整執行結果

scales 1:8   -> 6 個峰，峰寬中位數  22 點
scales 1:16  -> 6 個峰，峰寬中位數  22 點
scales 1:32  -> 6 個峰，峰寬中位數  24 點
scales 1:56  -> 6 個峰，峰寬中位數  24 點
Error : scale 64 is too large!
   SNR.Th ridgeLength n_peak
1       1           5      6
2       3           5      6
3       5           5      5
4      10           5      5
5       1          10      4
6       3          10      4
7       5          10      4
8      10          10      4
9       1          20      4
10      3          20      4
11      5          20      4
12     10          20      4
minNoiseLevel = 1e-03 ->  21 個峰
minNoiseLevel = 1e-04 ->  26 個峰
minNoiseLevel = 1e-05 ->  26 個峰
          1 
gap =  1 -> 4 個峰叢集
gap =  3 -> 4 個峰叢集
gap =  5 -> 4 個峰叢集
gap = 10 -> 4 個峰叢集
gap = 30 -> 4 個峰叢集
視窗 150-320 的位移估計 =  -41 點
視窗 440-560 的位移估計 =  -59 點
視窗 600-800 的位移估計 =  -70 點
  slack      r
1    20 0.1387
2    50 0.3190
3   100 0.7819
4   200 0.7819
8 次獨立執行：平均 0.7819，標準差 0.00000，範圍 0.7819-0.7819
          1 
選到的參考層析圖： base 
  sample cor_before cor_after
1    s30    -0.0575    0.9999
2   sm25    -0.0115    0.9999
3    s55    -0.0908    0.9999
4   sm40    -0.0790    0.9999
5    s18     0.1148    0.9999

峰頂位置在樣本之間的差距（點數）：
         range_before range_after
peak@210           95           0
peak@390           95           0
peak@480           98           0
peak@500           85           0
peak@515           95           0
peak@711           95           0
          1 

模擬資料的峰界高度比中位數 = 0.112 （正常）
CWT 峰寬估計：26 個峰，峰界高度比中位數 0.408
門檻法      ： 5 個峰，峰界高度比中位數 0.042
          1 

邊界峰：人工位移 k 點之後再對齊回來的相關係數
       k no_pad pad150
[1,] -33 0.4134 0.7432
[2,] -22 0.5508 0.8404
[3,]  18 0.9997 0.9999
[4,]  40 0.9998 0.9995
對齊前：PC1 23.93%  PC2 22.06%  Q2 = -0.472  (置換檢定 p = 0.710)
對齊後：PC1 62.30%  PC2 37.60%  Q2 =  0.914  (置換檢定 p = 0.120)

本機速度：cwt 0.015 s，WhittakerSmooth 0.003 s

完成。圖在 figs/ 資料夾。

附錄：參數速查表

函式	參數	建議值	調大的效果
cwt()	scales	1:56	能偵測到更寬的峰；上限受訊號長度限制
identifyMajorPeaks()	SNR.Th	0–3	峰變少、只留高訊噪比的峰
	ridgeLength	5–10	峰變少、只留跨越多尺度的穩定峰
	minNoiseLevel	1e-3 ~ 1e-5	調小才抓得到小峰；動態範圍大時是關鍵參數
baselineCorrectionCWT()	lambda	10–1000	基線更平滑
baselineCorrectionCWT()	threshold	0.3	更不容易判定成「坐在斜基線上的峰」
peakClustering()	n	3–5（點數）	更多相鄰峰被合併成一個叢集
alignDE()	slack	最大位移的 1.5–2 倍	搜尋範圍變大，但也更慢、更容易對錯峰
	n	4	一次同時最佳化更多峰，維度變高不易收斂
	control$NP	2 × slack	族群變大，較穩但較慢
	control$itermax	≥ 150	世代變多，收斂後就沒有幫助

附錄：三個最常見的坑

症狀	原因	解法
基線校正後主峰被削掉一大半	管柱過載的寬峰讓 `widthEstimationCWT()` 把上升沿判成基線	改用門檻法自訂 peakWidth（進階篇 5.1），並跳過基線校正
第一個峰對不回去，其他都對得很好	峰太靠近訊號開頭，前面沒有足夠基線可以伸縮	前後各補一段平坦基線再對齊（進階篇 5.2）
對齊後 PCA / PLS-DA 完全沒變好	用了 autoscaling，幾萬個雜訊點把峰的資訊稀釋掉	改成面積歸一化 + 中心化（進階篇 6）

參考文獻

Zhang, Z.-M., Chen, S. and Liang, Y.-Z. (2011) Peak alignment using wavelet pattern matching and differential evolution. Talanta 83(4), 1108–1117.

Zhang, Z.-M., Chen, S. and Liang, Y.-Z. (2010) Baseline correction using adaptive iteratively reweighted penalized least squares. Analyst 135, 1138–1146.

Du, P., Kibbe, W. A. and Lin, S. M. (2006) Improved peak detection in mass spectrum by incorporating continuous wavelet transform-based pattern matching. Bioinformatics 22, 2059–2065.