作者：老王
链接：https://laowang2023.cn/2023/03/22/20230322-pheCorDist/
来源：请叫我老王（个人博客）
本文采用 CC BY 4.0 许可协议 进行共享。转载请注明作者姓名和链接
考察一个群体的多个表型或者一个表型的多个重复，我们想展示其分布和他们之间的相关关系可以使用柱状图和散点图（如下图所示）。
这幅图主要有两部分组成，一个是对角线上的柱状图，使用柱状图展示了每一个表型重复的分布；另一个就是对角线下面的散点图，用散点图展示两两之间的相关关系，并且用不同颜色表示点的密度，在上面标注其相关性。下面我们将使用R语言完成这幅图。
对于这幅图我们可以先分别绘制其中每一个部分，然后使用图片组合、拼接函数进行整合：
分图绘制
首先导入数据，数据格式如下，每一行代表一个样本，每一列代表一个重复：
>data <- read.table("./data.txt", header = T, row.names = 1, sep = "\t")
>data
pheno16rep1 pheno16rep2 pheno17 pheno18rep1 pheno18rep2
s1 28.4 24.9 27.740 27.72500 29.30000
s2 26.6 25.3 NA NA NA
s3 27.8 27.0 24.660 NA 27.97500
s4 25.5 26.9 22.680 29.27500 25.95000
s5 26.5 28.7 24.760 31.97500 27.52000
s6 .... .... ....... ........ ........
使用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) {
#
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 = 1) +
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()
