今天测试另外一个经典的人类pbmc的数据，是一个Seurat的对象。
#加载相关的包
library(SCENIC)
library(SCopeLoomR)
library(Seurat)
===========载入pbmc seurat对象=========
seurat.obj <- readRDS("pbmc.rds")
exprMat <- as.matrix(seurat.obj@assays$RNA@data)
cellInfo <- seurat.obj@meta.data[,c("cell_type","nCount_RNA","nFeature_RNA")]
colnames(cellInfo) <- c("CellType","nGene","nUMI")
===========初始化数据库设置================
data(list="motifAnnotations_hgnc_v9", package="RcisTarget")
motifAnnotations_hgnc <- motifAnnotations_hgnc_v9
scenicOptions <- initializeScenic(org="hgnc", dbDir="cisTarget_databases", nCores=10)
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)
#这一步太慢了，R的程序，在这一步跑了大约15个小时
========推断共表达模块========
### Build and score the GRN
exprMat_log <- log2(exprMat+1)
scenicOptions@settings$dbs <- scenicOptions@settings$dbs["10kb"] # Toy run settings
scenicOptions <- runSCENIC_1_coexNetwork2modules(scenicOptions)
scenicOptions <- runSCENIC_2_createRegulons(scenicOptions, coexMethod=c("top5perTarget")) # Toy run settings
scenicOptions <- runSCENIC_3_scoreCells(scenicOptions, exprMat_log)
scenicOptions <- runSCENIC_4_aucell_binarize(scenicOptions)
tsneAUC(scenicOptions, aucType="AUC") # choose settings
===========推断细胞特异的regulon/module===========
# Cell-type specific regulators (RSS):
regulonAUC <- loadInt(scenicOptions, "aucell_regulonAUC")
rss <- calcRSS(AUC=getAUC(regulonAUC), cellAnnotation=cellInfo[colnames(regulonAUC), "CellType"], )
rssPlot <- plotRSS(rss)
rssPlot$plot
======默认结果展示=======
下面是富集的motif的信息。
所有regulon在细胞的AUCscore热图：
二进制形式的热图。
