Science同款｜个性化通路富集-GSEA排序气泡图

摘要LLM 归库摘要

• 核心内容： 这是一篇复现 Hallmark GSEA 排序气泡图的单细胞生信教程，围绕 treatment vs control 上皮细胞差异分析、GSEA 富集和 ggplot2 气泡图展示展开。
• 适用场景： 已有 Seurat 对象、细胞注释和分组信息时，用于比较两组目标细胞的通路整体变化，并制作类似 CNS/Science 风格的 NES 排序气泡图。
• 输入： 质控/注释后的 Seurat 对象、目标细胞类型字段、实验分组字段、Hallmark 基因集。
• 输出： GSEA 结果表、每个通路的 NES/pvalue/setSize/检出率，以及一张以 NES 为纵轴、通路排序为横轴的气泡图。
• 关键参数/步骤： FindMarkers() 构建 avg_log2FC 排序向量；msigdbr(category = "H") 获取 Hallmark；clusterProfiler::GSEA() 运行富集；点大小映射检出率，透明度映射 -log10(pvalue)。
• 常见错误： 原文代码块被微信剪藏压缩，部分 R 语句缺少换行；第 8 步 detect_rate 处疑似代码缺失，实际复现时需补为对通路基因 gene_detect 的均值计算。
• 与口腔生信连接： 可迁移到牙周炎、OSCC 或口腔免疫单细胞项目中，用于展示上皮细胞、免疫细胞或肿瘤细胞在不同分组间的 Hallmark 通路方向性变化。
• 关联页面建议： 后续可整合到 03 生信教程与工具库/02_wiki/workflows/单细胞教程 或可视化方法页。
• 下一步处理： 若准备复现，建议先整理为可运行 R 脚本，修复换行和检出率计算，再保存示例输出图。
• 归库说明： 分类到 03 生信教程与工具库，置信度 0.88；理由是其核心是单细胞 GSEA 可视化教程，而非单篇论文精读。

整理版可运行 R 代码#

说明使用说明

• 下面代码把原文中被微信压缩到一行的 R 代码重新整理成可读、可改、可运行版本。
• 你只需要优先修改 用户参数区：对象路径、目标细胞类型、分组列名、实验组/对照组名称、输出目录。
• 第 8 步 detect_rate 已补全：对每个 Hallmark 通路内、在当前 Seurat 对象中存在的基因，计算其在目标细胞中的平均检出率。
• 若你使用 Seurat v5，代码会优先尝试 layer = "counts"；若是 Seurat v4，则自动回退到 slot = "counts"。

## ============================================================
## Hallmark GSEA 排序气泡图｜整理版
## 目的：基于 Seurat 对象的 treatment vs control 差异分析结果，
##      绘制 Hallmark gene sets 的 NES 排序气泡图。
## 输入：已经完成质控、降维、聚类、细胞注释和分组标注的 Seurat 对象
## 输出：marker 表、GSEA 表、绘图数据表、GSEA bubble plot
## ============================================================

## 0. 加载 R 包 -------------------------------------------------
suppressPackageStartupMessages({
  library(Seurat)
  library(dplyr)
  library(tibble)
  library(ggplot2)
  library(ggrepel)
  library(clusterProfiler)
  library(msigdbr)
})

## 1. 用户参数区 ------------------------------------------------
## 根据自己的数据修改这里即可。
work_dir <- "E:/xy/2026/five/5.22/"
obj_path <- "E:/xy/public/scRNA.sub.rdata"
out_dir  <- file.path(work_dir, "gsea_bubble_output")

## load(obj_path) 后，环境里 Seurat 对象的名字。
## 原文对象名是 scRNA.sub。
seurat_object_name <- "scRNA.sub"

## 元数据列名。
celltype_col <- "celltype"
group_col    <- "group"

## 要分析的细胞类型和分组。
target_celltype <- "Epithelial_cells"
treatment_group <- "treatment"
control_group   <- "control"

## 差异分析参数。
min_pct       <- 0.5
logfc_cutoff  <- 0.5
p_adj_cutoff  <- 0.05

## GSEA 参数。
species_name <- "Homo sapiens"
min_gs_size  <- 10
max_gs_size  <- 5000
seed_value   <- 123

