R语言绘图那些事（四）

缘由

生信之路，道阻且长，记录经验，但渡有缘人。

碎碎念

及时当勉励，岁月不饶人 ——陶渊明

原图 & 文章

Zhou, J., Hu, F., Berhe, M. et al. Genome-wide identification, classification, and expression profiling of LAC gene family in sesame. BMC Plant Biol 24, 1254 (2024).

功能：同时展示了聚类tree和表达数据数据

炫图

代码

1. 加载必要的库

library(tidyverse)    # 数据处理和可视化
library(reshape2)     # 数据框格式转换
library(ggsci)        # 颜色主题
library(gground)      # 自定义可视化工具
library(patchwork)    # 图形拼接
library(ape)          # 树状图生成和处理

---

2. 创建模拟数据

set.seed(42)  # 设置随机数种子，保证结果可重复
gene_names <- paste0("Gene", 1:30)  # 模拟基因名称
samples1 <- paste0("CK_Sample", 1:20)  # 对照组样本名称
samples2 <- paste0("TR_Sample", 1:20)  # 处理组样本名称
fpkm1_values <- matrix(runif(30 * 40, min = 0, max = 100), nrow = 30, ncol = 40)

• 随机生成 30 个基因在 40 个样本中的表达量数据，范围为 [0, 100]。
• samples1 和 samples2 分别表示对照组和处理组的样本。

---

3. 数据框构造与预处理

data1 <- data.frame(Gene_Name = gene_names, fpkm1_values)
colnames(data1)[2:41] <- c(samples1, samples2)  # 设置列名
selected_columns1 <- data1[, c(1, 2:41)]  # 选择基因表达相关列
rownames(selected_columns1) <- selected_columns1[, 1]  # 设置行名为基因名称
selected_columns1 <- selected_columns1[, 2:41]  # 移除基因名称列
log_data1 <- log2(selected_columns1 + 1)  # log2 转换，防止零值问题

• 将矩阵数据转为数据框格式，并对列名和行名进行设置。
• 数据进行 log2(FPKM+1) 转换，压缩数据范围并避免取对数时的零值问题。

---

4. 计算距离矩阵与层次聚类

dist_matrix <- dist(log_data1, method = "euclidean")  # 计算基因表达的欧几里得距离
hc <- hclust(dist_matrix, method = "ward.D2")         # 使用 Ward.D2 方法聚类
phylo_tree <- ape::as.phylo(hc)                       # 转换为树状结构

• 通过欧几里得距离衡量基因表达模式的差异。
• ward.D2 是一种层次聚类方法，适合构建树状关系。
• 最终结果转换为树状图对象 phylo_tree。

---

5. 热图输入数据处理

log_data_ordered <- log_data1[hc$order, ]  # 按聚类结果对基因排序
draw1 <- melt(as.matrix(log_data_ordered))  # 转换为长格式，适合 ggplot

• 使用聚类结果的顺序对基因重新排序。
• 使用 melt 将矩阵数据转为长格式，便于热图绘制。

---

6. 绘制树状图

tree_plot <- ggtree(phylo_tree, branch.length = "none")

• ggtree 绘制聚类结果的树状图，设置 branch.length = "none" 使分支长度一致。

draw1$Var1 <- factor(draw1$Var1,
                     level = tree_plot$data %>% arrange(y) %>% na.omit() %>% pull(label) %>% unique())

• 根据树状图数据的顺序调整热图中基因名称的顺序，确保热图和树状图的基因对齐。

---

**7. 画图与拼图

