Single Cell——转录因子调控网络分析（SCENIC）

转录因子调控网络

SCENIC (Single-Cell rEgulatory Network Inference and Clustering) 是从单细胞RNA数据推断基因调控网络及其相关细胞状态的工具。

原作者将SCENIC应用于肿瘤和小鼠大脑单细胞图谱数据，证明了顺式调控网络分析能够有助于深入挖掘细胞异质性背后的生物学意义，并为疾病的诊断、治疗以及发育分化的研究提供有价值的线索。

SCENIC在2017年首先发表于nature methods，2020年又将流程整理后发表于nature protocls。有R版和python版本。

SCENIC 工作流程

官方手册

vignetteFile <- file.path(system.file('doc', package='SCENIC'), "SCENIC_Running.Rmd")
file.copy(vignetteFile, "SCENIC_myRun.Rmd")
# or: 
vignetteFile <- "https://raw.githubusercontent.com/aertslab/SCENIC/master/vignettes/SCENIC_Running.Rmd"
download.file(vignetteFile, "SCENIC_myRun.Rmd")

构建基因调控网络（GRN）：

• 根据共表达确定每个 TF 的潜在靶基因：
  • 过滤表达式矩阵并运行 GENIE3/GRNBoosrt
  • 将GENIE3/GRNBoost中的目标格式化为共表达模块。
• 基于DNA-motif分析选择潜在的直接结合靶标（调控子）（RcisTarget: TF motif analysis）

识别细胞状态及其调节因子：

• 分析每个单独单元（AUCell）中的网络活动
  • 细胞中的评分规则（计算AUC)
  • 将网络活动转换为二进制活动矩阵
• 基于基因调控网络活性（细胞聚类）鉴定稳定的细胞状态并探索结果

准备包和数据

装包

SCENIC的workflow主要依赖三个R包的支持：

• GENIE3：用于共表达网络的计算(也可以用GRNBoost2替代)

• RcisTarget：转录因子结合motif的计算

