单细胞调控网络SCENIC¶

SCENIC的三步核心流程：

步骤1：基因调控网络推断（GRN Inference）
  输入：单细胞表达矩阵
  方法：GRNBoost2（梯度提升，快）或GENIE3（随机森林，准）
  输出：TF → 靶基因的共表达模块（粗略，含间接靶基因）

  思路：如果基因X和TF_A在很多细胞中表达模式相似
  → 可能X是TF_A的靶基因

步骤2：Regulon优化（Motif Enrichment）
  输入：步骤1的共表达模块
  方法：cisTarget顺式调控元件分析
  输出：去除假阳性的精炼regulon

  思路：如果X是TF_A的真正靶基因
  → X的启动子区域应该有TF_A的结合基序（motif）
  → 有motif的保留，没有的去除

步骤3：细胞活性打分（AUCell Scoring）
  输入：步骤2的regulon + 表达矩阵
  方法：AUCell（Area Under recovery Curve）
  输出：每个细胞×每个regulon的活性矩阵

  思路：一个细胞中regulon的靶基因越多被高表达
  → 该regulon在该细胞中越活跃
  → AUC分数越高

SCENIC依赖的外部数据库：

1. cisTarget数据库（motif排名数据库）
   - 包含每个基因周围的TF motif富集信息
   - 人类：hg38 (~25GB)
   - 小鼠：mm10 (~25GB)
   - 需要预先下载

2. TF列表
   - 定义哪些基因是转录因子
   - 人类：约1,600个TF
   - 小鼠：约1,400个TF

下载地址：
  https://resources.aertslab.org/cistarget/
  包括：
  - ranking数据库（.feather格式，大文件）
  - motif注释（.tbl格式）
  - TF列表（.txt格式）

# 推荐用conda安装pySCENIC
conda create -n scenic python=3.10  # 创建环境
conda activate scenic                # 激活环境

pip install pyscenic                  # 安装pySCENIC
pip install scanpy loompy             # 安装依赖

# 下载cisTarget数据库（重要！文件很大）
# 人类hg38数据库
wget https://resources.aertslab.org/cistarget/databases/homo_sapiens/hg38/refseq_r80/mc_v10_clust/gene_based/hg38_10kbp_up_10kbp_down_full_tx_v10_clust.genes_vs_motifs.rankings.feather

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

# TF列表
wget https://resources.aertslab.org/cistarget/tf_lists/allTFs_hg38.txt

import scanpy as sc
import loompy as lp
import numpy as np
import pandas as pd

# 读取单细胞数据
adata = sc.read_h5ad("processed_scRNA.h5ad")  # 读取h5ad

# 基本过滤（SCENIC需要原始计数矩阵）
sc.pp.filter_genes(adata, min_cells=3)   # 基因至少在3个细胞表达
sc.pp.filter_cells(adata, min_genes=200)  # 细胞至少200个基因

# 导出为loom格式（pySCENIC的输入格式）
# loom格式包含表达矩阵和元数据
row_attrs = {
    "Gene": np.array(adata.var_names)      # 基因名
}
col_attrs = {
    "CellID": np.array(adata.obs_names),   # 细胞ID
    "CellType": np.array(adata.obs["cell_type"])  # 细胞类型
}

lp.create(
    "input.loom",                          # 输出文件名
    adata.X.T.toarray(),                   # 表达矩阵（转置）
    row_attrs,                              # 行属性（基因）
    col_attrs                               # 列属性（细胞）
)
print(f"导出{adata.n_obs}个细胞×{adata.n_vars}个基因")

# pySCENIC有3个命令行步骤

# ===== 步骤1：GRN推断（GRNBoost2） =====
pyscenic grn \
  input.loom \                  # 输入loom文件
  allTFs_hg38.txt \             # TF列表
  -o adjacencies.csv \          # 输出邻接矩阵
  --num_workers 8 \             # 并行线程数
  --sparse                      # 稀疏矩阵模式（节省内存）

