Chapter 16 SCENIC

library(Seurat)
library(tidyverse)
library(magrittr)
library(SCENIC)
library(SingleCellExperiment)
library(SCopeLoomR)

16.1 Load Seurat Obj

combined <- readRDS('data/Demo_CombinedSeurat_SCT_Preprocess_FilterLQCells.rds')
Idents(combined) <- 'cluster'
DefaultAssay(combined) <- 'RNA'
combined <- NormalizeData(combined)

16.2 Filtering genes for calculating time

exprMat <- combined@assays$RNA@data
cellInfo <- combined@meta.data

loci1 <- which(rowSums(exprMat) > 1*.01*ncol(exprMat))
exprMat_filter <- exprMat[loci1, ]

16.3 Tranfer into loom files

add_cell_annotation <- function(loom, cellAnnotation)
{
  cellAnnotation <- data.frame(cellAnnotation)
  if(any(c("nGene", "nUMI") %in% colnames(cellAnnotation)))
  {
    warning("Columns 'nGene' and 'nUMI' will not be added as annotations to the loom file.")
    cellAnnotation <- cellAnnotation[,colnames(cellAnnotation) != "nGene", drop=FALSE]
    cellAnnotation <- cellAnnotation[,colnames(cellAnnotation) != "nUMI", drop=FALSE]
  }
  
  if(ncol(cellAnnotation)<=0) stop("The cell annotation contains no columns")
  if(!all(get_cell_ids(loom) %in% rownames(cellAnnotation))) stop("Cell IDs are missing in the annotation")
  
  cellAnnotation <- cellAnnotation[get_cell_ids(loom),,drop=FALSE]
  # Add annotation
  for(cn in colnames(cellAnnotation))
  {
    add_col_attr(loom=loom, key=cn, value=cellAnnotation[,cn])
  }
  
  invisible(loom)
}

loom <- build_loom("data/Demo_CombinedSeurat_SCT_Preprocess_FilterLQCells.loom", dgem=exprMat_filter)
loom <- add_cell_annotation(loom, cellInfo)
close_loom(loom)

16.4 Run SCENIC vis pySCENIC

  • GRN pyscenic grn {loom_file} hs_hgnc_tfs.txt -o adj.csv --num_workers 20

  • CLI pyscenic ctx adj.csv *feather --annotations_fname motifs-v9-nr.hgnc-m0.001-o0.0.tbl --expression_mtx_fname {loom_file} --output reg.csv --mask_dropouts --num_workers

  • AUCell pyscenic aucell {loom_file} reg.csv --output {out_loom_file} --num_workers 20

16.5 Retrieve AUC scores per cell

  • AUC scores
loom <- open_loom('data/Demo_CombinedSeurat_SCT_Preprocess_FilterLQCells.loom')
regulons_incidMat <- get_regulons(loom, column.attr.name="Regulons")
regulons <- regulonsToGeneLists(regulons_incidMat)
regulonAUC <- get_regulons_AUC(loom, column.attr.name='RegulonsAUC')

saveRDS(list('regulons' = regulons, 'regulonAUC' = regulonAUC), 
        file = 'data/Demo_CombinedSeurat_SCT_Preprocess_FilterLQCells-pyscenic-output-AUC.rds')
  • Binary AUC scores
from binarization import *
import pandas as pd
auc_mtx = pd.read_csv('test.auc.mat', sep = '\t')

auc_mtx_binary = binarize(auc_mtx)
test = pd.DataFrame(auc_mtx_binary[0])