In aertslab/SCENIC: SCENIC (Single Cell rEgulatory Network Inference and Clustering)

# Suppress loading messages when building the HTML
suppressPackageStartupMessages({
  library(SCENIC)
  library(AUCell)
  library(RcisTarget)
  library(SingleCellExperiment)
})

This tutorial provides the detailed explanation of runSCENIC_4_aucell_binarize(): Binarize the AUC and re-cluster.

All the code below is the content of the function runSCENIC_4_aucell_binarize(). This tutorial is meant for advanced users, who want know the details about what this function does internally, or to modify the workflow. There is no need to follow this tutorial for a regular run of SCENIC (see vignette("SCENIC_Running")).

Input

(See in the main vignette how to adjust the AUCell thresholds)

setwd("SCENIC_MouseBrain")
scenicOptions <- readRDS("int/scenicOptions.Rds")

skipBoxplot=FALSE # Whether to plot the boxplots
skipHeatmaps=FALSE # Whether to plot the Binary heatmaps
skipTsne=FALSE # Whether to caculate the Binary t-SNE

runSCENIC_4_aucell_binarize() code:

Main code: Assing cells and save as binary matrix

nCores <- getSettings(scenicOptions, "nCores")
regulonAUC <- loadInt(scenicOptions, "aucell_regulonAUC")
thresholds <- loadInt(scenicOptions, "aucell_thresholds")
thresholds <- getThresholdSelected(thresholds)

# Assign cells
regulonsCells <- setNames(lapply(names(thresholds), 
                                 function(x) {
                                   trh <- thresholds[x]
                                   names(which(getAUC(regulonAUC)[x,]>trh))
                                 }),names(thresholds))
### Convert to matrix (regulons with zero assigned cells are lost)
regulonActivity <- reshape2::melt(regulonsCells)
binaryRegulonActivity <- t(table(regulonActivity[,1], regulonActivity[,2]))
class(binaryRegulonActivity) <- "matrix"
saveRDS(binaryRegulonActivity, file=getIntName(scenicOptions, "aucell_binary_full"))

The binaryRegulonActivity contains some duplicated regulons (e.g. for some TFs, there is a regulon based on direct annotation, and also the extended version). We will also save a version that only keeps the "extended" regulons if there is not a regulon based on direct annotation:

# Alternative version: Keep only non-duplicated thresholds
binaryRegulonActivity_nonDupl <- binaryRegulonActivity[which(rownames(binaryRegulonActivity) %in% onlyNonDuplicatedExtended(rownames(binaryRegulonActivity))),]
saveRDS(binaryRegulonActivity_nonDupl, file=getIntName(scenicOptions, "aucell_binary_nonDupl"))

Info:

minCells <- ncol(binaryRegulonActivity) * .01
msg <- paste0("Binary regulon activity: ",
              nrow(binaryRegulonActivity_nonDupl), " TF regulons x ",
              ncol(binaryRegulonActivity), " cells.\n(",
              nrow(binaryRegulonActivity), " regulons including 'extended' versions)\n",
              sum(rowSums(binaryRegulonActivity_nonDupl)>minCells),
              " regulons are active in more than 1% (", minCells, ") cells.")
if(getSettings(scenicOptions, "verbose")) message(msg)

Plot boxplots (optional):

if(!skipBoxplot)
{
  .openDev(fileName=getOutName(scenicOptions, "s4_boxplotBinaryActivity"),
           devType=getSettings(scenicOptions, "devType"))
  par(mfrow=c(1,2))
  boxplot(rowSums(binaryRegulonActivity_nonDupl), main="nCells per regulon",
          sub='number of cells \nthat have the regulon active',
          col="darkolivegreen1", border="#001100", lwd=2, frame=FALSE)
  boxplot(colSums(binaryRegulonActivity_nonDupl), main="nRegulons per Cell",
          sub='number of regulons \nactive per cell',
          col="darkolivegreen1", border="#001100", lwd=2, frame=FALSE)
  dev.off()
}

Plot binary activity heatmap (optional)

The binary activity matrix can be visualized as heatmap.

Since there are usually regulons detected in very few cells that can hide more relevant regulons, we will plot several versions of the heatmap (e.g. one including all the regulons, and others with several subsets):

if(!skipHeatmaps)
{
  regulonSelection <- loadInt(scenicOptions, "aucell_regulonSelection", ifNotExists="null", verbose=FALSE)
  if(is.null(regulonSelection)) 
    regulonSelection <- regulonSelections(binaryRegulonActivity, binaryRegulonActivity_nonDupl, minCells)

  cellInfo <- loadFile(scenicOptions, getDatasetInfo(scenicOptions, "cellInfo"), ifNotExists="null")
  cellInfo <- data.frame(cellInfo)
  colVars <- loadFile(scenicOptions, getDatasetInfo(scenicOptions, "colVars"), ifNotExists="null")


  ### Plot heatmap:
  for(selRegs in names(regulonSelection$labels))
  {
    if(length(regulonSelection[[selRegs]])>1)
    {
      regulonSelection[[selRegs]] <- regulonSelection[[selRegs]][which(regulonSelection[[selRegs]] %in% rownames(binaryRegulonActivity))]
      binaryMat <- binaryRegulonActivity[regulonSelection[[selRegs]],,drop=FALSE]

      fileName <- paste0(getOutName(scenicOptions, "s4_binaryActivityHeatmap"),selRegs)

      fileName <- .openDevHeatmap(fileName=fileName, devType=getSettings(scenicOptions, "devType"))
        NMF::aheatmap(binaryMat, scale="none", revC=TRUE, main=selRegs,   
                      annCol=cellInfo[colnames(binaryMat),, drop=FALSE],
                      annColor=colVars,
                      color = c("white", "black"),
                      filename=fileName)
      if(getSettings(scenicOptions, "devType")!="pdf") dev.off()
    }
  }
}  

t-SNE on binary activity (optional)

if(!skipTsne)
{
  tSNE_fileName <- tsneAUC(scenicOptions, aucType="Binary", filePrefix=getIntName(scenicOptions, "tsne_prefix"), onlyHighConf=FALSE) # default: nPcs, perpl, seed
  tSNE <- readRDS(tSNE_fileName)

  # AUCell (activity) as html: 
  fileName <- getOutName(scenicOptions, "s4_binarytSNE_colAct")
  plotTsne_regulonActivityHTML(scenicOptions, exprMat, fileName, tSNE) #open the resulting html locally

  # Plot cell properties:
  sub <- ""; if("type" %in% names(tSNE)) sub <- paste0("t-SNE on ", tSNE$type)
  cellInfo <- loadFile(scenicOptions, getDatasetInfo(scenicOptions, "cellInfo"), ifNotExists="null")
  colVars <- loadFile(scenicOptions, getDatasetInfo(scenicOptions, "colVars"), ifNotExists="null")
  pdf(paste0(getOutName(scenicOptions, "s4_binarytSNE_colProps"),".pdf"))
  plotTsne_cellProps(tSNE$Y, cellInfo=cellInfo, colVars=colVars, cex=1, sub=sub)
  dev.off()
}

# Finished. Update status.
object@status$current <- 4

aertslab/SCENIC documentation built on April 7, 2024, 10 a.m.