SCENIC 与转录因子活性：基因调控网络分析
基因调控网络分析回答「哪些转录因子在驱动当前的表达状态」。一条路线从数据中推断 TF 与靶基因关系，代表是 SCENIC：先用 GRNBoost2 找与 TF 共表达的基因，再用 cisTarget 保留含该 TF motif 的靶基因形成 regulon，最后用 AUCell 给每个细胞打分。另一条路线借用整理好的 TF-靶基因先验集合，如 DoRothEA 或 CollecTRI，经 decoupleR 等工具从靶基因表达推算 TF 活性，bulk 和单细胞数据都适用。
什么时候用
- 在单细胞数据中寻找各细胞类型或状态特异的核心转录因子
- 比较疾病与对照、处理前后哪些转录因子活性发生改变
- 转录因子自身表达低或变化小，想从靶基因层面推断其活性
- 解释细胞分化或轨迹中调控程序的切换
- 为 ChIP-seq、CUT&Tag 或敲低实验挑选候选转录因子
不适合的情况
- 未经结合或扰动实验验证，就把推断的 TF-靶基因关系当作直接结合
- 物种缺少 cisTarget 数据库或 TF 先验集合时，结果只能依赖同源推断
- 在只有几十个细胞的亚群上单独运行 SCENIC，共表达估计不稳定
需要准备的数据
- 表达矩阵
- 单细胞用计数或标准化矩阵并完成质控；bulk 用标准化表达或差异统计量
- 细胞或样本注释
- 细胞类型、聚类或分组标签，用于计算 regulon 特异性和组间比较
- TF 列表与 cisTarget 数据库
- SCENIC 需要与物种、基因组版本匹配的 motif 排序数据库和 motif 注释文件
- TF 先验集合
- DoRothEA（带 A–E 置信等级）或 CollecTRI，可通过 decoupleR 获取，注意物种
分析步骤
- 1
确定分析路线
有单细胞数据、想发现细胞类型特异的调控单元时选 SCENIC；想快速比较已知 TF 的活性变化时选先验集合加 decoupleR。两条路线可以先后使用、互相印证。
- 2
预处理表达矩阵
单细胞数据先完成质控和细胞注释；运行 SCENIC 前常过滤只在极少数细胞中检测到的基因。数据集很大时可按细胞类型均衡抽样，避免稀有类型被淹没。
- 3
推断共表达关系
pySCENIC 的 grn 步骤用 GRNBoost2 计算每个 TF 对每个基因的重要性得分。这一步带有随机性，常见做法是多次运行，只保留在多数运行中重复出现的 TF-靶基因关系。
- 4
motif 修剪
ctx 步骤用 cisTarget 数据库检查共表达靶基因是否富集该 TF 的 motif，剔除缺乏 motif 支持的间接靶点，得到 regulon。检查 regulon 数量和大小，过小的 regulon 不稳定。
- 5
regulon 活性打分
aucell 步骤为每个细胞计算各 regulon 的 AUC 分数，可再二值化为开或关。用 AUC 矩阵重新降维，若细胞类型能按 regulon 活性分开，说明网络捕捉到了真实的细胞身份信息。
- 6
推断 TF 活性
在 decoupleR 中用 ulm 或 mlm 等方法，结合 DoRothEA（常取 A、B、C 级）或 CollecTRI，计算每个细胞、样本或差异比较中的 TF 活性分数，正值表示靶基因整体被激活。
- 7
特异性与可视化
计算 regulon 特异性评分（RSS）找出各细胞类型的特异 TF，画热图和 UMAP 活性图。挑选候选时同时看 TF 自身表达和靶基因活性是否一致，并核对已知谱系因子。
关键参数与决策
路线选择
SCENIC 依赖本数据的共表达，能发现新的调控关系，但耗时且有随机性；先验集合方法速度快、跨数据集可比，但只能评估已收录的 TF。由研究目的决定。
DoRothEA 置信等级
A 级证据最强，E 级主要来自计算推断。常用 A、B、C 三级平衡覆盖度与可靠性；只用 A 级时可分析的 TF 数量很少。
靶基因最小数量
靶基因太少的 TF 活性分数波动大。decoupleR 的 minsize 参数默认是 5，会过滤掉在数据中匹配到的靶基因少于 5 个的 TF。
cisTarget 数据库版本
motif 排序数据库要与物种和基因组版本一致，并区分 TSS 周围的搜索范围。基因名与数据库不匹配时，大量 TF 会被直接丢弃。
重复运行次数
GRNBoost2 每次结果不同，常见做法是重复运行若干次，只保留在多数运行中出现的 regulon。次数越多结果越稳，计算成本也成倍增加。
AUC 二值化阈值
AUCell 可为每个 regulon 自动确定开或关的阈值，但在双峰不明显时阈值很不稳定。建议画出 AUC 分布检查阈值位置，必要时手动调整，并在方法中写明。
常见错误
- 把 TF 表达量等同于 TF 活性：许多 TF 受磷酸化、核转位等调控，mRNA 变化小但靶基因变化大，应以靶基因层面的活性为准。
- SCENIC 只运行一次就报告：GRNBoost2 的随机性会让处在边缘的 regulon 时有时无，单次结果难以复现。
- 基因名与 cisTarget 数据库不一致：如小鼠数据配人类数据库，或使用 Ensembl ID，大量基因匹配失败，regulon 数显著减少。
- 只用差异基因推断 TF 活性：先验集合方法需要完整的统计量分布，截取后活性分数会被系统性高估。
- 把 regulon 靶基因当作已证实的直接结合靶点：motif 富集只说明可能结合，仍需 ChIP-seq、CUT&Tag 或扰动实验确认。
结果怎么看
regulon 活性热图
行为 regulon，列为细胞类型的平均 AUC 或 RSS。理想情况下每种细胞类型有一组特异高活性的 regulon；若各类型活性都差不多，可能是 regulon 过大或细胞注释不准。
regulon 特异性评分图
每种细胞类型按 RSS 排序，排在前面的是最特异的 TF。检查其中是否包含已知谱系因子，例如 B 细胞中出现 PAX5，已知因子排名靠前说明分析结果可信。
AUC 活性 UMAP
把单个 regulon 的 AUC 投到 UMAP 上，看活性是否集中在特定亚群。与该 TF 自身表达分布对比，两者一致增强可信度；不一致时需要结合生物学背景解释。
TF 活性分数表
decoupleR 输出每个 TF 在每个细胞、样本或比较中的分数和 P 值。正值表示靶基因整体被激活，负值表示被抑制；结合 TF 自身表达方向，可判断调控是否发生在转录水平。
二值化 regulon 热图
每个细胞的 regulon 只有开和关两种状态，按细胞类型排列后可直观看出哪些 regulon 在特定亚群中被开启。若大部分 regulon 在所有细胞中都处于开启状态，说明二值化阈值偏低。
和相近方法怎么选
SCENIC 与 DoRothEA/decoupleR
SCENIC 从本数据构建新的 regulon，适合发现细胞类型特异的调控程序，但计算慢；DoRothEA 等先验集合配合 decoupleR 计算快，可在 bulk 和单细胞中统一比较，但只能评估已收录的 TF。
TF 活性推断与 GSEA
两者都基于基因集打分。TF 活性方法使用带调控方向（激活或抑制）和权重的靶基因集合，推断的是某个调控因子的输出；GSEA 使用无方向的通路基因集，回答通路整体是否变化。
调控网络与 WGCNA
WGCNA 找的是共表达模块，不区分调控者和被调控者；调控网络以转录因子为中心，并用 motif 或先验知识约束关系，更接近可以用实验检验的调控假设。
SCENIC 与 SCENIC+
SCENIC 只使用转录组；SCENIC+ 结合单细胞 ATAC-seq 的染色质可及性，把增强子纳入 TF 与靶基因之间，适合同时测了 RNA 和 ATAC 的单细胞多组学数据。
常见问题
SCENIC 和 pySCENIC 有什么区别？
SCENIC 最早以 R 实现，pySCENIC 是 Python 版本，计算更快，适合大规模单细胞数据，目前更常用。两者流程一致：共表达推断、motif 修剪、AUCell 打分，结果可导入 R 或 Python 继续可视化。
转录因子活性和表达量有什么不同？
表达量是 TF 自身 mRNA 的多少；活性是根据其靶基因整体表达变化推断出的调控输出。由于翻译后修饰等原因，两者常不一致，活性推断能发现表达变化不明显但功能已经改变的 TF。
DoRothEA 和 CollecTRI 选哪个？
两者都是 TF-靶基因先验集合。CollecTRI 整合的来源更多、覆盖的 TF 更多，decoupleR 的官方教程目前以它为默认示例；DoRothEA 的优势是有 A–E 置信等级可供筛选。两者可以在同一数据上对照。
SCENIC 需要多少细胞？
没有固定门槛，但每种细胞类型需要足够的细胞才能稳定估计共表达。大数据集常按细胞类型抽样以节省计算；稀有亚群细胞过少时，其 regulon 结果应谨慎解读，最好用其他方法交叉验证。
bulk RNA-seq 能做转录因子活性分析吗？
可以。用 decoupleR 配合 DoRothEA 或 CollecTRI，把差异分析的 t 值或 stat 作为输入，推算各 TF 在比较中的活性变化；也可以对每个样本的标准化表达打分，再做组间比较。
在 Codar 里怎么做
在 Codar 中，上传已注释的单细胞数据或 bulk 表达与差异结果，说明物种与目标。Codar 会先确认用调控网络推断还是基于先验集合的 TF 活性推断；单细胞数据可完成 regulon 构建、活性打分和细胞类型特异性分析。交付 regulon 列表、活性矩阵、特异性热图、TF 活性表和报告。
这是我已注释的小鼠肺纤维化单细胞数据（h5ad），想用 pySCENIC 找各成纤维细胞亚群的特异转录因子，并比较疾病组与对照的 TF 活性。