Nature图片复现 | 彩虹热图

参考文献：

https://www.nature.com/articles/s41586-025-09686-5#MOESM1

数据来源

加载数据

library(tidyverse)
library(stats)
library(RColorBrewer)
library(paletteer)
library(ComplexHeatmap)
library(circlize)
library(grid)
df<- readxl::read_excel("41586_2025_9686_MOESM8_ESM.xlsx", sheet ="Fig4b")
mat <- df |>select(starts_with("response_Score")) |> as.matrix()
row_annot <- df |>select(-starts_with("response_Score"))

绘制圈图

# ###########################################################################
# 前置检查+数据清洗（过滤NA，确保数据合法性）
# ###########################################################################
if (!exists("mat") || !is.matrix(mat) || ncol(mat) != 4) {
  stop("请先定义4列数值矩阵 'mat'（核心热图数据）")
}
if (!exists("row_annot") || !is.data.frame(row_annot)) {
  stop("请先定义数据框 'row_annot'（包含样本注释信息）")
}

# 必需列检查
required_cols <- c(
  "Flu_Year", "cohort.cohortGuid", "subject.biologicalSex", "CMV",
  "Plasma Cell Freq Changes", "Response_Group-Flu B/Phuket HA",
  "Response_Group-Flu B/Washington HA", "Response_Group-A/Guangdong",
  "Response_Group-A/Cambodia"
)
missing_cols <- setdiff(required_cols, colnames(row_annot))
if (length(missing_cols) > 0) {
  stop(paste("row_annot 缺少必需列：", paste(missing_cols, collapse = ", ")))
}

# 过滤Flu_Year为NA的样本（避免扇区数量异常）
valid_idx <- !is.na(row_annot$Flu_Year)
row_annot <- row_annot[valid_idx, , drop = FALSE]
mat <- mat[valid_idx, , drop = FALSE]
if (nrow(row_annot) == 0) stop("过滤NA后无有效样本！")

# 显式定义扇区列表（统一所有绘图函数的扇区来源）
sectors <- unique(row_annot$Flu_Year)
n_sectors <- length(sectors)
message(paste("实际扇区数量：", n_sectors))

# ###########################################################################
# 1. 初始化PDF输出设备（指定图形参数，避免默认值冲突）
# ###########################################################################
pdf("Figure 4B.pdf", width = 10, height = 8, useDingbats = FALSE)

# ###########################################################################
# 2. 设置circos参数（确保间隙数量与扇区完全匹配）
# ###########################################################################
circos.clear()

# 定义合法间隙（非负数值，长度=扇区数量）
gap_after <- if (n_sectors == 1) {
  c(180)
} else {
  c(rep(3.5, n_sectors - 1), 180)
}
message(paste("间隙设置：", paste(gap_after, collapse = ", "), "度"))

# 显式设置全局参数，避免继承异常
circos.par(
  track.height = 0.1,
  start.degree = 180,
  gap.after = gap_after,
  canvas.xlim = c(-1.2, 1.2),
  canvas.ylim = c(-1.2, 1.2)
)

# ###########################################################################
# 3. 颜色映射函数（无参数问题，保留）
# ###########################################################################
create_color_fun <- function(col_idx) {
  col_data <- mat[, col_idx, drop = TRUE]
  range_min <- median(col_data, na.rm = TRUE) - 2 * mad(col_data, na.rm = TRUE)
  range_max <- median(col_data, na.rm = TRUE) + 2 * mad(col_data, na.rm = TRUE)
  colorRamp2(
    seq(range_min, range_max, length.out = 11),
    rev(brewer.pal(11, "Spectral"))
  )
}
col_fun1 <- create_color_fun(1)
col_fun2 <- create_color_fun(2)
col_fun3 <- create_color_fun(3)
col_fun4 <- create_color_fun(4)

# ###########################################################################
# 4. 注释热图函数（规范cell.lwd参数，确保合法）
# ###########################################################################
draw_annotation_heatmap <- function(annot_col, annot_colors, track_h = 0.02, cell_lwd = 0.2) {
  # 检查cell.lwd合法性（强制转为非负数值）
  if (!is.numeric(cell_lwd) || cell_lwd < 0) {
    warning(paste("非法cell.lwd值：", cell_lwd, "，已自动修正为0.1"))
    cell_lwd <- 0.1
  }
  
  annot_vec <- row_annot[[annot_col]]
  if (all(is.na(annot_vec))) stop(paste(annot_col, "列全为NA！"))
  
  # 处理注释因子
  annot_factor <- factor(annot_vec, levels = unique(annot_vec[!is.na(annot_vec)]))
  valid_levels <- intersect(levels(annot_factor), names(annot_colors))
  if (length(valid_levels) == 0) stop(paste("No matching levels for", annot_col))
  
  annot_factor <- factor(annot_factor, levels = valid_levels)
  anno_mat <- as.matrix(as.numeric(annot_factor))
  col_fun_anno <- colorRamp2(seq_along(valid_levels), annot_colors[valid_levels])
  
  # 合法参数绘制热图
  circos.heatmap(
    anno_mat,
    col = col_fun_anno,
    cluster = FALSE,
    cell.border = "white",
    cell.lwd = cell_lwd,
    track.height = track_h,
    split = factor(row_annot$Flu_Year, levels = sectors)
  )
}