## 检出率阈值。
## 原文思路是“基因是否被检测到”。单细胞常用 counts > 0；
## 如果想更严格，可以改成 1 或 5。
detect_count_threshold <- 0

## 标签数量：NES 最高和最低各标几个通路。
n_label_each_side <- 3

## 图形保存参数。
plot_width  <- 9
plot_height <- 6
plot_dpi    <- 300

set.seed(seed_value)
dir.create(out_dir, showWarnings = FALSE, recursive = TRUE)
setwd(work_dir)

## 2. 读取 Seurat 对象 -----------------------------------------
load(obj_path)

if (!exists(seurat_object_name)) {
  stop("没有在 rdata 中找到对象：", seurat_object_name)
}

atlas_obj <- get(seurat_object_name)

if (!inherits(atlas_obj, "Seurat")) {
  stop(seurat_object_name, " 不是 Seurat 对象，请检查 seurat_object_name。")
}

## 检查元数据列是否存在。
meta_cols <- colnames(atlas_obj@meta.data)
required_cols <- c(celltype_col, group_col)
missing_cols <- setdiff(required_cols, meta_cols)
if (length(missing_cols) > 0) {
  stop("Seurat meta.data 缺少这些列：", paste(missing_cols, collapse = ", "))
}

## 3. 提取目标细胞 ---------------------------------------------
focus_obj <- subset(
  atlas_obj,
  subset = .data<a href="celltype_col">celltype_col</a> == target_celltype
)

message("目标细胞数量：", ncol(focus_obj))
print(table(focus_obj@meta.data<a href="group_col">group_col</a>, useNA = "ifany"))

if (!all(c(treatment_group, control_group) %in% focus_obj@meta.data<a href="group_col">group_col</a>)) {
  stop("目标细胞中没有同时包含 treatment_group 和 control_group，请检查分组名称。")
}

## 4. 差异分析 --------------------------------------------------
Idents(focus_obj) <- group_col

marker_tbl <- FindMarkers(
  object = focus_obj,
  ident.1 = treatment_group,
  ident.2 = control_group,
  min.pct = min_pct,
  logfc.threshold = logfc_cutoff
)

## 兼容 Seurat 不同版本的 logFC 列名。
fc_col <- dplyr::case_when(
  "avg_log2FC" %in% colnames(marker_tbl) ~ "avg_log2FC",
  "avg_logFC"  %in% colnames(marker_tbl) ~ "avg_logFC",
  TRUE ~ NA_character_
)

if (is.na(fc_col)) {
  stop("FindMarkers 结果中没有 avg_log2FC 或 avg_logFC 列。")
}

marker_tbl <- marker_tbl %>%
  rownames_to_column("gene") %>%
  arrange(desc(.data<a href="fc_col">fc_col</a>)) %>%
  mutate(
    change_type = case_when(
      .data<a href="fc_col">fc_col</a> >  logfc_cutoff & p_val_adj < p_adj_cutoff ~ "up",
      .data<a href="fc_col">fc_col</a> < -logfc_cutoff & p_val_adj < p_adj_cutoff ~ "down",
      TRUE ~ "stable"
    )
  )

write.csv(
  marker_tbl,
  file = file.path(out_dir, "01_FindMarkers_marker_tbl.csv"),
  row.names = FALSE
)

## 5. 构建 GSEA 排序向量 ---------------------------------------
rank_tbl <- marker_tbl %>%
  select(gene, rank_score = all_of(fc_col)) %>%
  filter(!is.na(rank_score), !is.na(gene), gene != "") %>%
  group_by(gene) %>%
  summarise(rank_score = mean(rank_score), .groups = "drop") %>%
  arrange(desc(rank_score))

rank_vec <- rank_tbl$rank_score
names(rank_vec) <- rank_tbl$gene
rank_vec <- sort(rank_vec, decreasing = TRUE)

if (length(rank_vec) < min_gs_size) {
  stop("排序基因数量过少，无法进行 GSEA。")
}

## 6. 准备 Hallmark 基因集 -------------------------------------
## msigdbr 不同版本参数可能略有差异：新版推荐 collection = "H"，旧版使用 category = "H"。
hallmark_raw <- tryCatch(
  msigdbr(species = species_name, collection = "H"),
  error = function(e) msigdbr(species = species_name, category = "H")
)

