2026版单细胞测序scRNA-seq入门实战（十二）：pySCENIC转录因子调控网络分析——从上游数据准备到下游可视化的完整闭环

本节概览：

pySCENIC简介：理解转录因子调控网络分析的意义和pySCENIC的核心优势。
环境配置：Conda环境创建与数据库文件准备。
数据准备：从Seurat对象导出pySCENIC所需的loom文件。
上游分析：GRNBoost2 → cisTarget → AUCell 三步完整Shell脚本。
下游可视化：RSS计算、Regulon活性热图、UMAP叠加图的完整R代码。
结果解读：从AUC矩阵到生物学故事的完整演绎。

在之前的推文中，我们完成了细胞类型注释、差异表达分析、功能富集和细胞通讯分析，但还有一个更深层的问题需要解答：是谁在驱动这些基因表达的变化？

转录因子（Transcription Factor, TF）是能够与特定DNA序列结合、调控下游基因转录的关键蛋白质。SCENIC（Single-Cell rEgulatory Network Inference and Clustering）正是为了解决这个问题而设计的工具。

本文将提供一个从上游数据准备到下游可视化的完整闭环流程，使用pySCENIC进行转录因子调控网络分析。

一、pySCENIC vs R SCENIC：为什么选择Python版？

对比项	R SCENIC	pySCENIC
GRN算法	GENIE3（随机森林）	GRNBoost2（梯度提升）
计算速度	慢（数天）	快（数小时）
内存效率	较低	较高
并行支持	有限	多核/分布式
适用规模	< 5000细胞	> 5000细胞，可达十万级

结论：对于大多数实际项目（细胞数 > 3000），强烈推荐使用pySCENIC。

二、核心概念速览

什么是Regulon？

一个 Regulon = 一个转录因子（TF） + 其调控的一组靶基因。

什么是AUC score？

AUCell算法评估每个regulon在每个细胞中的活跃程度。AUC值与转录因子本身的表达量是不同的概念——即使TF因dropout而为零，其regulon活性仍然可能很高。

什么是RSS？

Regulon Specificity Score衡量某个regulon在不同细胞类型中的分布是否集中，用于筛选细胞类型特异性转录因子。

三、环境配置：Conda一键安装

# 创建专属conda环境
conda create -y -n pyscenic python=3.7
conda activate pyscenic

# 安装pySCENIC及依赖
pip install pyscenic loompy scanpy matplotlib seaborn

# 验证安装
pyscenic -h

四、数据库准备

# 创建数据库目录
mkdir -p ~/reference/cisTarget_databases/mouse
cd ~/reference/cisTarget_databases/mouse

# 下载小鼠数据库文件（以mm10为例）
# 1. 转录因子列表
wget https://resources.aertslab.org/cistarget/tf_lists/allTFs_mm.txt

# 2. Motif排名数据库（约1 GB）
wget https://resources.aertslab.org/cistarget/databases/mus_musculus/mm10/refseq_r80/mc_v10_clust/gene_based/mm10_10kbp_up_10kbp_down_full_tx_v10_clust.genes_vs_motifs.rankings.feather

# 3. Motif注释文件
wget https://resources.aertslab.org/cistarget/motif2tf/motifs-v10nr_clust-nr.mgi-m0.001-o0.0.tbl

五、数据准备：从Seurat导出loom文件

在R中准备数据，随机抽取200-500个细胞进行探索性分析：

rm(list = ls())
options(stringsAsFactors = FALSE)
library(Seurat)
library(loomR)
library(tidyverse)
set.seed(12345)

# 读取 Seurat 对象（根据实际文件调整）
seurat_obj <- readRDS("4_scRNA_celltype.rds") 

# （可选）随机取子集，减少数据量
set.seed(12345)
n_cells <- 500   # 自行调整
cells_sub <- sample(colnames(seurat_obj), size = min(n_cells, ncol(seurat_obj)), replace = FALSE)
seurat_sub <- seurat_obj[, cells_sub]

# 提取 counts 矩阵（基因 × 细胞）
counts_mat <- GetAssayData(seurat_sub, assay = "RNA", layer = "counts")  
# 转置为 细胞 × 基因（scanpy 默认期望）
expr_df <- as.data.frame(t(as.matrix(counts_mat)))

# 检查行名（细胞名）和列名（基因名）是否完整
stopifnot(!is.null(rownames(expr_df)), !is.null(colnames(expr_df)))

# 导出 CSV（建议用 data.table 加速大文件）
data.table::fwrite(expr_df, file = "scrna_exp.csv", row.names = TRUE, quote = FALSE)

