前言

不是在帮人解决问题，就是在解决问题的路上~~~。

之前所介绍的热图，其每个颜色块都是一个矩形，而今天要介绍的是如何绘制对角线分割热图。也就是每个颜色块矩形被对角线分割为上下两个三角形，然后两个三角形分别用不同的变量来设置填充色。

这种图形重要用于展示行列变量配对值的不同维度信息，比如，对于相关性矩阵，上下两个三角形的填充色可以分别用来表示相关性大小和显著性。类似于下面这张图

看到这张图，第一反应便是可以使用 geom_polygon 函数来分别绘制上三角和下三角，两个图层叠加便可以实现这种效果。而其中的点的数量表示的是显著性大小，可以使用点图来实现。

实现细节

现在已经有思路了，重点是如何将配对变量值转换成坐标信息。

首先，让我们来看看，如何使用 geom_polygon 绘制一个上三角和下三角

library(tidyverse)

upper <- data.frame(
  x = c(0,0,1),
  y = c(0,1,1)
)

lower <- data.frame(
  x = c(0,1,1),
  y = c(0,0,1)
)

ggplot(upper, aes(x, y)) +
  geom_polygon(fill = "red") +
  geom_polygon(data = lower, fill = "blue")

上下三角形之间的区别只是一个坐标点的不同而已，对角线上的两个点是重叠的。

这个图形只是一个配对变量值的形状，这些坐标点属于同一个分组，我们需要指定 group 变量来进行区分

那如何扩展到所有变量对呢？我们只需将每个坐标进行横向和纵向平移即可扩展到整个矩阵。

假设有个变量的取值如下

var1 <- 1:3
var2 <- 4:6

那么它们的组合为

> pairs <- merge(var1, var2)
> pairs
  x y
1 1 4
2 2 4
3 3 4
4 1 5
5 2 5
6 3 5
7 1 6
8 2 6
9 3 6

而每个组合的值便是我们需要的平移量，我们可以对行应用函数生成一个上三角形矩阵

df <- do.call(rbind,
        apply(pairs, 1, function (x) {
          a = x[1]
          b = x[2]
          data.frame(
            x = c(0, 0, 1) + a,
            y = c(0, 1, 1) + b,
            group = paste(a, b, sep = "-")
          )
        }))

ggplot(df, aes(x, y, group = group)) +
  geom_polygon(fill = "red")

现在，我们可以读入准备好的相关分析的数据

data <- read_delim('Downloads/gene_sig.txt')

数据中，每行代表一个组合，基因与免疫细胞之间的相关系数（cor）及显著性（p）

> data
# A tibble: 140 × 4
   gene     cell                               p    cor
   <chr>    <chr>                          <dbl>  <dbl>
 1 SIGLEC16 Plasma cells               0.0304    -0.146
 2 SIGLEC16 T cells CD8                0.0000880  0.261
 3 SIGLEC16 T cells follicular helper  0.000250   0.244
 4 SIGLEC16 T cells regulatory (Tregs) 0.000183   0.249
 5 SIGLEC16 Macrophages M0             0.00763   -0.179
 6 SIGLEC16 Macrophages M1             0.000108   0.258
 7 SIGLEC16 Macrophages M2             0.0000851  0.261
 8 SIGLEC16 Dendritic cells activated  0.000596  -0.229
 9 SIGLEC16 Mast cells activated       0.00235   -0.204
10 SIGLEC16 Neutrophils                0.00153   -0.212
# … with 130 more rows

为了方便将字符转换为对应的数值，我们将前两列转换为 factor

data <- mutate_at(data, 1:2, ~ as.factor(.))

如果输入的是矩阵形式，即形如行为基因列为免疫细胞，值为相关系数，可以转换为这种形式

我们可以将提取上三角和下三角的操作封装成函数，方便使用