# ###########################################################################
# 5. 绘制注释轨道（使用合法cell.lwd值）
# ###########################################################################
set_track_gap(cm_h(0.1))

# 5.1 Flu_Year注释：cell.lwd=0.01（避免0值）
flu_colors <- if (n_sectors < 3) {
  setNames(brewer.pal(3, "Set2")[1:n_sectors], sectors)
} else {
  setNames(brewer.pal(n_sectors, "Set2"), sectors)
}
draw_annotation_heatmap("Flu_Year", flu_colors, track_h = 0.01, cell_lwd = 0.01)

# 5.2 其他注释：统一cell.lwd=0.2
draw_annotation_heatmap(
  annot_col = "cohort.cohortGuid",
  annot_colors = c(BR1 = "#35978f", BR2 = "#bf812d"),
  track_h = 0.02,
  cell_lwd = 0.2
)

draw_annotation_heatmap(
  annot_col = "subject.biologicalSex",
  annot_colors = c(Male = "#5aae61", Female = "#9970ab"),
  track_h = 0.02,
  cell_lwd = 0.2
)

draw_annotation_heatmap(
  annot_col = "CMV",
  annot_colors = c(Positive = "#d6604d", Negative = "#4393c3"),
  track_h = 0.02,
  cell_lwd = 0.2
)

# ###########################################################################
# 6. 终极修复：circos.barplot按位置传参（移除所有命名参数）
# ###########################################################################
set_track_gap(cm_h(0.1))
circos.trackPlotRegion(
  factors = factor(row_annot$Flu_Year, levels = sectors),
  ylim = c(0, max(row_annot$`Plasma Cell Freq Changes`, na.rm = TRUE)),
  track.height = 0.06,
  panel.fun = function(x, y) {
    sector.name <- get.cell.meta.data("sector.index")
    idx <- row_annot$Flu_Year == sector.name
    y_val <- na.omit(row_annot$`Plasma Cell Freq Changes`[idx])
    if (length(y_val) == 0) return()
    x_val <- seq_len(length(y_val)) - 0.5
    
    # 核心修复：按位置传递y_val和x_val，无任何命名参数
    circos.barplot(
      y_val,  # 第1个位置参数：条形高度（无需y=）
      x_val,  # 第2个位置参数：条形位置（无需x=）
      col = "black",
      border = "black"
    )
  }
)

# ###########################################################################
# 7. 绘制响应组+核心热图（参数均合法）
# ###########################################################################
set_track_gap(cm_h(0.3))
response_colors <- c(
  `low responder` = "#bdd7e7",
  `middle responder` = "#6baed6",
  `high responder` = "#08519c"
)

# 7.1 Flu B/Phuket HA
draw_annotation_heatmap(
  annot_col = "Response_Group-Flu B/Phuket HA",
  annot_colors = response_colors,
  track_h = 0.02,
  cell_lwd = 0.3
)
set_track_gap(cm_h(0.1))
circos.heatmap(
  mat[, 1],
  col = col_fun1,
  cluster = FALSE,
  split = factor(row_annot$Flu_Year, levels = sectors),
  track.height = 0.06,
  cell.border = NA
)

# 7.2 Flu B/Washington HA
set_track_gap(cm_h(0.2))
draw_annotation_heatmap(
  annot_col = "Response_Group-Flu B/Washington HA",
  annot_colors = response_colors,
  track_h = 0.02,
  cell_lwd = 0.3
)
set_track_gap(cm_h(0.1))
circos.heatmap(
  mat[, 2],
  col = col_fun2,
  cluster = FALSE,
  split = factor(row_annot$Flu_Year, levels = sectors),
  track.height = 0.06,
  cell.border = NA
)

# 7.3 A/Guangdong
set_track_gap(cm_h(0.2))
draw_annotation_heatmap(
  annot_col = "Response_Group-A/Guangdong",
  annot_colors = response_colors,
  track_h = 0.02,
  cell_lwd = 0.3
)
set_track_gap(cm_h(0.1))
circos.heatmap(
  mat[, 3],
  col = col_fun3,
  cluster = FALSE,
  split = factor(row_annot$Flu_Year, levels = sectors),
  track.height = 0.06,
  cell.border = NA
)