六、pySCENIC上游分析：

这是pySCENIC的核心分析流程，使用GRNBoost2、cisTarget和AUCell完成三步分析：

#!/bin/bash
# pySCENIC 完整分析脚本
# 功能：GRN推断 → Motif富集 → Regulon活性评分

date
echo "################### Import DATA ##################"

# ---- 运行change.py（如需对loom文件做额外处理） ----
python change.py
# $ cat change.py 
# import os, sys
# os.getcwd()
# os.listdir(os.getcwd()) 
# import loompy as lp;
# import numpy as np;
# import scanpy as sc;
# x=sc.read_csv("scrna_exp.csv");
# row_attrs = {"Gene": np.array(x.var_names),};
# col_attrs = {"CellID": np.array(x.obs_names)};
# lp.create("scrna.loom",x.X.transpose(),row_attrs,col_attrs);

date

# ---- 设置数据库路径 ----
dir=/home/lilab/reference/index/cisTarget_databases/mouse
tfs=$dir/allTFs_mm.txt
feather=$dir/mm10_10kbp_up_10kbp_down_full_tx_v10_clust.genes_vs_motifs.rankings.feather
tbl=$dir/motifs-v10nr_clust-nr.mgi-m0.001-o0.0.tbl

input_loom=./scrna.loom

# 检查文件是否存在
ls $tfs $feather $tbl

# Step 1: GRN推断（最耗时步骤）

echo "###################  1. GRN (Gene Regulatory Network)  ####################"

pyscenic grn \
    --num_workers 4 \
    --output grn.tsv \
    --method grnboost2 \
    $input_loom \
    $tfs

# Step 2: cisTarget Motif富集（修剪为可靠Regulon）

echo "#################### 2. cisTarget (Motif Enrichment) #################"

pyscenic ctx \
    grn.tsv \
    $feather \
    --output reg.csv \
    --annotations_fname $tbl \
    --expression_mtx_fname $input_loom \
    --mode "dask_multiprocessing" \
    --mask_dropouts \
    --num_workers 4

# Step 3: AUCell活性评分

echo "#################### 3. AUCell (Regulon Activity Scoring) ###################"

pyscenic aucell \
    $input_loom \
    reg.csv \
    --output out_SCENIC.loom \
    --num_workers 4

echo "##############  All pySCENIC Work Done!!!!! #######"
date

脚本参数详解

步骤	参数	说明
GRN	`--num_workers 4`	并行核心数，根据机器配置调整
	`--method grnboost2`	推荐使用GRNBoost2（梯度提升），也可用genie3
	`--output grn.tsv`	输出TF-靶基因共表达关系文件
cisTarget	`--mask_dropouts`	处理dropout事件，提高鲁棒性
	`--mode dask_multiprocessing`	并行模式
	`--output reg.csv`	输出Regulon列表
AUCell	`--output out_SCENIC.loom`	输出包含AUC得分的loom文件

七、下游可视化：R中加载结果并绘图

7.1 加载SCENIC结果

# 第二部分：R环境 - 结果可视化
rm(list = ls())
library(Seurat)
library(SCopeLoomR)
library(AUCell)
library(SCENIC)
library(tidyverse)
library(ComplexHeatmap)
library(circlize)
library(RColorBrewer)
library(patchwork)

# 读取原始Seurat对象（用于UMAP坐标）
seurat_obj <- readRDS("4_scRNA_celltype.rds")

# ---- 打开pySCENIC输出的loom文件 ----
loom <- open_loom('C:/Users/LX/Desktop/seq/scRNA-seq/GSE108761/SCENIC/out_SCENIC.loom')

# 读取Regulon矩阵
regulons_incidMat <- get_regulons(loom, column.attr.name = "Regulons")
regulons <- regulonsToGeneLists(regulons_incidMat)

# 读取AUC矩阵
regulonAUC <- get_regulons_AUC(loom, column.attr.name = 'RegulonsAUC')

# 关闭loom文件
close_loom(loom)

# 查看基本信息
dim(regulonAUC)  # Regulon数 × 细胞数
names(regulons)[1:5]

7.2 RSS（Regulon Specificity Score）计算

RSS用于筛选细胞类型特异性的转录因子：

# ---- 准备细胞注释 ----
# 从 seurat_obj 提取 cell_type 和 seurat_clusters
cellinfo <- seurat_obj@meta.data[, c('cell_type', 'seurat_clusters'), drop = FALSE]
# 确保细胞顺序与 regulonAUC 的列名（细胞名）一致
cellinfo <- cellinfo[colnames(regulonAUC), , drop = FALSE]

