> 作者 (author): LKP <kunpeng.liao@abiosciences.com>
> Skill 及脚本参数来源: Chen et al., 2025, Cancer Cell 43, 1656–1676, <https://doi.org/10.1016/j.ccell.2025.06.020>

SCENIC 转录调控网络推断（R 版 1.3.1）

GENIE3 共表达 → RcisTarget motif 富集 → AUCell 打分，产出 regulon 与细胞活性矩阵。

安装与数据库

# 包（SCENIC 1.3.1 + AUCell + RcisTarget）
install.packages(c("AUCell","RcisTarget"), repos = "https://mirrors.tuna.tsinghua.edu.cn/CRAN/")
remotes::install_github("aertslab/SCENIC", ref = "v1.3.1")
# 数据库（hg38 500bp 版，~1.2GB；官方resources.aertslab.org 的 old/ 目录）
# hg38__refseq-r80__500bp_up_and_100bp_down_tss.mc9nr.feather 放 dbDir 下

0. 已知 bug 修复（Windows 环境）

# RcisTarget 的 motifAnnotations_hgnc.RData 内部对象名是 motifAnnotations（不带后缀），
# initializeScenic 内部 data()+eval 找不到，需手动加载并重命名：
rn <- system.file("data", "motifAnnotations_hgnc.RData", package = "RcisTarget")
e <- new.env(); obj <- load(rn, envir = e)
assign("motifAnnotations_hgnc", get(obj, envir = e), envir = .GlobalEnv)

1. 输入准备（表达矩阵 + 细胞注释）

exprMat <- as.matrix(GetAssayData(sub_obj, assay = "RNA", layer = "counts"))  # genes × cells
# 输入前基因过滤：count>=3 且至少 3 个细胞表达
keep <- rowSums(exprMat >= 3) >= min(3, ncol(exprMat))
cellInfo <- data.frame(celltype = meta$celltype, row.names = colnames(exprMat))

2. 初始化（必须先 setwd 到工作目录！）

setwd(scenic_workdir)  # int/、output/ 相对路径都落在 cwd，先固定工作目录
scenicOptions <- initializeScenic(org = "hgnc", dbDir = dbDir,
                                  dbs = c("hg38__refseq-r80__500bp_up_and_100bp_down_tss.mc9nr.feather"),
                                  datasetTitle = "分析名", nCores = 1)
scenicOptions@settings$seed <- 123
saveRDS(scenicOptions, "scenicOptions.Rds")

3. 官方顺序五步（顺序不可乱）

# Step A: 基因过滤（官方 geneFiltering）
genesKept <- geneFiltering(exprMat, scenicOptions)
exprMat_filtered <- exprMat[genesKept, ]
# Step B: 相关矩阵
runCorrelation(exprMat_filtered, scenicOptions)
# Step C: GENIE3（输入必须 log 化！生成 int/1.4_GENIE3_linkList.Rds）
exprMat_filtered_log <- log2(exprMat_filtered + 1)
runGenie3(exprMat_filtered_log, scenicOptions, resumePreviousRun = FALSE)
# Step D: 共表达模块（注意函数名是 coexNetwork2modules，不带 exprMat 参数）
scenicOptions <- runSCENIC_1_coexNetwork2modules(scenicOptions)
# Step E: regulon（motif 富集；coexMethod 用 top5perTarget）
scenicOptions <- runSCENIC_2_createRegulons(scenicOptions, coexMethod = c("top5perTarget"))
# Step F: AUCell 打分（输入全基因 log 矩阵；跳过重可视化）
exprMat_log <- log2(exprMat + 1)
scenicOptions <- runSCENIC_3_scoreCells(scenicOptions, exprMat_log,
                                        skipBinaryThresholds = TRUE,
                                        skipHeatmap = TRUE, skipTsne = TRUE)

4. 汇总输出（作者对接格式）

regulons   <- readRDS("int/2.6_regulons_asGeneSet.Rds")
regulonAUC <- readRDS("int/3.4_regulonAUC.Rds")
saveRDS(list(regulons = regulons, regulonAUC = regulonAUC), "neuro_scenic_AUC.rds")
saveRDS(cellInfo, "cellInfo.rds")
# adj.csv：GENIE3 linkList（TF/Target/weight）直接 write.csv

5. RSS 分析与绘图

library(SCENIC)
AUCmat <- AUCell::getAUC(regulonAUC)
rownames(AUCmat) <- gsub("[(+)]", "", rownames(AUCmat))  # 去掉 extended (xxg) 后缀
rss <- calcRSS(AUC = AUCmat, cellAnnotation = cellInfo$celltype)
# plotRSS_oneSet 在 1.3.1 有越界 bug（n 超过 regulon 数），手动自绘：
thisRss <- data.frame(regulon = names(sort(rss[, setName], decreasing = TRUE)),
                      rank = seq_len(nrow(rss)))
n <- min(30, nrow(thisRss)); thisRss$regulon[(n+1):nrow(thisRss)] <- NA  # 需 if(n < nrow) 守卫
ggplot(thisRss, aes(rank, rss)) + geom_point(color="blue", size=1) +
  geom_label_repel(aes(label=regulon), na.rm=TRUE) + theme_classic()
# 热图 feature 名转换：regulon 名 "IRF9_extended 36g"/"IRF9 32g" 提取核心 TF：
reg_tf <- ifelse(grepl("_extended", regulon_names),
                 sub("_extended.*", "", regulon_names),
                 sub(" [0-9]+g.*", "", regulon_names))

常见坑

1. File 'int/1.4_GENIE3_linkList.Rds' does not exist：跳过了 runGenie3 或 cwd 不对；必须 setwd 后 initializeScenic，五步顺序执行
2. runGenie3 输入必须 log2(x+1)：raw counts 会产生畸形权重
3. runSCENIC_1 不接受 exprMat 参数（1.3.1 签名），数据通过 int/ 中间文件传递
4. 小样本（<500 细胞）regulon 偏少是正常现象：motif 注释率低（TF 命中率 18% 量级）
5. plotRSS_oneSet 越界：(n+1):nrow 当 n==nrow 产生倒序索引，必须加 if 守卫