hallmark_map <- hallmark_raw %>%
  transmute(
    pathway = gsub("^HALLMARK_", "", gs_name),
    symbol  = gene_symbol
  ) %>%
  distinct() %>%
  filter(!is.na(pathway), !is.na(symbol), pathway != "", symbol != "")

write.csv(
  hallmark_map,
  file = file.path(out_dir, "02_Hallmark_TERM2GENE.csv"),
  row.names = FALSE
)

## 7. 运行 GSEA -------------------------------------------------
gsea_obj <- GSEA(
  geneList      = rank_vec,
  TERM2GENE     = hallmark_map,
  pvalueCutoff  = 1,
  minGSSize     = min_gs_size,
  maxGSSize     = max_gs_size,
  verbose       = FALSE,
  seed          = TRUE
)

gsea_df <- as.data.frame(gsea_obj@result)

if (nrow(gsea_df) == 0) {
  stop("GSEA 没有返回结果。请检查基因命名是否与 msigdbr 的 gene_symbol 一致。")
}

write.csv(
  gsea_df,
  file = file.path(out_dir, "03_GSEA_result.csv"),
  row.names = FALSE
)

## 8. 计算每个通路的平均检出率 --------------------------------
## 兼容 Seurat v5 layer 和 Seurat v4 slot。
expr_mat <- tryCatch(
  GetAssayData(focus_obj, layer = "counts"),
  error = function(e) GetAssayData(focus_obj, slot = "counts")
)

## 每个基因在目标细胞中的检出比例。
gene_detect <- Matrix::rowSums(expr_mat > detect_count_threshold) / ncol(expr_mat)

plot_df <- gsea_df %>%
  select(ID, Description, NES, setSize, pvalue, p.adjust) %>%
  rowwise() %>%
  mutate(
    detect_rate = {
      genes_in_pathway <- hallmark_map$symbol[hallmark_map$pathway == ID]
      genes_in_pathway <- intersect(genes_in_pathway, names(gene_detect))

      if (length(genes_in_pathway) == 0) {
        NA_real_
      } else {
        mean(gene_detect[genes_in_pathway], na.rm = TRUE)
      }
    }
  ) %>%
  ungroup() %>%
  arrange(desc(NES)) %>%
  mutate(
    pathway_id     = factor(ID, levels = ID),
    x_pos          = row_number(),
    detect_percent = 100 * detect_rate,
    neg_log10_p    = -log10(pmax(pvalue, .Machine$double.xmin)),
    direction      = if_else(NES >= 0, treatment_group, control_group)
  )

write.csv(
  plot_df,
  file = file.path(out_dir, "04_GSEA_bubble_plot_data.csv"),
  row.names = FALSE
)

## 9. 标注 NES 最高和最低的通路 -------------------------------
label_terms <- c(
  head(plot_df$ID, n_label_each_side),
  tail(plot_df$ID, n_label_each_side)
) %>%
  unique()

label_palette <- c("#D11E00", "#9972B2", "#A99909", "#D39E44", "#0072B2", "#009E73")

label_df <- plot_df %>%
  filter(ID %in% label_terms) %>%
  mutate(label_col = rep(label_palette, length.out = n()))

## 10. 绘制 GSEA 排序气泡图 -----------------------------------
y_low  <- floor(min(plot_df$NES, na.rm = TRUE))
y_high <- ceiling(max(plot_df$NES, na.rm = TRUE))

x_break_step <- 10
x_max <- max(plot_df$x_pos, na.rm = TRUE)

## 如果检出率存在 NA，绘图前给出提示。
if (any(is.na(plot_df$detect_percent))) {
  warning("部分通路没有可匹配基因，detect_percent 为 NA；这些点的大小可能无法正常显示。")
}

