R | ggpairs -- 可视化相关性

最近想要可视化样本间的相关性，但又不满足于常规的相关性热图。因此，就注意到GGally包中的ggpairs函数，可以方便地实现多方面的相关性可视化。

本文仅介绍ggpairs 在连续型变量方面的应用。它也可以用到离散型变量的可视化上。

下面以airway数据集进行演示:
这里我们在前4个样本中随机选取1000个基因进行展示

library(GGally)
# airway example
library(airway)
data(airway)
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 gene
# ggpairs default
ggpairs(log2(df+1))

ggpairs将输出的图划分为三个区域，分别是左下角的lower, 对角线的diag, 以及右上角的upper. 对于连续性数值变量，默认在lower区画pairwise scatter plot，diag区画density plot，upper区展示相应的pairwise Pearson's correaltion coefficient.

进一步，我还希望在左下角的散点图中加入y=x的拟合线，并在对角线的加上直方图。我们可以通过自定义画图的函数实现这些操作。

ggscatter <- function(data, mapping, ...) {
  x <- GGally::eval_data_col(data, mapping$x)
  y <- GGally::eval_data_col(data, mapping$y)
  df <- data.frame(x = x, y = y)
  sp1 <- ggplot(df, aes(x=x, y=y)) +
    geom_point() +
    geom_abline(intercept = 0, slope = 1, col = 'darkred')
  return(sp1)
}

ggdehist <- function(data, mapping, ...) {
  x <- GGally::eval_data_col(data, mapping$x)
  df <- data.frame(x = x)
  dh1 <- ggplot(df, aes(x=x)) +
    geom_histogram(aes(y=..density..), bins = 50, fill = 'steelblue', color='black', alpha=.4) +
    geom_density(aes(y=..density..)) + 
    theme_minimal()
  return(dh1)
}

ggpairs(log2(df+1),
        lower = list(continuous = wrap(ggscatter)),
        diag = list(continuous = wrap(ggdehist))) + 
  theme_minimal() +
  theme(panel.grid = element_blank(),
        panel.border = element_rect(fill=NA),
        axis.text =  element_text(color='black'))

再放一个高度修改的版本

# https://pascal-martin.netlify.app/post/nicer-scatterplot-in-gggally/
GGscatterPlot <- function(data, mapping, ..., 
                          method = "pearson") {
  
  #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, use="pairwise.complete.obs")
  #Assemble data frame
  df <- data.frame(x = x, y = y)
  df <- na.omit(df)
  # 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, alpha = 1/density, color = cols)) +
    ggplot2::geom_point(shape=16, show.legend = FALSE) +
    ggplot2::scale_color_viridis_c(direction = colDirection) +
    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_bw()
  return(pp)
}

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, align_percent = 0.8), 
                             mapping = ggplot2::aes(color = MoreThan15Exons))) +
  theme_minimal() +
  theme(panel.grid = element_blank(),
        panel.border = element_rect(fill=NA),
        axis.text =  element_text(color='black'))

• 对散点图的透明度进行调整，点越密的区域透明度越高，反之亦然。
• 散点图的颜色代表了基因表达量。
• 在散点图左上角添加Pearson's 相关系数
• 右上角展示了对小于15个外显子的基因和大于15个外显子基因的相关性的统计。

在我看来ggpairs相当于是一个ggplot2的集成可视化方法，可以很方便的一次性展示多个方面的相关性信息。同时，它的可定制性也很高，可以满足许多额外的可视化需求。唯一的缺陷可能是需要耗费一定功夫写出包装的函数。

ref
ggpairs doc: https://ggobi.github.io/ggally/articles/ggpairs.html
Nicer scatter plots in GGgally ggpairs-ggduo: https://pascal-martin.netlify.app/post/nicer-scatterplot-in-gggally
Creating a density histogram in ggplot2: https://stackoverflow.com/questions/21061653/creating-a-density-histogram-in-ggplot2