# 根据配对列表生成上、下三角坐标
triangle <- function(pairs, type = "up") {
  # 默认的上三角坐标基
  x = c(0, 0, 1)
  y = c(0, 1, 1)
  # 下三角的坐标基
  if (type == "lower") {
    x = c(0, 1, 1)
    y = c(0, 0, 1)
  }
  # 生成三角矩阵
  mat = do.call(
    rbind,
    apply(pairs, 1, function (row) {
      a = row[1]
      b = row[2]
      data.frame(
        x = x + a,
        y = y + b,
        group = paste(a, b, sep = "-")
      )
    }))
  return(mat)
}

triangle_data <- function(data, row = 1, col = 2) {
  # 这里设置的 row 和 col 表示要指定的行列变量所在列
  # 生成所有组合
  rows = length(unique(data[[row]]))
  cols = length(unique(data[[col]]))
  pairs = merge(1:rows, 1:cols)
  # 获取上三角坐标
  upper <- triangle(pairs)
  colnames(upper) <- c(paste0("upper.", colnames(upper)[1:2]), "group")
  # 获取下三角坐标
  lower <- triangle(pairs, type = "lower")[1:2]
  colnames(lower) <- paste0("lower.", colnames(lower))
  # 合并坐标
  upper_lower = cbind(upper, lower)
  # 根据分组信息将坐标连接到数据中
  data %>% transmute(across(where(is.factor), ~ as.character(as.numeric(.)))) %>%
    unite("group", row:col, sep = "-") %>%
    cbind(data, .) %>%
    right_join(upper_lower, by = "group")
}

转换数据

> trian_data <- triangle_data(data)
> head(trian_data)
      gene         cell            p    cor group upper.x upper.y lower.x lower.y
1 SIGLEC16 Plasma cells 3.040666e-02 -0.146 14-10      14      10      14      10
2 SIGLEC16 Plasma cells 3.040666e-02 -0.146 14-10      14      11      15      10
3 SIGLEC16 Plasma cells 3.040666e-02 -0.146 14-10      15      11      15      11
4 SIGLEC16  T cells CD8 8.796373e-05  0.261 18-10      18      10      18      10
5 SIGLEC16  T cells CD8 8.796373e-05  0.261 18-10      18      11      19      10
6 SIGLEC16  T cells CD8 8.796373e-05  0.261 18-10      19      11      19      11

由于这份数据中包含 NA 值，即有些 gene 和 cell 组合被删掉了，所以在这里需要将 NA 值替换掉

df <- mutate(trian_data, 
             cor = replace_na(cor, 0),
             p = replace_na(p, 1)) 

最后，绘制图形

ggplot(df) +
  geom_polygon(aes(upper.x, upper.y, fill = abs(cor), group = group)) +
  geom_polygon(aes(lower.x, lower.y, fill = p, group = group))

虽然形状都是正确的，但是只有一个填充色，我们明明设置了两个填充色变量的。

其实，在 ggplot 中是不允许在一张图中对同一个 aes 参数的标度进行设置的，但是好在有人帮我们实现了这一功能

ggnewscale 包提供的 new_scale 函数可以允许我们设置多个颜色变量，也适用于其他 aes 变量，如 shape、linetype 等等，先安装包

install.packages("ggnewscale")

使用方式也很简单，只需添加到两个对象之间，可以看到出现了两个图例

library(ggnewscale)

ggplot(df) +
  geom_polygon(aes(upper.x, upper.y, fill = abs(cor), group = group)) +
  new_scale("fill") +
  geom_polygon(aes(lower.x, lower.y, fill = p, group = group))

配置一下好看的颜色

ggplot(df) +
  geom_polygon(aes(upper.x, upper.y, fill = abs(cor), group = group)) +
  # 相关性颜色
  scale_fill_gradientn(colors = colorRampPalette(c("#1E3163", "#00C1D4", "#FFED99","#FF7600"))(10)) +
  new_scale("fill") +
  # 显著性颜色
  geom_polygon(aes(lower.x, lower.y, fill = p, group = group)) +
  scale_fill_gradientn(colours = RColorBrewer::brewer.pal(5, "YlGnBu"))

