Chapter 17 SCENIC Differential Regulons

library(Seurat)
library(tidyverse)
library(magrittr)
library(viridis)
library(SCENIC)
library(network)
library(igraph)

17.1 Load Seurat Obj

combined <- readRDS('data/SCENIC-Demo-Seurat.rds')
combined$celltype.stim <- paste(combined$Try_Cluster2, combined$orig.ident, sep = "_")
combined$celltype.stim <- factor(combined$celltype.stim, levels = c('SC_Hyper', 'SC_AS', 'CyclingTA_Hyper', 'CyclingTA_AS',
                                                                    'TA_Hyper', 'TA_AS', 'EC_Hyper', 'EC_AS', 'GC_Hyper', 'GC_AS', 'EEC_Hyper'))
Idents(combined) <- "celltype.stim"

17.2 Load AUC socre and binary mat

## auc socres
scenic <- readRDS('data/SCENIC-Demo-SCENIC-AUC.rds')
regulons <- scenic$regulons
regulonAUC <-  scenic$regulonAUC
AUCmat <- AUCell::getAUC(regulonAUC)
rownames(AUCmat) <- gsub("[(+)]", "", rownames(AUCmat))

## binary auc scores
Binarymat <- read.csv('data/SCENIC-Demo-SCENIC-AUC-Binary.mat',  sep = '\t')
rownames(Binarymat) <- Binarymat$cells
Binarymat$cells <- NULL
Binarymat <- t(Binarymat)
rownames(Binarymat) <- gsub("[...]", "", rownames(Binarymat))

combined[['AUC']] <- CreateAssayObject(data = AUCmat)
combined[['AUCBinary']] <- CreateAssayObject(data = Binarymat)

17.3 Identify differential regulons based AUC scores

DefaultAssay(combined) <- 'AUC'

deg.ls <- unique(combined$Try_Cluster2)[1:5] %>% map(~{
  id1 <- paste0(.x, '_AS')
  id2 <- paste0(.x, '_Hyper')
  
  deg <- FindMarkers(combined, ident.1 = id2, ident.2 = id1, logfc.threshold = 0.005)
  deg.up <- deg[which(deg$avg_logFC > 0), ]
  deg.dn <- deg[which(deg$avg_logFC < 0), ]
  deg.ls <- list(deg.up, deg.dn)
  names(deg.ls) <- c(id2, id1)
  return(deg.ls)
  
})

names(deg.ls) <- unique(combined$Try_Cluster2)[1:5]
deg.ls <- list(deg.ls$SC$SC_Hyper, deg.ls$SC$SC_AS, 
               deg.ls$CyclingTA$CyclingTA_Hyper, deg.ls$CyclingTA$CyclingTA_AS, 
               deg.ls$TA$TA_Hyper, deg.ls$TA$TA_AS)
names(deg.ls) <- c('SC-Hyper', 'SC-AS', 'CyclingTA-Hyper', 'CyclingTA-AS', 'TA-Hyper', 'TA-AS')

gene.ls <- lapply(deg.ls, rownames)

17.4 Visu check differential regulons

• AUC scores

DefaultAssay(combined) <- 'AUC'
combined <- ScaleData(combined, assay = 'AUC', features = rownames(AUCmat))

## Centering and scaling data matrix

mg <- unique(Reduce('c', gene.ls[c(2)]))

DoHeatmap(combined, features = mg, slot = 'scale.data', raster = F) + scale_fill_viridis()

## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.

DoHeatmap(combined, features = mg, slot = 'scale.data', raster = F) + scale_fill_gradient2(
  low = rev(c('#d1e5f0', '#67a9cf', '#2166ac')),
  mid = "white",
  high = rev(c('#b2182b', '#ef8a62', '#fddbc7')),
  midpoint = 0,
  guide = "colourbar",
  aesthetics = "fill"
)

## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.

• Binary AUC scores

DefaultAssay(combined) <- 'AUCBinary'
mg <- unique(Reduce('c', gene.ls[c(1,3,5)]))

DoHeatmap(combined, features = mg, slot = 'data') + scale_fill_gradientn(colors = c( "white", "black"))

## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.

17.5 Visu check regulons via network

visuNetwork <- function(regulon.name){
  
  adj.ls <- regulon.name %>% map(~{
    tmp <- adj[which(adj$TF == .x), ]
    tmp <- tmp[order(tmp$importance, decreasing = T),]
    loci <- unique(c(1:50, grep('MUC2', tmp$target)))
    tmp <- tmp[loci,]
    return(tmp)
  })
  
  adj.sub <- Reduce('rbind', adj.ls)
  
  ## generate network
  edge.df <- adj.sub
  colnames(edge.df) <- c('from', 'to', 'weights')
  edge.df$from <- as.character(edge.df$from)
  edge.df$to <- as.character(edge.df$to)
  
  #net1 <- network(edge.df, directed = T, loops = T)
  #plot(net1)
  
  vertex <- unique(c(edge.df$from, edge.df$to))
  net <- graph_from_data_frame(d = edge.df, vertices = vertex, directed = T)
  
  
  ## vertex size scaled to log2FC
  #combined.tmp <- AverageExpression(combined, assays = 'RNA', features = vertex)
  #fc <- combined.tmp$RNA$SC_Hyper/combined.tmp$RNA$SC_AS
  #fc <- combined.tmp$RNA$SC_AS/combined.tmp$RNA$SC_Hyper
  #fc[which(fc > quantile(fc, 0.9))] <- quantile(fc, 0.9)
  #fc[which(fc < quantile(fc, 0.1))] <- quantile(fc, 0.1)
  #vsize <- scale(fc, center = min(fc), scale = max(fc) - min(fc))
  #vsize <- vsize+0.1
  
  ## vertex color: whether DEG in hyper vs as
 # vcl <- ifelse(vertex %in% dn.deg, 'blue', ifelse(vertex %in% up.deg, 'red', 'grey'))
#  vlabel <- rep("", length(vertex))
  # vertex label 
 # loci <- which(vcl != 'grey' | vertex %in% edge.df$from | vertex == 'MUC2')
#  vlabel[loci] = vertex[loci] 
  # vertex size scaled to weights
  vsize <- scale(edge.df$weights, center = min(quantile(edge.df$weights)), scale = max(quantile(edge.df$weights)) - min(quantile(edge.df$weights)))
  
  
  plot(
    net,
   # vertex.label = vlabel,
    vertex.label.cex = 1,
    vertex.label.dist = 1,
    vertex.label.color = 'black',
    vertex.size = c(1, vsize+0.1)*10,
    #vertex.size = vsize*10,
    vertex.frame.color = NA,
    edge.arrow.size = 0.5,
  #  vertex.color = vcl, 
    main = regulon.name
  )
  

}

adj <- read.csv('data/SCENIC-Demo-SCENIC-adj.csv')
names(regulons) <- gsub("[(+)]", "", names(regulons))
lapply(names(regulons)[grep('MUC2', regulons)][3], visuNetwork)

## [[1]]
## NULL