单细胞SCENIC分析实战从R代码到可视化热图一步步教你搞定如果你刚刚拿到单细胞转录组数据看着那些代表细胞命运的基因表达矩阵是不是总想再挖深一点想知道是哪些关键的转录因子在幕后“发号施令”决定了细胞的身份和命运SCENIC分析就是为你解答这个问题而生的利器。它不仅仅是一个分析流程更像是一个侦探工具能从海量的单细胞数据中系统地推断出活跃的基因调控网络和细胞状态。对于生物信息学新手或是希望将分析深度从“细胞类型注释”推进到“调控机制解析”的研究者来说掌握SCENIC意味着你能从“看到是什么”升级到“理解为什么”。这篇文章不会重复那些你已经能在官方文档里找到的理论概述。我们将直接切入实战假设你已经完成了基本的单细胞预处理比如使用Seurat得到了一个干净的细胞-基因表达矩阵或者一个标准的Seurat对象。接下来我们将手把手地走过从准备数据、运行核心分析到生成具有生物学意义的可视化热图的完整路径。整个过程会聚焦于R语言环境我会分享一些官方教程里未必会提到的参数调整细节和避坑指南确保你能复现并理解每一步背后的逻辑。1. 分析前的环境搭建与数据准备在开始任何分析之前一个稳定、兼容的环境是成功的基石。SCENIC分析涉及多个R包它们之间的版本依赖有时会比较“挑剔”。直接安装最新版不一定是最佳选择我更推荐建立一个可复现的R环境。1.1 创建与管理分析环境我强烈建议使用renv或conda来管理你的R环境。这里以renv为例它能为你当前的项目创建一个独立的库避免与其他项目的包版本冲突。# 在RStudio中初始化一个新的renv环境 install.packages(renv) renv::init() # 初始化后R会重启。然后安装SCENIC及相关核心包 # 首先设置Bioconductor镜像如果需要 if (!require(BiocManager, quietly TRUE)) install.packages(BiocManager) BiocManager::install(version 3.18) # 指定一个稳定的Bioconductor版本 # 安装SCENIC及其依赖 BiocManager::install(c(AUCell, RcisTarget, GENIE3)) # 注意SCENIC包本身在Bioconductor上但有时从GitHub安装能获得最新修复 # devtools::install_github(aertslab/SCENIC) install.packages(SCENIC) # 安装可视化相关的关键包 BiocManager::install(ComplexHeatmap) install.packages(c(Seurat, ggplot2, patchwork, dplyr, data.table))安装完成后使用renv::snapshot()将当前所有包的版本信息锁定下来。这样未来你或你的同事在任何机器上运行renv::restore()就能完全复现这个分析环境。1.2 理解输入数据从Seurat对象到SCENIC所需格式SCENIC的核心输入是一个基因表达矩阵行是基因列是细胞。如果你是从Seurat对象开始转换过程很简单。但有几个细节需要注意它们直接影响后续分析的质量。library(Seurat) library(SCENIC) # 假设你的Seurat对象叫 seurat_obj # 1. 提取表达矩阵。推荐使用RNA assay的counts或log-normalized数据。 # 使用log-normalized数据如data槽通常效果更好。 exprMat - as.matrix(GetAssayData(seurat_obj, assay RNA, slot data)) # 检查矩阵维度应该是 基因数 x 细胞数 dim(exprMat) # 2. 准备细胞注释信息可选但对后续分群可视化至关重要。 # 这通常来自你的Seurat聚类结果或已知的细胞类型注释。 cellInfo - data.frame( row.names colnames(seurat_obj), CellType seurat_obj$seurat_clusters, # 或用你定义的细胞类型列如 seurat_obj$celltype nGene colSums(GetAssayData(seurat_obj, slot counts) 0), # 每个细胞的检测基因数 nUMI colSums(GetAssayData(seurat_obj, slot counts)) # 每个细胞的UMI总数 ) head(cellInfo)注意矩阵exprMat中的基因名必须是官方基因符号如“TP53”而不是Ensembl ID如“ENSG00000141510”。如果你的数据是Ensembl ID需要使用biomaRt等工具进行转换。此外确保矩阵中没有全是零的行基因或列细胞极端低表达的基因应在Seurat预处理阶段已被过滤。2. 运行SCENIC核心三部曲GRN推断、顺式调控靶点富集与细胞活性评分SCENIC分析流程可以清晰地分为三个逻辑步骤对应三个核心R函数。理解每一步在做什么比机械地运行代码更重要。2.1 第一步共表达网络与候选调控子推断 (GENIE3/GRNBoost2)这一步的目标是找出转录因子TF和其潜在靶基因之间的共表达关系。SCENIC默认使用GENIE3算法但它也支持更快的GRNBoost2需要通过Python调用。我们这里使用纯R的GENIE3。# 初始化SCENIC对象这是管理整个分析流程的核心容器 scenicOptions - initializeScenic(orghgnc, # 物种hgnc(人), mgi(鼠), dmel(果蝇) dbDircisTarget_databases, # 存放物种特异的motif数据库的目录 datasetTitleMy_SCENIC_Analysis, nCores4) # 根据你的电脑核心数调整 # 将表达矩阵和细胞信息存入SCENIC对象 scenicOptionsinputDatasetInfo$cellInfo - cellInfo saveRDS(scenicOptions, fileint/scenicOptions.Rds) # 保存设置 # 1. 过滤低表达基因并筛选潜在的转录因子 # 这一步会基于表达量过滤基因并从一个内置列表中识别转录因子。 genesKept - geneFiltering(exprMat, scenicOptionsscenicOptions, minCountsPerGene3*.01*ncol(exprMat), # 至少在1%的细胞中表达 minSamplesncol(exprMat)*.01) exprMat_filtered - exprMat[genesKept, ] dim(exprMat_filtered) # 2. 运行GENIE3推断共表达网络 # 这步计算量最大耗时最长。对于超过5000个细胞的数据集考虑在服务器上运行。 runGenie3(exprMat_filtered, scenicOptions)运行runGenie3后会在int目录下生成一个加权共表达网络的文件。这个网络包含了所有TF-基因对的关联强度但其中很多是间接的或假阳性的关联。2.2 第二步基于motif的顺式调控靶点富集分析 (RcisTarget)这一步是SCENIC的精华所在它利用已知的转录因子结合motif数据库对上一步得到的共表达网络进行“精炼”只保留那些TF的靶基因上游存在相应TF结合motif的关联从而得到更可靠的调控子。# 3. 为每个TF构建共表达模块co-expression modules runSCENIC_1_coexNetwork2modules(scenicOptions) # 4. 对每个模块进行motif富集分析识别直接调控靶点 # 你需要提前下载对应物种的motif数据库如hg38__refseq-r80__500bp_up_and_100bp_down_tss.mc9nr.feather # 并放在 cisTarget_databases 目录下。 runSCENIC_2_createRegulons(scenicOptions) # 可选但推荐使用“扩展调控子”Extended regulons # 这会将高置信度的直接靶点motif支持和其共表达的基因都纳入调控子通常生物学意义更强。 runSCENIC_3_scoreCells(scenicOptions, exprMat_filtered, skipHeatmapFALSE, skipTsneFALSE)执行完第二步后你会得到一系列“调控子”的定义。每个调控子以一个转录因子命名如“STAT3()”包含一组被该TF直接或间接调控的基因。括号里的“”号表示这是一个扩展调控子。2.3 第三步使用AUCell评估每个细胞的调控子活性最后一步我们要计算每个细胞中每个调控子的活性分数。SCENIC使用AUCell算法其原理是看一个调控子中的基因在单个细胞表达量排序中的富集程度。# 5. 计算每个细胞中每个调控子的AUC值活性分数 # 这一步会生成一个矩阵行为调控子列为细胞值为AUC分数。 scenicOptionsfileNames$output[AUC] - int/3.4_regulonAUC.Rds # 可以指定输出文件名 runSCENIC_4_aucell_binarize(scenicOptions, exprMat_filtered)实际上runSCENIC_3_scoreCells已经包含了AUCell计算。第三步的核心输出是regulonAUC: 一个矩阵存储每个细胞每个调控子的连续AUC值。regulonActivity_binary: 一个二值化矩阵基于细胞特异的阈值判断调控子在每个细胞中是否“活跃”。至此核心计算部分已完成。你得到了每个细胞的调控子活性图谱这是后续所有可视化和生物学解释的基础。3. 结果解读与高级可视化从热图中挖掘生物学故事拿到AUC矩阵后如何把它变成一张能讲故事的图单纯画热图很容易但如何让热图突出关键信息则需要一些策略。3.1 计算与可视化Regulon特异性分数Regulon Specificity Score (RSS) 是一个非常有用的指标用于量化一个调控子对某一细胞类型的特异程度。RSS值越高说明该调控子在该细胞类型中越特异、越可能发挥关键作用。library(AUCell) library(ComplexHeatmap) library(circlize) # 加载之前计算好的AUC矩阵和细胞注释 regulonAUC - readRDS(int/3.4_regulonAUC.Rds) cellInfo - readRDS(int/cellInfo.Rds) # 假设已保存 # 计算RSS rss - calcRSS(AUCgetAUC(regulonAUC), cellAnnotationcellInfo[colnames(regulonAUC), CellType]) # 可视化RSSZ-score热图 # 我们通常只展示每个细胞类型中Top N的特异调控子 rss_plot - plotRSS(rss, zThreshold 1.5, # Z值阈值通常1.5或2表示显著特异 cluster_columns TRUE, # 对细胞类型聚类 order_rows TRUE, # 对调控子排序 thr 0.01, # 过滤掉所有细胞类型中RSS最大值小于此值的调控子 varName CellType, col.low navyblue, col.mid white, col.high darkred, verbose TRUE) print(rss_plot$plot)这张热图能一目了然地告诉你例如“NEUROD1”这个调控子特异性地在神经元细胞中高活性而“SPI1”在髓系细胞中活跃。这是从全局视角发现关键调控因子的最快方法。3.2 构建平均活性热图与差异活性分析除了看特异性我们还想看不同细胞类型间调控子活性的整体模式。计算每个细胞类型内调控子活性的平均值并进行标准化是另一种常见的可视化方式。# 按细胞类型分组计算平均AUC cellTypes - unique(cellInfo$CellType) regulonActivity_byGroup - sapply(cellTypes, function(ct) { cells - rownames(cellInfo)[cellInfo$CellType ct] rowMeans(getAUC(regulonAUC)[, cells, dropFALSE]) }) rownames(regulonActivity_byGroup) - rownames(regulonAUC) # 行方向标准化Z-score使每个调控子在不同组间的活性具有可比性 regulonActivity_byGroup_Scaled - t(scale(t(regulonActivity_byGroup))) # 使用ComplexHeatmap绘制 Heatmap(regulonActivity_byGroup_Scaled, name Regulon activity\n(z-score), col colorRamp2(c(-2, 0, 2), c(blue, white, red)), row_title Regulons, column_title Cell Types, row_names_gp gpar(fontsize 8), column_names_gp gpar(fontsize 10), clustering_distance_rows euclidean, clustering_method_rows ward.D2, cluster_columns TRUE, show_row_dend TRUE, show_column_dend TRUE, heatmap_legend_param list(title_gp gpar(fontsize 10)) )为了从这张热图中提取更有洞察力的信息我们可以进行简单的差异活性分析找出在每个细胞类型中活性最特异的调控子。# 一个简单的差异活性筛选方法寻找在某一类中活性远高于其他类平均活性的调控子 top_regulons_per_type - list() for (ct in colnames(regulonActivity_byGroup_Scaled)) { # 计算当前类型与其他类型平均值的差值 diff_score - regulonActivity_byGroup_Scaled[, ct] - rowMeans(regulonActivity_byGroup_Scaled[, setdiff(colnames(regulonActivity_byGroup_Scaled), ct), dropFALSE]) # 选取差值最大的前5个调控子 top5 - names(sort(diff_score, decreasing TRUE))[1:5] top_regulons_per_type[[ct]] - top5 } # 去重后用于绘制聚焦的热图 unique_top_regulons - unique(unlist(top_regulons_per_type)) focused_heatmap_data - regulonActivity_byGroup_Scaled[unique_top_regulons, ]3.3 在细胞嵌入图上可视化特定调控子活性热图展示了群体水平的信息但我们还需要在单细胞分辨率上查看活性分布。这可以通过将AUC值映射到你的UMAP或t-SNE图上实现。library(Seurat) library(ggplot2) library(patchwork) # 将调控子AUC值添加到Seurat对象的metadata中方便用Seurat函数绘图 # 选取几个感兴趣的调控子例如从RSS分析中得到的关键因子 regulons_to_plot - c(ASCL1(), OLIG2(), SPI1()) # 确保它们存在于AUC矩阵中 regulons_to_plot - regulons_to_plot[regulons_to_plot %in% rownames(regulonAUC)] # 提取这些调控子的AUC值列为细胞 auc_for_plot - t(getAUC(regulonAUC)[regulons_to_plot, , dropFALSE]) # 将AUC值添加到Seurat对象 seurat_obj - AddMetaData(seurat_obj, as.data.frame(auc_for_plot)) # 使用FeaturePlot在UMAP上可视化活性 p_list - list() for (reg in regulons_to_plot) { p - FeaturePlot(seurat_obj, features reg, reduction umap, cols c(lightgrey, darkblue), order TRUE) ggtitle(reg) theme(legend.position right) p_list[[reg]] - p } wrap_plots(p_list, ncol 2)此外VlnPlot和DotPlot也能很好地展示调控子活性在不同细胞群中的分布和平均强度。# 组合图点图展示平均表达和表达比例小提琴图展示分布 dot_p - DotPlot(seurat_obj, features regulons_to_plot, group.by CellType) RotatedAxis() scale_color_gradient2(lowgrey, midyellow, highred) vln_p - VlnPlot(seurat_obj, features regulons_to_plot, group.by CellType, pt.size 0, ncol 1) dot_p | vln_p4. 实战技巧、常见问题与下游整合思路掌握了标准流程后一些实战技巧能帮你提升分析效率和结果的可靠性。这里分享几个我踩过坑后才明白的关键点。4.1 参数调整与性能优化基因过滤的松紧度geneFiltering中的minCountsPerGene和minSamples参数至关重要。过滤太严会丢失稀有但重要的转录因子过滤太松则会极大增加GENIE3的计算负担和噪音。对于10x Genomics数据从3*.01*ncol(exprMat)开始是个不错的默认值但你可以根据细胞数和测序深度微调。GENIE3的并行计算对于超过1万个细胞的数据集GENIE3可能运行数天。确保正确设置scenicOptionssettings$nCores。如果资源允许可以考虑使用GRNBoost2通过SCENIC的Python接口它通常更快且内存效率更高。AUCell的二值化阈值runSCENIC_4_aucell_binarize中默认使用自动阈值。但有时对于活性分布双峰不明显的调控子自动阈值可能不理想。你可以通过AUCell_exploreThresholds函数手动检查并调整关键调控子的阈值。4.2 解读结果时的注意事项调控子命名的“()”带有“()”的扩展调控子包含直接靶点和高度共表达的基因通常更稳定生物学解释性更强。在分析时可以优先关注扩展调控子。AUC值不是表达量调控子AUC值高不代表其核心TF的表达量一定高。AUC反映的是整个调控子基因集的协同活性。有时一个TF自身表达量低但其靶基因集却显示出高度协调的活性这同样具有重要生物学意义。结合已知生物学知识SCENIC是计算预测工具其结果尤其是预测的靶基因需要与已有文献或ChIP-seq等实验数据交叉验证。不要盲目相信所有预测的调控关系。4.3 将SCENIC结果整合到更广泛的分析中SCENIC的输出可以成为下游分析的强大输入伪时间分析将关键调控子的AUC值作为输入与伪时间轨迹一起分析可以发现驱动细胞分化的动态调控程序。细胞通讯分析将高活性的转录因子与配体-受体对关联可以推测哪些调控子可能影响了细胞间的信号传递。与差异表达基因结合比较差异表达基因列表和关键调控子的靶基因集可以验证调控子是否确实驱动了观察到的表达变化。最后记得妥善保存你的中间结果和最终参数。SCENIC分析流程较长保存好scenicOptions对象和关键的RDS文件如regulonAUC能让你在未来轻松地重新生成图表或进行二次分析。整个流程看似代码繁多但一旦跑通一次你就会发现它其实是一条逻辑清晰的流水线。最重要的是它赋予了你从单细胞数据中窥视基因调控逻辑的能力这是迈向机制性理解的关键一步。