颜色搭配好了之后，需要将标签添加上去

ggplot(df) +
  geom_polygon(aes(upper.x, upper.y, fill = abs(cor), group = group)) +
  scale_fill_gradientn(colors = colorRampPalette(c("#1E3163", "#00C1D4", "#FFED99","#FF7600"))(10)) +
  new_scale("fill") +
  geom_polygon(aes(lower.x, lower.y, fill = p, group = group)) +
  scale_fill_gradientn(colours = RColorBrewer::brewer.pal(5, "YlGnBu")) +
  scale_x_continuous(breaks = c(1:length(unique(data[[2]]))) + 0.5, expand = c(0,0),
                     labels = sort(unique(data[[2]]))) +
  scale_y_continuous(expand = c(0, 0), breaks = c(1:length(unique(data[[1]]))) + 0.5,
                     labels = sort(unique(data[[1]])), sec.axis = dup_axis()) +
  theme(
    plot.margin = margin(0.5,0.01,0.5,0.01, "cm"),
    axis.title = element_blank(),
    axis.text.y.left = element_blank(),
    axis.ticks.y.left = element_blank(),
    axis.text.x = element_text(angle = 270, hjust = 0, vjust = 0.5)
  )

添加灰色边框

ggplot(df) +
  geom_polygon(aes(upper.x, upper.y, fill = abs(cor), group = group), colour = "grey") +
  ...

好了，万事俱备，只欠点图了。

这里，我的想法是提取出之前画三角形时的起始点位置，并添加偏移到下三角的最右侧，而根据 p 值的不同程度，再添加数值方向上的偏移点，就可以了。

首先，提取起始位置

tmp <- data %>% transmute(across(where(is.factor), as.numeric)) %>%
  `names<-`(c("y", "x")) %>%
  cbind(data, .) %>%
  as.data.frame()

添加偏移点

points <- do.call(rbind, apply(tmp, 1, function(row) {
  p = as.numeric(row['p'])
  x = as.numeric(row['x'])
  y = as.numeric(row['y'])
  df = data.frame()
  if (p < 0.001) {
    df = rbind(df, data.frame(x = x + 0.9, y = y + 0.5))
  }
  if (p < 0.01) {
    df = rbind(df, data.frame(x = x + 0.9, y = y + 0.3))
  }
  if (p < 0.05) {
    df = rbind(df, data.frame(x = x + 0.9, y = y + 0.1))
  }
  df
}))

最后，使用 geom_point 将点添加到图形中

ggplot(trian_data) +
  geom_polygon(aes(upper.x, upper.y, fill = abs(cor), group = group), colour = "grey") +
  scale_fill_gradientn(colors = colorRampPalette(c("#1E3163", "#00C1D4", "#FFED99","#FF7600"))(10)) +
  new_scale("fill") +
  geom_polygon(aes(lower.x, lower.y, fill = p, group = group)) +
  scale_fill_gradientn(colours = RColorBrewer::brewer.pal(5, "YlGnBu")) +
  geom_point(data = points, aes(x, y), size = 0.4) +
  scale_x_continuous(breaks = c(1:length(unique(data[[2]]))) + 0.5, expand = c(0,0),
                     labels = sort(unique(data[[2]]))) +
  scale_y_continuous(expand = c(0, 0), breaks = c(1:length(unique(data[[1]]))) + 0.5,
                     labels = sort(unique(data[[1]])), sec.axis = dup_axis()) +
  theme(
    plot.margin = margin(0.5,0.01,0.5,0.01, "cm"),
    axis.title = element_blank(),
    axis.text.y.left = element_blank(),
    axis.ticks.y.left = element_blank(),
    axis.text.x = element_text(angle = 270, hjust = 0, vjust = 0.5)
  )

由于数据的问题，没有 NA 点的话图像会好看点。

代码和文件已上传到 GitHub：
https://github.com/dxsbiocc/learn/blob/main/R/plot/triangle_heatmap.R