摘要LLM 归库摘要
- 核心内容： 用 R 复现《Nature Cell Biology》2026 灵长类 CS9 胚胎空间转录组论文的 Fig. 2d——HOX 基因沿前后轴（A–P）表达分布的山脊图+热图组合，附完整 ggplot2/ggridges 代码。
- 适用场景： 展示一组基因沿空间轴/解剖轴的表达梯度；单细胞、空间转录组"山脊图+热图"发表级组合排版。
- 输入： 基因 × 空间层表达矩阵（31 HOX 基因 × 100 A–P 层），源自 Tan_2026_PrimateOrganogenesis_NatCellBiol_source_data_fig2.xlsx 的 Panel d。
- 输出： Fig2d_xinye_plot.png/pdf——HOXA/B/C/D 四组山脊图 + 右侧 z-score 行标准化热图，A→P 方向，含配色与图例。
- 关键参数/步骤： 每基因 min-max 归一化定山脊高度；smooth.spline(spar=0.45) 平滑；Cubic Hermite 尾巴延伸；热图 z-score 与蓝-白-红渐变（limits -4~4）；cowplot 拆行对齐山脊与热图。
- 常见错误： 直接用原始表达量会压扁图形；AP 第 100 层不是右边界需留白；图例/标题不拆行会导致两图错位。
- 与口腔生信连接： 可套用此模板展示牙周炎/OSCC 关键基因沿空间轴向（上皮-间质、炎症梯度）的表达分布，替换数据与基因列表即可。
- 下一步处理： 暂未试用；复现需先下载源数据 xlsx，确认 readxl 路径与包版本后试跑。
- 归库说明： 分类到 03 生信教程与工具库，置信度 0.88，理由：可视化复现教程、含完整 R 代码、非单篇论文剪藏。
新叶 生信技能树 2026年8月5日 19:56
今天学习2026年5月7日发表在《Nature Cell Biology》杂志上的一篇文献，标题为《A three-dimensional spatial transcriptome atlas reconstructs early organogenesis in primate Carnegie stages 9 and 10 embryos》。
使用AI复现这个文献中的Figure2D，是一个山脊图+热图的组合模式！
山脊图我们前面有过介绍： 顶刊Nature杂志同款配色一流的山峦图
原图如下：
example07-Fig2d
图的含义#
Fig. 2d 的核心是：把 HOX 基因表达沿胚胎前后轴（A–P，Anterior→Posterior）展开，观察其空间分布，而不是做 GO 富集。
图怎么读
- 横轴：胚胎前端 A → 后端 P。
- 每条 ridge：一个 HOX 基因在 A–P 轴上的表达曲线。
- 峰的位置：该基因主要在哪个轴向区域表达。
- 峰越靠右，说明越偏后部；峰越窄，说明空间限制越强。
- 右侧热图：对应基因在不同 A–P 位置的平均表达，蓝色低、红色高。
- 四组分别是 HOXA、HOXB、HOXC 和 HOXD 基因簇。
作者图注称其为：
“Ridgeline plots showing the distribution of HOX gene expression along the A–P axis…”
主要观察结果
总体呈现典型的 HOX 前后轴分区规律：
- 较前端的 HOX 基因，例如 HOXA1、HOXB1、HOXB2，表达范围相对更宽，覆盖较长的 A–P 区域。
- 随着 HOX 基因位置偏向中部和后部，表达峰逐渐向 P 端移动并变窄。
- HOXA7、HOXA9、HOXB9、HOXC10 等基因主要集中在胚胎后部。
- HOXC10、HOXD13 等后部 HOX 基因的表达峰尤其尖锐，提示其对应的轴向身份已经比较局限。
因此，图 2d 支持一个结论：随着胚胎沿前后轴向后延伸，不同区域逐渐获得更加明确、受限制的 HOX 组合，也就是前后轴区域化正在建立。
生物学意义
HOX 基因并不是简单地“越往后表达越高”，而是形成一种组合式的 positional code：
- 前部 HOX：表达较早、范围较宽；
- 中部 HOX：偏向躯干和中后部；
- 后部 HOX：集中在尾侧、脊髓、脊索、NMP、PSM、后肠等区域。
这说明 CS9 猕猴胚胎已经出现明显的轴向身份划分。作者据此认为，中央和后部 HOX 基因逐渐向后端限制，反映了早期胚胎的 progressive axial regionalization，即“渐进式前后轴区域化”。
此外，正文还指出 HOXA10、HOXA11、HOXA13 在尿囊中也有表达。由于尿囊位于胚胎后部，这提示尿囊在早期仍然具有类似后部中胚层的分子特征。
需要注意的地方
Fig. 2d 展示的是 CS9 胚胎的连续空间分布，不是时间序列，也不能单独证明 HOX 基因直接造成了这些组织命运。
示例数据#
Fig. 2d 的数据来源分为三层：
生物学样本:作者自己的 1 个 CS9 猕猴胚胎：
- 10 张近中线矢状切片；
- Stereo-seq 空间转录组；
- bin50 分辨率，即 50×50 DNBs，约 25 μm；
- 质控后 CS9 共 35,568 个 spots。
因此，Fig. 2d 不是来自公开 scRNA-seq 数据 GSE193007，也不是外部数据库的 HOX 表达数据。
图中实际使用的处理数据
作者将 CS9 空间转录组矩阵进行 3D 重建，排除非胚胎区域后：
- 以胚胎前端建立参考点；
- 计算每个 bin50 到参考点的欧氏距离；
- 按距离排序，形成 A→P 梯度；
- 划分为 100 个 AP 层；
- 计算每层中各 HOX 基因的平均表达；
- 以这些表达值绘制 Fig. 2d 的 ridgeline 和热图。
对应的作者源数据文件是：
Tan_2026_PrimateOrganogenesis_NatCellBiol_source_data_fig2.xlsx
其中 Panel d 工作表包含 31 个 HOX 基因 × 100 个 A–P 空间层（ AP_1 – AP_100 ）。
原始测序数据
该研究原始测序数据存放于 CNGBdb，登录号为 CNP0007017 。
所以最准确地说：Fig. 2d 是基于作者自测的 CS9 猕猴胚胎 Stereo-seq 空间转录组数据，经 3D 重建和 A–P 分层后得到的 HOX 表达空间分布。
解释的非常清楚！！！
我们去把附表下载下来就好了！
AI复现#
提示词比较简单：当然一次沟通并不能达到想要的效果，需要多沟通一些，我这里的xinye-plot是使用技能树的两个手搓代码的绘图专辑蒸馏得到的一个专属绘图skill。如果你想进群交流AI绘图和chatgpt的使用，添加新叶微信：jzhang0527。请备注来意~
根据fig2d的数据绘制一个同款图，图片的输出尺寸使用 $xinye-plot 。
复现结果如下：
本次的绘图代码如下：
suppressPackageStartupMessages({
library(readxl)
library(ggplot2)
library(ggridges)
library(patchwork)
library(cowplot)
library(scales)
})
input_xlsx <- file.path(
"references", "Tan_2026_PrimateOrganogenesis_NatCellBiol",
"Tan_2026_PrimateOrganogenesis_NatCellBiol_source_data_fig2.xlsx"
)
output_png <- "Fig2d_xinye_plot.png"
output_pdf <- "Fig2d_xinye_plot.pdf"
figure_width_in <- 7.0
figure_height_in <- figure_width_in * 918 / 1410
figure_dpi <- 1410 / figure_width_in
ridge_x_max <- 118
ridge_smooth_spar <- 0.45
ridge_tail_end <- 112
panel_d <- as.data.frame(read_excel(input_xlsx, sheet = "Panel d"))
gene_col <- names(panel_d)[1]
genes <- as.character(panel_d<a href="gene_col">gene_col</a>)
ap_names <- names(panel_d)[-1]
ap <- seq_along(ap_names)
expr <- as.matrix(panel_d[, -1, drop = FALSE])
rownames(expr) <- genes
colnames(expr) <- ap_names
groups <- list(
HOXA = grep("^HOXA", genes),
HOXB = grep("^HOXB", genes),
HOXC = grep("^HOXC", genes),
HOXD = grep("^HOXD", genes)
)
ridge_palettes <- list(
HOXA = c("#c95665", "#b18fa6", "#c66e60", "#d0a16e", "#888b90",
"#7e98bb", "#90aaa5", "#d09d68", "#8fa38d", "#eb5d60"),
HOXB = c("#af91a9", "#b56d60", "#d2a570", "#888e92", "#7895b8",
"#91a59a", "#d19c6c", "#8da382", "#eb5d60"),
HOXC = c("#8b9091", "#7894b6", "#91aaa3", "#d7945d", "#8da17e", "#ed5b5e"),
HOXD = c("#8b9091", "#7894b6", "#91aaa3", "#d7945d", "#8da17e", "#ed5b5e")
)
# Xinye-plot adaptation:
# 2025-065 “顶刊Nature杂志同款配色一流的山峦图” (code-03/04)
# supplies the ggridges + manual palette + scale/theme-classic pattern.
# Here geom_ridgeline is used because Panel d already contains a continuous
# expression profile at 100 AP layers, rather than one observation per cell.
row_z <- function(x) {
z <- t(apply(x, 1, function(v) {
s <- sd(v)
if (is.na(s) || s == 0) return(rep(0, length(v)))
as.numeric(scale(v))
}))
rownames(z) <- rownames(x)
colnames(z) <- colnames(x)
z
}
make_panel <- function(group_name, idx) {
group_genes <- genes[idx]
group_expr <- expr[idx, , drop = FALSE]
# Normalize each gene only for ridge height; this preserves the spatial shape
# while preventing high-abundance genes from dominating the panel geometry.
ridge_height <- t(apply(group_expr, 1, function(v) {
v <- v - min(v)
if (max(v) == 0) return(rep(0, length(v)))
v / max(v)
}))
rownames(ridge_height) <- group_genes
colnames(ridge_height) <- ap_names
ridge_df <- do.call(rbind, lapply(seq_along(group_genes), function(i) {
# Lightly smooth the 100-layer profile for visual rendering. The post-100
# tail is a continuous display extension, not additional AP measurements.
smooth_fit <- smooth.spline(
ap, as.numeric(ridge_height[i, ]), spar = ridge_smooth_spar
)
smooth_x <- seq(min(ap), ridge_x_max, length.out = 360)
smooth_y <- predict(smooth_fit, x = pmin(smooth_x, max(ap)))$y
tail_idx <- smooth_x > max(ap)
if (any(tail_idx)) {
# Cubic Hermite tail: preserve the fitted endpoint slope, then join
# smoothly to zero at the end of the visual tail.
x0 <- max(ap)
x1 <- ridge_tail_end
y0 <- predict(smooth_fit, x = x0)$y
y1 <- 0
slope0 <- (predict(smooth_fit, x = x0)$y -
predict(smooth_fit, x = x0 - 5)$y) / 5
tail_x <- pmin(smooth_x[tail_idx], x1)
tail_fraction <- (tail_x - x0) / (x1 - x0)
tail_length <- x1 - x0
smooth_y[tail_idx] <-
(2 * tail_fraction^3 - 3 * tail_fraction^2 + 1) * y0 +
(tail_fraction^3 - 2 * tail_fraction^2 + tail_fraction) *
tail_length * slope0 +
(-2 * tail_fraction^3 + 3 * tail_fraction^2) * y1
smooth_y[smooth_x > x1] <- 0
}
data.frame(
AP = smooth_x,
gene = factor(group_genes[i], levels = rev(group_genes)),
y = length(group_genes) - i + 1,
height = pmax(0, smooth_y)
)
}))
ridge <- ggplot(ridge_df, aes(x = AP, y = y, height = height, fill = gene)) +
geom_ridgeline(
scale = 1.15, min_height = 0, color = "#1f1f1f", linewidth = 0.18,
alpha = 1
) +
scale_fill_manual(values = setNames(ridge_palettes<a href="group_name">group_name</a>[seq_along(group_genes)], group_genes)) +
# In the reference, 100 is the last tick, not the right panel boundary.
# Keep a small post-100 margin so the terminal ridge is not clipped.
scale_x_continuous(
breaks = c(1, 21, 41, 61, 81, 100),
labels = c("0", "20", "40", "60", "80", "100"),
limits = c(0, ridge_x_max),
expand = c(0, 0)
) +
scale_y_continuous(
breaks = seq_along(group_genes),
labels = rev(group_genes),
expand = expansion(mult = c(0.02, 0.03))
) +
coord_cartesian(xlim = c(0, ridge_x_max), expand = FALSE, clip = "off") +
labs(x = NULL, y = NULL) +
theme_classic(base_size = 6) +
theme(
legend.position = "none",
axis.text.y = element_text(face = "italic", size = 6.2),
axis.text.x = element_text(size = 5.2),
axis.title.x = element_blank(),
axis.ticks.y = element_blank(),
panel.grid.major.x = element_line(color = "#e5e5e5", linewidth = 0.15),
plot.margin = margin(1.5, 4, 1, 1)
)
heat_z <- row_z(group_expr)
heat_df <- do.call(rbind, lapply(seq_along(group_genes), function(i) {
data.frame(
gene = factor(group_genes[i], levels = rev(group_genes)),
AP = ap,
z = as.numeric(heat_z[i, ])
)
}))
heat <- ggplot(heat_df, aes(x = AP, y = gene, fill = z)) +
geom_tile(width = 1, height = 1) +
scale_x_continuous(
breaks = c(1, 100), labels = c("A", "P"),
expand = expansion(mult = c(0, 0))
) +
scale_y_discrete(position = "right", drop = FALSE) +
scale_fill_gradient2(
low = "#425794", mid = "#f3f3f3", high = "#e4515a", midpoint = 0,
limits = c(-4, 4), oob = squish,
breaks = c(-4, 0, 4),
name = "Exp."
) +
labs(x = NULL, y = NULL) +
theme_classic(base_size = 6) +
theme(
axis.text.y = element_text(
face = "italic", size = 6.2, hjust = 0, margin = margin(l = 2)
),
axis.ticks.y = element_blank(),
axis.text.x = element_text(size = 5.2),
axis.title.x = element_blank(),
legend.position = "bottom",
legend.direction = "horizontal",
legend.title = element_text(size = 5.2),
legend.text = element_text(size = 4.8),
legend.key.width = unit(0.35, "cm"),
legend.key.height = unit(0.10, "cm"),
legend.box.spacing = unit(0, "pt"),
plot.margin = margin(1.5, 8, 1, 1)
) +
guides(fill = guide_colorbar(
title.position = "top", title.hjust = 0.5,
barwidth = unit(0.85, "cm"), barheight = unit(0.10, "cm"),
ticks = FALSE
))
ridge_title <- ggdraw() + draw_label(
paste0(group_name, " clusters"), x = 0.5, y = 0.5,
hjust = 0.5, vjust = 0.5, size = 9
)
heat_title <- ggdraw() + draw_label(
"Average expression level", x = 0.5, y = 0.5,
hjust = 0.5, vjust = 0.5, size = 7
)
# Split titles, equal-height main plots, and the heatmap colorbar into rows.
# This keeps the ridge/heatmap bodies aligned and places the colorbar on a
# shared horizontal guide row below the two main plots.
heat_legend <- get_legend(heat)
heat_no_legend <- heat + theme(legend.position = "none")
# Reference-like balance: narrower ridge plot, wider heatmap + gene labels.
panel_widths <- c(1.20, 0.08, 0.92)
blank <- ggplot() +
theme_void() +
theme(
plot.background = element_rect(fill = "white", colour = NA),
panel.background = element_rect(fill = "white", colour = NA)
)
title_row <- plot_grid(
ridge_title, blank, heat_title,
ncol = 3, rel_widths = panel_widths
)
body_row <- plot_grid(
ridge, blank, heat_no_legend,
ncol = 3, rel_widths = panel_widths
)
ridge_footer <- ggdraw() + draw_label(
"A -> P", x = 0.5, y = 0.5,
hjust = 0.5, vjust = 0.5, size = 6.0
)
heat_ap_footer <- ggdraw() + draw_label(
"A -> P", x = 0.38, y = 0.5,
hjust = 0.5, vjust = 0.5, size = 6.0
)
heat_footer <- plot_grid(
heat_ap_footer, heat_legend,
ncol = 2, rel_widths = c(0.78, 1.22)
)
footer_row <- plot_grid(
ridge_footer, blank, heat_footer,
ncol = 3, rel_widths = panel_widths
)
plot_grid(
title_row, body_row, footer_row,
ncol = 1, rel_heights = c(0.15, 1, 0.16), align = "v"
)
}
fig <- wrap_plots(
make_panel("HOXA", groups$HOXA),
make_panel("HOXB", groups$HOXB),
make_panel("HOXC", groups$HOXC),
make_panel("HOXD", groups$HOXD),
ncol = 2
) +
plot_annotation(
title = "d",
theme = theme(
plot.title = element_text(hjust = 0, face = "bold", size = 11),
plot.margin = margin(1, 2, 1, 2)
)
)
# Explicit output block following the Xinye-plot recipe convention.
# Flexible dimensions: preserve the reference aspect ratio while allowing the
# physical export width to be changed without using the Xinye example's fixed
# 5.5 x 5 inch output.
ggsave(filename = output_png, width = figure_width_in, height = figure_height_in,
dpi = figure_dpi, plot = fig, bg = "white")
ggsave(filename = output_pdf, width = figure_width_in, height = figure_height_in,
device = "pdf", plot = fig, bg = "white")
数据代码都给你了，但是如果是你自己做复现，可能得不到这样的结果，试试看！
画了这么久了，现在根本就再也不想回去手搓代码了怎么办？
生信技能树在8月份有一个长沙的线下培训班，感兴趣的瞄一瞄：
暑期特训长沙站：Codex赋能单细胞+空间转录组线下培训班来袭！
此外，虚拟细胞开了一个复现群，来吧： 2026年最火的双虚拟敲除：CellOracle 虚拟 KO + scTenifoldKnk 验证 → 锁定目标基因
绘图小技巧2026 · 目录
闪记
复制 LaTeX 公式
