在ggplot2散点图中使用伪彩色来表示密度

Question

有人知道如何创建像屏幕截图中那样的图形吗？ 我尝试通过调整alpha来获得类似的效果，但这将使异常值几乎不可见。 我只从一个叫做FlowJo的软件中了解到这种类型的图表，他们称之为“伪彩色点图”。 不确定这是否是正式术语。

我想在ggplot2中特别完成它，因为我需要分面选项。 我附上了另一张我的图表的屏幕快照。 垂直线表示某些基因组区域的突变簇。 其中一些簇比其他簇密集得多。 我想使用密度颜色来说明这一点。

数据很大且难以模拟，但这里有一次尝试。 它看起来与实际数据不同，但数据格式相同。

chr <- c(rep(1:10,1000))
position <- runif(10000, min=0, max=5e8)
distance <- runif(10000, min=1, max=1e5)
log10dist <- log10(distance)

df1 <- data.frame(chr, position, distance, log10dist)

ggplot(df1, aes(position, log10dist)) + 
  geom_point(shape=16, size=0.25, alpha=0.5, show.legend = FALSE) +
  facet_wrap(~chr, ncol = 5, nrow = 2, scales = "free_x")

非常感谢您的帮助。

- Peer Wünsche

---

感谢您快速的回复。看起来正是我想要的东西。然而，问题在于我需要 ggplot2 的分面选项。我将编辑帖子以提供更精确的示例。 - Peer Wünsche

那看起来像是一个六边形图。请参见 geom_hex。 - Roland

smoothScatter() 调用（通过 grDevices:::.smoothScatterCalcDensity()） KernSmooth::bkde2D() 然后过滤掉非异常值。您可以使用 ggalt::geom_bkde2d() 进行密度图绘制，并在其下方绘制 geom_point()。您没有提供任何数据供其他人模拟。 - hrbrmstr

这不是一个十六进制二元图。 - hrbrmstr

@hrbrmstr：不确定我是否理解你的意思，但是geom_bkde2D()似乎并没有给我想要的结果。也许我必须尝试smoothScatter，并将单个染色体粘贴在一页上以获得分面效果。我还会尝试使用hexbin。让我们看看。 - Peer Wünsche

Answer 1

library(ggplot2)
library(ggalt)
library(viridis)

chr <- c(rep(1:10,1000))
position <- runif(10000, min=0, max=5e8)
distance <- runif(10000, min=1, max=1e5)
log10dist <- log10(distance)

df1 <- data.frame(chr, position, distance, log10dist)

ggplot(df1, aes(position, log10dist)) + 
  geom_point(shape=16, size=0.25, show.legend = FALSE) +
  stat_bkde2d(aes(fill=..level..), geom="polygon") +
  scale_fill_viridis() +
  facet_wrap(~chr, ncol = 5, nrow = 2, scales = "free_x")

实际上，我会先进行初始带宽猜测，然后找出最佳带宽。除了采用懒惰的方法，直接绘制未经过滤的数据点（smoothScatter() 根据 npoints 过滤所有但离群值之外的数据）生成类似于您发布的示例的“平滑散点图”。

smoothScatter() 使用不同的默认值，因此结果略有不同：

par(mfrow=c(nr=2, nc=5))
for (chr in unique(df1$chr)) {
  plt_df <- dplyr::filter(df1, chr==chr)
  smoothScatter(df1$position, df1$log10dist, colramp=viridis)
}

geom_hex()会显示异常值，但并不是作为明显的点：

ggplot(df1, aes(position, log10dist)) + 
  geom_point(shape=16, size=0.25, show.legend = FALSE, color="red") +
  scale_fill_viridis() +
  facet_wrap(~chr, ncol = 5, nrow = 2, scales = "free_x")

这个：