# 这一步最耗时（小时级别），输出TF-靶基因对及其权重
# adjacencies.csv: TF, target, importance

# ===== 步骤2：Regulon优化（cisTarget） =====
pyscenic ctx \
  adjacencies.csv \             # 步骤1的输出
  hg38_10kbp_up_10kbp_down_full_tx_v10_clust.genes_vs_motifs.rankings.feather \  # 排名数据库
  --annotations_fname motifs-v10nr_clust-nr.hgnc-m0.001-o0.0.tbl \  # motif注释
  --expression_mtx_fname input.loom \  # 表达矩阵
  --output regulons.csv \       # 输出regulon
  --num_workers 8 \             # 线程数
  --mask_dropouts               # 忽略dropout的零值

# 输出regulon.csv：每个regulon包含TF名和靶基因列表
# 只保留启动子区有TF motif的靶基因

# ===== 步骤3：AUCell打分 =====
pyscenic aucell \
  input.loom \                  # 输入loom文件
  regulons.csv \                # 步骤2的regulon
  -o output.loom \              # 输出loom（含AUCell分数）
  --num_workers 8               # 线程数

# 输出output.loom包含：
# 原始表达矩阵 + 每个细胞×每个regulon的AUCell分数

# 也可以在Python脚本中直接调用
from pyscenic.rnkdb import FeatherRankingDatabase as RankingDatabase
from pyscenic.utils import modules_from_adjacencies
from pyscenic.prune import prune2df, df2regulons
from pyscenic.aucell import aucell
from arboreto.algo import grnboost2

import pandas as pd

# 步骤1：GRNBoost2
adjacencies = grnboost2(
    expression_data=adata.to_df(),     # 表达矩阵
    tf_names=tf_names,                  # TF列表
    verbose=True
)

# 步骤2：cisTarget pruning
dbs = [RankingDatabase(fname=ranking_db_path, name="hg38")]
modules = modules_from_adjacencies(adjacencies, adata.to_df())

df = prune2df(dbs, modules, motif_annotations_fname)
regulons = df2regulons(df)

# 步骤3：AUCell
auc_mtx = aucell(adata.to_df(), regulons, num_workers=8)

# auc_mtx: DataFrame, 行=细胞, 列=regulon
print(f"检测到{len(regulons)}个regulon")
print(f"AUCell矩阵: {auc_mtx.shape}")

import scanpy as sc
import matplotlib.pyplot as plt

# 将AUCell结果添加到adata
adata.obsm["X_aucell"] = auc_mtx.loc[adata.obs_names, :].values

# 在UMAP上可视化特定regulon的活性
for tf in ["PAX6(+)", "FOXP1(+)", "MYC(+)"]:
    if tf in auc_mtx.columns:
        adata.obs[tf] = auc_mtx.loc[adata.obs_names, tf].values
        sc.pl.umap(adata, color=tf, cmap="viridis",
                   title=f"{tf} Activity")

# 二值化——判断regulon在每个细胞中是"开"还是"关"
from pyscenic.binarization import binarize

binary_mtx, thresholds = binarize(auc_mtx)
# binary_mtx: 0/1矩阵
# thresholds: 每个regulon的二值化阈值

# 统计每个细胞类型中活跃的regulon
cell_types = adata.obs["cell_type"]
regulon_activity_by_type = binary_mtx.groupby(cell_types).mean()
# 值代表该细胞类型中regulon活跃的细胞比例

import seaborn as sns

# 1. Regulon活性热图（按细胞类型汇总）
mean_auc = auc_mtx.groupby(adata.obs["cell_type"]).mean()

# 选择变异最大的Top30 regulon
top_regulons = mean_auc.var().nlargest(30).index
plt.figure(figsize=(15, 8))
sns.clustermap(mean_auc[top_regulons].T,
               cmap="YlOrRd",                # 黄-橙-红配色
               standard_scale=0,               # 按行标准化
               figsize=(12, 10),
               yticklabels=True)              # 显示regulon名
