R-SCENIC pipline

SCENIC 功能说明：基于物种转录因子数据库，用以识别 单细胞转录组 数据中的转录因子(TFs)和对应的靶标基因(Target gene)，建立共表达网络。

SCENIC 软件分析的一般用途：
1、识别TFs与潜在靶基因 的共表达模块，鉴定regulon，计算每个细胞的regulon活性评分（RAS）
2、regulon-specific 得分（RSS），用以判断regulon与细胞类型的富集关系；
3、regulon关联特异性指数（CSI），判断不同regulon的关联性，同时具有较高CSI的regulon，也即相关性较强的regulon可能共同调节下游基因，负责细胞功能。

0、环境配置

BiocManager::install(c("AUCell", "RcisTarget"))
BiocManager::install(c("GENIE3"))
BiocManager::install(c("zoo", "mixtools", "rbokeh"))
BiocManager::install(c("DT", "NMF", "ComplexHeatmap", "R2HTML", "Rtsne"))
BiocManager::install(c("doMC", "doRNG"))
remotes::install_github("aertslab/SCopeLoomR", build_vignettes = TRUE)
remotes::install_github("aertslab/SCENIC") 
remotes::install_github("FloWuenne/scFunctions")

cis_database链接 : https://pan.baidu.com/s/1ivtABzpxkWyesvs3Omswpw 提取码: fq8y

1、SCENIC R语言流程脚本

library(Seurat);

### load scRNA data
sc <- readRDS("seurat.object.path")

### Featch data for SCENIC
exprMat <- as.matrix(sc@assays$SCT@data)
cellInfo <- sc@meta.data

### Initialize SCENIC settings
library(SCENIC)
dir.create("int")
saveRDS(cellInfo, file="int/cellInfo.Rds")

scenicOptions <- initializeScenic(org="hgnc", dbDir="cisTarget_databases", nCores=20)
scenicOptions@inputDatasetInfo$cellInfo <- "int/cellInfo.Rds"
saveRDS(scenicOptions, file="int/scenicOptions.Rds") 

### Co-expression network
exprMat_filtered <- exprMat ### has done QC in snRNA
runCorrelation(exprMat_filtered, scenicOptions)
runGenie3(exprMat_filtered, scenicOptions)

### Build and score the GRN
scenicOptions@settings$verbose <- TRUE
scenicOptions@settings$seed <- 123
scenicOptions@settings$dbs <- scenicOptions@settings$dbs["10kb"]

exprMat_log <- exprMat_filtered
scenicOptions <- runSCENIC_1_coexNetwork2modules(scenicOptions)
scenicOptions <- runSCENIC_2_createRegulons(scenicOptions, coexMethod=c("top10perTarget")) 
library(doParallel)
scenicOptions <- runSCENIC_3_scoreCells(scenicOptions, exprMat_log)
scenicOptions <- runSCENIC_4_aucell_binarize(scenicOptions)
tsneAUC(scenicOptions, aucType="AUC") # choose settings

### save main results ...
saveRDS(scenicOptions, file="int/scenicOptions.Rds")

2、SCENIC regulons 结果导出

###  Result exploring
library(AUCell);library(SCENIC);library(BiocParallel)

cellInfo <- readRDS("./int/cellInfo.Rds")
scenicOptions <- readRDS("int/scenicOptions.Rds")
Class <- "celltype"  ### cellInfo 表中记录细胞类型的列名，后续RSS，Regulators相关计算用到

### load regulon data，返回一个不计算置信度低的结果的相关性，其中 tf_extend 结果是在数据库中置信度较低的 target-gene 注释结果
regulons <- loadInt(scenicOptions, "aucell_regulons")  ### 提取 regulons list记录，记录了每个regulon和对应的target-gene
dat <- tibble::enframe(regulons) %>% data.frame() %>% rename("regulon"=1,"targets"=2)  ### 将list转换成dataframe 导出
dat$targets <- vapply(dat$targets, paste, collapse = ", ", character(1L))
write.table(dat,"regulon-targets.list.tsv",quote = F,row.names = F,sep = '\t',col.names = T) 

### AUCell score 结果
regulonAUC <- loadInt(scenicOptions, "aucell_regulonAUC") ### regulon 和在每个细胞的 AUC 评分
regulonAUC <- regulonAUC[onlyNonDuplicatedExtended(rownames(regulonAUC)),]  ### TFs * Cell value== AUC data.frame ，这一部过滤了低置信度的regulon结果

### calculate Average Regulon Activity for cell type
regulonActivity_byCellType <- sapply(split(rownames(cellInfo), cellInfo[[Class]]),function(cells) rowMeans(getAUC(regulonAUC)[,cells]))
regulonActivity_byCellType_Scaled <- t(scale(t(regulonActivity_byCellType), center = T, scale=T)) ### 标准化细胞类型的AUC评分
topRegulators <- reshape2::melt(regulonActivity_byCellType_Scaled)
colnames(topRegulators) <- c("Regulon", "CellType", "RelativeActivity")
write.table(topRegulators,"regulon_celltype_AUC.tsv",col.names = T,row.names = F,quote = F,sep = "\t")

### 计算每个细胞类型正相关regulon的 AUC评分
topRegulators <- topRegulators[which(topRegulators$RelativeActivity>0),]

### Cell-type specific regulators  RSS
regulonAUC <- loadInt(scenicOptions, "aucell_regulonAUC")
rss <- calcRSS(AUC=getAUC(regulonAUC), cellAnnotation=cellInfo[colnames(regulonAUC), Class]) ###  TFs * Cell value== RSS data.frame
write.table(rss,"regulon_RSS.tsv",col.names = T,row.names = T,quote = F,sep = "\t")

### Calculate connection specificity index (CSI) for all regulons
library(scFunctions)
csi <- calculate_csi(regulonAUC, calc_extended = FALSE, verbose = FALSE) ### calc_extended 不计算tf_extend的结果，数据框
write.table(csi,"regulon_CSI.tsv",col.names = T,row.names = F,quote = F,sep = "\t")