Commit 1eca02fc authored by smorabit's avatar smorabit
Browse files

module score functions and enrichr functions

parent 83d24a9f

R/WGCNA_functions.R

+201 −9
Original line number Diff line number Diff line
@@ -76,11 +76,19 @@ SelectNetworkGenes <- function(seurat_obj, gene_select="variable", fraction=0.05 
#' @export

#' @examples

#' SetupForWGCNA(pbmc)

SetupForWGCNA <- function(seurat_obj, wgcna_name, ...){

SetupForWGCNA <- function(seurat_obj, wgcna_name, group=NULL, ...){



  # set the active WGCNA variable

  seurat_obj <- SetActiveWGCNA(seurat_obj, wgcna_name)



  # set the active WGCNA group, this is used for actually

  # making the co-expression network

  if(is.null(group)){

    seurat_obj <- SetWGCNAGroup(seurat_obj, group='all', wgcna_name)

  } else{

    seurat_obj <- SetWGCNAGroup(seurat_obj, group, wgcna_name)

  }



  # select genes for WGCNA:

  seurat_obj <- SelectNetworkGenes(seurat_obj, ...)
@@ -102,15 +110,15 @@ SetupForWGCNA <- function(seurat_obj, wgcna_name, ...){ 
TestSoftPowers <- function(

  seurat_obj,

  use_metacells = TRUE,

  group.by=NULL, group_name=NULL,

  powers=c(seq(1,10,by=1), seq(12,30, by=2)),

  make_plot=TRUE, outfile="softpower.pdf", figsize=c(7,7)

){



  # add datExpr if not already added:

  if(!("datExpr" %in% names(GetActiveWGCNA(seurat_obj)))){

    seurat_obj <- SetDatExpr(seurat_obj, use_metacells)

    seurat_obj <- SetDatExpr(seurat_obj, use_metacells, group.by=group.by, group_name=group_name)

  }



  # get datExpr

  datExpr <- GetDatExpr(seurat_obj)
@@ -180,7 +188,8 @@ TestSoftPowers <- function( 
#' @examples

#' ConstructNetwork(pbmc)

ConstructNetwork <- function(

  seurat_obj, soft_power=NULL, use_metacells=TRUE, cur_celltype='net', tom_outdir="TOM",

  seurat_obj, soft_power=NULL, use_metacells=TRUE, group.by=NULL, group_name=NULL,

  tom_outdir="TOM",

  blocks=NULL, maxBlockSize=30000, randomSeed=12345, corType="pearson",

  consensusQuantile=0.3, networkType = "signed", TOMType = "unsigned",

  TOMDenom = "min", scaleTOMs = TRUE, scaleQuantile = 0.8,
@@ -192,20 +201,24 @@ ConstructNetwork <- function( 


  # add datExpr if not already added:

  if(!("datExpr" %in% names(GetActiveWGCNA(seurat_obj)))){

    seurat_obj <- SetDatExpr(seurat_obj, use_metacells)

    seurat_obj <- SetDatExpr(seurat_obj, use_metacells, group.by=group.by, group_name=group_name)

  }



  # get datExpr

  datExpr <- GetDatExpr(seurat_obj)



  if(is.null(group_name)){

    group_name <- 'all'

  }



  # Add functionality to accept multiExpr and perform consensus WGCNA

  # TODO



  nSets = 1

  setLabels = gsub(' ', '_', cur_celltype)

  setLabels = gsub(' ', '_', group_name)

  shortLabels = setLabels

  multiExpr <- list()

  multiExpr[[cur_celltype]] <- list(data=datExpr)

  multiExpr[[group_name]] <- list(data=datExpr)

  checkSets(multiExpr) # check data size



  if(!dir.exists(tom_outdir)){
@@ -235,7 +248,7 @@ ConstructNetwork <- function( 
    consensusTOMFilePattern = "ConsensusTOM-block.%b.rda", ...)



  # rename consensusTOM file:

  file.rename('ConsensusTOM-block.1.rda', paste0('TOM/', gsub(' ', '_',cur_celltype), '_ConsensusTOM-block.1.rda'))

  file.rename('ConsensusTOM-block.1.rda', paste0('TOM/', gsub(' ', '_',group_name), '_ConsensusTOM-block.1.rda'))



  # add network parameters to the Seurat object:
@@ -267,6 +280,7 @@ ConstructNetwork <- function( 


  # set the modules df in the Seurat object

  mods <- GetNetworkData(seurat_obj)$colors

  print(mods)

  seurat_obj <- SetModules(

    seurat_obj, mod_df = data.frame(

      "gene_name" = names(mods),
@@ -392,7 +406,8 @@ ModuleEigengenes <- function(seurat_obj, group.by.vars=NULL, verbose=TRUE, wgcna 
  }



  # set module factor levels based on order

  MEs <- GetMEs(cur_seurat, harmonized, wgcna_name)

  MEs <- GetMEs(seurat_obj, harmonized, wgcna_name)

  print(colnames(MEs))

  modules$module <- factor(

    as.character(modules$module),

    levels=colnames(MEs)
@@ -407,7 +422,129 @@ ModuleEigengenes <- function(seurat_obj, group.by.vars=NULL, verbose=TRUE, wgcna 
  seurat_obj

}



#' ModuleExprScores

#'

#' Computes module eigengenes for all WGCNA co-expression modules

#'

#' @param seurat_obj A Seurat object

#' @param method Seurat or UCell?

#' @keywords scRNA-seq

#' @export

#' @examples

#'  ModuleExprScores (pbmc)

ModuleExprScore <- function(

  seurat_obj, n_genes = 25,

   wgcna_name=NULL, method='Seurat', ...

 ){



  # set as active assay if wgcna_name is not given

  if(is.null(wgcna_name)){wgcna_name <- seurat_obj@misc$active_wgcna}



  # get modules:

  modules <- GetModules(seurat_obj, wgcna_name)

  mods <- as.character(unique(modules$module))

  mods <- mods[mods != 'grey']



  # use all genes in the module?

  if(n_genes == "all"){

    gene_list <- lapply(mods, function(cur_mod){

      subset(modules, module == cur_mod) %>% .$gene_name

    })

  } else{

    gene_list <- lapply(mods, function(cur_mod){

      cur <- subset(modules, module == cur_mod)

      cur[,c('gene_name', paste0('kME_', cur_mod))] %>%

        top_n(n_genes) %>% .$gene_name

    })

  }

  names(gene_list) <- mods



  # compute Module Scores with Seurat or UCell:

  if(method == "Seurat"){

    mod_scores <- Seurat::AddModuleScore(

      seurat_obj, features=gene_list, ...

    )@meta.data

  } else if(method == "UCell"){

    mod_scores <- UCell::AddModuleScore_UCell(

      seurat_obj, features=gene_list, ...

    )@meta.data

  } else{

    stop("Invalid method selection. Valid choices are Seurat, UCell")

  }



  mod_scores <- mod_scores[,(ncol(mod_scores)-length(mods)+1):ncol(mod_scores)]



  # rename module scores:

  colnames(mod_scores) <- mods



  # order cols

  col_order <- levels(modules$module)

  col_order <- col_order[col_order != 'grey']

  mod_scores <- mod_scores[,col_order]



  # add module scores to seurat object

  seurat_obj <- SetModuleScores(seurat_obj, mod_scores, wgcna_name)



  seurat_obj



}





#' AverageModuleExpr

#'

#' Computes module eigengenes for all WGCNA co-expression modules

#'

#' @param seurat_obj A Seurat object

#' @keywords scRNA-seq

#' @export

#' @examples

#'  AverageModuleExpr (pbmc)

AvgModuleExpr <- function(seurat_obj, n_genes = 25, wgcna_name=NULL, ...){



  # set as active assay if wgcna_name is not given

  if(is.null(wgcna_name)){wgcna_name <- seurat_obj@misc$active_wgcna}



  # get modules:

  modules <- GetModules(seurat_obj, wgcna_name)

  mods <- as.character(unique(modules$module))

  mods <- mods[mods != 'grey']



  # get wgcna params

  params <- GetWGCNAParams(seurat_obj, wgcna_name)



  # datExpr for full expression dataset

  datExpr <- GetAssayData(

    seurat_obj,

    assay=params$metacell_assay,

    slot=params$metacell_slot

  )



  # use all genes in the module?

  if(n_genes == "all"){

    gene_list <- lapply(mods, function(cur_mod){

      subset(modules, module == cur_mod) %>% .$gene_name

    })

  } else{

    gene_list <- lapply(mods, function(cur_mod){

      cur <- subset(modules, module == cur_mod)

      cur[,c('gene_name', paste0('kME_', cur_mod))] %>%

        top_n(n_genes) %>% .$gene_name

    })

  }

  names(gene_list) <- mods



  # for each module, compute average expression of genes:

  avg_mods <- t(do.call(rbind,lapply(gene_list, function(cur_genes){colMeans(datExpr[cur_genes,])})))



  # order cols

  col_order <- levels(modules$module)

  col_order <- col_order[col_order != 'grey']

  avg_mods <- avg_mods[,col_order]



  # add avg module expression to Seurat object

  seurat_obj <- SetAvgModuleExpr(seurat_obj, avg_mods, wgcna_name)

  seurat_obj

}



#' ModuleConnectivity

#'
@@ -459,3 +596,58 @@ ModuleConnectivity <- function(seurat_obj, harmonized=TRUE, wgcna_name=NULL, ... 
  seurat_obj



}



#' RunEnrichr

#'

#' Computes intramodular connectivity (kME) based on module eigengenes.

#'

#' @param seurat_obj A Seurat object

#' @param dbs List of EnrichR databases

#' @param max_genes Max number of genes to include per module, ranked by kME.

#' @param wgcna_name

#' @keywords scRNA-seq

#' @export

#' @examples

#' RunEnrichr

RunEnrichr <- function(

  seurat_obj,

  dbs = c('GO_Biological_Process_2021','GO_Cellular_Component_2021','GO_Molecular_Function_2021'),

  max_genes = 100,

  wgcna_name=NULL, ...

){



  # get data from active assay if wgcna_name is not given

  if(is.null(wgcna_name)){wgcna_name <- seurat_obj@misc$active_wgcna}



  # get modules:

  modules <- GetModules(seurat_obj, wgcna_name)

  mods <- as.character(unique(modules$module))



  # exclude grey

  mods <- mods[mods != 'grey']



  # run EnrichR for loop:

  combined_output <- data.frame()

  for(i in 1:length(mods)){

  	cur_mod <- mods[i]

    cur_info <- subset(modules, module == cur_mod)

    cur_info <- cur_info[,c('gene_name', paste0('kME_', cur_mod))]

    cur_genes <- top_n(cur_info, max_genes) %>% .$gene_name %>% as.character

    enriched <- enrichR::enrichr(cur_genes, dbs)



    # collapse into one db

    for(db in names(enriched)){

      cur_df <- enriched[[db]]

      if (nrow(cur_df) > 1){

        cur_df$db <- db

        cur_df$module <- cur_mod

        combined_output <- rbind(combined_output, cur_df)

      }

    }

  }



  # add GO term table to seurat object

  seurat_obj <- SetEnrichrTable(seurat_obj, combined_output, wgcna_name)

  seurat_obj



}

R/getters_and_setters.R

+88 −5
Original line number Diff line number Diff line
@@ -30,6 +30,25 @@ GetWGCNA <- function(seurat_obj, wgcna_name=NULL){ 
  seurat_obj@misc[[wgcna_name]]

}



############################

# WGCNA Group

###########################



SetWGCNAGroup <- function(seurat_obj, group, wgcna_name){



  if(is.null(wgcna_name)){wgcna_name <- seurat_obj@misc$active_wgcna}



  # add gene list to Seurat obj

  seurat_obj@misc[[wgcna_name]]$wgcna_group <- group

  seurat_obj

}



GetWGCNAGroup <- function(seurat_obj, wgcna_name){

  if(is.null(wgcna_name)){wgcna_name <- seurat_obj@misc$active_wgcna}

  seurat_obj@misc[[wgcna_name]]$wgcna_group

}





############################

# metacell object

###########################
@@ -88,7 +107,7 @@ GetWGCNAGenes <- function(seurat_obj, wgcna_name=NULL){ 
#' @export

#' @examples

#' SetDatExpr(pbmc)

SetDatExpr <- function(seurat_obj, use_metacells=TRUE, wgcna_name=NULL, ...){

SetDatExpr <- function(seurat_obj, use_metacells=TRUE, wgcna_name=NULL, group.by=NULL, group_name=NULL, ...){



  # get data from active assay if wgcna_name is not given

  if(is.null(wgcna_name)){wgcna_name <- seurat_obj@misc$active_wgcna}
@@ -105,25 +124,32 @@ SetDatExpr <- function(seurat_obj, use_metacells=TRUE, wgcna_name=NULL, ...){ 
    s_obj <- seurat_obj

  }



  # columns to group by

  if(!is.null(group.by)){

    cells <- s_obj@meta.data %>% subset(get(group.by) == group_name) %>% rownames

  } else{

    cells <- colnames(s_obj)

  }



  # get expression data from seurat obj

  datExpr <- as.data.frame(

    Seurat::GetAssayData(

      s_obj,

      assay=assay,

      slot='data'

    )[genes_use,]

    )[genes_use,cells]

  )



  # transpose data

  datExpr <- as.data.frame(t(datExpr))



  # only get good genes:

  good_genes <- WGCNA::goodGenes(datExpr, ...)

  gene_list = GetWGCNAGenes(seurat_obj)[WGCNA::goodGenes(datExpr, ...)]



  # update the WGCNA gene list:

  seurat_obj <- SetWGCNAGenes(seurat_obj, good_genes, wgcna_name)

  seurat_obj <- SetWGCNAGenes(seurat_obj, gene_list, wgcna_name)



  datExpr <- datExpr[,good_genes]

  datExpr <- datExpr[,gene_list]



  # set the datExpr in the Seurat object

  seurat_obj@misc[[wgcna_name]]$datExpr <- datExpr
@@ -277,3 +303,60 @@ GetMEs <- function(seurat_obj, harmonized=TRUE, wgcna_name=NULL){ 
  }

  MEs

}



############################

# GO term table

###########################



SetEnrichrTable <- function(seurat_obj, enrich_table, wgcna_name=NULL){



  # get data from active assay if wgcna_name is not given

  if(is.null(wgcna_name)){wgcna_name <- seurat_obj@misc$active_wgcna}



  # set enrichr table

  seurat_obj@misc[[wgcna_name]]$enrichr_table <- enrich_table

  seurat_obj

}





GetEnrichrTable <- function(seurat_obj,  wgcna_name=NULL){

  if(is.null(wgcna_name)){wgcna_name <- seurat_obj@misc$active_wgcna}

  seurat_obj@misc[[wgcna_name]]$enrichr_table

}





############################

# Module Scores

###########################



SetModuleScores <- function(seurat_obj, mod_scores, wgcna_name=NULL){



  # get data from active assay if wgcna_name is not given

  if(is.null(wgcna_name)){wgcna_name <- seurat_obj@misc$active_wgcna}

  seurat_obj@misc[[wgcna_name]]$module_scores <- mod_scores

  seurat_obj

}





GetModuleScores <- function(seurat_obj,  wgcna_name=NULL){

  if(is.null(wgcna_name)){wgcna_name <- seurat_obj@misc$active_wgcna}

  seurat_obj@misc[[wgcna_name]]$module_scores

}



############################

# Average Module Expression

###########################



SetAvgModuleExpr <- function(seurat_obj, avg_mods, wgcna_name=NULL){



  # get data from active assay if wgcna_name is not given

  if(is.null(wgcna_name)){wgcna_name <- seurat_obj@misc$active_wgcna}

  seurat_obj@misc[[wgcna_name]]$avg_modules <- avg_mods

  seurat_obj

}





GetAvgModuleExpr <- function(seurat_obj,  wgcna_name=NULL){

  if(is.null(wgcna_name)){wgcna_name <- seurat_obj@misc$active_wgcna}

  seurat_obj@misc[[wgcna_name]]$avg_modules

}

R/metacells.R

+12 −2
Original line number Diff line number Diff line
@@ -13,7 +13,11 @@ 
#' @export

#' @examples

#' ConstructMetacells(pbmc)

ConstructMetacells <- function(seurat_obj, name='agg', k=50, reduction='umap', assay='RNA', slot='counts',  meta=NULL, return_metacell=FALSE){

ConstructMetacells <- function(

  seurat_obj, name='agg', ident.group='seurat_clusters', k=50,

  reduction='umap', assay='RNA',

  slot='counts',  meta=NULL, return_metacell=FALSE

){



  # check reduction

  if(!(reduction %in% names(seurat_obj@reductions))){
@@ -137,7 +141,10 @@ ConstructMetacells <- function(seurat_obj, name='agg', k=50, reduction='umap', a 
#' @export

#' @examples

#' MetacellsByGroups(pbmc)

MetacellsByGroups <- function(seurat_obj, group.by=c('seurat_clusters'), k=50, reduction='umap', assay='RNA', slot='counts'){

MetacellsByGroups <- function(

  seurat_obj, group.by=c('seurat_clusters'), ident.group='seurat_clusters',

  k=50, reduction='umap', assay='RNA', slot='counts'

){



  # setup grouping variables

  if(length(group.by) > 1){
@@ -179,6 +186,9 @@ MetacellsByGroups <- function(seurat_obj, group.by=c('seurat_clusters'), k=50, r 
  # combine metacell objects

  metacell_obj <- merge(metacell_list[[1]], metacell_list[2:length(metacell_list)])



  # set idents for metacell object:

  Idents(metacell_obj) <- metacell_obj@meta.data[[ident.group]]



  # add seurat metacell object to the main seurat object:

  seurat_obj <- SetMetacellObject(seurat_obj, metacell_obj)

R/plotting.R

+180 −33
Original line number Diff line number Diff line
@@ -48,7 +48,9 @@ MECorrelogram <- function( 
  type = 'upper', order='original', method='ellipse',

  tl.col = 'black', tl.srt=45,

  sig.level = c(0.0001, 0.001, 0.01, 0.05),

  insig='label_sig', pch.cex=0.7, col=NULL, ncolors=200, ...){

#  insig='label_sig',

  pch.cex=0.7, col=NULL, ncolors=200, ...

){



  if(is.null(wgcna_name)){wgcna_name <- seurat_obj@misc$active_wgcna}
@@ -76,7 +78,8 @@ MECorrelogram <- function( 
    type=type, order=order,

    method=method, tl.col=tl.col,

    tl.srt=tl.srt, sig.level=sig.level,

    insig=insig, pch.cex=pch.cex,

    # insig=insig,

    pch.cex=pch.cex,

    col = col, ...

  )
@@ -114,11 +117,12 @@ ModuleCorrNetwork <- function( 
  mods  <-  colnames(MEs)



  # what clusters are we using?

  if(is.null(clusters)){

    cluster_col <- Idents(seurat_obj)

  if(is.null(cluster_col)){

    clusters <- Idents(seurat_obj)

  } else{

    clusters <- droplevels(seurat_obj@meta.data[[cluster_col]])

  }

  clusters <- droplevels(seurat_obj$annotation)



  MEs$cluster <- clusters



  # compute average MEs for each cluster
@@ -183,12 +187,24 @@ ModuleCorrNetwork <- function( 


  # igraph layout

  e <-  get.edgelist(g, name=FALSE)

  l <- igraph::layout_with_fr(

    g, weights=E(g)$value,

    coords = as.matrix(data.frame(x=V(g)$x, y=V(g)$y)),

    niter=niter

  # l <- igraph::layout_with_fr(

  #   g,

  #   weights=E(g)$value,

  #   coords = as.matrix(data.frame(x=V(g)$x, y=V(g)$y)),

  #   niter=niter

  # )

  #

  l <- qgraph::qgraph.layout.fruchtermanreingold(

    e, vcount = vcount(g),

    weights=E(g)$value,

    repulse.rad=(vcount(g)),

    #cool.exp = 0.5,

    # init = as.matrix(data.frame(x=V(g)$x, y=V(g)$y)),

    niter=niter,

    #max.delta = vcount(g)/2

  )





  # ggplot just to get the colors

  plot_df <- rbind(cor_df, data.frame(Var1=c('x', 'y'), Var2=c('y', 'x'), value=c(-1,1)))

  temp <- ggplot(plot_df, aes(x=value, y=value, color=value)) +
@@ -218,36 +234,55 @@ ModuleCorrNetwork <- function( 
  )



  # plot colorbar:

  colfunc <- colorRampPalette(c( "seagreen", "white", "darkorchid1"))

  image.plot(legend.only=T, zlim=c(-1,1), col=colfunc(256) )



  # TODO

  # Possibly get rid of this because it depends on a super random package

  # (fields) and kinda looks ugly?

  # colfunc <- colorRampPalette(c( "seagreen", "white", "darkorchid1"))

  # image.plot(legend.only=T, zlim=c(-1,1), col=colfunc(256) )

  #



}





#' MEFeaturePlot

#' ModuleFeaturePlot

#'

#' Plot module eigengenes as a FeaturePlot

#'

#' @param seurat_obj A Seurat object

#' @param features What to plot? Can select hMEs, MEs, scores, or average

#' @param order TRUE, FALSE, or "shuffle" are valid options

#' @keywords scRNA-seq

#' @export

#' @examples

#' MEFeaturePlot

MEFeaturePlot<- function(

#' ModuleFeaturePlot

ModuleFeaturePlot<- function(

  seurat_obj, module_names=NULL, wgcna_name = NULL,

  harmonized=TRUE, reduction='umap',

  order=TRUE, restrict_range=TRUE, point_size = 0.5,

  label_legend = FALSE

  reduction='umap', features = 'hMEs',

  order=TRUE, restrict_range=TRUE, point_size = 0.5, alpha=1,

  label_legend = FALSE, ucell = FALSE

){



  if(is.null(wgcna_name)){wgcna_name <- seurat_obj@misc$active_wgcna}



  # get MEs, module data from seurat object

  MEs <- GetMEs(seurat_obj, harmonized, wgcna_name)

  modules <- GetModules(seurat_obj, wgcna_name)

  if(features == 'hMEs'){

    MEs <- GetMEs(seurat_obj, TRUE, wgcna_name)

  } else if(features == 'MEs'){

    MEs <- GetMEs(seurat_obj, FALSE, wgcna_name)

  } else if(features == 'scores'){

    MEs <- GetModuleScores(seurat_obj, wgcna_name)

  } else if(features == 'average'){

    MEs <- GetAvgModuleExpr(seurat_obj, wgcna_name)

    restrict_range <- FALSE

  } else(

    stop('Invalid feature selection. Valid choices: hMEs, MEs, scores, average')

  )



  # override restrict_range if ucell

  if(ucell){restrict_range <- FALSE}



  # use all modules except gray if not specified by the user

  modules <- GetModules(seurat_obj, wgcna_name)

  if(is.null(module_names)){

    module_names <- colnames(MEs)

    module_names <- module_names[module_names != 'grey']
@@ -280,28 +315,41 @@ MEFeaturePlot<- function( 
    }



    # order points:

    if(order){

    if(order == TRUE){

      plot_df <- plot_df %>% dplyr::arrange_(cur_mod)

    } else if(order == "shuffle"){

      plot_df <- plot_df[sample(nrow(plot_df)),]

    }



    # label for plot:

    if(harmonized){

      label <- paste0('hME_', cur_mod)

    } else{

      label <- paste0('ME_', cur_mod)

    }

    label <- cur_mod



    # plot with ggplot

    plot_list[[cur_mod]] <- plot_df %>%

     p <- plot_df %>%

      ggplot(aes_string(x=x_name, y=y_name, color=cur_mod)) +

      geom_point(size=point_size) +

      scale_color_gradient2(

        low='grey75', mid='grey95', high=cur_color,

        breaks = c(plot_range[1], 0, plot_range[2]),

        labels = c('-', '0', '+'),

      ) +

      geom_point(size=point_size, alpha=alpha) +

      ggtitle(label) + umap_theme +

      labs(color="")



    # UCell?

    if(!ucell){

      p <- p + scale_color_gradient2(

        low='grey75', mid='grey95', high=cur_color,

        breaks = c(plot_range[1], plot_range[2]),

        labels = c('-', '+'),

        guide = guide_colorbar(ticks=FALSE, barwidth=0.5, barheight=4)

      )

    } else{

      p <- p + scale_color_gradient(

        low='grey95', high=cur_color,

        breaks = c(plot_range[1], plot_range[2]),

        labels = c('0', '+'),

        guide = guide_colorbar(ticks=FALSE, barwidth=0.5, barheight=4)

      )

    }



    plot_list[[cur_mod]] <- p



  }



  # return plot
@@ -314,3 +362,102 @@ MEFeaturePlot<- function( 
  p



}











#' PlotEnrichr

#'

#' Makes barplots from Enrichr data

#'

#' @param seurat_obj A Seurat object

#' @param dbs List of EnrichR databases

#' @param max_genes Max number of genes to include per module, ranked by kME.

#' @param wgcna_name

#' @keywords scRNA-seq

#' @export

#' @examples

#' PlotEnrichr

PlotEnrichr <- function(

  seurat_obj, outdir = "enrichr_plots",

  n_terms = 25, plot_size = c(6,15),

  wgcna_name=NULL, logscale=FALSE, ...

){



  # get data from active assay if wgcna_name is not given

  if(is.null(wgcna_name)){wgcna_name <- seurat_obj@misc$active_wgcna}



  # get modules:

  modules <- GetModules(seurat_obj, wgcna_name)

  mods <- as.character(unique(modules$module))

  mods <- mods[mods != 'grey']



  # get Enrichr table

  enrichr_df <- GetEnrichrTable(seurat_obj, wgcna_name)



  # helper function to wrap text

  wrapText <- function(x, len) {

      sapply(x, function(y) paste(strwrap(y, len), collapse = "\n"), USE.NAMES = FALSE)

  }



  # make output dir if it doesn't exist:

  if(!dir.exists(outdir)){dir.create(outdir)}



  # loop through modules:

  for(i in 1:length(mods)){



    cur_mod <- mods[i]

    cur_terms <- subset(enrichr_df, module == cur_mod)

    print(cur_mod)



    # get color for this module:

    cur_color <- modules %>% subset(module == cur_mod) %>% .$color %>% unique %>% as.character



    # skip if there are not any terms for this module:

    if(nrow(cur_terms) == 0){next}

    cur_terms$wrap <- wrapText(cur_terms$Term, 45)



    # plot top n_terms as barplot

    plot_list <- list()

    for(cur_db in dbs){



      plot_df <- subset(cur_terms, db==cur_db) %>% top_n(n_terms, wt=Combined.Score)



      # text color:

      if(cur_mod == 'black'){

        text_color = 'grey'

      } else {

        text_color = 'black'

      }



      # logscale?

      if(logscale){

        plot_df$Combined.Score <- log2(plot_df$Combined.Score)

        lab <- 'Enrichment log2(combined score)'

        x <- 0.2

      } else{lab <- 'Enrichment (combined score)'; x <- 5}



      # make bar plot:

      plot_list[[cur_db]] <- ggplot(plot_df, aes(x=Combined.Score, y=reorder(wrap, Combined.Score)))+

        geom_bar(stat='identity', position='identity', color='white', fill=cur_color) +

        geom_text(aes(label=wrap), x=x, color=text_color, size=3.5, hjust='left') +

        ylab('Term') + xlab(lab) + ggtitle(cur_db) +

        theme(

          panel.grid.major=element_blank(),

          panel.grid.minor=element_blank(),

          legend.title = element_blank(),

          axis.ticks.y=element_blank(),

          axis.text.y=element_blank(),

          plot.title = element_text(hjust = 0.5)

        )

    }



    # make pdfs in output dir

    pdf(paste0(outdir, '/', cur_mod, '.pdf'), width=plot_size[1], height=plot_size[2])

    for(plot in plot_list){

      print(plot)

    }

    dev.off()

  }

}

CRA Git | Maintained and supported by SUSTech CRA and CCSE