挑选出符合模式的峰

Question

我得到了 x 轴上的时间（秒）和 y 轴上的强度（以相对荧光单位或 rfu 表示）的数据。它是通过观察 DNA 片段通过相机而生成的——DNA 片段越大，时间就越长。有 23 个已知大小的片段（以 DNA 碱基对单位 bp 为单位），因此应该有 23 个峰。由于我知道 DNA 片段的大小（以 bp 为单位），因此我想使用线性模型重新校准 x 轴，从时间（秒）到碱基对 (bp)。

不幸的是，数据中存在大量噪声，会产生虚假峰值。自信地区分真假的唯一方法是假碱基与 DNA 碱基对中的预期模式不符。

我已在此链接中的名为 demo 的数据框中提供了一个示例的数据。不幸的是它太大了，无法粘贴到下面。

https://1drv.ms/t/s!AvBi5ipmBYfrhf0v_kvWuN2foLyBgg?e=RWfdXZ

我可以如下选出所有峰值。

library(ggplot2)
library(ggpmisc)
ggplot(demo, aes(x=time, y=rfu)) +
geom_line(size=0.1) +
stat_peaks(col = "red", span=5, ignore_threshold=0.1) +
theme_bw()

然而，这些是预期的 DNA 片段大小（以 bp 为单位）：

ladder <- c(50, 75, 100, 125, 150, 200, 250, 300, 350, 400, 450, 475, 500, 550, 600, 650, 700, 750, 800, 850, 900, 950, 1000)

这是一个选出了正确峰值的图：

有没有办法让峰值查找器（例如 stat_peak 或 ggpmisc:::find_peaks）在选择峰值以实现第二个图时考虑峰值模式？

附录：假峰可能与真峰处于相同的 rfu 范围内，因此峰高不能用于排除它们（请参阅下面不同示例中的示例）

Answer 1

对于给定的数据，我们可以通过查找 1000 到 6000 之间的局部最大值来隔离峰值。这会识别 23 个峰值。

我使用该slider包来识别 21 bp（时间 - 10 到时间 +10）范围内的最大值，然后排除rfu1000-6000 范围之外的点或与之前的点匹配的点。

library(tidyverse); library(slider)
demo %>%
  mutate(peak = rfu == slider::slide(
    rfu, max, .before = 10, .after = 10) & 
      rfu > 1000 & rfu < 6000 & rfu != lag(rfu)) %>%
  ggplot(aes(time, rfu, color = peak, alpha = peak)) +
  geom_point() 

我探索了是否有可能找到一个强力的“最佳拟合”，这可能会使它更加自动化，并且对与底层信号处于相同范围内的噪声更加鲁棒。我确信这是可能的，但它比我预期的更困难，因为“时间拉伸”在整个范围内变化，因此具有两个参数的简单偏移+缩放模型是不够的。梯子和数据之间的“时间延伸”逐渐变化，从开始时略高于 10 倍（即梯子上从 50 到 75 的 25 diff 转换为数据 1578 到 1831 中的 253 diff）到结束时的 5 倍以下。

在这种情况下，看起来二次拟合可能效果很好，但这可能无法转化为时间扭曲不同的其他数据。

如果运行之间的畸变完全一致，则使用内置畸变定义梯子可能更有用，如下面的“ladder_scaled”列。那么问题将简化为找到最适合数据的单个偏移值，在您的示例中为+1528。

ladder_scaled <- tibble::tribble(
      ~ladder, ~ladder_scaled,
   50 ,          50,
   75 ,         303,
  100 ,         557,
  125 ,         811,
  150 ,        1068,
  200 ,        1580,
  250 ,        2104,
  300 ,        2630,
  350 ,        3159,
  400 ,        3684,
  450 ,        4204,
  475 ,        4451,
  500 ,        4705,
  550 ,        5203,
  600 ,        5683,
  650 ,        6143,
  700 ,        6574,
  750 ,        6983,
  800 ,        7357,
  850 ,        7701,
  900 ,        8010,
  950 ,        8283,
  1000,         8523
  )

如果我们可以依赖运行之间一致的时间对齐，我们就可以解决时序偏移问题，从而提供与数据的最佳对齐，即使我们的数据中存在重叠的噪声峰值。关于这个特定示例的一个棘手的事情是，梯子“签名”并不是非常独特 - 梯子上大多数步骤之间的间距相似，因此即使您偏移了一个或一个，您也可以很好地适应大多数丢失的功能。梯子的两个“台阶”。

这是一种方法，我可以强制拟合一堆偏移值。我依赖于这样的假设（也许没有根据）：梯子具有一致的时间长度，并且唯一需要解决的变量是时间偏移。我相信，如果不能假设这一点，问题就会变得更加困难。

首先，我使用了更广泛的合理值，其中除了 23 个真实峰值之外，还包括 6 个噪声峰值：

demo %>%
  mutate(peak = rfu == slider::slide(
    rfu, max, .before = 10, .after = 10) & 
      rfu > 500 & rfu < 30000 & rfu != lag(rfu)) -> demo_labels
# includes 29 "peaks": 23 real + 6 noisy

然后，我从 0:5000 开始交叉一系列可能的偏移，并将每个峰值模糊连接到 200 个时间单位内可能的阶梯峰值。对于每个可能的偏移和峰值，我选择最接近的拟合，然后对于每个可能的偏移，我选择 23 个最接近的拟合。最后，我绘制了每个偏移值的最差对齐情况。这表明在原始数据中，大约 1528 的偏移量提供了最佳拟合。但这个特殊的“阶梯签名”的一个困难是，1976 年的偏移量大约高出 450，通过这种衡量标准看起来并没有特别糟糕，尽管它肯定是错误的。因此，可能需要更多的领域知识来识别比“最差拟合”更好的函数来选择良好的匹配。

library(fuzzyjoin)
demo_labels %>%
  filter(peak) %>%
  crossing(offset = seq(0, 5000, by = 2)) %>%
  mutate(time_adj = time - offset) %>%
  distance_left_join(ladder_scaled, 
                  by = c("time_adj" = "ladder_scaled"),
                  max_dist = 200) %>%
  mutate(time_error = time_adj - ladder_scaled) %>% 
  group_by(offset, time_adj) %>%
  slice_min(abs(time_error)) %>%
  group_by(offset) %>%
  slice_min(abs(time_error), n = 23) %>%  ### pick the 23 best fits %>%
  summarise(worst_error = max(abs(time_error)),
            matching_peaks = n()) %>%
  arrange(worst_error) %>% 
  ggplot(aes(offset, worst_error, alpha = matching_peaks == 23)) +
  geom_point()