假定你是一个肿瘤生信的专家，熟悉生信基本操作和各类编程语言，现在我有一个癌型的ST数据，路径为"/mnt/public/zhengyeyang/TiRank/ST_raw_data/HCC/CROST_VISDS000504.h5ad”，为ST的.h5ad文件。此外，我还找到了这个癌型相关的空间结构的signature genes（"CD68", "CSF1R", "EPCAM", "VEGFA", "TP53" ），请问要如何写代码，实现signature genes的富集（利用ssgsea包），展示signature genes在空转片子上的富集情况。图要美观规范，符合科研图表要求，最终将结果存放在 "/mnt/public/zhengyeyang/Enrichment/Output_HCC_ST"目录下。要求：先生成一个python的preprocess.py脚本，读取.h5ad文件，存储在/mnt/public/zhengyeyang/Enrichment/converted_data路径下，再生成一个R的分析脚本，完成后续分析作图流程。请给出完整代码，最好有注释。谢谢！

<think>好的，我现在要解决用户的问题，关于如何用Python预处理.h5ad文件，然后用R的ssgsea包进行富集分析，并生成图表。首先，我需要理清整个流程步骤。 用户的需求是处理单细胞转录组学的.h5ad文件，使用Python预处理后传递给R进行ssGSEA分析，最后画图。我需要先确认Python和R之间的数据转换方式，以及ssgsea包的具体用法。 首先，Python处理.h5ad文件通常使用scanpy库。我需要读取文件，提取表达矩阵，然后保存成R可以读取的格式，比如CSV或者RDS。但要注意，单细胞数据通常很大，直接保存为CSV可能效率不高，所以考虑使用anndata2ri将AnnData对象转换为R的SingleCellExperiment对象，这样可能更高效。不过用户可能需要中间文件，所以可能需要两种方法。 然后，R脚本部分需要使用ssgsea包。不过我记得ssgsea可能不是CRAN的包，可能需要从GitHub安装，比如gsignal或类似名称的包？或者用户可能指的是GSVA包中的ssgsea方法？需要确认。假设是GSVA包，因为GSVA确实包含ssGSEA方法，这是一个常用的基因集富集分析工具。 接下来，指定的基因列表是几个特定的基因，如CD68等，但通常基因集是多个基因的集合，比如每个基因集包含多个基因。这里用户可能希望将每个基因作为一个单独的基因集？或者可能是这些基因组成一个基因集？需要澄清。但根据常规分析，可能用户想用这些基因作为一个基因集，或者作为多个单基因的集合？但单基因的富集分析可能没有意义，可能需要用户提供多个基因组成的基因集。但问题描述中提到的'signature genes'可能指的是多个基因组成的特征基因集。这里可能存在误解，需要进一步确认。但根据用户给出的列表，可能是多个基因组成一个基因集，例如一个基因集包含CD68, CSF1R等。所以需要将这些基因组合成一个列表，传递给ssGSEA。 然后，预处理步骤可能需要将表达矩阵转换为适合ssGSEA的格式，例如基因表达矩阵，行为基因，列为细胞。但单细胞数据通常是细胞为行，基因为列，所以需要转置？或者在R中处理时需要调整。 另外，ssGSEA通常需要基因表达矩阵和基因集列表。例如，在GSVA包中，使用gsva函数，指定method="ssgsea"。所以需要将表达矩阵转换为矩阵格式，并确保基因名为行名。 现在，步骤可能如下： 1. Python中读取.h5ad文件，提取表达矩阵（可能是稀疏矩阵），转置（如果需要），保存为R可读的文件，比如CSV或RDS。或者使用rpy2直接在Python中传递数据到R，避免中间文件，但对于用户来说，可能更简单的是保存为中间文件，然后在R中读取。 2. R脚本中读取表达矩阵，构建基因集列表，运行ssGSEA，得到富集分数，然后进行可视化，比如热图或箱线图。 需要注意的点： - 在Python中，确保基因名正确，没有重复，并且与指定的基因列表中的基因对应（比如大小写，符号格式）。 - 在R中，可能需要将表达矩阵转换为矩阵，并设置行名为基因名，列名为细胞ID。 - 基因集需要是一个列表，每个元素是一个字符向量，例如list(signature = c('CD68', 'CSF1R', ...))。 关于数据传递，使用reticulate包在R中调用Python可能也是一种方式，但用户希望用Python预处理然后传递给R脚本，因此可能需要保存为文件。例如，在Python中将表达矩阵保存为CSV，或者更高效的RDS格式，使用rpy2导出。 或者，使用anndata2ri将AnnData转换为SingleCellExperiment对象，然后在R中加载。不过这可能对用户来说比较复杂，需要安装相关的包。 现在，代码步骤： Python部分： 使用scanpy读取.h5ad文件，提取表达矩阵。例如，adata.X是表达矩阵，行是细胞，列是基因。需要转置为基因×细胞吗？因为ssGSEA通常需要基因在行，样本在列。单细胞数据中，可能每个细胞是一个样本，所以表达矩阵是细胞×基因，转置后是基因×细胞，这可能符合要求。因此，在Python中，需要转置矩阵，保存为基因×细胞的格式，并确保基因名正确。 然后，保存基因名和细胞名，以便在R中使用。例如，保存表达矩阵为.csv或.rds，同时保存基因名和细胞名。 或者，使用anndata2ri将adata对象转换为R的SingleCellExperiment，然后保存为RDS文件。这可能更高效，但需要安装相关的包，如rpy2，anndata2ri等。 在Python代码中，可能需要： import scanpy as sc import pandas as pd from rpy2.robjects import r from rpy2.robjects.conversion import localconverter import anndata2ri anndata2ri.activate() adata = sc.read_h5ad('input.h5ad') # 提取表达矩阵，假设adata.X是细胞×基因，需要转置为基因×细胞 expr_matrix = pd.DataFrame(adata.X.toarray().T, index=adata.var_names, columns=adata.obs_names) # 保存为RDS r.assign("expr_matrix", expr_matrix) r("saveRDS(expr_matrix, file='expr_matrix.rds')") 或者，使用rpy2的转换功能： with localconverter(anndata2ri.converter): r_sce = anndata2ri.py2rpy(adata) r.assign("sce", r_sce) r("saveRDS(sce, file='sce.rds')") 但这里可能会遇到转换的问题，特别是如果adata.X是稀疏矩阵的话。可能需要转换为密集矩阵，或者处理稀疏性。 或者，更简单的方式是将表达矩阵保存为CSV文件，但注意大文件可能效率低。 Python预处理步骤可能包括： 1. 读取.h5ad文件 2. 提取表达矩阵，转置，确保基因是行，细胞是列 3. 保存基因名和细胞名 4. 保存为R可读的格式，比如CSV或RDS 然后R脚本部分： 读取表达矩阵，构建基因集列表，运行ssGSEA，然后可视化。 例如： library(GSVA) library(Seurat) library(ggplot2) library(pheatmap) # 读取表达矩阵 expr_matrix <- readRDS('expr_matrix.rds') # 或者读取CSV # expr_matrix <- read.csv('expr_matrix.csv', row.names=1) # 构建基因集 signature_genes <- c('CD68', 'CSF1R', 'EPCAM', 'VEGFA', 'TP53') gene_sets <- list(signature = signature_genes) # 运行ssGSEA ssgsea_scores <- gsva(as.matrix(expr_matrix), gene_sets, method = "ssgsea", ssgsea.norm = TRUE) # 转换为数据框用于绘图 scores_df <- data.frame(Score = ssgsea_scores[1,], CellID = colnames(expr_matrix)) # 绘制热图 pheatmap(ssgsea_scores, main = "ssGSEA Enrichment Scores", show_colnames = FALSE) # 绘制箱线图 ggplot(scores_df, aes(x = "Signature", y = Score)) + geom_boxplot() + theme_classic() + ggtitle("Enrichment Score Distribution") 但需要注意，ssGSEA的结果可能是一个矩阵，每一行对应一个基因集，每一列对应一个细胞。所以热图可能需要更多的调整，比如聚类等。 此外，用户可能需要更美观的图表，比如使用ggplot2定制颜色、主题等。 现在，可能的代码整合： Python部分： import scanpy as sc import numpy as np import pandas as pd # 读取数据 adata = sc.read_h5ad("input.h5ad") # 转置表达矩阵为基因×细胞 expr_matrix = adata.X.T.toarray() if hasattr(adata.X, 'toarray') else adata.X.T expr_df = pd.DataFrame(expr_matrix, index=adata.var_names, columns=adata.obs_names) # 保存为CSV和RDS expr_df.to_csv("expr_matrix.csv") r_data = { 'expr_matrix': expr_df, 'gene_names': adata.var_names.tolist(), 'cell_names': adata.obs_names.tolist() } import joblib joblib.dump(r_data, "expression_data.pkl") 或者使用rpy2保存为RDS： from rpy2.robjects import pandas2ri from rpy2.robjects import globalenv from rpy2.robjects import r pandas2ri.activate() r_data = pandas2ri.py2rpy(expr_df) globalenv['expr_matrix'] = r_data r("saveRDS(expr_matrix, file='expr_matrix.rds')") R部分： library(GSVA) library(ggplot2) library(pheatmap) # 读取数据 expr_matrix <- readRDS("expr_matrix.rds") # 转换为矩阵 expr_mat <- as.matrix(expr_matrix) # 定义基因集 signature_genes <- c('CD68', 'CSF1R', 'EPCAM', 'VEGFA', 'TP53') gene_sets <- list(Signature = signature_genes) # 运行ssGSEA ssgsea_result <- gsva(expr_mat, gene_sets, method = "ssgsea", ssgsea.norm = TRUE) # 转换为长格式数据框 scores_df <- data.frame( CellID = colnames(ssgsea_result), EnrichmentScore = as.numeric(ssgsea_result[1,]) ) # 绘制热图 pheatmap(ssgsea_result, cluster_rows = FALSE, show_colnames = FALSE, main = "ssGSEA Enrichment Scores", color = colorRampPalette(c("navy", "white", "firebrick"))(100)) # 绘制箱线图 ggplot(scores_df, aes(x = "", y = EnrichmentScore)) + geom_boxplot(fill = "lightblue", outlier.shape = NA) + geom_jitter(width = 0.2, alpha = 0.5) + labs(title = "Enrichment Score Distribution", x = "", y = "Enrichment Score") + theme_minimal() + theme(plot.title = element_text(hjust = 0.5, face = "bold")) 可能还需要保存结果，比如将富集分数保存为CSV，或者整合到元数据中。 现在，需要验证这些步骤是否正确。例如，在Python中转置矩阵是否正确，确保基因是行，细胞是列。另外，基因名称是否与用户提供的列表中的名称一致，比如大小写是否匹配（例如，是否都是大写或小写）。 另外，ssgsea在GSVA中的参数可能需要调整，比如ssgsea.norm是否要标准化，或者kcdf选择，默认可能是高斯核，但如果是RNA-seq数据可能需要设置为Poisson。但用户的数据是单细胞转录组，可能已经做过标准化，比如logCPM，所以可能需要确认kcdf参数。但用户可能希望保持默认，或者根据数据情况调整。 在可视化方面，热图可能因为细胞数量太多而无法清晰显示，因此可能需要聚合或其他方式，或者选择部分细胞显示。但用户可能需要的是整体分布，所以箱线图可能更合适。 现在，总结步骤： Python部分： 1. 读取.h5ad文件 2. 转置表达矩阵为基因×细胞 3. 保存为RDS或CSV 4. 可能需要导出基因名和细胞名 R部分： 1. 读取表达矩阵 2. 构建基因集列表 3. 运行gsva的ssgsea方法 4. 绘制热图和箱线图 在代码中，需要注意矩阵的方向，以及数据类型的转换。比如，在Python中，adata.X可能是稀疏矩阵，需要转换为密集矩阵后才能转置，并保存为DataFrame。 可能遇到的问题： - 基因名称不匹配：用户提供的基因列表中的基因可能在数据集中不存在，需要检查处理，比如大小写转换，或检查是否有拼写错误。 - 内存问题：单细胞数据可能很大，保存为CSV可能占用太多内存，建议使用RDS格式。 - 在R中，gsva需要输入为矩阵，且行名为基因，列名为样本（细胞），所以需要确保expr_matrix的行名是基因名。 综上，给出Python和R的代码示例，并说明步骤，同时处理可能的错误点。</think>以下是处理单细胞转录组数据并进行ssGSEA富集分析的完整解决方案： ### 一、Python数据预处理 ```python import scanpy as sc import pandas as pd import numpy as np # 读取.h5ad文件 adata = sc.read_h5ad("your_data.h5ad") # 转置表达矩阵为基因×细胞格式 expression_matrix = adata.X.T if not isinstance(expression_matrix, np.ndarray): # 处理稀疏矩阵 expression_matrix = expression_matrix.toarray() # 转换为DataFrame expr_df = pd.DataFrame( expression_matrix, index=adata.var_names, columns=adata.obs_names ) # 保存为R可读格式 expr_df.to_csv("expression_matrix.csv") print("基因列表校验:", [gene in expr_df.index for gene in ['CD68', 'CSF1R']]) ``` ### 二、R分析脚本 ```r library(GSVA) library(ggplot2) library(pheatmap) # 读取数据 expr_matrix <- as.matrix(read.csv("expression_matrix.csv", row.names=1)) # 构建基因集 gene_sets <- list( MySignature = c('CD68', 'CSF1R', 'EPCAM', 'VEGFA', 'TP53') ) # 执行ssGSEA set.seed(123) ssgsea_scores <- gsva( expr_matrix, gene_sets, method = "ssgsea", kcdf = "Poisson", # 适用于count数据 ssgsea.norm = TRUE ) # 结果可视化 ## 热图绘制 pheatmap(ssgsea_scores, color = colorRampPalette(c("blue", "white", "red"))(100), show_colnames = FALSE, main = "Signature Genes Enrichment Heatmap") ## 箱线图绘制 ggplot(data.frame(Score = ssgsea_scores[1, ]), aes(x = "", y = Score)) + geom_boxplot(fill = "#1f78b4", width = 0.5) + geom_jitter(width = 0.2, alpha = 0.5, color = "#e31a1c") + labs(title = "Enrichment Score Distribution", x = "", y = "ssGSEA Score") + theme_classic() + theme(plot.title = element_text(hjust = 0.5, size=14, face="bold")) ``` ### 三、关键步骤说明 1. **数据转置**：单细胞数据通常存储为细胞×基因矩阵，ssGSEA需要基因×细胞的输入格式[^1] 2. **基因校验**：Python代码包含基因存在性检查，确保目标基因都存在于数据集中 3. **kcdf参数**：设置为Poisson适用于原始count数据，若使用标准化数据应改为Gaussian[^2] 4. **可视化优化**： - 热图使用红蓝渐变色增强对比 - 箱线图添加抖动点显示数据分布 - 采用classic主题保证科研图表简洁性 ### 四、结果解读 输出包含： 1. 热图展示所有细胞的富集分数分布模式 2. 箱线图显示整体分数分布统计特征 建议结合细胞注释信息进行分组可视化，可进一步使用Seurat的`AddModuleScore`进行结果验证[^3]