在 ggplot2 散点图中使用伪彩色来指示密度
Using pseudocolour in ggplot2 scatter plot to indicate density
有人知道如何创建像屏幕截图中那样的图表吗？我试图通过调整 alpha 获得类似的效果，但这会使离群值几乎不可见。我只从一个叫做 FlowJo 的软件中知道这种类型的图表，在这里他们将其称为 "pseudocolored dot plot"。不确定这是否是官方术语。
我想专门在 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")
非常感谢任何帮助。
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")
实际上，我会先猜测初始带宽，然后计算出最佳带宽。除了采用懒惰的方法并仅绘制点 w/o 过滤（smoothScatter() 过滤除基于 npoints 的异常值之外的所有内容）这会像您发布的示例一样生成 "smoothed scatterplot" .
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 轴限制就可以巧妙地完成 nrpoints 过滤代码在 smoothScatter() 函数中所做的大部分工作。
叫我老派，但为什么不使用包 latticeExtra 中的 panel.smoothScatter。它提供对 smoothScatter 的直接访问，但鉴于它是一个面板功能，它会自动将其应用于已定义面板的每个子集。你说你需要 "facetting" 所以 lattice 是一个明显的选择，因为它被明确设计为产生小的倍数（即小平面，或者在格子中，面板）。使用 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)
这样您就可以完全控制平滑功能，无需修改。
虽然生成可能有数百万个点的图可能需要大量计算，但这里有一个解决方案，可以根据每个点的局部密度（即“伪彩色”点图）为每个点着色。
计算局部密度的通用函数（比较快）。
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")
(pseudo-colored dot plot)
有人知道如何创建像屏幕截图中那样的图表吗？我试图通过调整 alpha 获得类似的效果，但这会使离群值几乎不可见。我只从一个叫做 FlowJo 的软件中知道这种类型的图表，在这里他们将其称为 "pseudocolored dot plot"。不确定这是否是官方术语。
我想专门在 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")
非常感谢任何帮助。
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")
实际上，我会先猜测初始带宽，然后计算出最佳带宽。除了采用懒惰的方法并仅绘制点 w/o 过滤（smoothScatter() 过滤除基于 npoints 的异常值之外的所有内容）这会像您发布的示例一样生成 "smoothed scatterplot" .
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 轴限制就可以巧妙地完成 nrpoints 过滤代码在 smoothScatter() 函数中所做的大部分工作。
叫我老派，但为什么不使用包 latticeExtra 中的 panel.smoothScatter。它提供对 smoothScatter 的直接访问，但鉴于它是一个面板功能，它会自动将其应用于已定义面板的每个子集。你说你需要 "facetting" 所以 lattice 是一个明显的选择，因为它被明确设计为产生小的倍数（即小平面，或者在格子中，面板）。使用 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)
这样您就可以完全控制平滑功能，无需修改。
虽然生成可能有数百万个点的图可能需要大量计算，但这里有一个解决方案，可以根据每个点的局部密度（即“伪彩色”点图）为每个点着色。
计算局部密度的通用函数（比较快）。
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")
(pseudo-colored dot plot)