p <-  ggplot(draw1 )+
    geom_round_rect(data=data.frame(xmin=0,xmax=20.4,ymin=0, ymax=12.3),
                    aes(xmin=xmin,xmax=xmax, ymin=ymin, ymax=ymax),radius = unit(0.1,"npc"),fill="#ecf7dd")+
    geom_round_rect(data=data.frame(xmin=0,xmax=20.4,ymin=12.5, ymax=15.3),
                    aes(xmin=xmin,xmax=xmax, ymin=ymin, ymax=ymax),radius = unit(0.2,"npc"),fill="#fae8f3")+
    geom_round_rect(data=data.frame(xmin=0,xmax=20.4,ymin=15.5, ymax=18.3),
                    aes(xmin=xmin,xmax=xmax, ymin=ymin, ymax=ymax),radius = unit(0.2,"npc"),fill="#e0f3ec")+
    geom_round_rect(data=data.frame(xmin=0, xmax=20.4,ymin=18.5, ymax=Inf),
                    aes(xmin=xmin,xmax=xmax, ymin=ymin, ymax=ymax),radius = unit(0.1,"npc"),fill="#e9ecf8")+
    geom_round_rect(data=data.frame(xmin=20.6,xmax=41,ymin=0, ymax=12.3),
                    aes(xmin=xmin,xmax=xmax, ymin=ymin, ymax=ymax),radius = unit(0.1,"npc"),fill="#ecf7dd")+
    geom_round_rect(data=data.frame(xmin=20.6,xmax=41,ymin=12.5, ymax=15.3),
                    aes(xmin=xmin,xmax=xmax, ymin=ymin, ymax=ymax),radius = unit(0.2,"npc"),fill="#fae8f3")+
    geom_round_rect(data=data.frame(xmin=20.6,xmax=41,ymin=15.5, ymax=18.3),
                    aes(xmin=xmin,xmax=xmax, ymin=ymin, ymax=ymax),radius = unit(0.2,"npc"),fill="#e0f3ec")+
    geom_round_rect(data=data.frame(xmin=20.6, xmax=41,ymin=18.5, ymax=Inf),
                    aes(xmin=xmin,xmax=xmax, ymin=ymin, ymax=ymax),radius = unit(0.1,"npc"),fill="#e9ecf8")+
    geom_point(data=draw1,aes(Var2, Var1, fill = value, size = value), stroke =1, shape=21, color="black",inherit.aes = F,alpha=0.75)
 
p <- p + scale_size(c(1,6))+scale_x_discrete(position="top")+
    scale_fill_gradient2(
        low = "white", mid = "#e58537", high = "red",
        midpoint = median(draw1$value),  # 中间值的定义
        na.value = "grey"  # NA 值使用灰色
    ) + theme_minimal() + # 调整圆的大小范围
    labs(x = NULL, y = NULL, fill = "Expression\n(log2(FPKM+1))") +
    scale_y_discrete(position = "right") +  # 将基因名放右边
    guides(
        fill = guide_legend(title = "(log2(FPKM+1))", override.aes = list(size = 5)), 
        size = guide_legend(title = "(log2(FPKM+1))", override.aes = list(shape = 21, color = "black"))  
    )+
    theme(axis.text.x=element_text(angle=90),
          panel.grid = element_blank()
          )
          
  tree_plot+p+plot_layout(ncol=2,widths = c(1,2))   

ALL