ggplot(df1, aes(position, log10dist)) + 
  geom_point(shape=16, size=0.25) +
  stat_bkde2d(bandwidth=c(18036446, 0.05014539), 
              grid_size=c(128, 128), geom="polygon", aes(fill=..level..)) +
  scale_y_continuous(limits=c(3.5, 5.1)) +
  scale_fill_viridis() +
  facet_wrap(~chr, ncol = 5, nrow = 2, scales = "free_x") +
  theme_bw() +
  theme(panel.grid=element_blank())

该方法让你接近smoothScatter()的默认设置，但通过限制y轴的极限值，hackish地完成了smoothScatter()中的nrpoints过滤代码几乎所有功能。

- hrbrmstr

---

哇，看起来很不错。对于我的实际数据，geom_hex似乎效果最好。（但是我认为你发布的代码是错误的？！）另外两个解决方案不够好，会导致分辨率低下。是否可以增加对ggalt解决方案的灵敏度，以便单独的聚类变得可见？还有：如何将其更改为log刻度或..levels..。...我尝试了..count..但它不起作用。 - Peer Wünsche

Answer 2

虽然我有点老派，但为什么不使用来自latticeExtra包的panel.smoothScatter呢？它直接访问smoothScatter，但由于它是一个面板函数，因此会自动将其应用于定义的每个子集。你说你需要"facetting"，所以显然选择lattice，因为它专门设计用于生成小多个（即facet或在lattice中称为panel）。可以使用y~x|g轻松创建面板，其中g是用于定义小多个的变量。对于你的例子，这只需简单地写成：

library(latticeExtra)

chr <- c(rep(1:10,1000))
position <- runif(10000, min=0, max=5e8)
distance <- runif(10000, min=1, max=1e5)
log10dist <- log10(distance)

df1 <- data.frame(chr, position, distance, log10dist)

clrs <- colorRampPalette(brewer.pal(9, "Reds"))

xyplot(log10dist ~ position | chr, data = df1,
       panel = panel.smoothScatter, layout = c(5, 2),
       as.table = TRUE)

通过这种方式，您可以完全控制平滑函数，无需任何黑客手段。

- TimSalabim

Answer 3

虽然生成可能有数百万个点的图表可能需要大量计算，但这里有一种解决方案，可以根据每个点的局部密度（即“伪彩色”点图）对每个点进行着色。

通用函数用于计算局部密度（相对较快）。

densVals <- function(x, y = NULL, nbin = 128, bandwidth, range.x) {
  dat <- cbind(x, y)
  # limit dat to strictly finite values
  sel <- is.finite(x) & is.finite(y)
  dat.sel <- dat[sel, ]
  # density map with arbitrary graining along x and y
  map   <- grDevices:::.smoothScatterCalcDensity(dat.sel, nbin, bandwidth)
  map.x <- findInterval(dat.sel[, 1], map$x1)
  map.y <- findInterval(dat.sel[, 2], map$x2)
  # weighted mean of the fitted density map according to how close x and y are
  # to the arbitrary grain of the map
  den <- mapply(function(x, y) weighted.mean(x = c(
    map$fhat[x, y], map$fhat[x + 1, y + 1],
    map$fhat[x + 1, y], map$fhat[x, y + 1]), w = 1 / c(
    map$x1[x] + map$x2[y], map$x1[x + 1] + map$x2[y + 1],
    map$x1[x + 1] + map$x2[y], map$x1[x] + map$x2[y + 1])),
    map.x, map.y)
  # replace missing density estimates with NaN
  res <- rep(NaN, length(sel))
  res[sel] <- den
  res
}

对于给定染色体的分组，将此应用于每个点。

library(dplyr)
library(ggplot2)

df1 %>% group_by(chr) %>% mutate(point_density = densVals(position, log10dist)) %>% 
  arrange(chr, point_density) %>% 
  ggplot(aes(x = position, y = log10dist, color = point_density)) +
  geom_point(size = .5) +
  scale_color_viridis_c() +
  facet_wrap(vars(chr), ncol = 5, scales = "free_x")

（伪彩色点图）

- ninjaminb