R绘图 | 群体表型分布和相关关系

考察一个群体的多个表型或者一个表型的多个重复，我们想展示其分布和他们之间的相关关系可以使用柱状图和散点图（如下图所示）。

这幅图主要有两部分组成，一个是对角线上的柱状图，使用柱状图展示了每一个表型重复的分布；另一个就是对角线下面的散点图，用散点图展示两两之间的相关关系，并且用不同颜色表示点的密度，在上面标注其相关性。下面我们将使用R语言完成这幅图。

对于这幅图我们可以先分别绘制其中每一个部分，然后使用图片组合、拼接函数进行整合：

分图绘制

首先导入数据，数据格式如下，每一行代表一个样本，每一列代表一个重复：

data <- read.table("./data.txt", header = T, row.names = 1, sep = "\t")
head(data)
#>    pheno16rep1 pheno16rep2 pheno17 pheno18rep1 pheno18rep2
#> s1        28.4        24.9   27.74      27.725    29.30000
#> s2        26.6        25.3      NA          NA          NA
#> s3        27.8        27.0   24.66          NA    27.97500
#> s4        25.5        26.9   22.68      29.275    25.95000
#> s5        26.5        28.7   24.76      31.975    27.52000
#> s6        28.9        27.5   23.20      28.500    32.03333

使用 ggplot2 扩展包绘制每一个分图。柱状图使用 geom_histogram() 绘制，散点图使用 ggpointdensity 包的 geom_pointdensity() 函数绘制，使用 cor() 函数计算两个重复之间的相关系数，并将其放在图片标题位置，并使用 ggtext 包的 element_markdown() 函数设置标题的主题，同时使用 cowplot 包的 theme_half_open() 函数设置整体主题。

library(tidyverse)
library(ggpointdensity)
library(cowplot)
library(ggtext)

# 柱状图
p1 <- ggplot(data, aes(x = pheno16rep1)) + 
        geom_histogram(binwidth = 1) +
        labs(x = NULL, y = NULL, title = cn[i]) +
        theme_half_open() +
        theme(plot.title = element_text(hjust = 0.5))
# 散点密度图
p2 <- ggplot(data, aes(x = pheno16rep2, y = pheno16rep1)) +
        geom_pointdensity() +
        scale_color_continuous(type = "viridis") +                # 设置点密度颜色梯度
        labs(x = NULL, y = NULL, title = paste("*R*: ", round(cor(df$col1, df$col2, use = "na.or.complete"), 2), sep = "")) +
        theme_half_open() +
        theme(legend.position = "NA", 
           plot.title = element_markdown(hjust = 0.5,
                           face = "plain"))

组合图片

使用 customLayout 包进行图片组合，这个包可以对 base 绘图和 ggplot2 绘图进行整合，而且比较灵活。首先需要 lay_new() 函数创建一个拼接画布，然后使用 lay_grid() 函数组合各个图片。因为总共有 5 个重复，因此需要一个 5×5 的画图，如下图所示，各个分图从左上角开始往下排列走”之“字形排列。

lay <- lay_new(mat = matrix(1:25, nrow = 5), widths = rep(1, 1), heights = rep(1, 1))
lay_show(lay)

现在出现了一个问题，我们并没有在对角线上方安排图片，而 lay_new() 产生的是一个矩形排列画布，因此我们需要在右上角填充空白图片，并将空白图和柱状图、散点密度图整合。

p <- ggplot() + theme_nothing()
lay_grid(list(p1, p2, p3, ...), lay)

整理以上过程

在一个 5×5 的组合中我们总共需要绘制 25 个分图，其中有多次重复的过程，并且最终图片是矩形有规律分布，因此为了减少代码长度我们可以使用循环来处理每个分图。根据 lay_new() 的组合形式可以设置两层循环分别处理行和列，并且因为组合图是从左上角开始向下排布，因此外层循环用来处理行，内层分布处理列。最后一点就是可以把这一系列代码写成一个 function，方便以后使用。

最终代码如下所示：

library(tidyverse)
library(ggpointdensity)
library(cowplot)
library(ggtext)
library(customLayout)

## 定义pheCorDist函数
pheCorDist <- function(data, binwidth = NULL, bins = 30) {
  #
  n <- ncol(data)
  cn <- colnames(data)
  Pall <- list()
  index <- 1
  #
  for (j in 1:n) {
    for (i in 1:n) {
      df <- data[, c(i,j)]
      colnames(df) <- c("col1", "col2")
      if (i == j) {
        p <- ggplot(df, aes(x = col1)) + 
          geom_histogram(binwidth = binwidth, bins = bins) +
          labs(x = NULL, y = NULL, title = cn[i]) +
          theme_half_open() +
          theme(plot.title = element_text(hjust = 0.5))
      } else if (i > j) {
        p <- ggplot(df, aes(x = col2, y = col1)) +
          geom_pointdensity() +
          scale_color_continuous(type = "viridis") +
          labs(x = NULL, y = NULL, title = paste("*R*: ", round(cor(df$col1, df$col2, use = "na.or.complete"), 2), sep = "")) +
          theme_half_open() +
          theme(legend.position = "NA", 
                plot.title = element_markdown(hjust = 0.5,
                                              face = "plain"))
      } else if(i < j) {
        p <- ggplot() + theme_nothing()
      }
      Pall[index][[1]] <- p
      index = index + 1
    }
  }
  lay <- lay_new(mat = matrix(1:n^2, nrow = n), widths = rep(1, n), heights = rep(1, n))
  lay_grid(Pall, lay)
}

# 导入数据并绘图
data <- read.table("./data.txt", header = T, row.names = 1, sep = "\t")
png(filename = "test.png", width = 10, height = 8, units = "in", res = 500)
pheCorDist(data)
dev.off()