从零构建单细胞调控网络：pySCENIC保姆级教程+Seurat联动可视化

从零构建单细胞调控网络：pySCENIC保姆级教程与Seurat联动可视化实战

如果你已经完成了单细胞转录组的基础分析，比如用Seurat做了细胞分群和差异表达，那么下一步，你可能会好奇：驱动这些细胞类型分化和功能差异的“幕后导演”是谁？答案很可能就藏在转录因子（Transcription Factors, TFs）及其构建的基因调控网络（Gene Regulatory Network, GRN）中。传统的差异表达分析能告诉我们哪些基因在不同细胞间有表达差异，但它无法揭示这些差异背后的调控逻辑——是哪些关键的转录因子“按下开关”，激活或抑制了特定的基因程序，从而塑造了独特的细胞状态？

这正是SCENIC（Single-Cell rEgulatory Network Inference and Clustering）这类工具大显身手的地方。它不再局限于静态的基因表达量比较，而是通过整合共表达模式与DNA顺式作用元件（motif） 信息，从单细胞数据中直接推断出活跃的转录因子-靶基因调控关系（即Regulon），并量化每个Regulon在每个细胞中的活性。最终，我们能得到一张动态的“调控电路图”，直观地看到哪些转录因子在哪些细胞亚群中扮演着核心调控角色。

对于熟悉Python生态的分析者而言，pySCENIC因其更快的计算速度（尤其是GRNBoost2算法）和与Python数据分析栈（如pandas, scanpy）的无缝衔接而备受青睐。然而，许多单细胞分析流程始于R语言的Seurat，这就引出了一个常见的痛点：如何将Seurat对象中精心注释的细胞信息，与pySCENIC强大的网络推断能力结合起来，并进行高效、直观的可视化？

本文将手把手带你走通这条融合之路。我们将从一个真实的膀胱癌单细胞数据集出发，详细演示如何将Seurat对象转换为pySCENIC所需的输入格式，完整执行其三步分析流程（GRNBoost2 -> cisTarget -> AUCell），并最终将计算得到的转录因子活性（AUC值）完美映射回原始的Seurat UMAP图中进行解读。整个过程，你将不仅学会工具的使用，更能理解每一步背后的生物学意义和计算逻辑。

1. 环境准备与数据基础

在开始任何计算之前，搭建一个稳定、可复现的分析环境至关重要。pySCENIC的依赖相对复杂，强烈建议使用Conda进行环境管理，这能有效避免包版本冲突。

1.1 创建并激活Conda环境

首先，我们创建一个名为pyscenic_env的独立Python环境，并安装指定版本的pySCENIC及其核心依赖。

# 创建新环境，指定Python版本为3.9（经测试兼容性较好）
conda create -n pyscenic_env python=3.9 -y

# 激活环境
conda activate pyscenic_env

# 通过conda-forge频道安装pySCENIC及其主要依赖
conda install -c conda-forge pyscenic scanpy loompy -y

# 安装一些额外的实用数据处理包
conda install -c conda-forge pandas numpy scikit-learn matplotlib seaborn -y

# 安装用于加速计算的dask（可选，但对于大数据集强烈推荐）
conda install -c conda-forge dask distributed -y

注意：pySCENIC对pandas和numpy的版本较为敏感。如果后续运行报错，可以尝试固定版本，例如 conda install pandas==1.5.3 numpy==1.23.5。

1.2 准备参考数据库与转录因子列表

pySCENIC的核心步骤——cisTarget分析，需要物种特异的参考数据库（包含基因启动子区域的motif富集信息）和转录因子列表。对于人类（hg19）和小鼠（mm9/mm10），官方提供了预编译的数据库。

# 创建一个专门存放参考数据的目录
mkdir -p scenic_references
cd scenic_references

# 下载人类（hg19）的参考数据库（约1GB每个）
# TSS上游500bp区域数据库
wget https://resources.aertslab.org/cistarget/databases/homo_sapiens/hg19/refseq_r45/mc9nr/gene_based/hg19-500bp-upstream-7species.mc9nr.feather
# TSS中心±10kb区域数据库（更常用，覆盖范围更广）
wget https://resources.aertslab.org/cistarget/databases/homo_sapiens/hg19/refseq_r45/mc9nr/gene_based/hg19-tss-centered-10kb-7species.mc9nr.feather

# 下载人类转录因子列表
wget https://raw.githubusercontent.com/aertslab/pySCENIC/master/resources/hs_hgnc_tfs.txt

# 下载Motif到TF的注释文件
wget https://resources.aertslab.org/cistarget/motif2tf/motifs-v9-nr.hgnc-m0.001-o0.0.tbl

# 返回工作目录
cd ..

关键文件说明：

• *.feather文件：是cisTarget分析的基石，存储了全基因组范围内每个基因附近区域（如TSS±10kb）的motif富集得分。mc9nr代表第9版Motif合集，包含约2.4万个motif。
• hs_hgnc_tfs.txt：一个简单的文本文件，每行列出一个已知的人类转录因子基因符号（HGNC格式）。pySCENIC将只对这些基因进行调控网络推断。
• motifs-v9-nr.hgnc-m0.001-o0.0.tbl：建立了motif ID与转录因子之间的对应关系，用于将cisTarget富集的motif注释到具体的TF上。

1.3 准备单细胞表达矩阵：从Seurat到Loom

假设你已经有一个处理好的Seurat对象（例如名为seurat_obj），其中包含了归一化后的表达矩阵和细胞注释信息。pySCENIC的输入需要是细胞为行、基因为列的计数矩阵（推荐使用归一化后的计数，如SCTransform或LogNormalize的结果），并保存为loom格式。

在R环境中操作：

# 加载必要的R包
library(Seurat)
library(loo