gsea_bubble <- ggplot(plot_df, aes(x = x_pos, y = NES)) +
  geom_hline(yintercept = 0, linetype = "dashed", linewidth = 0.4, color = "grey50") +
  geom_point(
    aes(size = detect_percent, alpha = neg_log10_p),
    shape = 21,
    stroke = 0.7,
    fill = "#0000ff80",
    colour = "black"
  ) +
  geom_text_repel(
    data = label_df,
    aes(label = ID),
    size = 3,
    color = label_df$label_col,
    force = 20,
    point.padding = 0.5,
    min.segment.length = 0,
    hjust = 1.2,
    segment.color = "grey20",
    segment.size = 0.3,
    segment.alpha = 0.8,
    nudge_y = -0.1,
    max.overlaps = Inf
  ) +
  scale_size_continuous(
    range = c(0.5, 8),
    name = "Detection rate (%)"
  ) +
  scale_alpha_continuous(
    range = c(0.15, 0.9),
    name = "Significance\n(-log10 p-value)"
  ) +
  scale_x_continuous(
    breaks = seq(0, ceiling(x_max / x_break_step) * x_break_step, by = x_break_step),
    labels = seq(0, ceiling(x_max / x_break_step) * x_break_step, by = x_break_step),
    expand = expansion(mult = c(0.03, 0.05))
  ) +
  scale_y_continuous(
    breaks = seq(y_low, y_high, by = 1),
    limits = c(y_low, y_high),
    expand = expansion(mult = c(0.08, 0.08))
  ) +
  labs(
    title = "Epithelial Cell Pathway Enrichment Analysis",
    subtitle = paste0("Hallmark gene sets based on ", treatment_group, " versus ", control_group),
    x = "Hallmark gene sets ranked by NES",
    y = "Normalized Enrichment Score (NES)"
  ) +
  guides(
    size  = guide_legend(title = "Detection rate (%)", order = 1),
    alpha = guide_legend(title = "Significance\n(-log10 p-value)", order = 2)
  ) +
  theme_classic(base_size = 15) +
  theme(
    plot.title      = element_text(face = "bold", size = 16, hjust = 0.5),
    plot.subtitle   = element_text(size = 12, hjust = 0.5, color = "gray40"),
    axis.line       = element_line(color = "black", linewidth = 0.6),
    axis.text       = element_text(face = "bold"),
    axis.title      = element_text(size = 14, face = "bold"),
    legend.title    = element_text(face = "bold"),
    legend.position = "right",
    legend.box      = "vertical",
    panel.grid      = element_blank(),
    panel.background = element_rect(fill = "white", color = NA)
  )

print(gsea_bubble)

ggsave(
  filename = file.path(out_dir, "05_GSEA_bubble_plot.pdf"),
  plot = gsea_bubble,
  width = plot_width,
  height = plot_height
)

ggsave(
  filename = file.path(out_dir, "05_GSEA_bubble_plot.png"),
  plot = gsea_bubble,
  width = plot_width,
  height = plot_height,
  dpi = plot_dpi
)

message("分析完成。结果已输出到：", out_dir)

关键修复点#

1. 加载包部分：原文 library(ggrepel) library(clusterProfiler) 中间缺少分隔；整理版已拆成多行。
2. 读取对象部分：原文把 work_dir、obj_path、setwd()、load()、对象赋值压到一行；整理版加入对象名检查。
3. 差异分析部分：整理版兼容 avg_log2FC / avg_logFC 两种列名，并保存 marker 表。
4. GSEA 排序向量：整理版去除 NA、处理重复基因名，并保证按 logFC 从高到低排序。
5. Hallmark 基因集：整理版兼容新版 msigdbr(collection = "H") 和旧版 msigdbr(category = "H")。
6. 检出率计算：补全原文缺失逻辑：mean(gene_detect[genes_in_pathway], na.rm = TRUE)。
7. 显著性映射：使用 -log10(pmax(pvalue, .Machine$double.xmin))，避免 p 值为 0 时出现无限大。
8. 输出文件：自动保存 marker、Hallmark 映射、GSEA 结果、绘图数据、PDF 和 PNG。

华哥生信 2026年5月23日 20:51

CNS | 参考文献

文章 Fig1 c图

(C) GSEA for 50 Hallmark gene sets within DEGs after KRAS siRNA treatment. Detection is based on the presence of a gene in all eight PDAC cell lines at >5 reads.

解读