• AUCell：识别单细胞测序数据中的激活基因集(gene-network)

  dir.create(“./SCENIC”) #构建SCENIC分析文件夹
  setwd(“SCENIC”)
  dir.create(“D:/R/seurat/SCENIC/cisTarget_databases”)#构建数据库文件夹
  setwd(“D:/R/seurat/SCENIC/cisTarget_databases”)

  #导入相关R包
  if(!require(SCENIC))
  BiocManager::install(c(“AUCell”, “RcisTarget”),ask = F,update = F);
  BiocManager::install(c(“GENIE3”),ask = F,update = F)#这三个包显然是必须安装的
  if(!require(SCopeLoomR))
  devtools::install_github(“aertslab/SCopeLoomR”, build_vignettes = TRUE)#SCopeLoomR用于获取测试数据

  if (!requireNamespace(“arrow”, quietly = TRUE))
  BiocManager::install(‘arrow’)

  if(!require(SCENIC))
  devtools::install_github(“aertslab/SCENIC”)

  if (!requireNamespace(“BiocManager”, quietly = TRUE))

  install.packages(“BiocManager”)

  options(“repos” = c(CRAN=“https://mirrors.tuna.tsinghua.edu.cn/CRAN/”))

  options(BioC_mirror=“https://mirrors.tuna.tsinghua.edu.cn/bioconductor”)

  install.packages(“feather”)
  install.packages(“doRNG”)
  install.packages(“KernSmooth”)
  install.packages(“doMC”, repos=“http://R-Forge.R-project.org”)
  install.packages(“DT”)

  install.packages(“shiny”)

  BiocManager::install(‘shiny’)

  #bioconductor版本小于4.0或你的R（3.6）的版本也比较老，你可能得试试下面的方法进行安装
  devtools::install_github(“aertslab/SCENIC@v1.1.2”)
  devtools::install_github(“aertslab/AUCell”)
  devtools::install_github(“aertslab/RcisTarget”)
  devtools::install_github(“aertslab/GENIE3”)

加载数据

输入数据的表达矩阵为gene-summarized counts，不同于其他进阶分析的是，输入的数据最好是raw counts或normalized counts。计算后的TPM或FPKM/RPKM之类也是允许的, 但是这些可能会导致人为的共变量引入，避免节外生枝，推荐用raw count。但总的来说，** raw (logged) UMI counts**, normalized UMI counts, 和TPM都能产生值得信赖的效果

### Load data
loomPath <- system.file(package="SCENIC", "examples/mouseBrain_toy.loom")
library(SCopeLoomR)
loom <- open_loom(loomPath)#获取loom文件
exprMat <- get_dgem(loom)#取出表达信息
cellInfo <- get_cell_annotation(loom)#获取注释信息
close_loom(loom)

而Python版的SCENIC是推荐使用loom作为输入数据的,在参考文章SCENIC单细胞转录因子预测|3.软件安装与数据准备中会教大家如何利用表达矩阵和注释信息创建loom文件(to loom)

SCENIC工作流的主要命令的概述

在后面会有更详细的步骤

###初始化 SCENIC 设置,设置分析环境
library(SCENIC)
#options(max.connections = 256)#增加最大连接数
data(list="motifAnnotations_mgi_v9", package="RcisTarget")
motifAnnotations_mgi <- motifAnnotations_mgi_v9
scenicOptions <- initializeScenic(org="mgi", dbDir="D:/R/seurat/SCENIC/cisTarget_databases", nCores=10)
# scenicOptions@inputDatasetInfo$cellInfo <- "int/cellInfo.Rds"
saveRDS(scenicOptions, file="int/scenicOptions.Rds") 

##### 构建Co-expression network
#基因过滤/选择，去除最有可能是噪音的基因
genesKept <- geneFiltering(exprMat, scenicOptions,minCountsPerGene=3*.01*ncol(exprMat),minSamples=ncol(exprMat)*.01)
exprMat_filtered <- exprMat[genesKept, ]
runCorrelation(exprMat_filtered, scenicOptions)##计算相关性矩阵，1.2_corrMat.Rds：基因之间的相关性矩阵
exprMat_filtered_log <- log2(exprMat_filtered+1) 
runGenie3(exprMat_filtered_log, scenicOptions)

### 建立GRN并评分
exprMat_log <- log2(exprMat+1)
scenicOptions@settings$dbs <- scenicOptions@settings$dbs["10kb"]
scenicOptions <- runSCENIC_1_coexNetwork2modules(scenicOptions)
scenicOptions <- runSCENIC_2_createRegulons(scenicOptions, coexMethod=c("top5perTarget"))
scenicOptions <- runSCENIC_3_scoreCells(scenicOptions, exprMat_log)
scenicOptions <- runSCENIC_4_aucell_binarize(scenicOptions)
tsneAUC(scenicOptions, aucType="AUC") 

# 输出:
# saveRDS(cellInfo, file=getDatasetInfo(scenicOptions, "cellInfo")) # Temporary, to add to loom
export2loom(scenicOptions, exprMat)

# 若要保存当前状态或任何设置更改，请再次保存对象：
saveRDS(scenicOptions, file="int/scenicOptions.Rds") 

# /Step2_MotifEnrichment_preview.html in detail/子集:
motifEnrichment_selfMotifs_wGenes <- loadInt(scenicOptions, "motifEnrichment_selfMotifs_wGenes")
tableSubset <- motifEnrichment_selfMotifs_wGenes[highlightedTFs=="Sox8"]
viewMotifs(tableSubset) 

# output/Step2_regulonTargetsInfo.tsv in detail: 
regulonTargetsInfo <- loadInt(scenicOptions, "regulonTargetsInfo")
tableSubset <- regulonTargetsInfo[TF=="Stat6" & highConfAnnot==TRUE]
viewMotifs(tableSubset)

# 细胞类型特异性调节因子 (RSS): 
regulonAUC <- loadInt(scenicOptions, "aucell_regulonAUC")
rss <- calcRSS(AUC=getAUC(regulonAUC), cellAnnotation=cellInfo[colnames(regulonAUC), "CellType"], )
rssPlot <- plotRSS(rss)
plotly::ggplotly(rssPlot$plot) 

详细步骤从这开始

构建目录

在此流程中，我们将保存多个文件。为了保持它们的整洁，建议将工作目录设置为一个新文件夹。

dir.create("SCENIC_MouseBrain")
setwd("SCENIC_MouseBrain")

SCENIC 的主要输出存储在output文件夹中，其中还包括一些自动生成的图和报告，可以使用这些图和报告对结果进行概述。

此外，一些中间/临时文件将被保存到int文件夹中，并带有编号前缀以保持它们的顺序。可以使用这些文件来检查每个步骤的详细信息，或者使用不同的设置重新运行部分分析。

输入

SCENIC的输入是单细胞RNA-seq表达矩阵（以基因符号作为行名）。第一步是加载这个矩阵。

在本教程中，使用了一个示例，只有来自小鼠大脑的200个细胞和<1000个基因：

loomPath <- system.file(package="SCENIC", "examples/mouseBrain_toy.loom")

打开loom文件并加载表达式矩阵（如果可用的话，还有单元格注释）

library(SCopeLoomR)
#可直接下载文件
#download.file("http://loom.linnarssonlab.org/clone/Previously%20Published/Cortex.loom", "Cortex.loom")
#loomPath <- "Cortex.loom"
loom <- open_loom(loomPath)
exprMat <- get_dgem(loom)
cellInfo <- get_cell_annotation(loom)
close_loom(loom)

dim(exprMat)

细胞信息/ phenodata

在步骤3-4 （GRN评分和聚类）中，将结果与已知的细胞信息进行比较是很有趣的。您已经可以指示要绘制的变量，并为它们分配特定的颜色（否则将自动分配一个）。

head(cellInfo)
cellInfo <- data.frame(cellInfo)
cbind(table(cellInfo$CellType))
dir.create("int")
saveRDS(cellInfo, file="int/cellInfo.Rds")

# 要分配给变量的颜色（与NMF::aheatmap格式相同）
colVars <- list(CellType=c("microglia"="forestgreen", 
                           "endothelial-mural"="darkorange", 
                           "astrocytes_ependymal"="magenta4", 
                           "oligodendrocytes"="hotpink", 
                           "interneurons"="red3", 
                           "pyramidal CA1"="skyblue", 
                           "pyramidal SS"="darkblue"))
colVars$CellType <- colVars$CellType[intersect(names(colVars$CellType), cellInfo$CellType)]
saveRDS(colVars, file="int/colVars.Rds")
plot.new(); legend(0,1, fill=colVars$CellType, legend=names(colVars$CellType))

初始化SCENIC设置

SCENIC包中的大多数函数使用一个公共对象来存储当前运行的选项。该对象取代了大多数函数的“参数”，应该在运行SCENIC开始时使用initializeScenic（）函数创建。

默认设置对于大多数分析都是有效的。在所有运行中需要指定的参数是生物体（mgi用于鼠标，hgnc用于人类，dmel用于苍蝇），以及RcisTarget数据库存储的目录

#初始化SCENIC设置
library(SCENIC)
org <- "mgi" # or hgnc, or dmel
dbDir <- "D:/R/seurat/SCENIC/cisTarget_databases" # RcisTarget数据库位置
myDatasetTitle <- "SCENIC example on Mouse brain" # 为您的分析选择一个名称
data(defaultDbNames)
dbs <- defaultDbNames[[org]]
scenicOptions <- initializeScenic(org=org, dbDir=dbDir, dbs=dbs, datasetTitle=myDatasetTitle, nCores=10) 

# 按需修改
scenicOptions@inputDatasetInfo$cellInfo <- "int/cellInfo.Rds"
scenicOptions@inputDatasetInfo$colVars <- "int/colVars.Rds"
# 数据库:
# scenicOptions@settings$dbs <- c("mm9-5kb-mc8nr"="mm9-tss-centered-5kb-10species.mc8nr.feather")
# scenicOptions@settings$db_mcVersion <- "v8"

# 保存以供后续使用....
saveRDS(scenicOptions, file="int/scenicOptions.Rds") 

共表达网络

SCENIC工作流程的第一步是根据表达数据推断潜在的转录因子靶标。

GENIE3或GRNBoost。这两种工具的输入都是表达矩阵（过滤后的）和转录因子列表（潜在调节因子）。GENIE3/GRBBoost的输出和一个相关矩阵将用于创建共同表达式模块（runSCENIC_1_coexNetwork2modules()）。

• GENIE3 (Huynh-Thu et al.(2010))允许识别非线性关系，即使它们只存在于样本子集中，并且它在网络推理DREAM5挑战中表现最佳（Marbach et al.(2012)）。GENIE3可以很容易地在R中运行。然而，GENIE3非常耗费时间和计算量（在3-5k单元的数据集上需要几个小时或几天）。

• GRNBoost在很短的时间内提供了与GENIE3相似的结果

基因筛选/选择

使用软基因过滤器，以去除表达水平非常低或在太少细胞中的基因。在这里，根据基因计数的总数和检测到的细胞数量应用过滤。

• 过滤每个基因的总读取数
  这种过滤器旨在去除最可能是噪音的基因。默认情况下，它只保留所有样本中至少有6个UMI计数的基因（例如，如果基因在1%的细胞中以3的值表达，则该基因的总数）。根据数据集调整此值（minCountsPerGene）

• 通过检测到基因的细胞数量进行筛选（例如>0 （UMI），或>log2 (TPM)）
  默认情况下（minSamples），在至少1%的细胞中检测到的基因被保留。这种过滤是为了去除来自少数“嘈杂”细胞的基因（只在一个或很少的细胞中表达的基因，如果它们碰巧在一个给定的细胞中重合，就会增加很多重量）。建议设置比要检测的最小细胞群低的百分比。

• 保留在RcisTarget数据库中可用的基因
  为了节省GENIE3/GRNBoost的运行时间，因为数据库中不可用的基因将不会在接下来的步骤中使用。

  genesKept <- geneFiltering(exprMat, scenicOptions=scenicOptions,
  minCountsPerGene=3*.01ncol(exprMat),
  minSamples=ncol(exprMat).01)

  #检查是否有已知的相关基因被过滤掉
  interestingGenes <- c(“Sox9”, “Sox10”, “Dlx5”)
  interestingGenes[which(!interestingGenes %in% genesKept)]

  exprMat_filtered <- exprMat[genesKept, ]
  dim(exprMat_filtered)
  #去除exprMat不会用于共表达式分析
  rm(exprMat)

相关性

GENIE3/GRNBoost可以检测正相关和负相关目标，以区分潜在激活和抑制（默认TF与潜在目标之间的Spearman相关）

runCorrelation(exprMat_filtered, scenicOptions)
# Optional: add log (if it is not logged/normalized already)
exprMat_filtered <- log2(exprMat_filtered+1)

GENIE3

runGenie3函数将使用默认设置运行GENIE3，这对于大多数数据集来说通常是足够的，使用RcisTarget数据库中可用的转录因子作为候选调节因子。

由于GENIE3基于随机森林方法，每次运行结果都会略有不同。使用的树数（n棵树）越多，可变性就越低。

GENIE3通常需要几个小时（或几天）才能运行。

#跑GENIE3
runGenie3(exprMat_filtered, scenicOptions)

建立GRN并评分

主要步骤：

• 构建基因调控网络：
  • 获取共表达式模块
  • 获取调控子使用RcisTarge： TF基序分析
• 识别细胞状态：
  • 用AUCell对细胞中的GRN（调控子）进行评分。根据GRN活性将细胞聚集

如果需要，重新加载表达式矩阵：

loom <- open_loom(loomPath)
exprMat <- get_dgem(loom)
close_loom(loom)
# Optional:log表达式（对于TF表达式图，不影响其他任何计算）
exprMat_log <- log2(exprMat+1)
dim(exprMat)


library(SCENIC)
scenicOptions <- readRDS("D:/R/seurat/SCENIC/cisTarget_databases/int/scenicOptions.Rds")
scenicOptions@settings$verbose <- TRUE
scenicOptions@settings$nCores <- 10
scenicOptions@settings$seed <- 123

# 缩短时间: 
# coexMethod=c("top5perTarget")
scenicOptions@settings$dbs <- scenicOptions@settings$dbs["10kb"] 

scenicOptions <- runSCENIC_1_coexNetwork2modules(scenicOptions)
#closeAllConnections()
scenicOptions <- runSCENIC_2_createRegulons(scenicOptions, coexMethod=c("top5perTarget")) 
scenicOptions <- runSCENIC_3_scoreCells(scenicOptions, exprMat_log)

saveRDS(scenicOptions, file="int/scenicOptions.Rds") 

将网络活动二值化(可选做）

#将网络活动二值化
library(AUCell)

aucellApp <- plotTsne_AUCellApp(scenicOptions, exprMat_log)#该函数里的内置函数could not find function "AUCell_createViewerApp"
savedSelections <- shiny::runApp(aucellApp)

# Save the modified thresholds:
newThresholds <- savedSelections$thresholds
scenicOptions@fileNames$int["aucell_thresholds",1] <- "int/newThresholds.Rds"
saveRDS(newThresholds, file=getIntName(scenicOptions, "aucell_thresholds"))
saveRDS(scenicOptions, file="int/scenicOptions.Rds")
scenicOptions <- runSCENIC_4_aucell_binarize(scenicOptions)

调节子活动的聚类/降维

可以根据规则活动对细胞进行分组/聚类，可以是连续的，也可以是二值化的

如果使用t-SNE作为可视化，建议尝试不同的设置来评估状态/集群的稳定性。

nPcs <- c(5) 
scenicOptions@settings$seed <- 123 # 设置相同的种子数
# 运行不同设置的t-SNE:
fileNames <- tsneAUC(scenicOptions, aucType="AUC", nPcs=nPcs, perpl=c(5,15,50))
fileNames <- tsneAUC(scenicOptions, aucType="AUC", nPcs=nPcs, perpl=c(5,15,50), onlyHighConf=TRUE, filePrefix="int/tSNE_oHC")
#pdf格式的绘图（单个文件在int/中）:
fileNames <- paste0("int/",grep(".Rds", grep("tSNE_", list.files("int"), value=T), value=T))
par(mfrow=c(length(nPcs), 3))
fileNames <- paste0("int/",grep(".Rds", grep("tSNE_AUC", list.files("int"), value=T, perl = T), value=T))
plotTsne_compareSettings(fileNames, scenicOptions, showLegend=FALSE, varName="CellType", cex=.5)

par(mfrow=c(3,3))
fileNames <- paste0("int/",grep(".Rds", grep("tSNE_oHC_AUC", list.files("int"), value=T, perl = T), value=T))
plotTsne_compareSettings(fileNames, scenicOptions, showLegend=FALSE, varName="CellType", cex=.5)

选择的t-SNE可以保存为默认值以用于绘图（也可以是“二进制”）：

scenicOptions@settings$defaultTsne$aucType <- "AUC"
scenicOptions@settings$defaultTsne$dims <- 5
scenicOptions@settings$defaultTsne$perpl <- 15
saveRDS(scenicOptions, file="int/scenicOptions.Rds")

导出到loom/范围

可以使用export2loom（）函数创建.loom文件（需要使用SCopeLoomR包）。这个函数将SCENIC的主要结果保存到一个.loom文件中：

• 调节子

• 调节活性（AUC矩阵和阈值）

• 嵌入（例如规则活动上的t-SNE和UMAP）
  motif富集分析和共表达模块（例如GRNBoost/GENIE3输出）存储在独立的文本文件中（主要是因为它们的大小更大）

  exprMat <- get_dgem(open_loom(loomPath))
  dgem <- exprMat
  head(colnames(dgem)) #should contain the Cell ID/name

  Export:

  scenicOptions@fileNames$output[“loomFile”,] <- “output/mouseBrain_SCENIC.loom”
  export2loom(scenicOptions, exprMat)

从.loom文件加载结果

library(SCopeLoomR)
scenicLoomPath <- getOutName(scenicOptions, "loomFile")
loom <- open_loom(scenicLoomPath)

# Read information from loom file:
regulons_incidMat <- get_regulons(loom)
regulons <- regulonsToGeneLists(regulons_incidMat)
regulonsAUC <- get_regulons_AUC(loom)
regulonsAucThresholds <- get_regulon_thresholds(loom)
embeddings <- get_embeddings(loom)

探索/解释结果

输出文件夹包含几个文件，这些文件提供了每个步骤的结果概述。可以通过中间文件（保存在int文件夹中，可以通过loadInt（scenicOptions）列出）更详细地研究这些结果。下面是几个研究方向

细胞状态

AUCell提供跨细胞的调控活性。通过基于这种调节子活性的细胞聚类（无论是连续的还是二元的AUC矩阵），我们可以看到是否有一组细胞倾向于具有相同的调节子活性，并揭示在多个细胞中反复出现的网络状态。

将这些聚类与不同的可视化方法相结合，我们可以探索细胞状态与特定规则的关联。

将AUC和TF表达投影到t-SNEs上

t-SNE是细胞的二维投影，如果细胞（点）具有相似的输入特征（调节活动），则它们彼此靠近。

t-SNE可以很好地识别不同的类，但它不适用于动态/连续过程（例如轨迹类可视化）

exprMat_log <- exprMat # 最好是有先验知识/规范化
aucellApp <- plotTsne_AUCellApp(scenicOptions, exprMat_log) # 默认是t-SNE
savedSelections <- shiny::runApp(aucellApp)

print(tsneFileName(scenicOptions))

tSNE_scenic <- readRDS(tsneFileName(scenicOptions))
aucell_regulonAUC <- loadInt(scenicOptions, "aucell_regulonAUC")

# Show TF expression:
par(mfrow=c(2,3))
AUCell::AUCell_plotTSNE(tSNE_scenic$Y, exprMat, aucell_regulonAUC[onlyNonDuplicatedExtended(rownames(aucell_regulonAUC))[c("Dlx5", "Sox10", "Sox9","Irf1", "Stat6")],], plots="Expression")

#install.packages("Cairo")
library(Cairo)
Cairo::CairoPDF("output/Step4_BinaryRegulonActivity_tSNE_colByAUC.pdf", width=20, height=15)
par(mfrow=c(4,6))
AUCell::AUCell_plotTSNE(tSNE_scenic$Y, cellsAUC=aucell_regulonAUC, plots="AUC")
dev.off()

密度图用于检测最可能的稳定状态（t-SNE中的高密度区域）：

library(KernSmooth)
library(RColorBrewer)
dens2d <- bkde2D(tSNE_scenic$Y, 1)$fhat
image(dens2d, col=brewer.pal(9, "YlOrBr"), axes=FALSE)
contour(dens2d, add=TRUE, nlevels=5, drawlabels=FALSE)

同时显示几个调节子：

par(mfrow=c(1,2))

regulonNames <- c( "Dlx5","Sox10")
cellCol <- SCENIC::plotEmb_rgb(scenicOptions, regulonNames, aucType="AUC", aucMaxContrast=0.6)

regulonNames <- list(red=c("Sox10", "Sox8"),
                     green=c("Irf1"),
                     blue=c( "Tef"))
cellCol <- SCENIC::plotEmb_rgb(scenicOptions, regulonNames, aucType="Binary")

GRN：调控靶标和基序

包含在调节子中的基因：

regulons <- loadInt(scenicOptions, "regulons")
regulons[c("Dlx5", "Irf1")]

只有含有10个或更多基因的调节子才会被AUCell评分：

regulons <- loadInt(scenicOptions, "aucell_regulons")
head(cbind(onlyNonDuplicatedExtended(names(regulons))))

关于TF-target链接的详细信息：
对于每个TF-target对，中间步骤的统计信息总结在loadInt（scenicOptions, “regulonTargetsInfo“）中（保存为文本：getOutName(scenicOptions, ”s2_regulonTargetsInfo”): output/Step2_regulonTargetsInfo.tsv）。该表可用于研究对特定链接的支持。因为它通常包含几千行（在这次运行中：1276行），所以在大多数情况下，建议在将其导出为HTML之前对其进行子集处理。

regulonTargetsInfo <- loadInt(scenicOptions, "regulonTargetsInfo")
tableSubset <- regulonTargetsInfo[TF=="Stat6" & highConfAnnot==TRUE]
viewMotifs(tableSubset, options=list(pageLength=5)) 

支持这些调控的TF基序的完整列表可以在RcisTarget基序富集结果（共表达模块）中看到。这些都保存在motifEnrichment_selfMotifs_wGenes中。这些结果的预览在output/ step2_motifenrichment .html中导出为html（并在output/Step2_MotifEnrichment.tsv中导出为文本）

motifEnrichment_selfMotifs_wGenes <- loadInt(scenicOptions, "motifEnrichment_selfMotifs_wGenes")
tableSubset <- motifEnrichment_selfMotifs_wGenes[highlightedTFs=="Dlx5"]
viewMotifs(tableSubset) 

对已知细胞类型或细胞簇的调节因子

SCENIC的调控分析可以与其他分析相结合（通常是聚类分析），或者专注于特定细胞类型的调控。

• 按簇划分的平均规则活动

  #cellInfo <- data.frame(seuratCluster=Idents(seuratObject))
  regulonAUC <- loadInt(scenicOptions, “aucell_regulonAUC”)
  regulonAUC <- regulonAUC[onlyNonDuplicatedExtended(rownames(regulonAUC)),]
  regulonActivity_byCellType <- sapply(split(rownames(cellInfo), cellInfo$CellType),
  function(cells) rowMeans(getAUC(regulonAUC)[,cells]))
  regulonActivity_byCellType_Scaled <- t(scale(t(regulonActivity_byCellType), center = T, scale=T))

  ComplexHeatmap::Heatmap(regulonActivity_byCellType_Scaled, name=“Regulon activity”)

topRegulators <- reshape2::melt(regulonActivity_byCellType_Scaled)
colnames(topRegulators) <- c("Regulon", "CellType", "RelativeActivity")
topRegulators <- topRegulators[which(topRegulators$RelativeActivity>0),]
viewTable(topRegulators)

• 二值化方式（该细胞类型/集群中具有活性调控子的细胞百分比）

  minPerc <- .7
  binaryRegulonActivity <- loadInt(scenicOptions, “aucell_binary_nonDupl”)
  cellInfo_binarizedCells <- cellInfo[which(rownames(cellInfo)%in% colnames(binaryRegulonActivity)), drop=FALSE]
  regulonActivity_byCellType_Binarized <- sapply(split(rownames(cellInfo_binarizedCells), cellInfo_binarizedCells$CellType),
  function(cells) rowMeans(binaryRegulonActivity[,cells, drop=FALSE]))
  binaryActPerc_subset <- regulonActivity_byCellType_Binarized[which(rowSums(regulonActivity_byCellType_Binarized>minPerc)>0),]
  ComplexHeatmap::Heatmap(binaryActPerc_subset, name=“Regulon activity (%)”, col = c(“white”,“pink”,“red”))

topRegulators <- reshape2::melt(regulonActivity_byCellType_Binarized)
colnames(topRegulators) <- c("Regulon", "CellType", "RelativeActivity")
topRegulators <- topRegulators[which(topRegulators$RelativeActivity>minPerc),]
viewTable(topRegulators)

• 细胞类型特异性调节因子

  • 基于Suo等人在2018年为小鼠细胞图谱提出的Regulon特异性评分（RSS）。
  • 可用于多种细胞类型的大分析，以确定细胞类型的特定规则。

  rss <- calcRSS(AUC=getAUC(regulonAUC), cellAnnotation=cellInfo[colnames(regulonAUC), “CellType”])
  rssPlot <- plotRSS(rss)
  plotly::ggplotly(rssPlot$plot)

plotRSS_oneSet(rss, setName = "interneurons")

参考文章

单细胞转录调控网络分析工具（TF，motifs）-SCENIC

单细胞转录组高级分析二：转录调控网络分析

Running SCENIC

单细胞个性化分析之转录因子篇

SCENIC单细胞转录因子预测|3.软件安装与数据准备

SCENIC单细胞转录因子预测|4.精简版流程

SCENIC单细胞转录因子预测|5.step1+step2构建共表达网络与regulon

SCENIC单细胞转录因子预测|7.Step4 二元矩阵的计算与可视化