2019-08-28第七天：关于SCENIC的使用教程
文章来源：scenic: single-cell regulatory network inference and clustering
一、软件安装
安装教程
1、安装依赖包和本包
2、下载Rcistarget使用的物种特异性的数据集
例如人TSS上游500bp和人TSS周围20kb及鼠TSS上游500bp和鼠TSS周围20kb
默认的，SCENIC使用的数据库是对基因启动子（TSS上游500bp和TSS周围20kb）上的motif打分的。
同时，该包目前只支持人、小鼠和果蝇三个模式物种。
二、SCENIC的运行
1、SCENIC的运行流程
建立基因调控网络（GRN）需要以下步骤
- 以共表达为基础确定每个转录因子的潜在靶点
选出基因表达矩阵并运行GENIE3/GRNBOOST
从GENIE3/GRNBoost中将靶点格式化为共表达模块
- 根据DNA-motif分析选择潜在的直接相关的靶点（RcisTarget：TF motif analysis）
确定细胞状态和它们的调控子
- 分析每个细胞的网络活性（AUCell）
对每个细胞中的调控子进行打分（计算AUC）
选做：将网络活性转化为二位矩阵（设定一个阈值）
- 根据细胞的基因调控网路确定稳定的细胞状态（细胞聚类）并且对结果进行解读
2、命令列表
以下是SCENIC流程的主要的命令行
library(foreach)#注意一定要加载这个包，不然会报错！
library(SCENIC)
### Load data
loomPath <- system.file(package="SCENIC", "examples/mouseBrain_toy.loom")
library(SCopeLoomR)
loom <- open_loom(loomPath, mode="r")
exprMat <- get_dgem(loom)
cellInfo <- get_cellAnnotation(loom)
close_loom(loom)
### Initialize settings
library(SCENIC)
scenicOptions <- initializeScenic(org="mgi", dbDir="cisTarget_databases", nCores=10)
# scenicOptions@inputDatasetInfo$cellInfo <- "int/cellInfo.Rds"
saveRDS(scenicOptions, file="int/scenicOptions.Rds")
### Co-expression network
genesKept <- geneFiltering(exprMat, scenicOptions)
exprMat_filtered <- exprMat[genesKept, ]
runCorrelation(exprMat_filtered, scenicOptions)
exprMat_filtered_log <- log2(exprMat_filtered+1)
runGenie3(exprMat_filtered_log, scenicOptions)
### Build and score the GRN
exprMat_log <- log2(exprMat+1)
scenicOptions@settings$dbs <- scenicOptions@settings$dbs["10kb"] # Toy run settings
runSCENIC_1_coexNetwork2modules(scenicOptions)
runSCENIC_2_createRegulons(scenicOptions, coexMethod=c("top5perTarget")) # Toy run settings
runSCENIC_3_scoreCells(scenicOptions, exprMat_log)
# Export:
export2scope(scenicOptions, exprMat)
# Binarize activity?
# aucellApp <- plotTsne_AUCellApp(scenicOptions, exprMat_log)
# savedSelections <- shiny::runApp(aucellApp)
# 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")
runSCENIC_4_aucell_binarize(scenicOptions)
### Exploring output
# Check files in folder 'output'
# .loom file @ http://scope.aertslab.org
# output/Step2_MotifEnrichment_preview.html in detail/subset:
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)
3、目录
在这个工作流程中我们将要保存多个文件，为了保证他们的整洁性，我们需要建立一个新的工作目录
比如说：
dir.create("SCENIC_MouseBrain")
setwd("SCENIC_MouseBrain") # Or `knitr::opts_knit$set(root.dir = 'example_results/SCENIC_MouseBrain')` in the first chunk if running a notebook
默认的，中间文件和图文件将会保存到int文件夹中，并且会按顺序整齐排列好。
你可能使用这些文件来检查之前做过的每个步骤，或者是以一些不同的设置重新运行其中的分析。重要的输出结果会保存到output文件夹中。
1）数据输入
- 表达矩阵
SCENIC的输入是一个单细胞的表达矩阵，其中基因缩写名是行名，第一步是去加载这个矩阵。
在这个教程中，我们提供了一个只有200个细胞以及少于1000个基因的小鼠脑数据集
loomPath <- system.file(package="SCENIC", "examples/mouseBrain_toy.loom")
打开loom文件并加载这个表达矩阵（如果有的话也可以加载细胞注释文件）
library(SCopeLoomR)
loom <- open_loom(loomPath, mode="r")
exprMat <- get_dgem(loom)
cellInfo <- get_cellAnnotation(loom)
close_loom(loom)
dim(exprMat)
- 细胞信息或者表行数据
在3-4步（对GRN进行打分并聚类）中，我们或许会对分析的结果与这些细胞已知的信息之间的差异感兴趣。你可以指定使用那些变量作图，同时给他们设定不同的颜色。
# cellInfo$nGene <- colSums(exprMat>0)
head(cellInfo)
cellInfo <- data.frame(cellInfo)
cellTypeColumn <- "Class"
colnames(cellInfo)[which(colnames(cellInfo)==cellTypeColumn)] <- "CellType"
cbind(table(cellInfo$CellType))
dir.create("int")
saveRDS(cellInfo, file="int/cellInfo.Rds")
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))
2）初始化SCENIC设置
为了保持SCENIC的多步中的设置参与的一致性，SCENIC包中的大部分的函数使用一个共同的对象，这个对象中储存了当前环境下运行参数的信息。这个对象会替换很多函数的参数，并且应该在SCENIC的开始时就使用initializescenic（）这个函数来创建它。
默认的设置对大多数分析过程都有效。在所有分析中比较特殊的参数是物种这个选项，其中mgi表示小鼠，hgnc表示人，dmel表示果蝇。同时其中的dir应该是RcisTarget数据库保存的位置。
这个设置的详细信息可以参考帮助文档
library(SCENIC)
org="mgi" # or hgnc, or dmel
dbDir="cisTarget_databases" # RcisTarget databases location
myDatasetTitle="SCENIC example on Mouse brain" # choose a name for your analysis
data(defaultDbNames)
dbs <- defaultDbNames[[org]]
scenicOptions <- initializeScenic(org=org, dbDir=dbDir, dbs=dbs, datasetTitle=myDatasetTitle, nCores=10)
## Motif databases selected:
## mm9-500bp-upstream-7species.mc9nr.feather
## mm9-tss-centered-10kb-7species.mc9nr.feather
# Modify if needed
scenicOptions@inputDatasetInfo$cellInfo <- "int/cellInfo.Rds"
scenicOptions@inputDatasetInfo$colVars <- "int/colVars.Rds"
# Databases:
# scenicOptions@settings$dbs <- c("mm9-5kb-mc8nr"="mm9-tss-centered-5kb-10species.mc8nr.feather")
# scenicOptions@settings$db_mcVersion <- "v8"
# Save to use at a later time...
saveRDS(scenicOptions, file="int/scenicOptions.Rds")
3）共表达网路建立
SCENIC流程的第一步是根据表达数据推测潜在的转录因子靶点。因为我们使用了GENIE3或者GRNBoost这两种方法。这两个工具的输入文件都是经过筛选的表达矩阵，以及一个转录因子列表。其输出文件及相关性矩阵将会被用来建立共表达模块（使用函数runscenic——1——coexnetwork2modules（））
选择这两个工具的依据为：原则上，推测共表达网络的目前大多方法都可以使用这个过程的输出文件，只要他的输出文件是一个包含每个转录因子的潜在靶点的列表即可。我们选择GENIE3是因为它可以确定非线性的关系，即使他们只存在于一部分样品中。同时它在Network inference DREAM5挑战中具有很好的表现。同时这个方法可以很好的运行在R中。
然而，GENIE3是非常耗时的以及对计算力有严重的消耗，比如在处理3k到5k个细胞的时候可能会花费几个小时。因此为了处理大数据集，我们开发了GRNBoost算法，这个算法提供与GENIE3相似的结果，但是时间上大大节省了，因此我们强烈推荐大家在处理大数据集的时候使用它。
二次取样细胞：当在你的数据集中存在大量的低质量的细胞时，也会消耗大量的时间，因此我们可以使用这种方法去挑选其中一部分的细胞来推断调控网络，比如随机挑选或者使用高质量的细胞的表达矩阵作为输入对象来进行共表达分析。其中调控网络的活性会根据这个小细胞集进行训练，但是在之后的第三步中可以使用全部的细胞在AUCell中重新评估。值得注意的是，为了避免丢失分辨率，这个小细胞集应该尽可能代表所有细胞，比如可以包含所有的细胞种类。
基因筛选和选择
为了运行GENIE3/GRNBoost，我们推荐使用软基因筛选方式，来去除在细胞中低表达或者仅仅在很少的细胞中表达的基因。在这里，我们应用了基于基因表达总数和其表达的细胞数目的双筛选标准进行筛选基因。
-
通过基因表达总数进行筛选。这个筛选是为了去除那些看起来像噪音的基因。默认这个过程只保留在所有细胞中至少有6个umi检测到的基因（这里使用了200细胞），比如基因的总数应该为31%细胞数目。可以使用mincountspergene这个参数来调节这个值。
-
通过该基因表达的细胞数目筛选，默认的，minsample参数设置为至少在1%的细胞中表达的基因会被保留。这个筛选的目的是去除那些所有的read都在同一个或极少细胞中检测到的基因，保证了留下来的基因都是在很多细胞中有较高的表达量。为了避免去除一些小的细胞群体，我们推荐设置的值一定要小于可以检测到的最小的cluster中的细胞数目。
-
最后，只有在RcisTatget数据库中存在的基因才会被保留下来，这个过程的是为了节省在GENIE3运行的时间，因为在数据库中不存在的基因在之后的步骤中没有用。
# (Adjust minimum values according to your dataset)
genesKept <- geneFiltering(exprMat, scenicOptions=scenicOptions,
minCountsPerGene=3*.01*ncol(exprMat),
minSamples=ncol(exprMat)*.01)
在进行network inference之前，检查一下是否一些已知的相关基因被筛掉了，通过这个方式可以检查一下这个筛选过程是否合理。
interestingGenes <- c("Sox9", "Sox10", "Dlx5")
# any missing?
interestingGenes[which(!interestingGenes %in% genesKept)]
我们现在可以使用筛选得到的基因集与表达矩阵做一个交集，最终得到的表达矩阵可以作为共表达过程的的输入文件
exprMat_filtered <- exprMat[genesKept, ]
dim(exprMat_filtered)
rm(exprMat)
相关性分析
GENIE3/GRNBoost可以检测到正相关和负相关。为了区分潜在的转录激活和转录抑制，我们将得到的靶点分为正靶点和负靶点，在这里使用的相关性分析的方法为spearman法。
runCorrelation(exprMat_filtered, scenicOptions)
GENIE3
如果是运行GRNBoost（python），那需要查看？exportsforgrnboost
GENIE3的输入是典型的表达矩阵和待选的转录因子的列表。函数rungenie3会按照默认的设置进行GENIE3，这些设置对大多数数据集都足够用了，同时使用了RcisTarget数据库中的转录因子作为候选的转录因子。
因为GENIE3是根据随机森林的方法进行的，每次运行的时候结果都稍微有点不一样。树的数目越多，变异性越低。我们推荐使用set.seed来保证多次运行的重复性。
GENIE3会消耗几个小时去运行，如果你在Rsturio中运行这个过程，我们推荐你停止并且在新的独立的R控制台中运行这个代码块，或者在服务器中运行。随后的代码块将会通过加载GENIE3的输出结果来恢复这个流程。
## If launched in a new session, you will need to reload...
# setwd("...")
# loomPath <- "..."
# loom <- open_loom(loomPath, mode="r")
# exprMat <- get_dgem(loom)
# close_loom(loom)
# genesKept <- loadInt(scenicOptions, "genesKept")
# exprMat_filtered <- exprMat[genesKept,]
# library(SCENIC)
# scenicOptions <- readRDS("int/scenicOptions.Rds")
# Optional: add log (if it is not logged/normalized already)
exprMat_filtered <- log2(exprMat_filtered+1)
# Run GENIE3
runGenie3(exprMat_filtered, scenicOptions)
4）建立和计分GRN
在完成了上述过程之后，剩下的流程可以使用预先包装的函数进行。其中每个函数都与SCENIC流程一一对应
建立共表达网络：1⃣️ 得到共表达模块 2⃣️得到调控子（使用RcisTarget）：TFmotif分析
确定细胞状态：3⃣️ 使用AUCell进行打分 4⃣️根据GRN活性来聚类细胞
如果需要的话，重新加载表达数据
loom <- open_loom(loomPath, mode="r")
exprMat <- get_dgem(loom)
close_loom(loom)
# Optional: log expression (for TF expression plot, it does not affect any other calculation)
logMat <- log2(exprMat+1)
dim(exprMat)
使用包装函数运行之后的步骤
library(SCENIC)
scenicOptions <- readRDS("int/scenicOptions.Rds")
scenicOptions@settings$verbose <- TRUE
scenicOptions@settings$nCores <- 10
scenicOptions@settings$seed <- 123
# For a very quick run:
# coexMethod=c("top5perTarget")
scenicOptions@settings$dbs <- scenicOptions@settings$dbs["10kb"] # For toy run
# save...
runSCENIC_1_coexNetwork2modules(scenicOptions)
runSCENIC_2_createRegulons(scenicOptions, coexMethod=c("top5perTarget")) #** Only for toy run!!
runSCENIC_3_scoreCells(scenicOptions, logMat)
#这一步可能会报
#Error in foreach(param = allParams) %dopar% { :
#could not find function "%dopar%"
#使用library(foreach)这个包解决，下面runSCENIC_4也是一样！
5）可选步骤
对网络活性的结果进行二维化（on/off）
建立GRN和使用AUCell进行打分对于非常清楚的细胞类型数据集来说是足够了，然而，在很多情况下，将活性的结果转化为二维矩阵是有用的，比如可以更容易理解或者最大化细胞类型之间的差异。我们一般只对某些特殊的调控子进行二维化处理，但是，当我们对所有的调控子处理时可以获得一个完整的二维调控子活性矩阵，这个可以用来做随后的聚类分析，这种方法也可以用来降低技术误差，或者用于物种之间的比较等。
为了确定每个调控子在细胞中的激活情况，我们使用了一个AUC的阈值。AUCell可以自动计算一个可能阈值，但是一般都过于保守。我们推荐用户手动确定这些阈值。这会是一个迭代的过程，这些阈值会在运行过程中重新调整。一旦最后的阈值被选择，细胞调控子的活性将会被总结为一个二维矩阵，其中行代表调控子，列代表细胞。矩阵中的坐标表示在细胞中激活的使用1表示，在其他细胞中用0表示。
我们可以在以前的输出的文件中看到最后设置的阈值。
aucellApp <- plotTsne_AUCellApp(scenicOptions, logMat)
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")
一旦你确定了阈值，那就可以运行函数runscenic_4_aucell_binarize来对AUC进行二维化，并且产生一些额外的数值和类别。
# scenicOptions@settings$devType="png"
runSCENIC_4_aucell_binarize(scenicOptions)
tsne图也可以使用二维矩阵去创作，只需要运行tsneauc函数即可。
在调节子活性的矩阵上对细胞降维并聚类
细胞可以根据调节子活性进行聚类，这个调节子活性矩阵可以是二维化后的也可以是连续的矩阵。
如果使用tsne做可视化，我们推荐尝试不同的设置来评估簇的稳定性，也可以是用UMAP及其他聚类方法或者是轨迹推断方法
SCENIC中的函数可以以不同的设置运行多重tsne。它会在选择的number of PCs和perpiexity中创建所有的组合，消耗的时间取决于细胞的数目
nPcs <- c(5) # For toy dataset
# nPcs <- c(5,15,50)
scenicOptions@settings$seed <- 123 # same seed for all of them
# Run t-SNE with different settings:
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")
# Plot as pdf (individual files in int/):
fileNames <- paste0("int/",grep(".Rds", grep("tSNE_", list.files("int"), value=T), value=T))
其中被选择的tsne则会做为默认的选项用于之后的作图
scenicOptions@settings$defaultTsne$aucType <- "AUC"
scenicOptions@settings$defaultTsne$dims <- 5
scenicOptions@settings$defaultTsne$perpl <- 15
saveRDS(scenicOptions, file="int/scenicOptions.Rds")
以loom或者scope的格式导出文件
将结果加载到loom文件中
四、对于结果的理解和挖掘
在output文件夹中包含了一些文件，这是每一步结果的输出文件的概览。这些结果可以通过多个中间文件（储存在int文件夹中）进行更细致的探索。
一些探索的例子如下所述
1）细胞状态
AUCell提供每个细胞的调节子活性状态。通过这些数据对细胞进行聚类（可以是连续的矩阵或者二维化的矩阵），我们可以看到是否存在一些类群的细胞有相同的调控子活性以及揭示了多细胞间存在的网络状态。这些状态等同于网路的吸引子的状态（attractor）。结合不同的可视化方法，我们可以探索一些特殊的调控子与细胞状态之间的关联。
SCENIC提供了一些包装的函数去快速地运行。比如说，将AUC的结果和TF的表达投射到tsne上，并且将AUC以热图的形式可视化，但是用户也可以根据自己的需求去探索别的聚类和可视化的工具
- 将AUC和TF的表达投射到tsne上
简单来说，tsne是一个2D的细胞投影，在这个图中如果细胞基因表达模式相近的话细胞就会彼此靠近。tsne通常可以让我们快速的在一个数据集中看到细胞状态的概览。值得注意的是，tsne可以很好的对细胞类群进行区分，但是它并不适合去探索一些连续的或者多变的过程，比如说轨迹的可视化
logMat <- exprMat # Better if it is logged/normalized
aucellApp <- plotTsne_AUCellApp(scenicOptions, logMat) # default t-SNE
savedSelections <- shiny::runApp(aucellApp)
#AUCell_plotTSNE() to save static plots:
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")
# Save AUC as PDF:
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()
密度图可以检测最稳定的细胞状态
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(bg = "black")
par(mfrow=c(1,2))
regulonNames <- c( "Dlx5","Sox10")
cellCol <- plotTsne_rgb(scenicOptions, regulonNames, aucType="AUC", aucMaxContrast=0.6)
text(0, 10, attr(cellCol,"red"), col="red", cex=.7, pos=4)
text(-20,-10, attr(cellCol,"green"), col="green3", cex=.7, pos=4)
regulonNames <- list(red=c("Sox10", "Sox8"),
green=c("Irf1"),
blue=c( "Tef"))
cellCol <- plotTsne_rgb(scenicOptions, regulonNames, aucType="Binary")
text(5, 15, attr(cellCol,"red"), col="red", cex=.7, pos=4)
text(5, 15-4, attr(cellCol,"green"), col="green3", cex=.7, pos=4)
text(5, 15-8, attr(cellCol,"blue"), col="blue", cex=.7, pos=4)
- GRN:调控子靶点和motif
展示调控子中的基因
regulons <- loadInt(scenicOptions, "regulons")
regulons[c("Dlx5", "Irf1")]
注意只有在具有十个及以上的调控子才会被AUCell打分
regulons <- loadInt(scenicOptions, "aucell_regulons")
head(cbind(onlyNonDuplicatedExtended(names(regulons))))
关于转录因子-靶点对的详细信息：对于每个转录因子靶点对，中间步骤的统计统计数据结果总结在loadInt(scenicOptions, "regulonTargetsInfo") (saved as text in: getOutName(scenicOptions, "s2_regulonTargetsInfo"): output/Step2_regulonTargetsInfo.tsv).这个表格可以被用来探索一些特殊的连接。因为它包含了几千行，所以在对特殊连接进行探索时需要先筛选。
regulonTargetsInfo <- loadInt(scenicOptions, "regulonTargetsInfo")
tableSubset <- regulonTargetsInfo[TF=="Stat6" & highConfAnnot==TRUE]
viewMotifs(tableSubset)
转录因子的motif 的全部列表可以从结果文件RcisTarget motif enrichment results中查找。这些保存在motifEnrichment_selfMotifs_wGenes中。这些结果的可视化可以在html文件中查看output/Step2_MotifEnrichment_preview.html
- regulators for clusters or known cell types
使用SCENIC做的调控分析可以与其他分析结合使用，比如说聚类分析或者是关注特殊细胞类型的调控子等。有很多选项可以做这些分析，这取决于你的想象力。以下是一些可以用来考虑的分析过程
对于每个cluster的平均调节子活性分析
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))
pheatmap::pheatmap(regulonActivity_byCellType_Scaled, #fontsize_row=3,
color=colorRampPalette(c("blue","white","red"))(100), breaks=seq(-3, 3, length.out = 100),
treeheight_row=10, treeheight_col=10, border_color=NA)
# filename="regulonActivity_byCellType.pdf", width=10, height=20)
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),]
pheatmap::pheatmap(binaryActPerc_subset, # fontsize_row=5,
color = colorRampPalette(c("white","pink","red"))(100), breaks=seq(0, 1, length.out = 100),
treeheight_row=10, treeheight_col=10, border_color=NA)
#filename="regulonActivityBinary_byCellType.pdf", width=10, height=20)
topRegulators <- reshape2::melt(regulonActivity_byCellType_Binarized)
colnames(topRegulators) <- c("Regulon", "CellType", "RelativeActivity")
topRegulators <- topRegulators[which(topRegulators$RelativeActivity>minPerc),]
viewTable(topRegulators)
library(Seurat)
dr_coords <- Embeddings(seuratObject, reduction="tsne")
tfs <- c("Sox10","Irf1","Sox9", "Dlx5")
par(mfrow=c(2,2))
AUCell::AUCell_plotTSNE(dr_coords, cellsAUC=selectRegulons(regulonAUC, tfs), plots = "AUC")
五、对于输出结果的调用
可以使用getIntName(scenicOptions)查看中间文件，使用loadInt(scenicOptions,"文件")加载到环境中，
如加载regulons信息：regulons<-loadInt(scenicOptions,"regulons"),就可以把int/2.6_regulons_asGeneSet.Rds 这个文件加载到环境中
以下是int文件夹中的文件列表
fileName
genesKept "int/1.1_genesKept.Rds" #通过了筛选的基因列表
corrMat "int/1.2_corrMat.Rds" #相关性分析的结果矩阵
genie3wm "int/1.3_GENIE3_weightMatrix.Rds"
genie3ll "int/1.4_GENIE3_linkList.Rds"
genie3weighPlot "int/1.5_weightPlot"
tfModules_asDF "int/1.6_tfModules_asDF.Rds"
tfModules_forEnrichment "int/2.1_tfModules_forMotifEnrichmet.Rds"
motifs_AUC "int/2.2_motifs_AUC.Rds"
motifEnrichment_full "int/2.3_motifEnrichment.Rds"
motifEnrichment_selfMotifs_wGenes "int/2.4_motifEnrichment_selfMotifs_wGenes.Rds"
regulonTargetsInfo "int/2.5_regulonTargetsInfo.Rds" #所有转录因子及其靶基因的信息，包含motif信息
regulons "int/2.6_regulons_asGeneSet.Rds" #regulon及其调控的基因
regulons_incidMat "int/2.6_regulons_asIncidMat.Rds"
aucell_regulons "int/3.1_regulons_forAUCell.Rds" #与regulons中的信息一样
aucell_genesStatsPlot "int/3.2_aucellGenesStats"
aucell_rankings "int/3.3_aucellRankings.Rds"
aucell_regulonAUC "int/3.4_regulonAUC.Rds" #储存每个细胞的每个regulons的AUC数值，可以使用getAUC(loadInt(scenicOptions,"aucell_regulonAUC")) 提取
aucell_thresholds "int/3.5_AUCellThresholds.Rds" #储存用于进行二进制化的阈值，以及通过阈值的细胞名称
aucell_thresholdsTxt "int/3.5_AUCellThresholds_Info.tsv"
aucell_binary_full "int/4.1_binaryRegulonActivity.Rds" #储存 二进制化的矩阵
aucell_binary_nonDupl "int/4.2_binaryRegulonActivity_nonDupl.Rds"
aucell_regulonSelection "int/4.3_regulonSelections.Rds"
aucell_binaryRegulonOrder "int/4.4_binaryRegulonOrder.Rds"