# 转换为命名因子向量（calcRSS 要求）
annotation <- cellinfo$cell_type
names(annotation) <- rownames(cellinfo)   # 细胞名作为名字

# 检查有无缺失值
stopifnot(!anyNA(annotation))

# ---- 计算 RSS（直接传入 aucellResults 对象） ----
library(AUCell)   # 确保已加载
rss <- calcRSS(
  AUC = regulonAUC,       # 直接传入对象
  cellAnnotation = annotation
)

# 查看结果维度：应为 198（调控子） × 细胞类型数
dim(rss)
colnames(rss)   # 细胞类型名称

# ---- 绘制 RSS 热图 ----
pdf("SCENIC/RSS_heatmap.pdf", width = 10, height = 20)
rssPlot <- plotRSS(
  rss,
  zThreshold = 1.5,
  thr = 0.01,
  cluster_columns = FALSE,
  order_rows = TRUE,
  varName = "cell_type",
  col.low = '#330066',
  col.mid = '#66CC66',
  col.high = '#FFCC33'
)
print(rssPlot$plot)
dev.off()

# ---- 保存RSS结果 ----
write.csv(rss, "SCENIC/RSS_celltype.csv")

# ---- 查看Top特异性Regulon ----
rss_avg <- rowMeans(rss, na.rm = TRUE)
top_regulons <- names(sort(rss_avg, decreasing = TRUE))[1:10]
cat("Top 10 特异性Regulon:\n", paste(top_regulons, collapse = "\n"))

7.3 Regulon活性热图

展示各细胞类型中Regulon的平均活性：

library(ComplexHeatmap)
library(circlize)

# ---- 1. 提取 AUC 矩阵并转置 ----
auc_matrix <- getAUC(regulonAUC)          # 198 × 500（调控子 × 细胞）
auc_t <- t(auc_matrix)                    # 500 × 198（细胞 × 调控子）

# ---- 2. 准备细胞类型注释（确保顺序一致） ----
cell_types <- seurat_obj$cell_type[rownames(auc_t)]   # 顺序与 auc_t 行名一致
# 检查有无缺失
stopifnot(!anyNA(cell_types))

# ---- 3. 按细胞类型计算平均 AUC ----
# 方法：将 auc_t 转为数据框，用 aggregate
df_auc <- as.data.frame(auc_t)
df_auc$cell_type <- cell_types
auc_avg_df <- aggregate(. ~ cell_type, data = df_auc, FUN = mean)
rownames(auc_avg_df) <- auc_avg_df$cell_type
auc_avg_df <- auc_avg_df[, -1]                     # 去掉 cell_type 列
auc_avg <- t(auc_avg_df)                            # 转置为 调控子 × 细胞类型

# 查看维度
dim(auc_avg)   # 198 × 细胞类型数

# ---- 4. 筛选 Top 30 调控子（按方差） ----
top_regulons <- names(sort(apply(auc_avg, 1, var), decreasing = TRUE))[1:30]
auc_avg_top <- auc_avg[top_regulons, , drop = FALSE]

# ---- 5. Z-score 标准化（按行） ----
auc_avg_scaled <- t(scale(t(auc_avg_top)))

# ---- 6. 绘制热图 ----
pdf("SCENIC/regulon_activity_heatmap.pdf", width = 12, height = 10)
Heatmap(
  auc_avg_scaled,
  name = "AUC (scaled)",
  col = colorRamp2(c(-2, 0, 2), c("blue", "white", "red")),
  cluster_rows = TRUE,
  cluster_columns = TRUE,
  show_row_names = TRUE,
  show_column_names = TRUE,
  row_names_gp = gpar(fontsize = 8),
  column_names_gp = gpar(fontsize = 10)
)
dev.off()

7.4 在UMAP上展示特定Regulon活性

# ---- 1. 提取 AUC 矩阵 ----
auc_mat <- getAUC(regulonAUC)   # 普通矩阵：调控子 × 细胞

# ---- 2. 选择 Top 特异性 Regulon ----
# 假设 rss_avg 已经计算（例如 rowMeans(rss)）
top_regulon <- names(sort(rss_avg, decreasing = TRUE))[1]

# ---- 3. 提取该调控子的 AUC 向量 ----
# auc_mat[top_regulon, ] 返回命名向量（细胞名→AUC值）
auc_vector <- auc_mat[top_regulon, ]

# ---- 4. 只保留与 seurat_obj 共有的细胞 ----
common_cells <- intersect(colnames(seurat_obj), names(auc_vector))

