Using pseudocolour in ggplot2 scatter plot to indicate density

Question

Does someone know how to create a graph like the one in the screenshot? I've tried to get a similar effect adjusting alpha, but this renders outliers to be almost invisible. I know this type of graph only from a software called FlowJo, here they refer to it as "pseudocolored dot plot". Not sure if this an official term.

I'd like to do it specifically in ggplot2, as I need the faceting option. I attached another screenshot of one of my graphs. The vertical lines depict clusters of mutations at certain genomic regions. Some of these clusters are much denser than others. I'd like to illustrate this using the density colors.

The data is quite big and hard to simulate, but here's a try. I doesn't look like the actual data, but the data format is the same.

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")

Any help is highly appreciated.

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")

In practice, I'd take the initial bandwidth guess and then figure out an optimal bandwidth. Apart from taking the lazy approach and just plotting the points w/o filtering (smoothScatter() filters everything but the outliers based on npoints) this is generating the "smoothed scatterplot" like the example you posted.

smoothScatter() uses different defaults, so it comes out a bit differently:

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() is going to show the outliers, but not as distinct points:

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")

This:

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())

gets you very close to the defaults smoothScatter() uses, but hackishly accomplishes most of what the nrpoints filtering code does in the smoothScatter() function solely by restricting the y axis limits.

Answer 2

Call me oldschool, but why not use panel.smoothScatter from package latticeExtra. It provides direct access to smoothScatter but given that it is a panel function it automatically applies it to each subset of the defined panels. You say you need "facetting" so lattice is an obvious choice as it is explicitly designed to produce small multiples (i.e. facets or, in lattice-speak, panels). Panels can easily be created with y ~ x | g, where g is the variable used to define the small multiples. For your example, this would simply be:

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)

This way you get full control over the smoothing function, no hacking needed.

Answer 3

Although it might be computational intensive to generate a plot with probably millions of points, here is a solution to color each point based on its local density (i.e. a ‘pseudo-colored’ dot plot).

Generic function to calculate the local density (comparably fast).

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
}

Apply this on each point given the grouping of chromosomes.

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)