Hallmark 通路 GSEA 结果气泡图。纵轴是 NES 标准化富集分数 ，NES 越大，说明该通路越偏向上调方向；NES 越小，说明该通路越偏向下调方向。横轴不是连续变量，而是把不同 Hallmark gene sets 按某种顺序排列出来。每一个圆点代表一个 Hallmark 通路；圆点越大，说明这个通路中被检测到的基因比例越高；颜色越深，说明显著性越强，也就是 -log10(p-value) 越大。

图中上方的 KRAS_SIGNALING_DN、INTERFERON_ALPHA_RESPONSE 是正向富集；下方的 E2F_TARGETS、G2M_CHECKPOINT、MYC_TARGETS、KRAS_SIGNALING_UP、EMT、TNFA_SIGNALING_VIA_NFKB 是负向富集，说明这些通路在比较组中整体呈下降趋势，尤其是 E2F、G2M、MYC 这类增殖/细胞周期相关通路，显著性和检测比例都很高。

欢迎加入CNS生信交流群一起讨论与研究时下热门CNS文章各类生信分析方法（若二维码过期可以添加微信：huagesx）

跟着二十篇CNS文章学 虚拟细胞大模型

新课发布：虚拟细胞（跟着二十篇顶刊文章学虚拟细胞）

复现#

1 加载环境#

suppressPackageStartupMessages({  library(Seurat); library(dplyr);   library(ggplot2); library(ggrepel)  library(clusterProfiler); library(msigdbr)})

2 读取对象#

先设置工作路径，然后读取已经保存好的 Seurat 对象。这里的 scRNA.sub.rdata 是前面完成质控、降维、聚类和细胞注释后的对象。读取完成后，将对象赋值给 atlas_obj ，方便后面继续分析。

work_dir <- "E:/xy/2026/five/5.22/"obj_path <- "E:/xy/public/scRNA.sub.rdata"setwd(work_dir); load(obj_path)atlas_obj <- scRNA.sub

需要示例数据请联系原作者或查看原文渠道。

3 提取上皮#

目的是从完整单细胞对象中单独提取上皮细胞。因为后续分析关注的是 treatment 和 control 条件下上皮细胞的通路变化，所以不能把所有细胞类型混在一起做差异分析。 subset() 根据 celltype == "Epithelial_cells" 筛选目标细胞。提取后再查看细胞数量和分组情况，确认 treatment 与 control 两组都存在。

focus_obj <- subset(atlas_obj,                     subset = celltype == "Epithelial_cells")

4 差异分析#

用 treatment 和 control 两组上皮细胞进行差异基因分析。先把 Seurat 对象的身份设置为 group ，这样 FindMarkers() 才能按照分组比较。这里设置 ident.1 = "treatment" ，表示正的 log2FC 代表 treatment 组表达更高。分析完成后根据 log2FC 和校正 P 值给基因标记为上调、下调或稳定，方便后面理解整体变化方向。

Idents(focus_obj) <- "group"marker_tbl <- FindMarkers(  focus_obj, ident.1 = "treatment", ident.2 = "control",  min.pct = 0.5, logfc.threshold = 0.5)
marker_tbl <- marker_tbl %>%  arrange(desc(avg_log2FC)) %>%  mutate(change_type = case_when(    avg_log2FC > 0.5 & p_val_adj < 0.05 ~ "up",    avg_log2FC < -0.5 & p_val_adj < 0.05 ~ "down",    TRUE ~ "stable"))

' fill='%23FFFFFF'%3E%3Crect x='249' y='126' width='1' height='1'%3E%3C/rect%3E%3C/g%3E%3C/g%3E%3C/svg%3E)

5 构建排序#

GSEA 不是只看显著差异基因，而是需要把所有参与排序的基因按照某个统计量从高到低排列。这里使用 avg_log2FC 作为排序指标，正值代表 treatment 组更高，负值代表 control 组更高。

rank_vec <- marker_tbl$avg_log2FCnames(rank_vec) <- rownames(marker_tbl)rank_vec <- sort(rank_vec, decreasing = TRUE)

6 准备基因集#

