在R中针对DNA序列窗口识别计数数据的显著峰值

步骤1：数据预处理 #

先标准化两个数据框的列名，方便后续操作：

library(tidyverse)

# 重命名df2的列，让含义更清晰
df2 <- df2 %>%
  rename(seqnames = V1, position = V2, count = V3)

# 确保基因名字段格式一致，避免匹配出错
df$seqnames <- as.character(df$seqnames)
df2$seqnames <- as.character(df2$seqnames)

步骤2：定义显著峰值的识别规则 #

这里提供结合局部最大值和统计显著性的双规则判断逻辑，你可以根据需求调整：

• 局部最大值：位点计数大于左右相邻位点（边缘位点只需单侧大于）
• 统计显著性：计数高于当前窗口内所有位点的平均值+2倍标准差（倍数可自定义）

编写识别函数：

identify_peaks <- function(count_vec, pos_vec) {
  # 标记局部最大值
  is_local_max <- logical(length(count_vec))
  # 中间位点：同时大于左右
  is_local_max[2:(length(count_vec)-1)] <- 
    count_vec[2:(length(count_vec)-1)] > count_vec[1:(length(count_vec)-2)] &
    count_vec[2:(length(count_vec)-1)] > count_vec[3:length(count_vec)]
  # 第一个位点：大于右侧
  if(length(count_vec) >= 2) is_local_max[1] <- count_vec[1] > count_vec[2]
  # 最后一个位点：大于左侧
  if(length(count_vec) >= 2) is_local_max[length(count_vec)] <- count_vec[length(count_vec)] > count_vec[length(count_vec)-1]
  
  # 计算显著性阈值
  mean_count <- mean(count_vec)
  sd_count <- sd(count_vec)
  sig_threshold <- mean_count + 2*sd_count # 可改为3倍标准差提升严格度
  
  # 筛选同时满足两个条件的峰值
  peak_indices <- which(is_local_max & count_vec > sig_threshold)
  
  # 返回峰值的位置和计数
  data.frame(
    position = pos_vec[peak_indices],
    count = count_vec[peak_indices]
  )
}

步骤3：按窗口匹配并识别峰值 #

将df2的数据按df中的窗口拆分，逐个窗口运行峰值识别：

# 关联窗口与对应计数数据，批量识别峰值
peak_results <- df %>%
  rowwise() %>%
  mutate(
    # 筛选当前窗口内的计数数据（position落在start和end之间）
    window_data = list(df2 %>% filter(seqnames == !!seqnames, position >= start, position <= end)),
    # 调用峰值识别函数
    peaks = list(if(nrow(window_data) > 0) identify_peaks(window_data$count, window_data$position) else data.frame())
  ) %>%
  ungroup() %>%
  # 展开峰值结果
  unnest(peaks, keep_empty = TRUE) %>%
  # 保留核心字段
  select(seqnames, start, end, strand, position, count)

最终输出 #

peak_results即为每个窗口内的显著峰值列表，结构示例：

seqnames start  end strand position count
gene     65    1237 +      120      15
gene     65    1237 +      450      22
gene     1262  2134 +      1800     18

可选调整 #

• 若需考虑负链(strand == "-")，可在处理窗口数据时反转position顺序，保证局部最大值判断符合链的方向
• 可改用中位数+四分位距作为显著性阈值，避免极端值干扰
• 大数据量场景下，建议用data.table替代tidyverse提升处理速度

内容的提问来源于stack exchange，提问作者margo