保姆级教程：从Seurat对象到热图，用DecoupleR完整跑通单细胞转录因子活性分析

单细胞转录因子活性分析实战从Seurat到DecoupleR的完整解决方案在单细胞转录组数据分析中我们常常需要超越基因表达层面探索调控这些表达模式的幕后推手——转录因子(TF)的活性状态。传统方法局限于观察TF基因本身的表达水平而现代计算生物学提供了更精密的工具来推断TF的功能活性。本文将手把手带你使用R语言生态中的Seurat和decoupleR工具包构建一套可重复的转录因子活性分析流程。1. 环境准备与数据加载1.1 软件包安装与配置工欲善其事必先利其器。我们需要确保所有必要的R包都已正确安装。建议在RStudio中创建一个新项目专门用于本次分析# 安装CRAN来源的依赖包 install.packages(c(Seurat, dplyr, tibble, tidyr, patchwork, ggplot2, pheatmap)) # 安装Bioconductor依赖 if (!require(BiocManager, quietly TRUE)) install.packages(BiocManager) BiocManager::install(GSVA) # 安装decoupleR开发版 devtools::install_github(saezlab/decoupleR)安装完成后加载所有必要的库library(Seurat) library(decoupleR) library(dplyr) library(tibble) library(tidyr) library(patchwork) library(ggplot2) library(pheatmap)1.2 示例数据加载我们使用GSE148190数据集中的SM4455931-10X样本作为演示。假设数据已经完成了基本的质控、标准化和聚类步骤load(step1.final.Rdata) sce - step1.final检查细胞聚类情况table(Idents(sce)) DimPlot(sce, label TRUE, repel TRUE)2. 转录因子调控网络获取与理解2.1 CollecTRI网络简介CollecTRI是目前最全面的转录因子-靶基因互作数据库之一整合了多个来源的调控证据包含1,186个转录因子涵盖43,178个TF-基因互作关系每个互作都有置信度评分(mor值)net - get_collectri(organism human, split_complexes FALSE) head(net)提示split_complexes参数决定是否将蛋白质复合物拆分为单个组分。设为FALSE时保留复合物的原始形式TRUE则分解为亚基。2.2 网络质量评估在实际分析前建议检查网络的基本统计特征network_stats - net %% group_by(source) %% summarise(targets n_distinct(target), avg_mor mean(mor)) summary(network_stats$targets) summary(network_stats$avg_mor)3. 单细胞转录因子活性计算3.1 数据矩阵准备从Seurat对象中提取标准化后的表达矩阵mat - as.matrix(sceassays$RNAdata) dim(mat)3.2 并行计算设置ULM(单变量线性模型)计算较为耗时建议启用多线程library(future) plan(multisession, workers 4) # 根据CPU核心数调整3.3 运行ULM分析核心计算步骤通常需要5-30分钟不等acts - run_ulm(mat, net, minsize 5, .source source, .target target, .mor mor)理解输出结构source: 转录因子名称condition: 细胞barcodescore: 活性分数(t值)p_value: 统计显著性4. 结果整合与可视化4.1 结果回填到Seurat对象将计算结果整合回Seurat对象便于后续分析sce[[tfs_ulm]] - acts %% pivot_wider(id_cols source, names_from condition, values_from score) %% column_to_rownames(source) %% CreateAssayObject() DefaultAssay(sce) - tfs_ulm sce - ScaleData(sce) sceassays$tfs_ulmdata - sceassays$tfs_ulmscale.data4.2 单转录因子活性可视化比较转录因子活性与表达模式的差异DefaultAssay(sce) - tfs_ulm p_activity - FeaturePlot(sce, features E2F1) scale_colour_gradient2(low blue, mid white, high red) ggtitle(E2F1 activity) DefaultAssay(sce) - RNA p_expression - FeaturePlot(sce, features E2F1) ggtitle(E2F1 expression) p_activity | p_expression4.3 聚类水平活性热图识别各细胞群中特异性活跃的转录因子# 提取前25个变异最大的TF n_tfs - 25 df - t(as.matrix(sceassays$tfs_ulmdata)) %% as.data.frame() %% mutate(cluster Idents(sce)) %% pivot_longer(cols -cluster, names_to source, values_to score) %% group_by(cluster, source) %% summarise(mean mean(score)) tfs - df %% group_by(source) %% summarise(std sd(mean)) %% arrange(-abs(std)) %% head(n_tfs) %% pull(source) # 构建热图矩阵 top_acts_mat - df %% filter(source %in% tfs) %% pivot_wider(id_cols cluster, names_from source, values_from mean) %% column_to_rownames(cluster) %% as.matrix() # 自定义颜色标度 palette_length - 100 my_color - colorRampPalette(c(Darkblue, white, red))(palette_length) my_breaks - c(seq(-3, 0, length.out ceiling(palette_length/2) 1), seq(0.05, 3, length.out floor(palette_length/2))) pheatmap(top_acts_mat, border_color NA, color my_color, breaks my_breaks, cluster_rows TRUE, cluster_cols TRUE, show_colnames TRUE, fontsize_row 8, fontsize_col 8)5. 高级分析与结果解读5.1 活性-表达相关性分析探索转录因子活性与表达水平的关系tf - E2F1 # 提取活性分数 activity_scores - sceassays$tfs_ulmdata[tf, ] # 提取表达水平 expression_levels - sceassays$RNAdata[tf, ] # 创建数据框 cor_df - data.frame(activity activity_scores, expression expression_levels, cluster Idents(sce)) # 计算全局相关系数 global_cor - cor(cor_df$activity, cor_df$expression) # 分群计算相关系数 cluster_cors - cor_df %% group_by(cluster) %% summarise(correlation cor(activity, expression)) ggplot(cor_df, aes(x expression, y activity, color cluster)) geom_point(alpha 0.5) geom_smooth(method lm, se FALSE) labs(title paste0(Activity vs Expression: , tf), subtitle paste0(Global Pearson r , round(global_cor, 2)), x Expression level (log-normalized), y Activity score (t-value)) theme_minimal()5.2 关键转录因子识别基于活性分数识别各细胞群的关键调控因子# 计算各群特异的TF活性 cluster_specific_tfs - df %% group_by(source) %% mutate(z_score scale(mean)) %% group_by(cluster) %% arrange(desc(z_score)) %% slice_head(n 5) %% ungroup() # 可视化 ggplot(cluster_specific_tfs, aes(x reorder(source, z_score), y z_score, fill cluster)) geom_col() coord_flip() facet_wrap(~cluster, scales free_y) labs(x Transcription Factor, y Z-scored Activity, title Top Active TFs by Cluster) theme_minimal() theme(legend.position none)5.3 功能富集分析将高活性转录因子与通路数据库关联library(clusterProfiler) library(org.Hs.eg.db) # 获取所有显著活跃的TF significant_tfs - df %% group_by(source) %% summarise(max_activity max(abs(mean))) %% filter(max_activity quantile(max_activity, 0.9)) %% pull(source) # 转换基因符号到ENTREZID tf_entrez - bitr(significant_tfs, fromType SYMBOL, toType ENTREZID, OrgDb org.Hs.eg.db) # GO富集分析 go_results - enrichGO(gene tf_entrez$ENTREZID, OrgDb org.Hs.eg.db, keyType ENTREZID, ont BP, pvalueCutoff 0.05, qvalueCutoff 0.1) dotplot(go_results, showCategory 15) ggtitle(GO Biological Processes for Active TFs)