plt.title("Top 30 Regulon Activity by Cell Type")
plt.savefig("regulon_heatmap.pdf", dpi=300)
plt.show()

# 2. 特定细胞类型的关键regulon
# 找每个细胞类型特异性活跃的regulon
from scipy.stats import mannwhitneyu

specific_regulons = {}
for ct in cell_types.unique():
    ct_cells = cell_types == ct
    pvalues = []
    for reg in auc_mtx.columns:
        _, p = mannwhitneyu(
            auc_mtx.loc[ct_cells, reg],      # 该类型细胞
            auc_mtx.loc[~ct_cells, reg]        # 其他细胞
        )
        pvalues.append(p)
    sig_regs = auc_mtx.columns[np.array(pvalues) < 0.001]
    specific_regulons[ct] = sig_regs.tolist()
    print(f"{ct}: {len(sig_regs)}个特异regulon")

# RSS衡量每个regulon对某个细胞类型的特异性
# 基于Jensen-Shannon散度

from pyscenic.rss import regulon_specificity_scores

rss = regulon_specificity_scores(auc_mtx, adata.obs["cell_type"])

# 找每个细胞类型最特异的Top5 regulon
for ct in rss.columns:
    top5 = rss[ct].nlargest(5)
    print(f"\n{ct} Top5 regulon:")
    for reg, score in top5.items():
        print(f"  {reg}: RSS={score:.3f}")

# 将pySCENIC结果导入R/Seurat
library(Seurat)
library(SCopeLoomR)

# 读取pySCENIC输出的loom文件
loom <- open_loom("output.loom")
regulon_auc <- get_regulons_AUC(loom)  # 获取AUCell矩阵
close_loom(loom)

# 添加到Seurat对象
# 转换为矩阵
auc_matrix <- as.matrix(regulon_auc)

# 添加为metadata
for (reg in colnames(auc_matrix)) {
  seurat_obj[[reg]] <- auc_matrix[colnames(seurat_obj), reg]
}

# 在Seurat中可视化
FeaturePlot(seurat_obj,
            features = c("PAX6(+)", "FOXP1(+)"),
            cols = c("grey", "red"),
            ncol = 2)

# 按cluster汇总regulon活性
avg_auc <- AverageExpression(seurat_obj,
                              features = colnames(auc_matrix),
                              group.by = "seurat_clusters")

# SCENIC+整合scATAC-seq和scRNA-seq数据
# 从开放染色质区域发现增强子-TF-靶基因关系

# 安装
# pip install scenicplus

# SCENIC+的额外功能：
# 1. 整合scATAC-seq鉴定增强子区域
# 2. TF-增强子-靶基因三元调控关系
# 3. 预测TF扰动的效果
# 4. GRN velocity（调控网络的动态变化）

# SCENIC+需要：
# - 配对的scRNA-seq和scATAC-seq数据
# - 或多模态数据（10x Multiome）

# pySCENIC三步流程
pyscenic grn → adjacencies.csv（GRN推断）
pyscenic ctx → regulons.csv（Motif验证）
pyscenic aucell → output.loom（活性打分）

# 所需文件
表达矩阵: input.loom（原始计数）
TF列表: allTFs_hg38.txt
排名数据库: *.feather（~25GB）
Motif注释: *.tbl

# 关键参数
grn --num_workers: 线程数（越多越快）
grn --sparse: 稀疏模式（节省内存）
ctx --mask_dropouts: 忽略dropout零值
aucell: 默认参数通常足够

# 结果解读
AUCell分数高: regulon在该细胞中活跃
二值化"开": regulon分数超过阈值
RSS高: regulon对该细胞类型特异

# SCENIC版本选择
R版SCENIC: 已弃用，不再维护
pySCENIC: 推荐使用，Python实现，快10倍
SCENIC+: 整合scATAC-seq，最新版本

# 计算资源建议
GRNBoost2: 8-16核，16-64GB内存
cisTarget: I/O密集，推荐SSD
AUCell: 相对快速
总时间: 5000细胞约2-4小时