如何在R中高效绘制基因区间事件的密度图

你有如下格式的基因区间数据集：

gene     start       end
1 RP1-178F15.4 153633035 153633399
2 RP1-178F15.4 153633035 153634110
3        CRTC2 153633035 153948664
4        CRTC2 153633035 153951365
5        CRTC2 153633035 153951365
6        CRTC2 153633035 153951891

其中start和end是基因事件的位置区间，需要绘制x轴为位置、y轴为事件发生频率的密度图。

你原本的代码通过展开区间内所有位置生成数据，导致数据量暴增、运行缓慢：

List <- NULL
for (i in (1:nrow(pan))){ 
  vec <- c(seq(pan$start[i],pan$end[i]))
  List[[i]] <- vec
}
vveec <- unlist(List)
df <- data.frame(cat=c(rep("pan",length(vveec) )), weight=vveec)

# 绘制密度图
d <- density(df$weight) 
plot(d)

高效解决方案 #

核心思路是利用区间断点的计数变化，无需展开所有位置，仅处理区间的起始和结束点，大幅降低数据量：

# 1. 生成事件点：每个start位置事件计数+1，每个end+1位置计数-1
events <- data.frame(
  position = c(pan$start, pan$end + 1),
  delta = c(rep(1, nrow(pan)), rep(-1, nrow(pan)))
)

# 2. 按位置排序事件点
events <- events[order(events$position), ]

# 3. 计算每个区间的累积覆盖度（即该区间内的事件发生次数）
events$coverage <- cumsum(events$delta)

# 4. 提取区间断点、长度及对应覆盖度，用于密度计算
breaks <- unique(events$position)
interval_lengths <- diff(breaks)
coverage_vec <- events$coverage[-nrow(events)] # 每个区间的覆盖度

# 5. 生成密度计算所需的加权数据
density_data <- data.frame(
  x = breaks[-length(breaks)] + interval_lengths / 2, # 区间中点
  weight = interval_lengths * coverage_vec # 权重为区间长度×覆盖度
)

# 6. 计算并绘制密度图
d <- density(density_data$x, weights = density_data$weight, 
             from = min(pan$start), to = max(pan$end))
plot(d, xlab = "位置", ylab = "密度", main = "基因事件位置密度图")

如果使用ggplot2绘制更美观的图形：

library(ggplot2)

ggplot(density_data) +
  stat_density(aes(x = x, weight = weight), 
               geom = "line", linewidth = 1, color = "#2c3e50") +
  geom_area(stat = "density", aes(x = x, weight = weight), 
            fill = "#3498db", alpha = 0.3) +
  labs(x = "位置", y = "密度", title = "基因事件位置密度图") +
  theme_minimal()

效率说明 #

这种方法仅处理2×nrow(pan)个数据点，避免了展开百万/千万级的位置数据，运行速度提升显著，同时完全保留了原始区间的事件频率信息。

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