```r
#devtools::install_github("dxsbiocc/gground")
library(tidyverse)
library(reshape2)
library(ggsci)
library(gground)
library(patchwork)
library(ape)


# 创建一个模拟数据框
set.seed(42)
gene_names <- paste0("Gene", 1:30) # 30个基因名称
samples1 <- paste0("CK_Sample", 1:20) # 20个样本
samples2 <- paste0("TR_Sample", 1:20) # 20个样本
fpkm1_values <- matrix(runif(30 * 40, min = 0, max = 100), nrow = 30, ncol = 40)

# 构建数据框，并添加基因名称列
data1 <- data.frame(Gene_Name = gene_names, fpkm1_values)
colnames(data1)[2:41] <- c(samples1,samples2)
selected_columns1<- data1[, c(1, 2:41)]
rownames(selected_columns1) <- selected_columns1[, 1]
selected_columns1 <- selected_columns1[, 2:41]
log_data1 <- log2(selected_columns1 + 1)

# 计算距离矩阵
dist_matrix <- dist(log_data1, method = "euclidean")
# 进行层次聚类
hc <- hclust(dist_matrix, method = "ward.D2")

# 转换聚类结果为树结构
phylo_tree <- ape::as.phylo(hc)

# 绘制热图数据的矩阵，按聚类顺序排列
log_data_ordered <- log_data1[hc$order, ]
draw1 <- melt(as.matrix(log_data_ordered))

#tree
tree_plot <- ggtree(phylo_tree, branch.length = "none",)#+geom_tiplab(size=3, color="purple")

draw1$Var1 <- factor(draw1$Var1,
                     level=tree_plot$data %>% arrange(y)
                     %>%na.omit() %>%  pull(label) %>% unique()
                     )
#geom_point 
p <-  ggplot(draw1 )+
    geom_round_rect(data=data.frame(xmin=0,xmax=20.4,ymin=0, ymax=12.3),
                    aes(xmin=xmin,xmax=xmax, ymin=ymin, ymax=ymax),radius = unit(0.1,"npc"),fill="#ecf7dd")+
    geom_round_rect(data=data.frame(xmin=0,xmax=20.4,ymin=12.5, ymax=15.3),
                    aes(xmin=xmin,xmax=xmax, ymin=ymin, ymax=ymax),radius = unit(0.2,"npc"),fill="#fae8f3")+
    geom_round_rect(data=data.frame(xmin=0,xmax=20.4,ymin=15.5, ymax=18.3),
                    aes(xmin=xmin,xmax=xmax, ymin=ymin, ymax=ymax),radius = unit(0.2,"npc"),fill="#e0f3ec")+
    geom_round_rect(data=data.frame(xmin=0, xmax=20.4,ymin=18.5, ymax=Inf),
                    aes(xmin=xmin,xmax=xmax, ymin=ymin, ymax=ymax),radius = unit(0.1,"npc"),fill="#e9ecf8")+
    geom_round_rect(data=data.frame(xmin=20.6,xmax=41,ymin=0, ymax=12.3),
                    aes(xmin=xmin,xmax=xmax, ymin=ymin, ymax=ymax),radius = unit(0.1,"npc"),fill="#ecf7dd")+
    geom_round_rect(data=data.frame(xmin=20.6,xmax=41,ymin=12.5, ymax=15.3),
                    aes(xmin=xmin,xmax=xmax, ymin=ymin, ymax=ymax),radius = unit(0.2,"npc"),fill="#fae8f3")+
    geom_round_rect(data=data.frame(xmin=20.6,xmax=41,ymin=15.5, ymax=18.3),
                    aes(xmin=xmin,xmax=xmax, ymin=ymin, ymax=ymax),radius = unit(0.2,"npc"),fill="#e0f3ec")+
    geom_round_rect(data=data.frame(xmin=20.6, xmax=41,ymin=18.5, ymax=Inf),
                    aes(xmin=xmin,xmax=xmax, ymin=ymin, ymax=ymax),radius = unit(0.1,"npc"),fill="#e9ecf8")+
    geom_point(data=draw1,aes(Var2, Var1, fill = value, size = value), stroke =1, shape=21, color="black",inherit.aes = F,alpha=0.75)

p <- p + scale_size(c(1,6))+scale_x_discrete(position="top")+
    scale_fill_gradient2(
        low = "white", mid = "#e58537", high = "red",
        midpoint = median(draw1$value),  # 中间值的定义
        na.value = "grey"  # NA 值使用灰色
    ) + theme_minimal() + # 调整圆的大小范围
    labs(x = NULL, y = NULL, fill = "Expression\n(log2(FPKM+1))") +
    scale_y_discrete(position = "right") +  # 将基因名放右边
    guides(
        fill = guide_legend(title = "(log2(FPKM+1))", override.aes = list(size = 5)), 
        size = guide_legend(title = "(log2(FPKM+1))", override.aes = list(shape = 21, color = "black"))  
    )+
    theme(axis.text.x=element_text(angle=90),
          panel.grid = element_blank()
          )
          
  tree_plot+p+plot_layout(ncol=2,widths = c(1,2))