# Python版pySCENIC实战：从单细胞数据到转录因子网络可视化（附避坑指南） 如果你正在单细胞转录组领域深耕，想从海量的基因表达数据中挖掘出驱动细胞命运的“幕后导演”——转录因子，那么pySCENIC绝对是你绕不开的利器。这个基于Python的版本，相比其R语言前身，在分析速度上实现了质的飞跃，尤其适合处理成千上万个细胞的庞大数据集。然而，从环境配置、数据准备到最终的可视化，这条路上布满了依赖冲突、版本不兼容、路径错误等“暗坑”。这篇文章，我将结合自己多次在实验室服务器上部署和运行pySCENIC的实战经验，为你梳理一套清晰、可复现的完整流程，并重点分享那些官方文档里不会写的“避坑”细节。无论你是第一次接触pySCENIC，还是曾被其复杂的安装过程劝退，这篇指南都将助你顺利通关。 ## 1. 环境搭建：用Conda构筑稳固的分析基石 在Linux服务器上，一个独立、纯净的Python环境是成功运行pySCENIC的第一步。Conda是管理此类复杂生物信息学软件依赖的黄金标准。但请注意，pySCENIC对Python版本有特定要求，直接使用`pip install`很容易导致依赖地狱。 ### 1.1 创建并激活专用环境 首先，我们创建一个名为`pyscenic`的Conda环境，并指定Python版本为3.9（经测试，3.7至3.9版本兼容性较好，3.10及以上版本可能存在部分库的兼容性问题）。 ```bash # 创建环境 conda create -n pyscenic python=3.9 -y # 激活环境 conda activate pyscenic ``` > **注意**：有些教程会建议使用`mamba`来加速包解决过程，如果你的服务器已经安装了mamba，可以用`mamba create`替代`conda create`，速度会快很多。 ### 1.2 安装核心依赖与pySCENIC 接下来，我们需要安装一系列科学计算和单细胞分析的基础包，最后安装pySCENIC本身。这里的关键是**按顺序安装**，并优先使用Conda渠道的预编译包。 ```bash # 1. 安装基础科学计算栈 conda install -c conda-forge numpy pandas scipy scanpy -y # 2. 安装一些必要的工具库 conda install -c anaconda cytoolz dask -y # 3. 安装pySCENIC（使用pip从PyPI安装） pip install pyscenic ``` **常见坑点与解决方案：** * **`pandas`版本冲突**：这是最常见的错误。如果你遇到类似“`ImportError: cannot import name 'DtypeArg' from 'pandas._typing'`”的报错，说明pandas版本过高。解决方法是指定安装兼容版本： ```bash pip install "pandas<2.0" "pyscenic" ``` * **网络超时**：从PyPI下载包时可能因网络问题失败。可以尝试使用国内镜像源，或分步安装： ```bash pip install pyscenic -i https://pypi.tuna.tsinghua.edu.cn/simple ``` * **内存不足**：在资源有限的服务器上，Conda解决依赖关系时可能内存溢出。可以尝试增加交换空间，或使用更轻量的`mamba`。 安装完成后，强烈建议进行一个快速验证： ```bash python -c "import pyscenic; print(f'pySCENIC version: {pyscenic.__version__}')" ``` 如果成功输出版本号（如`0.12.1`），则说明核心环境已就绪。 ## 2. 关键数据准备：数据库与输入矩阵 pySCENIC的分析依赖于两个外部数据库：**转录因子-靶基因调控关系数据库（Feather格式）**和**转录因子注释文件（TBL格式）**。同时，你需要将单细胞表达矩阵转换为pySCENIC接受的Loom格式。 ### 2.1 下载并配置物种数据库 数据库需要根据你的研究物种（人或小鼠）从官方源下载。**务必确认基因组版本（如hg19/hg38, mm9/mm10）与你的单细胞数据比对所用版本一致**。 以下以人类hg19为例，在服务器上创建目录并下载： ```bash # 创建数据库存储目录 mkdir -p /path/to/your/db/cisTarget_databases cd /path/to/your/db/cisTarget_databases # 下载人类（hg19）的数据库文件 # 1. 基因启动子区域（TSS上游500bp）的调控信息 wget -c https://resources.aertslab.org/cistarget/databases/homo_sapiens/hg19/refseq_r45/mc9nr/gene_based/hg19-500bp-upstream-10species.mc9nr.genes_vs_motifs.rankings.feather # 2. 基因TSS上下游10kb区域的调控信息（更全面） wget -c https://resources.aertslab.org/cistarget/databases/homo_sapiens/hg19/refseq_r45/mc9nr/gene_based/hg19-tss-centered-10kb-10species.mc9nr.genes_vs_motifs.rankings.feather # 3. 下载Motif到转录因子的注释文件 wget -c https://resources.aertslab.org/cistarget/motif2tf/motifs-v9-nr.hgnc-m0.001-o0.0.tbl # 4. 下载人类转录因子列表 wget -c https://github.com/aertslab/pySCENIC/raw/master/resources/hs_hgnc_tfs.txt ``` 下载完成后，检查文件大小和完整性。两个Feather文件通常都很大（几百MB到几GB），TBL文件约100MB，TFs列表文件很小。 ### 2.2 准备单细胞输入数据：从Seurat到Loom pySCENIC的输入需要是Loom格式的文件，其中行为基因，列为细胞。大多数单细胞分析始于R语言的Seurat对象。转换过程分为两步：从Seurat导出矩阵，再用Python转为Loom。 **第一步：在R中导出表达矩阵** 假设你的Seurat对象名为`seurat_obj`，使用原始计数矩阵（`RNA@counts`）进行分析效果更佳。**关键点是矩阵需要转置**，因为pySCENIC期望列是细胞。 ```r # 激活包含Seurat的R环境 library(Seurat) # 导出转置后的计数矩阵为CSV # 注意：细胞数过多时，CSV文件会非常大，可以考虑导出为`.h5ad`或`.mtx`格式，但Loom转换脚本需相应调整。 write.csv(t(as.matrix(seurat_obj@assays$RNA@counts)), file = "sc_counts_for_pyscenic.csv", row.names = TRUE) # 行名是基因名 ``` **第二步：在Python中将CSV转换为Loom** 将上一步生成的`sc_counts_for_pyscenic.csv`上传到服务器你的工作目录。然后创建一个Python脚本`csv_to_loom.py`： ```python import loompy as lp import numpy as np import scanpy as sc import pandas as pd # 读取CSV文件，注意设置第一列为行名（基因名） print("Reading CSV file...") adata = sc.read_csv("sc_counts_for_pyscenic.csv", first_column_names=True) # scanpy读取 # 检查维度：行应为基因，列应为细胞 print(f"Data shape: {adata.shape} (genes, cells)") # 创建Loom文件所需的属性 row_attrs = {"Gene": np.array(adata.var_names)} col_attrs = {"CellID": np.array(adata.obs_names)} # 创建Loom文件。注意：adata.X通常是稀疏矩阵，需要转置并转换为稠密数组（或保持稀疏存储） # 对于大型数据集，直接转换为稠密数组可能内存爆炸。loompy支持稀疏矩阵输入。 print("Creating loom file...") lp.create("sample_input.loom", adata.X.T, # 关键转置：loom文件要求行为基因，列为细胞，而scanpy读取后行为细胞，列为基因。 row_attrs=row_attrs, col_attrs=col_attrs) print("Loom file 'sample_input.loom' created successfully.") ``` 在激活的`pyscenic`环境中运行此脚本： ```bash python csv_to_loom.py ``` 如果一切顺利，你将得到`sample_input.loom`文件，这就是pySCENIC三大步骤的输入起点。 > **避坑提示**：如果细胞数量巨大（>10万），CSV文件的读写和转换会非常缓慢且耗内存。此时，可以考虑在R中使用`SeuratDisk`或`sceasy`包直接将Seurat对象转换为H5AD格式，然后使用`scanpy`读取H5AD再转换为Loom，效率更高。 ## 3. 核心三步分析流程：GRNBoost2、cisTarget与AUCell pySCENIC的分析逻辑清晰分为三步，每一步都对应一个独立的命令行工具。我们将在一个Shell脚本中串联执行，并合理设置计算资源。 ### 3.1 第一步：GRNBoost2推断共表达网络 这一步利用GRNBoost2算法，基于基因间的共表达关系，推断每个转录因子（TF）的潜在靶基因。 ```bash #!/bin/bash # scenic_pipeline.sh # 设置路径（根据你的实际路径修改） TF_LIST="/path/to/your/db/cisTarget_databases/hs_hgnc_tfs.txt" INPUT_LOOM="sample_input.loom" OUTPUT_PREFIX="pyscenic_results" NUM_WORKERS=20 # 根据服务器CPU核心数调整，通常设置为可用核心数-2 echo "Step 1: Running GRNBoost2 to infer co-expression network..." pyscenic grn \ --num_workers $NUM_WORKERS \ --output ${OUTPUT_PREFIX}_adj.tsv \ --method grnboost2 \ $INPUT_LOOM \ $TF_LIST # 检查上一步是否成功 if [ $? -eq 0 ]; then echo "GRNBoost2 completed successfully." else echo "GRNBoost2 failed. Exiting." exit 1 fi ``` **参数解读与避坑：** * `--num_workers`：并行进程数。设置过高可能导致内存不足，建议从10开始尝试，观察内存使用情况。 * `--method grnboost2`：默认算法，比原始的GENIE3更快更节省内存。 * 输出文件`*_adj.tsv`是一个三列（TF， Target， Importance）的表格，记录了所有推断出的调控关系及其权重。 ### 3.2 第二步：cisTarget进行Motif富集分析 这一步对上一步得到的共表达网络进行“精炼”，通过扫描靶基因启动子区域的DNA motif，筛选出那些可能被转录因子直接结合的、生物学意义更强的调控关系（Regulon）。 ```bash # 接上一步的脚本 FEATHER_DB="/path/to/your/db/cisTarget_databases/hg19-tss-centered-10kb-10species.mc9nr.genes_vs_motifs.rankings.feather" MOTIF_ANNOT="/path/to/your/db/cisTarget_databases/motifs-v9-nr.hgnc-m0.001-o0.0.tbl" echo "Step 2: Running cisTarget for motif enrichment..." pyscenic ctx \ ${OUTPUT_PREFIX}_adj.tsv \ $FEATHER_DB \ --annotations_fname $MOTIF_ANNOT \ --expression_mtx_fname $INPUT_LOOM \ --mode "dask_multiprocessing" \ --output ${OUTPUT_PREFIX}_regulons.csv \ --num_workers $NUM_WORKERS \ --mask_dropouts if [ $? -eq 0 ]; then echo "cisTarget completed successfully." else echo "cisTarget failed. Exiting." exit 1 fi ``` **关键点与避坑：** * `--mask_dropouts`：这个参数非常重要。它会忽略表达矩阵中的零值（dropouts），避免因技术噪音导致的假阳性motif富集。 * 输出文件`*_regulons.csv`包含了鉴定出的Regulon。每个Regulon以`TF(+)`或`TF(±)`命名，括号内的符号表示该Regulon被预测为激活(+)或双功能(±)状态。文件内容是该TF调控的靶基因列表。 ### 3.3 第三步：AUCell计算Regulon活性 最后一步，计算每个Regulon在每个细胞中的活性分数（AUC值）。这个分数代表了该调控网络在特定细胞中的“活跃程度”。 ```bash # 接上一步脚本 echo "Step 3: Running AUCell to calculate regulon activity..." pyscenic aucell \ $INPUT_LOOM \ ${OUTPUT_PREFIX}_regulons.csv \ --output ${OUTPUT_PREFIX}_aucell.loom \ --num_workers $NUM_WORKERS if [ $? -eq 0 ]; then echo "AUCell completed successfully." echo "All pySCENIC steps finished. Final output: ${OUTPUT_PREFIX}_aucell.loom" else echo "AUCell failed." exit 1 fi ``` 至此，核心分析完成。你得到的`*_aucell.loom`文件包含了原始的基因表达矩阵、细胞和基因的元数据，以及最重要的——一个名为`RegulonsAUC`的层（layer），存储了每个细胞中每个Regulon的AUC值矩阵。这个文件是后续所有可视化和下游分析的基石。 ## 4. 结果可视化与生物学解读：在R中释放洞察力 虽然分析在Python中完成，但强大的单细胞可视化生态仍在R中（特别是Seurat和ggplot2）。我们需要将pySCENIC的结果与原始的Seurat对象整合。 ### 4.1 在R中加载并整合AUC矩阵 首先，确保安装了必要的R包：`SCopeLoomR`, `AUCell`, `Seurat`, `pheatmap`, `RColorBrewer`, `dplyr`等。 ```r # 加载R包 library(SCopeLoomR) library(AUCell) library(Seurat) library(pheatmap) library(dplyr) # 1. 读取pySCENIC输出的loom文件 scenic_loom <- open_loom("pyscenic_results_aucell.loom") # 2. 提取Regulon的AUC矩阵 regulon_auc <- get_regulons_AUC(scenic_loom, column.attr.name = "RegulonsAUC") # regulon_auc 是一个 SummarizedExperiment 对象，assay部分存储了AUC矩阵 auc_matrix <- assay(regulon_auc) # 矩阵维度：Regulons x Cells # 3. 提取Regulon的具体基因列表（可选，用于后续分析） regulons <- get_regulons(scenic_loom, column.attr.name = "Regulons") regulon_gene_list <- regulonsToGeneLists(regulons) # 转换为TF->靶基因列表 close_loom(scenic_loom) # 4. 加载你原始的Seurat对象（包含细胞聚类、UMAP坐标、细胞类型注释等） seurat_obj <- readRDS("your_original_seurat_object.rds") # 5. 确保AUC矩阵的细胞与Seurat对象完全匹配 # 通常，在创建loom文件时，细胞顺序是一致的。但务必检查： all(colnames(auc_matrix) == colnames(seurat_obj)) # 应该返回TRUE # 6. 将AUC矩阵作为一个新的Assay添加到Seurat对象中 seurat_obj[["SCENIC"]] <- CreateAssayObject(counts = auc_matrix) # 默认情况下，`counts`槽位存储原始值。我们可以将其复制到`data`槽位用于标准化和可视化。 seurat_obj <- SetAssayData(seurat_obj, slot = "data", assay = "SCENIC", new.data = auc_matrix) # 设置默认assay为SCENIC，方便后续函数调用 DefaultAssay(seurat_obj) <- "SCENIC" ``` ### 4.2 可视化策略：从整体到特异 整合完成后，就可以利用Seurat丰富的绘图函数进行可视化了。 **（1）Regulon活性热图（按细胞聚类）** 这是最直观的展示方式，可以一眼看出哪些Regulon在哪些细胞类群中特异性活跃。 ```r # 计算每个细胞簇（cluster）的Regulon平均活性 cluster.averages <- AverageExpression(seurat_obj, assays = "SCENIC", group.by = "seurat_clusters") # 按聚类分组 avg_auc <- cluster.averages$SCENIC # 对行（Regulon）进行Z-score标准化，便于比较 avg_auc_scaled <- t(scale(t(avg_auc))) # 绘制热图 pheatmap(avg_auc_scaled, color = colorRampPalette(rev(brewer.pal(11, "RdBu")))(100), cluster_rows = TRUE, cluster_cols = TRUE, show_rownames = FALSE, # Regulon太多，通常不显示全部名称 fontsize_col = 8, main = "Regulon Activity (AUC) by Cluster") ``` 为了更清晰，我们通常只展示在至少一个簇中高活性的Regulon。 ```r # 找出在任一簇中Z-score > 1.5的Regulon high_activity_regulons <- rownames(avg_auc_scaled)[apply(avg_auc_scaled, 1, max) > 1.5] # 绘制筛选后的热图，并显示Regulon名称 pheatmap(avg_auc_scaled[high_activity_regulons, ], color = colorRampPalette(rev(brewer.pal(11, "RdBu")))(100), cluster_rows = TRUE, cluster_cols = TRUE, show_rownames = TRUE, fontsize_row = 7, angle_col = 45) ``` **（2）在UMAP图上展示特定Regulon的活性** 类似于绘制基因表达量，我们可以将Regulon的AUC值映射到细胞的UMAP坐标上，观察其空间分布。 ```r # 切换回SCENIC assay DefaultAssay(seurat_obj) <- "SCENIC" # 选择几个感兴趣的Regulon进行可视化，例如与特定功能相关的TF regulons_to_plot <- c("STAT1(+)", "IRF7(+)", "NFKB1(+)") # 示例：炎症相关TF FeaturePlot(seurat_obj, features = regulons_to_plot, cols = c("lightgrey", "red"), ncol = 2, order = TRUE) # order=TRUE使高活性细胞点绘制在上层 ``` **（3）Regulon特异性评分（RSS）与气泡图** RSS可以帮助量化一个Regulon对某个细胞类型的特异性。我们可以用`AUCell`包自带的函数计算，并用气泡图展示。 ```r # 计算RSS library(AUCell) # 需要细胞类型注释信息，假设存储在seurat_obj$cell_type中 cell_types <- seurat_obj$cell_type names(cell_types) <- colnames(seurat_obj) rss <- calcRSS(AUC = getAUC(regulon_auc), # 使用之前提取的regulon_auc对象 cellAnnotation = cell_types[colnames(regulon_auc)]) # 提取RSS矩阵 rss_mat <- rss$RSS # 筛选每个细胞类型中RSS最高的前N个Regulon top_n <- 5 top_regulons <- lapply(colnames(rss_mat), function(ct){ head(rownames(rss_mat)[order(rss_mat[, ct], decreasing = TRUE)], top_n) }) top_regulons <- unique(unlist(top_regulons)) # 准备绘图数据 plot_data <- rss_mat[top_regulons, ] plot_data <- reshape2::melt(plot_data) colnames(plot_data) <- c("Regulon", "CellType", "RSS") # 使用ggplot2绘制气泡图 library(ggplot2) ggplot(plot_data, aes(x = CellType, y = Regulon)) + geom_point(aes(size = RSS, color = RSS)) + scale_color_gradient(low = "blue", high = "red") + theme_minimal() + theme(axis.text.x = element_text(angle = 45, hjust = 1)) + labs(title = "Top Regulon Specificity Score (RSS) per Cell Type") ``` 这套从环境搭建、数据分析到结果可视化的完整流程，结合了Python的高效计算和R的卓越可视化能力。在实际操作中，最耗时的往往是第一步GRNBoost2，尤其是当细胞数和基因数很大时。因此，在服务器上使用足够的CPU核心并合理分配内存至关重要。记住，成功的单细胞调控网络分析，一半在于精准的计算，另一半在于对结果的生物学解读。当你看到那些在特定细胞亚群中高活性的转录因子时，不妨结合已有的生物学知识，去推测它们可能主导的细胞功能或状态，这才是分析的价值所在。