Commit dd7fdbe7 authored by smorabit's avatar smorabit
Browse files

added enrichment tutorial

parent bef3e8d0
Loading
Loading
Loading
Loading
+27 −17
Original line number Diff line number Diff line
@@ -919,12 +919,12 @@ ModuleConnectivity <- function(seurat_obj, harmonized=TRUE, wgcna_name=NULL, ...

#' RunEnrichr
#'
#' Computes intramodular connectivity (kME) based on module eigengenes.
#' Run Enrichr gene set enrichment tests on scWGCNA modules
#'
#' @param seurat_obj A Seurat object
#' @param dbs List of EnrichR databases
#' @param dbs character vector of EnrichR databases
#' @param max_genes Max number of genes to include per module, ranked by kME.
#' @param wgcna_name
#' @param wgcna_name The name of the scWGCNA experiment in the seurat_obj@misc slot
#' @keywords scRNA-seq
#' @export
#' @examples
@@ -975,23 +975,30 @@ RunEnrichr <- function(

#' OverlapModulesDEGs
#'
#' Computes intramodular connectivity (kME) based on module eigengenes.
#' Performs Fisher's Exact Test for overlap between DEGs and scWGCNA modules.
#'
#' @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
#' @param deg_df DEG table formatted like the output from Seurat's FindMarkers
#' @param fc_cutoff log fold change cutoff for DEGs to be included in the overlap test
#' @param group_col the name of the column in deg_df containing the cell grouping information
#' @param wgcna_name The name of the scWGCNA experiment in the seurat_obj@misc slot
#' @keywords scRNA-seq
#' @export
#' @examples
#' OverlapModulesDEGs
OverlapModulesDEGs <- function(
  seurat_obj, deg_df, wgcna_name = NULL, fc_cutoff = 0.5, ...
  seurat_obj,
  deg_df,
  fc_cutoff = 0.5,
  group_col = 'cluster',
  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}

  deg_df$group <- deg_df[,group_col]
  cell_groups <- deg_df$group %>% unique

  # get modules,
@@ -1015,8 +1022,10 @@ OverlapModulesDEGs <- function(

  # run overlaps between module gene lists and DEG lists:
  overlap_df <- do.call(rbind, lapply(mods, function(cur_mod){
    #print(cur_mod)
    cur_module_genes <- modules %>% subset(module == cur_mod) %>% .$gene_name
    cur_overlap_df <- do.call(rbind, lapply(cell_groups, function(cur_group){
      #print(cur_group)
      cur_DEGs <- deg_df %>% subset(group == cur_group & p_val_adj <= 0.05 & avg_log2FC > fc_cutoff) %>% .$gene
      cur_overlap <- testGeneOverlap(newGeneOverlap(
          cur_module_genes,
@@ -1034,6 +1043,7 @@ OverlapModulesDEGs <- function(
    cur_overlap_df
  }))


  # adjust for multiple comparisons:
  overlap_df$fdr <- p.adjust(overlap_df$pval, method='fdr')

@@ -1061,7 +1071,7 @@ OverlapModulesDEGs <- function(
#' @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
#' @param wgcna_name The name of the scWGCNA experiment in the seurat_obj@misc slot
#' @keywords scRNA-seq
#' @export
#' @examples
@@ -1477,7 +1487,7 @@ MotifScan <- function(
#' Overlap modules with TF target genes
#'
#' @param seurat_obj A Seurat object
#' @param wgcna_name
#' @param wgcna_name The name of the scWGCNA experiment in the seurat_obj@misc slot
#' @keywords scRNA-seq
#' @export
#' @examples
@@ -1604,22 +1614,23 @@ MotifTargetScore <- function(
 }



 #' RunModuleUMAP
 #'
#' Run UMAP on co-expression matrix using hub genes as features.
#'
#' @param seurat_obj A Seurat object
#' @param n_hubs
#' @param exclude_grey
#' @param wgcna_name
#' @param
#' @param wgcna_name The name of the scWGCNA experiment in the seurat_obj@misc slot
#' @param ... Additional parameters supplied to uwot::umap
#' @keywords scRNA-seq
#' @export
#' @examples
#' RunModuleUMAP
RunModuleUMAP <- function(
  seurat_obj,
  n_hubs = 50,
  exclude_grey = TRUE,
  harmonized=TRUE,
  wgcna_name = NULL,
  n_neighbors= 25,
  metric = "cosine",
@@ -1634,8 +1645,7 @@ RunModuleUMAP <- function(
  # get the TOM
  TOM <- GetTOM(seurat_obj, wgcna_name)

  # get modules, MEs:
  MEs <- GetMEs(seurat_obj, harmonized, wgcna_name)
  # get modules,
  modules <- GetModules(seurat_obj, wgcna_name)
  mods <- levels(modules$module)

@@ -1725,7 +1735,7 @@ RunModuleUMAP <- function(
#' Traits must be a categorical variable (not a character vector), or a numeric variable.
#' @param features Which features to use to summarize each modules? Valid choices are hMEs, MEs, or scores
#' @param cor_meth Which method to use for correlation? Valid choices are pearson, spearman, kendall.
#' @param wgcna_name
#' @param wgcna_name The name of the scWGCNA experiment in the seurat_obj@misc slot
#' @keywords scRNA-seq
#' @export
#' @examples
+53 −32
Original line number Diff line number Diff line
@@ -414,11 +414,13 @@ ModuleFeaturePlot<- function(

#' EnrichrBarPlot
#'
#' Makes barplots from Enrichr data
#' Makes barplots from RunEnrhcr output.
#'
#' @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 outdir directory to place output .pdf files
#' @param n_terms the number of terms to plot in each barplot
#' @param plot_size the size of the output .pdf files (width, height)
#' @param logscale logical controlling whether to plot the enrichment on a log scale
#' @param wgcna_name The name of the scWGCNA experiment in the seurat_obj@misc slot
#' @keywords scRNA-seq
#' @export
@@ -427,7 +429,7 @@ ModuleFeaturePlot<- function(
EnrichrBarPlot <- function(
  seurat_obj, outdir = "enrichr_plots",
  n_terms = 25, plot_size = c(6,15),
  wgcna_name=NULL, logscale=FALSE, ...
  logscale=FALSE, wgcna_name=NULL,  ...
){

  # get data from active assay if wgcna_name is not given
@@ -513,17 +515,20 @@ EnrichrBarPlot <- function(
#' 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 database name of the enrichr database to plot.
#' @param mods names of modules to plot. All modules are plotted if mods='all' (default)
#' @param n_terms number of enriched terms to plot for each module
#' @param break_ties logical controlling whether or not to randomly select terms with equal enrichments to precisely enforce n_terms.
#' @param logscale logical controlling whether to plot the enrichment on a log scale.
#' @param wgcna_name The name of the scWGCNA experiment in the seurat_obj@misc slot
#' @keywords scRNA-seq
#' @export
#' @examples
#' EnrichrBarPlot
#' EnrichrDotPlot
EnrichrDotPlot <- function(
  seurat_obj, database, mods="all", outdir = "enrichr_plots",
  seurat_obj, database, mods="all",
  n_terms = 3, break_ties=TRUE,
  wgcna_name=NULL, logscale=TRUE, ...
   logscale=TRUE, wgcna_name=NULL, ...
){

  # get data from active assay if wgcna_name is not given
@@ -758,7 +763,7 @@ HubGeneNetworkPlot <- function(
  return_graph=FALSE,
  edge.alpha=0.25,
  vertex.label.cex=0.5,
  hub.vertex.size=4,Ω
  hub.vertex.size=4,
  other.vertex.size=1,
  wgcna_name=NULL,
  ...
@@ -899,11 +904,13 @@ HubGeneNetworkPlot <- function(
#' Makes a igraph network plot using the module UMAP
#'
#' @param seurat_obj A Seurat object
#' @param
#' @param
#' @param
#' @param
#' @param
#' @param sample_edges logical determining whether we downsample edges for plotting (TRUE), or take the strongst edges.
#' @param edge_prop proportion of edges to plot. If sample_edges=FALSE, the strongest edges are selected.
#' @param label_hubs the number of hub genes to label in each module
#' @param edge.alpha scaling factor for edge opacity
#' @param vertex.label.cex font size for labeled genes
#' @param return_graph logical determining whether to plot thr graph (FALSE) or return the igraph object (TRUE)
#' @param wgcna_name The name of the scWGCNA experiment in the seurat_obj@misc slot
#' @keywords scRNA-seq
#' @export
#' @examples
@@ -915,8 +922,6 @@ ModuleUMAPPlot <- function(
  label_hubs = 5, # how many hub genes to label?
  edge.alpha=0.25,
  vertex.label.cex=0.5,
  hub.vertex.size=4,
  other.vertex.size=1,
  return_graph = FALSE, # this returns the igraph object instead of plotting
  wgcna_name=NULL,
  ...
@@ -1060,17 +1065,21 @@ ModuleUMAPPlot <- function(
#'
#' 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 The name of the scWGCNA experiment in the seurat_obj@misc slot
#' @param overlap_df the Module/DEG overlap table from OverlapModulesDEGs
#' @param plot_var the name of the overlap statistic to plot
#' @param logscale logical controlling whether to plot the result on a log scale, useful for odds ratio
#' @param neglog logical controlling wehether to plot the result as a negative log, useful for p-value / FDR
#' @param plot_significance logical controlling whether to plot the significance levels on top of the dots
#' @keywords scRNA-seq
#' @export
#' @examples
#' OverlapDotPlot
OverlapDotPlot <- function(
  overlap_df, plot_var = 'odds_ratio',
  logscale=TRUE, neglog=FALSE, plot_significance=TRUE, ...
  logscale=TRUE,
  neglog=FALSE,
  plot_significance=TRUE,
  ...
){

  label <- plot_var
@@ -1106,19 +1115,25 @@ OverlapDotPlot <- function(
  p
}

#' OverlapBarPlot
#'
#' Plots the results from OverlapModulesDEGs as a bar plot
#'
#' @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 The name of the scWGCNA experiment in the seurat_obj@misc slot
#' @param overlap_df the Module/DEG overlap table from OverlapModulesDEGs
#' @param plot_var the name of the overlap statistic to plot
#' @param logscale logical controlling whether to plot the result on a log scale, useful for odds ratio
#' @param neglog logical controlling wehether to plot the result as a negative log, useful for p-value / FDR
#' @param label_size the size of the module labels in the bar plot
#' @keywords scRNA-seq
#' @export
#' @examples
#' OverlapBarPlot
OverlapBarPlot <- function(
  overlap_df, plot_var = 'odds_ratio',
  logscale=TRUE, neglog=FALSE, ...
  overlap_df,
  plot_var = 'odds_ratio',
  logscale=FALSE, neglog=FALSE,
  label_size=2,
  ...
){

  label <- plot_var
@@ -1135,9 +1150,9 @@ OverlapBarPlot <- function(
    yint = log(yint)
  }
  if(neglog){
    overlap_df[[plot_var]] <- -1 * overlap_df[[plot_var]]
    label <- paste0('-', label)
    yint = -1 * yint
    overlap_df[[plot_var]] <- -1 * log(overlap_df[[plot_var]])
    label <- paste0('-log(', label, ')')
    yint = -1 * log(yint)
  }

  groups <- overlap_df$group %>% as.character %>% unique
@@ -1165,6 +1180,12 @@ OverlapBarPlot <- function(
        p <- p + geom_hline(yintercept=yint, linetype='dashed', color='gray')
      }

      # add the labels:
      p <- p +
        geom_text(
          aes(label=module, x=module, y=get(plot_var)), color='black', size=label_size, hjust='inward'
        )

      plot_list[[cur_group]] <- p
  }

+16 −4
Original line number Diff line number Diff line
@@ -40,19 +40,31 @@ navbar:
    - text: All vignettes
      href: articles/index.html



reference:
- title: Metacells
  desc: Functions for constructing metacells from single-cell data
  contents:
  - '`ConstructMetacells`'
  - '`MetacellsByGroups`'
- title: Plotting
  desc: Functions for generating plots with scWGCNA
- title: Network Visualization
  desc: Functions for visualizing the scWGCNA co-expression network
  contents:
    - '`ModuleNetworkPlot`'
    - '`HubGeneNetworkPlot`'
    - '`RunModuleUMAP`'
    - '`ModuleUMAPPlot`'
- title: Enrichment Analysis
  desc: Functions for Enrichr analysis and DEG overlap analysis
  contents:
    - '`RunEnrichr`'
    - '`EnrichrBarPlot`'
    - '`EnrichrDotPlot`'
    - '`OverlapModulesDEGs`'
    - '`OverlapBarPlot`'
    - '`OverlapDotPlot`'
- title: Plotting
  desc: Functions for generating plots with scWGCNA
  contents:
  - '`ModuleFeaturePlot`'
- title: Seurat wrappers
  desc: Wrapper functions to run Seurat commands on the metacell data
+187 −7

File changed.

Preview size limit exceeded, changes collapsed.

+330 KiB
Loading image diff...
Loading