# ---- 5. 将 AUC 添加到 Seurat 对象 ----
# 先创建空列（所有细胞为 NA）
seurat_obj[[top_regulon]] <- NA_real_
# 给共有细胞赋值（注意顺序）
seurat_obj[[top_regulon]][common_cells, 1] <- auc_vector[common_cells]

# ---- 6. 绘制 FeaturePlot ----
p <- FeaturePlot(
  seurat_obj,
  features = top_regulon,
  order = TRUE,
  cols = c("grey90", "red")
) + ggtitle(paste("Regulon:", top_regulon))

# 保存
ggsave(paste0("SCENIC/umap_", top_regulon, ".pdf"), p, width = 8, height = 6)

7.5 Regulon活性箱线图

展示Top Regulon在不同细胞类型中的活性分布：

library(tidyverse)
library(ggplot2)

# ---- 1. 提取 AUC 矩阵 ----
auc_mat <- getAUC(regulonAUC)   # 调控子 × 细胞，普通矩阵

# ---- 2. 选择 Top 6 调控子 ----
# 假设 rss_avg 已计算（例如 rowMeans(rss)）
top6_regulons <- names(sort(rss_avg, decreasing = TRUE))[1:6]

# ---- 3. 提取并转置数据 ----
# 先子集再转置，得到 细胞 × 调控子
auc_sub <- t(auc_mat[top6_regulons, ])
plot_df <- as.data.frame(auc_sub)
# 添加细胞类型，按行名匹配
plot_df$cell_type <- seurat_obj$cell_type[rownames(plot_df)]

# ---- 4. 转为长格式 ----
plot_long <- plot_df %>%
  pivot_longer(
    cols = all_of(top6_regulons),
    names_to = "regulon",
    values_to = "auc"
  )

# ---- 5. 绘制箱线图 ----
p_box <- ggplot(plot_long, aes(x = regulon, y = auc, fill = cell_type)) +
  geom_boxplot(outlier.size = 0.5) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(x = "", y = "AUC Score", title = "Top 6 Regulon Activity by Cell Type")

# ---- 6. 保存 ----
ggsave("SCENIC/regulon_boxplot.pdf", p_box, width = 10, height = 6)

八、结果解读指南

8.1 RSS热图怎么看？

行：不同的Regulon（转录因子调控程序）
列：不同的细胞类型
颜色：Z-score（标准化后的RSS值）
红色方块：该Regulon在该细胞类型中高度特异激活

生物学含义：如果一个Regulon在神经元中为红色，说明该转录因子调控程序在神经元中被特异激活。

8.2 Regulon活性热图怎么看？

横轴为细胞类型，纵轴为Regulon
颜色表示AUC值的高低
可以发现哪些转录因子在特定细胞类型中整体活跃

8.3 如何从热图筛选候选TF？

在RSS热图中找到红色方块 → 该Regulon具有细胞类型特异性
检查该Regulon在活性热图中的表现 → 确认其在该细胞类型中确实高表达
查看该Regulon的靶基因 → 是否与该细胞类型的已知功能相关

九、常见问题与解决方案

9.1 数据库文件不匹配

症状：步骤报错，提示找不到motif
解决：

确认feather、tbl、tfs三个文件来自同一物种
检查文件路径是否正确，文件名是否一致

9.2 GRNBoost2运行时间过长

症状：GRN步骤运行超过24小时
解决：

减少输入细胞数（随机抽样200-500个）
减少--num_workers避免内存溢出
适当过滤低表达基因

9.3 AUC矩阵全为零

症状：所有AUC值都为0
解决：

检查表达矩阵是否为counts格式（pySCENIC推荐使用counts而非标准化数据）
确认基因名格式一致（gene-symbol）
检查数据库物种是否与数据匹配

十、总结

至此，我们完成了一套从上游数据准备到下游可视化的pySCENIC完整闭环流程：

阶段	工具/语言	核心任务
数据准备	R + loomR	Seurat → loom文件导出
GRN推断	pySCENIC (GRNBoost2)	TF-靶基因共表达关系
Motif富集	pySCENIC (cisTarget)	Regulon构建与剪枝
活性评分	pySCENIC (AUCell)	每个细胞的Regulon活性
下游分析	R	RSS计算、热图、UMAP叠加

核心优势：

速度：GRNBoost2比GENIE3快数十倍
可扩展：支持多核并行处理
可重复：Shell脚本一键运行
完整闭环：从上游到下游无缝衔接

这里是两栖生物手册，中科院生物医学博士，持续记录生物医学实验和生物信息学笔记，干湿结合两不误～正努力成为最贴心负责的生信数据分析者和实验技术分享者，欢迎大家关注～