Nicer scatter plots in GGgally ggpairs-ggduo

GGally is an extension of ggplot2. I use it mainly as a replacement/refinement of pairs in base R.

Here, I want to replace the default scatter plot used to compare 2 continuous variables with something nicer, especially adapted to plotting large number of points.

Load the library and some RNA-seq data:

library(GGally)
library(airway)
data(airway)

The airway dataset contains data for 4 cell lines treated or not with dexamethasone:

colData(airway)[,1:3]
SampleName cell dex
SRR1039508 GSM1275862 N61311 untrt
SRR1039509 GSM1275863 N61311 trt
SRR1039512 GSM1275866 N052611 untrt
SRR1039513 GSM1275867 N052611 trt
SRR1039516 GSM1275870 N080611 untrt
SRR1039517 GSM1275871 N080611 trt
SRR1039520 GSM1275874 N061011 untrt
SRR1039521 GSM1275875 N061011 trt

Let’s keep a fraction of the data for simplicity

df <- as.data.frame(assays(airway)$counts[,1:4]) #first 4 columns
df <- df[rowSums(df)>4,] #keep genes with some counts
set.seed(123)
df <- df[sample.int(nrow(df),1e3),] #random 1K genes

default behavior of GGally::ggpair:

GGally::ggpairs(log2(df+1))

Here is a better function to make the scatterplots in the lower triangle. It is largely based on this blog post and implements:

  • viridis coloring of points based on expression level (practically based on effect size captured by the first principal component). I added some code to get consistent coloring even if orientation of the first PC differs between pairs of columns.
  • alpha based on inverse of 2D density, i.e. low alpha when there are lots of points and higher when points are more scattered

Also, the use of GGally:eval_data_col is described here:

GGscatterPlot <- function(data, mapping, ..., 
                        method = "spearman") {

#Get correlation coefficient
    x <- GGally::eval_data_col(data, mapping$x)
    y <- GGally::eval_data_col(data, mapping$y)

    cor <- cor(x, y, method = method)
#Assemble data frame
    df <- data.frame(x = x, y = y)
# PCA
    nonNull <- x!=0 & y!=0
    dfpc <- prcomp(~x+y, df[nonNull,])
    df$cols <- predict(dfpc, df)[,1]
# Define the direction of color range based on PC1 orientation:
    dfsum <- x+y
    colDirection <- ifelse(dfsum[which.max(df$cols)] < 
                               dfsum[which.min(df$cols)],
                           1,
                           -1)
#Get 2D density for alpha
    dens2D <- MASS::kde2d(df$x, df$y)
    df$density <- fields::interp.surface(dens2D , 
                                         df[,c("x", "y")])

if (any(df$density==0)) {
    mini2D = min(df$density[df$density!=0]) #smallest non zero value
    df$density[df$density==0] <- mini2D
}
#Prepare plot
    pp <- ggplot(df, aes(x=x, y=y, color = cols, alpha = 1/density)) +
                ggplot2::geom_point(shape=16, show.legend = FALSE) +
                ggplot2::scale_color_viridis_c(direction = colDirection) +
#                scale_color_gradient(low = "#0091ff", high = "#f0650e") +
                ggplot2::scale_alpha(range = c(.05, .6)) +
                ggplot2::geom_abline(intercept = 0, slope = 1, col="darkred") +
                ggplot2::geom_label(
                        data = data.frame(
                                        xlabel = min(x, na.rm = TRUE),
                                        ylabel = max(y, na.rm = TRUE),
                                        lab = round(cor, digits = 3)),
                        mapping = ggplot2::aes(x = xlabel, 
                                               y = ylabel, 
                                               label = lab),
                        hjust = 0, vjust = 1,
                        size = 3, fontface = "bold",
                        inherit.aes = FALSE # do not inherit anything from the ...
                        ) +
                theme_minimal()

return(pp)
}

Now use in ggpairs (adjust correlation method to be consistent with defaults in GGscatterPlot):

GGally::ggpairs(log2(df+1),
                1:4,
                lower = list(continuous = GGscatterPlot),
                upper = list(continuous = wrap("cor", method= "spearman")))

Now with Pearson correlation coefficient (and removing the correlation coefficient in the upper triangle since we already have them in the lower triangle):

GGally::ggpairs(log2(df+1),
                1:4,
                lower = list(continuous = wrap(GGscatterPlot, method="pearson")),
                upper = "blank") +
    theme_minimal() +
    theme(panel.grid.minor = element_blank())

We can use the upper triangle to plot other info since the scatter plots include the correlation coefficient.
I illustrate it here with the correlation coefficient for genes with >15 exons vs those with <15 exons:

exonNumber <- elementNROWS(rowRanges(airway[rownames(df),]))
df$MoreThan15Exons <- ifelse(exonNumber>15,
                             ">15ex", "<15ex")
df[,1:4] <- log2(df[,1:4]+1)

GGally::ggpairs(df,
                1:4,
        lower = list(continuous = wrap(GGscatterPlot, method="pearson")),
        upper = list(continuous = wrap(ggally_cor, alignPercent = 0.8), 
                     mapping = ggplot2::aes(color = MoreThan15Exons))) +
    theme_minimal() +
    theme(panel.grid.minor = element_blank())

The GGscatterPlot function can also be used with GGally::ggduo to study the link between 2 sets of columns/samples. Here an example comparing control vs DEX samples:

GGally::ggduo(df,
              c(1,3), 
              c(2,4),
              types = list(continuous = GGscatterPlot)) +
    theme_minimal() +
    theme(panel.grid.minor = element_blank())

Related