【SCENIC】R版本例子测试1---pbmc

今天测试另外一个经典的人类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的信息。

富集motif

所有regulon在细胞的AUCscore热图：

二进制形式的热图。