# 7.4 A/Cambodia
set_track_gap(cm_h(0.2))
draw_annotation_heatmap(
  annot_col = "Response_Group-A/Cambodia",
  annot_colors = response_colors,
  track_h = 0.02,
  cell_lwd = 0.3
)
set_track_gap(cm_h(0.1))
circos.heatmap(
  mat[, 4],
  col = col_fun4,
  cluster = FALSE,
  split = factor(row_annot$Flu_Year, levels = sectors),
  track.height = 0.06,
  cell.border = NA
)

# ###########################################################################
# 8. 清理资源
# ###########################################################################
set_track_gap(cm_h(0.2))
circos.clear()
dev.off()

添加图例

# ###########################################################################
# 注意：请确保已提前定义 row_annot 数据框和 col_fun1 颜色映射函数
# （若未定义，可参考下方“快速测试数据”部分，取消注释后运行）
# ###########################################################################

# ---------------------- 核心：图例绘制主代码（独立运行）----------------------
# 1. 队列名称替换（BR1→Young，BR2→Older）
row_annot <- row_annot |>
  mutate(cohort.cohortGuid = case_match(
    cohort.cohortGuid,
    "BR1" ~ "Young",
    "BR2" ~ "Older"
  ))

# 2. 定义所有图例颜色映射（完整独立定义，不依赖外部变量）
# 2.1 Flu_Year颜色
annotation_colors_Flu_Year <- if (length(unique(row_annot$Flu_Year)) < 3) {
  setNames(brewer.pal(3, "Set2")[1:length(unique(row_annot$Flu_Year))], unique(row_annot$Flu_Year))
} else {
  setNames(brewer.pal(length(unique(row_annot$Flu_Year)), "Set2"), unique(row_annot$Flu_Year))
}

# 2.2 年龄组颜色（Young/Older）
annotation_colors_cohort <- c(Young = "#35978f", Older = "#bf812d")[levels(factor(row_annot$cohort.cohortGuid))]

# 2.3 性别颜色
annotation_colors_sex <- c(Male = "#5aae61", Female = "#9970ab")[levels(factor(row_annot$subject.biologicalSex))]

# 2.4 CMV状态颜色
annotation_colors_CMV <- c(Positive = "#d6604d", Negative = "#4393c3")[levels(factor(row_annot$CMV))]

# 2.5 响应组颜色（关键：补充用户代码中缺失的 annotation_colors 定义）
annotation_colors <- c(
  `low responder` = "#bdd7e7",
  `middle responder` = "#6baed6",
  `high responder` = "#08519c"
)

# 2.6 各病毒株响应组颜色（过滤有效水平）
annotation_colors_fluPhuket <- annotation_colors[intersect(
  unique(row_annot$`Response_Group-Flu B/Phuket HA`),
  names(annotation_colors)
)]

# 3. 创建单个图例
lgd_response_score <- Legend(
  title = "Response score (HAI)",
  col_fun = col_fun1,
  direction = "horizontal"
)

lgd_Flu_Year <- Legend(
  title = "Flu_Year",
  at = names(annotation_colors_Flu_Year),
  legend_gp = gpar(fill = annotation_colors_Flu_Year)
)

lgd_cohort <- Legend(
  title = "Age group",
  at = names(annotation_colors_cohort),
  legend_gp = gpar(fill = annotation_colors_cohort)
)

lgd_sex <- Legend(
  title = "Sex",
  at = names(annotation_colors_sex),
  legend_gp = gpar(fill = annotation_colors_sex)
)

lgd_CMV <- Legend(
  title = "CMV",
  at = names(annotation_colors_CMV),
  legend_gp = gpar(fill = annotation_colors_CMV)
)

lgd_responder_groups <- Legend(
  title = "Responder groups",
  at = names(annotation_colors_fluPhuket),
  legend_gp = gpar(fill = annotation_colors_fluPhuket)
)

# 4. 打包图例（横向排列，按用户要求位置绘制）
all_lgd <- packLegend(
  lgd_response_score, lgd_Flu_Year, lgd_cohort, lgd_sex, lgd_CMV, lgd_responder_groups,
  direction = "horizontal",
  gap = unit(5, "mm")  # 图例间距，可调整
)

# 5. 初始化PDF设备（单独生成图例文件）
pdf("Figure_4B_Legend.pdf", width = 15, height = 3, useDingbats = FALSE)  # 宽画布适配横向图例

# 6. 绘制图例（按用户指定位置：画布中心偏下）
draw(
  all_lgd,
  x = unit(0.5, "npc"),
  y = unit(0.4, "npc"),
  just = "center"
)

# 7. 关闭设备，保存文件
dev.off()

合并图片