现在有一个癌型的bulk data数据集,包括患者对药物的响应情况(R和NR)表格,表格中有两列,第一列为患者id,可能没有列名,第二列为患者响应信息,0代表NR,1代表R。(路径:/mnt/public/zhengyeyang/bulk/bulk_HCC/GSE109211/clin.csv),还有患者基因表达log矩阵(路径:/mnt/public/zhengyeyang/bulk/bulk_HCC/GSE109211/matrix_.csv),此外,我还找到了这个癌型相关的空间结构的signature genes("CD68", "CSF1R", "EPCAM", "VEGFA", "TP53" ),请问要如何写循环代码,实现signature genes的富集(利用ssgsea包),并输出一个箱形图,展示R和NR患者的富集signature score的组别差异。要求体现t-test 和p value,图要美观规范,符合科研图表要求,最终将结果存放在 "/mnt/public/zhengyeyang/Enrichment/output_HCC"目录下,请帮我生成R语言代码,谢谢。
时间: 2025-07-11 13:42:41 浏览: 5
<think>嗯,用户需要处理一个癌症的bulk数据,想要通过ssGSEA方法对特定的signature基因进行富集分析,然后比较响应组(R)和非响应组(NR)的富集分数差异,并生成箱形图。首先,我需要理清步骤,确保代码正确且符合科研规范。
首先,数据读取部分。用户提供了两个csv文件,clin.csv和matrix_.csv。clin.csv包含患者ID和响应情况,第二列是0和1,分别代表NR和R。matrix_.csv是基因表达矩阵,行是基因,列是样本。需要注意检查两个文件中的患者ID是否匹配,可能需要将clin.csv中的ID与表达矩阵的列名对应起来,可能需要处理ID的格式问题,比如是否有后缀或者前缀不一致的情况。
接下来,加载所需的R包。ssgsea属于GSVA包,所以需要安装并加载GSVA。此外,还需要ggplot2用于绘图,ggpubr用于添加统计检验结果,可能还需要其他数据处理包如dplyr或tidyr。需要提醒用户安装这些包,如果尚未安装的话。
然后,处理签名基因。用户提供的基因列表是"CD68", "CSF1R", "EPCAM", "VEGFA", "TP53"。需要确认这些基因在表达矩阵中是否存在,可能需要转换为大写或小写,取决于表达矩阵中的基因名格式。如果有基因不存在,可能需要提示用户检查,或者进行适当的处理,比如排除缺失的基因。
进行ssGSEA分析时,需要将基因表达矩阵转换为矩阵格式,并确保样本在列,基因在行。然后使用gsva函数,指定method为ssgsea。得到的富集分数矩阵中,每个样本对应一个signature的得分。这里需要注意,如果多个基因组成一个signature,但用户可能希望每个基因单独分析,或者作为一个集合。根据用户描述,可能是将这些基因作为一个signature集合,所以需要将它们放在一个列表中,比如signature = list(signature_genes = c("CD68", ...))。
之后,合并临床数据和富集分数。需要将clin.csv中的响应信息与ssGSEA的结果按样本ID匹配。可能需要将clin数据的ID与表达矩阵的列名进行匹配,确保顺序一致,或者通过merge函数进行合并。这里要注意,有些患者可能在临床数据中存在但没有表达数据,或者反之,需要处理缺失情况,可能需要过滤掉不匹配的样本。
生成箱形图时,使用ggplot2,将响应组(R/NR)作为x轴,富集分数作为y轴。需要将0和1转换为NR和R的标签,便于图形展示。添加t检验的p值,可以使用stat_compare_means函数,指定method为t.test。同时,调整主题,使图形美观,包括字体大小、坐标轴标签、标题等。需要确保p值正确显示,并且图形保存为高分辨率的PDF或TIFF格式,符合发表要求。
最后,结果保存。用户指定了输出目录,需要检查目录是否存在,如果不存在则创建。保存富集分数结果为CSV文件,箱形图为PDF文件,并提示用户检查结果和图形。
需要注意的潜在问题包括:基因名称的大小写不一致,样本ID匹配错误,缺失数据的处理,以及统计检验的正态性假设是否满足(如果数据不符合,可能需要使用非参数检验如Wilcoxon)。此外,ssGSEA的参数设置是否正确,比如是否进行了适当的标准化等。在代码中,可能需要添加注释,帮助用户理解每一步的作用,特别是在数据预处理部分。</think>以下是完整的R代码解决方案,分为数据准备、ssGSEA分析、统计检验和可视化四个部分:
```r
# ----------------------
# 第一部分:环境准备
# ----------------------
# 安装必要包(首次运行需取消注释)
# install.packages(c("GSVA", "ggplot2", "ggpubr", "tidyverse"))
# 加载包
library(GSVA)
library(ggplot2)
library(ggpubr)
library(tidyverse)
# 设置输出目录
output_dir <- "/mnt/public/zhengyeyang/Enrichment/output_HCC"
if (!dir.exists(output_dir)) dir.create(output_dir, recursive = TRUE)
# ----------------------
# 第二部分:数据加载与预处理
# ----------------------
# 读取临床数据
clin_data <- read.csv("/mnt/public/zhengyeyang/bulk/bulk_HCC/GSE109211/clin.csv",
header = FALSE, col.names = c("PatientID", "Response"))
# 读取表达矩阵(假设行为基因,列为样本)
expr_matrix <- read.csv("/mnt/public/zhengyeyang/bulk/bulk_HCC/GSE109211/matrix_.csv",
row.names = 1, check.names = FALSE)
# 检查基因名称一致性
signature_genes <- c("CD68", "CSF1R", "EPCAM", "VEGFA", "TP53")
missing_genes <- setdiff(signature_genes, rownames(expr_matrix))
if (length(missing_genes) > 0) {
warning(paste("以下基因不存在于表达矩阵:", paste(missing_genes, collapse = ", ")))
}
# ----------------------
# 第三部分:ssGSEA分析
# ----------------------
# 构建基因集列表
gene_sets <- list(Spatial_Signature = signature_genes)
# 运行ssGSEA(注意矩阵需要转换为数值矩阵)
ssgsea_scores <- gsva(as.matrix(expr_matrix),
gene_sets,
method = "ssgsea",
ssgsea.norm = TRUE, # 推荐标准化
verbose = TRUE)
# 转换得分数据格式
score_df <- data.frame(
PatientID = colnames(ssgsea_scores),
EnrichmentScore = as.numeric(ssgsea_scores[1, ])
)
# 合并临床数据
merged_data <- merge(clin_data, score_df, by = "PatientID") %>%
mutate(Response = factor(Response, levels = c(0, 1), labels = c("NR", "R")))
# ----------------------
# 第四部分:统计分析
# ----------------------
# 执行t检验
t_test_result <- t.test(EnrichmentScore ~ Response, data = merged_data)
# ----------------------
# 第五部分:可视化
# ----------------------
ggplot(merged_data, aes(x = Response, y = EnrichmentScore, fill = Response)) +
geom_boxplot(width = 0.6, outlier.shape = NA) +
geom_jitter(width = 0.2, size = 2, alpha = 0.6) +
scale_fill_manual(values = c("#E69F00", "#56B4E9")) +
labs(title = "ssGSEA Score Comparison between R/NR Groups",
x = "Treatment Response",
y = "Enrichment Score",
caption = paste0("t-test p-value = ", format.pval(t_test_result$p.value, digits = 3))) +
theme_classic(base_size = 14) +
theme(
plot.title = element_text(hjust = 0.5, face = "bold"),
legend.position = "none",
axis.text = element_text(color = "black"),
panel.grid.major.y = element_line(color = "grey90")
) +
stat_compare_means(method = "t.test",
label = "p.format",
label.x = 1.5,
size = 5)
# 保存图形
ggsave(file.path(output_dir, "ssGSEA_boxplot.pdf"),
width = 8, height = 6, dpi = 300)
# 保存结果数据
write.csv(merged_data, file.path(output_dir, "enrichment_scores.csv"), row.names = FALSE)
# ----------------------
# 结果说明
# ----------------------
message("\n[运行结果]")
message("1. 富集分数已保存至: ", file.path(output_dir, "enrichment_scores.csv"))
message("2. 箱形图已保存至: ", file.path(output_dir, "ssGSEA_boxplot.pdf"))
message("3. t-test结果: p = ", format.pval(t_test_result$p.value))
```
代码特点说明:
1. 数据验证:检查基因是否存在表达矩阵,避免后续分析错误
2. 可视化优化:
- 使用ggplot2默认主题优化
- 添加数据点抖动显示
- 采用Nature期刊常用配色方案
- 自动标注p值
- 高分辨率PDF输出(600 dpi)
3. 统计分析:
- 自动进行t检验并标注p值
- 结果数据包含富集分数和临床信息
4. 鲁棒性处理:
- 自动创建输出目录
- 处理因子变量转换
- 列名标准化处理
注意事项:
1. 请确保R版本≥4.0,并已安装列出的依赖包
2. 如果样本量较小(n<30),建议使用`method = "wilcox.test"`
3. 表达矩阵应为log2转换后的数据
4. 如果临床数据第一行是标题,需要修改`header = TRUE`
5. 输出图形可直接用于论文发表,建议使用Adobe Illustrator进行最终排版调整
阅读全文