从 msigdbr 中读取 Hallmark 基因集。Hallmark 是 MSigDB 中比较常用的一套通路集合，数量相对精简，适合做整体生物学过程解释。这里选择人类基因集，所以 species 设置为 Homo sapiens 。随后只保留通路名称和基因符号两列，并去掉 HALLMARK_ 前缀，让后面的图更简洁。

hallmark_raw <- msigdbr(species = "Homo sapiens", category = "H")hallmark_map <- hallmark_raw %>%  transmute(pathway = gsub("HALLMARK_", "", gs_name),             symbol = gene_symbol)

' fill='%23FFFFFF'%3E%3Crect x='249' y='126' width='1' height='1'%3E%3C/rect%3E%3C/g%3E%3C/g%3E%3C/svg%3E)

7 运行 GSEA#

把前面构建好的排序向量和 Hallmark 基因集放进 GSEA() 。这里 pvalueCutoff = 1 表示先保留全部结果，后面可以根据 NES 和 P 值自己筛选或展示。 minGSSize 和 maxGSSize 设置得比较宽，是为了避免通路因为基因数量限制被提前过滤。运行完成后，将 GSEA 结果转成普通数据框，方便后续整理和画图。

gsea_obj <- GSEA(  geneList = rank_vec, TERM2GENE = hallmark_map,  pvalueCutoff = 1, minGSSize = 1, maxGSSize = 500000,  verbose = FALSE)
gsea_df <- as.data.frame(gsea_obj@result)

' fill='%23FFFFFF'%3E%3Crect x='249' y='126' width='1' height='1'%3E%3C/rect%3E%3C/g%3E%3C/g%3E%3C/svg%3E)

需要示例数据请联系原作者或查看原文渠道。

8 计算检出率#

每个 Hallmark 通路增加一个新的展示维度：平均检出率。代码先从上皮细胞对象中提取 counts 矩阵，然后计算每个基因在多少比例的细胞中被检测到。接着对每个通路内部的基因取平均检出率，得到该通路在目标细胞中的整体表达覆盖程度。最后把 NES、P 值和检出率整理到一个表格中，作为气泡图的输入数据。

expr_mat <- GetAssayData(focus_obj, slot = "counts")gene_detect <- rowSums(expr_mat > 1) / ncol(expr_mat)
plot_df <- gsea_df %>%  select(ID, NES, setSize, pvalue) %>%  rowwise() %>%  mutate(detect_rate = {    gs <- hallmark_map$symbol[hallmark_map$pathway == ID]    gs <- intersect(gs, names(gene_detect))0
  }) %>%  ungroup() %>%  arrange(desc(NES)) %>%  mutate(pathway_id = factor(ID, levels = ID), x_pos = row_number(), detect_percent = detect_rate * 100)

' fill='%23FFFFFF'%3E%3Crect x='249' y='126' width='1' height='1'%3E%3C/rect%3E%3C/g%3E%3C/g%3E%3C/svg%3E)

9 标注通路#

NES 排名前 3 和排名后 3 的通路作为标注对象。这样图中既能展示 treatment 方向富集的通路，也能展示 control 方向富集的通路。 label_col 给不同标签设置颜色，方便在最终图中区分标注文字。

label_terms <- c(head(as.character(plot_df$ID), 3), tail(as.character(plot_df$ID), 3))
label_df <- plot_df %>%  filter(ID %in% label_terms) %>%  mutate(label_col = rep(c("#D11E00", "#9972B2", "#A99909", "#D39E44"), length.out = n()))

' fill='%23FFFFFF'%3E%3Crect x='249' y='126' width='1' height='1'%3E%3C/rect%3E%3C/g%3E%3C/g%3E%3C/svg%3E)

10 绘制气泡#

把 GSEA 结果整合成一张通路气泡图。横坐标表示 Hallmark 通路的排序位置，纵坐标表示 NES，点越靠上说明越偏 treatment 组富集，越靠下说明越偏 control 组富集。气泡大小表示通路基因在上皮细胞中的平均检出率，透明度表示统计显著性。最后通过 geom_text_repel() 标注关键通路，避免标签互相遮挡。

