今天測試另外一個經(jīng)典的人類pbmc的數(shù)據(jù)泼舱，是一個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")

===========初始化數(shù)據(jù)庫設置================

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")

===========構建共表達網(wǎng)絡===========

### 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熱圖：

二進制形式的熱圖年缎。