Network inference based on normalized mutual information scores

### Define the network inference function
infer_network <- function(sce, interest_celltype = 'all', pct = 0.2, log2FC = 1, n_gene = 250){
  cat('1.Finding significant expressed genes...\n')
  marker <- FindAllMarkers(sce, min.pct = pct, only.pos = T, verbose = F) %>% 
    filter(avg_log2FC > log2FC & p_val_adj < 0.05)
  if (interest_celltype == 'all') {
    expr <- GetAssayData(sce[unique(marker$gene),], layer = 'counts') %>% t() %>% data.frame()
  } else {
    sce <- subset(sce, subset = cellType == interest_celltype)
    marker <- marker %>% filter(cluster == interest_celltype)
    expr <- GetAssayData(sce[unique(marker$gene),], layer = 'counts') %>% t() %>% data.frame()
  }
  cat('2.Inferring raw network with normalized mutual information scores...\n')
  net_raw <- nmiMatrix(expr)
  cat('3.Finding the threshold for network inferring...\n')
  ratio <- 10 * n_gene / (ncol(expr) ^ 2)
  net <- ifelse(net_raw > quantile(net_raw, (1-ratio)), 1, 0)
  while(nrow(net) > n_gene) {
    net <- ifelse(net_raw > quantile(net_raw, (1-ratio)), 1, 0)
    net <- net[rowSums(net)>1, colSums(net)>1]
    ratio <- ratio / 1.2
  }
  rownames(net) <- gsub('\\.','-',rownames(net))
  colnames(net) <- gsub('\\.','-',colnames(net))
  return(net)
}
### Define the network inference plot function
plot_network <- function(net, interest_celltype = 'all'){
  p1 <- GGally::ggnet2(net, label = T, size = 'degree', legend.position = 'none', 
                       label.size = 2.5, color = 'pink3', 
                       edge.alpha = 0.1, label.color = 'blue4',
                       layout.par = list(repulse.rad = 5e3)) + 
    ggtitle(paste0('Network of the ', interest_celltype))
  top_genes <- rowSums(net) %>% sort(decreasing = T) %>% head(30)
  top_genes <- data.frame(id = names(top_genes), degree = top_genes)
  p2 <- top_genes %>% ggplot(aes(x = degree, y = reorder(id, X = degree))) + 
    geom_col(fill='pink3') + theme_minimal() + ylab('') + theme(legend.position = 'none') + 
    ggtitle('Top 30 most-connected genes') + 
    geom_vline(xintercept = 10, linetype = 2, colour = 'grey')
  p1 + p2 + patchwork::plot_layout(design = 'AAAAB')
}

### Load the data and process the analysis
p_load(praznik, tidyverse, qs, Seurat)
sce <- qread('sce.qs', nthreads = 8)
interest_celltype <- 'all'
net <- infer_network(sce, interest_celltype)
plot_network(net, interest_celltype)

Example network inferred from GSE120839