y_low <- floor(min(plot_df$NES, na.rm = TRUE))y_high <- ceiling(max(plot_df$NES, na.rm = TRUE))gsea_bubble <- ggplot(plot_df, aes(x = x_pos, y = NES)) +  geom_point(aes(size = detect_percent, alpha = -log10(pvalue)),             shape = 21, stroke = 0.7, fill = "#0000ff80", colour = "black") +  geom_text_repel(data = label_df, aes(label = ID), size = 3,                  color = label_df$label_col, force = 20, point.padding = 0.5,                  min.segment.length = 0, hjust = 1.2, segment.color = "grey20",                  segment.size = 0.3, segment.alpha = 0.8, nudge_y = -0.1) +  scale_size_continuous(range = c(0.2, 8), name = "Detection rate (%)") +  scale_alpha_continuous(range = c(0.1, 0.9), name = "Significance\n(-log10 p-value)") +  scale_x_continuous(breaks = seq(0, 50, by = 10), labels = seq(0, 50, by = 10)) +  scale_y_continuous(breaks = seq(y_low, y_high, by = 1),                     labels = seq(y_low, y_high, by = 1),                     limits = c(y_low, y_high), expand = c(0.1, 0.1)) +  labs(title = "Epithelial Cell Pathway Enrichment Analysis",       subtitle = "Hallmark gene sets based on treatment versus control",       x = "Hallmark gene sets", y = "Normalized Enrichment Score (NES)") +  guides(size = guide_legend(title = "Detection rate (%)", order = 1),         alpha = guide_legend(title = "Significance\n(-log10 p-value)", order = 2)) +  theme_classic(base_size = 15) +  theme(plot.title = element_text(face = "bold", size = 16, hjust = 0.5),        plot.subtitle = element_text(size = 12, hjust = 0.5, color = "gray40"),        axis.line = element_line(color = "black", linewidth = 0.6),        axis.text = element_text(face = "bold"),        axis.title = element_text(size = 14, face = "bold"),        legend.title = element_text(face = "bold"),        legend.position = "right", legend.box = "vertical",        panel.grid = element_blank(),        panel.background = element_rect(fill = "white", color = NA))print(gsea_bubble)

' fill='%23FFFFFF'%3E%3Crect x='249' y='126' width='1' height='1'%3E%3C/rect%3E%3C/g%3E%3C/g%3E%3C/svg%3E)

需要示例数据请联系原作者或查看原文渠道。

图形解读#

这张图的核心是把差异分析、通路富集和表达检出率放在同一个坐标系中。NES 大于 0 的通路更偏 treatment 组，NES 小于 0 的通路更偏 control 组。气泡越大，说明该通路相关基因在上皮细胞中整体检出越充分，结果解释时更有表达基础。透明度越高，代表 P 值越小，也就是该通路的统计信号更强。

' fill='%23FFFFFF'%3E%3Crect x='249' y='126' width='1' height='1'%3E%3C/rect%3E%3C/g%3E%3C/g%3E%3C/svg%3E)

解读重点：
1. NES > 0：treatment 方向富集
2. NES < 0：control 方向富集
3. 气泡大小：通路基因平均检出率
4. 透明度：-log10(pvalue)，越高越显著

AI大模型：当别人还在使用模型的时候，我们已经开始构建大模型了

AI大模型基于病理切片预测空间蛋白，CELL和NM文章复现及原理全解析

孟德尔＋单细胞全文复现

孟德尔+单细胞全文复现（纯生信公共数据挖掘）

零基础想学生信，看 👇 这个文章 ：

零基础也能学的AI生信课（AI助力生信入门班即将开始）

想学机器学习，看 👇 这个文章 ：

AI助力多组学与机器学习联合分析（机器学习分析代谢组、蛋白组、宏基因组、网络药理学、转录组）

想系统学单细胞，看 👇 这个文章 ：

CNS文章单细胞数据分析内容汇总（基础--进阶--python分析）

想系统学空间转录组，看 👇 这个文章 ：

Nature单细胞空间组学1:1复现训练营

跟着CNS文章学空间转录组分析线上直播课（Visium HD、Xenium、Stereo-seq）

欢迎关注 | 华哥科研平台

亲，写的这么辛苦，记得关注、点赞、转发哟！

继续滑动看下一个

Bioinformation

向上滑动看下一个

闪记

复制 LaTeX 公式