Commit 687007cd authored by smorabit's avatar smorabit
Browse files

fix kME issue

parent 8f438000

DESCRIPTION

+2 −2
Original line number Diff line number Diff line 
Package: scWGCNA

Package: hdWGCNA

Title: WGCNA

Version: 0.0.0.9000

Version: 0.0.1

Authors@R: c(

    person("Sam", "Morabito", , "smorabit@uci.edu", role = c("aut", "cre"),

           comment = c(ORCID = "0000-0002-7768-4856")),

NAMESPACE

+2 −0
Original line number Diff line number Diff line
@@ -13,6 +13,7 @@ export(GetActiveWGCNA) 
export(GetAvgModuleExpr)

export(GetDatExpr)

export(GetEnrichrTable)

export(GetMELoadings)

export(GetMEs)

export(GetMetacellObject)

export(GetModulePreservation)
@@ -74,6 +75,7 @@ export(SetActiveWGCNA) 
export(SetAvgModuleExpr)

export(SetDatExpr)

export(SetEnrichrTable)

export(SetMELoadings)

export(SetMEs)

export(SetMetacellObject)

export(SetModulePreservation)

R/WGCNA_functions.R

+143 −171
Original line number Diff line number Diff line
@@ -274,7 +274,7 @@ ConstructNetwork <- function( 
  multi.group.by = NULL,

  multi_groups = NULL,

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

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

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

  TOMDenom = "min", scaleTOMs = TRUE, scaleQuantile = 0.8,

  sampleForScaling = TRUE, sampleForScalingFactor = 1000,

  useDiskCache = TRUE, chunkSize = NULL,
@@ -417,128 +417,6 @@ ConstructNetwork <- function( 


}



#

# #' ConstructNetwork

# #'

# #' This function constructs a co-expression network from a Seurat object

# #'

# #' @param seurat_obj A Seurat object

# #' @param soft_power

# #' @keywords scRNA-seq

# #' @export

# #' @examples

# #' ConstructNetwork(pbmc)

# ConstructNetwork <- function(

#   seurat_obj, soft_power=NULL, use_metacells=TRUE,

#   setDatExpr=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,

#   sampleForScaling = TRUE, sampleForScalingFactor = 1000,

#   useDiskCache = TRUE, chunkSize = NULL,

#   deepSplit = 4, pamStage=FALSE, detectCutHeight = 0.995, minModuleSize = 50,

#   mergeCutHeight = 0.2, saveConsensusTOMs = TRUE, ...

# ){

#

#   # add datExpr if not already added:

#   if(!("datExpr" %in% names(GetActiveWGCNA(seurat_obj))) | setDatExpr == TRUE){

#     seurat_obj <- SetDatExpr(

#       seurat_obj,

#       group_name = group_name,

#       group.by=group.by,

#       use_metacells=use_metacells,

#       return_seurat=TRUE

#      )

#   }

#

#   # get datExpr

#   datExpr <- GetDatExpr(seurat_obj)

#

#   if(is.null(group_name)){

#     group_name <- 'all'

#   }

#

#   nSets = 1

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

#   shortLabels = setLabels

#   multiExpr <- list()

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

#   checkSets(multiExpr) # check data size

#

#   if(!dir.exists(tom_outdir)){

#     dir.create(tom_outdir)

#   }

#

#

#   net <- WGCNA::blockwiseConsensusModules(

#     multiExpr,

#     power = soft_power,

#     blocks = blocks,

#     maxBlockSize = maxBlockSize, ## This should be set to a smaller size if the user has limited RAM

#     randomSeed = randomSeed,

#     corType = corType,

#     consensusQuantile = consensusQuantile,

#     networkType = networkType,

#     TOMType = TOMType,

#     TOMDenom = TOMDenom,

#     scaleTOMs = scaleTOMs, scaleQuantile = scaleQuantile,

#     sampleForScaling = sampleForScaling, sampleForScalingFactor = sampleForScalingFactor,

#     useDiskCache = useDiskCache, chunkSize = chunkSize,

#     deepSplit = deepSplit,

#     pamStage=pamStage,

#     detectCutHeight = detectCutHeight, minModuleSize = minModuleSize,

#     mergeCutHeight = mergeCutHeight,

#     saveConsensusTOMs = saveConsensusTOMs,

#     consensusTOMFilePattern = "ConsensusTOM-block.%b.rda", ...)

#

#   # rename consensusTOM file:

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

#

#   # add network parameters to the Seurat object:

#

#   params <- list(

#     power = soft_power,

#     blocks = blocks,

#     maxBlockSize = maxBlockSize, ## This should be set to a smaller size if the user has limited RAM

#     randomSeed = randomSeed,

#     corType = corType,

#     consensusQuantile = consensusQuantile,

#     networkType = networkType,

#     TOMType = TOMType,

#     TOMDenom = TOMDenom,

#     scaleTOMs = scaleTOMs, scaleQuantile = scaleQuantile,

#     sampleForScaling = sampleForScaling, sampleForScalingFactor = sampleForScalingFactor,

#     useDiskCache = useDiskCache, chunkSize = chunkSize,

#     deepSplit = deepSplit,

#     pamStage=pamStage,

#     detectCutHeight = detectCutHeight, minModuleSize = minModuleSize,

#     mergeCutHeight = mergeCutHeight,

#     saveConsensusTOMs = saveConsensusTOMs

#   )

#

#   # add parameters:

#   seurat_obj <- SetWGCNAParams(seurat_obj, params)

#

#   # append working directory to the TOM file so it has the full path:

#   net$TOMFiles <- paste0(getwd(), '/TOM/', group_name, '_', net$TOMFiles)

#

#   # add network to seurat obj

#   seurat_obj <- SetNetworkData(seurat_obj, net)

#

#   # set the modules df in the Seurat object

#   mods <- GetNetworkData(seurat_obj)$colors

#   seurat_obj <- SetModules(

#     seurat_obj, mod_df = data.frame(

#       "gene_name" = names(mods),

#       "module" = factor(mods, levels=unique(mods)),

#       "color" = mods

#     )

#   )

#

#   seurat_obj

#

# }



#' ComputeModuleEigengene

#'
@@ -546,13 +424,24 @@ ConstructNetwork <- function( 
#'

#' @param seurat_obj A Seurat object

#' @param cur_mod name of a module found in seurat_obj@misc[[seurat_obj@misc$active_wgcna]]$wgcna_net$colors

#' @param modules table containing module / gene assignments, as in GetModules(seurat_obj).

#' @param group.by.vars groups to harmonize by

#' @param verbose logical indicating whether to print messages

#' @param vars.to.regress character vector of variables in seurat_obj@meta.data to regress when running ScaleData

#' @param scale.model.use model to scale data when running ScaleData choices are "linear", "poisson", or "negbinom"

#' @param wgcna_name name of the WGCNA experiment

#' @keywords scRNA-seq

#' @export

#' @examples

#' ConstructNetwork(pbmc)

#' ComputeModuleEigengene(pbmc)

ComputeModuleEigengene <- function(

  seurat_obj, cur_mod, modules,

  group.by.vars=NULL, verbose=TRUE,

  seurat_obj,

  cur_mod,

  modules,

  group.by.vars=NULL,

  verbose=TRUE,

  vars.to.regress = NULL,

  scale.model.use = 'linear',

  wgcna_name=NULL, ...

){
@@ -562,35 +451,91 @@ ComputeModuleEigengene <- function( 
  # get genes in this module:

  cur_genes <- modules %>% subset(module == cur_mod) %>% .$gene_name



  # subset seurat object by these genes only:

  cur_seurat <- seurat_obj[cur_genes,]



  # scale the subsetted expression dataset:

  if(is.null(vars.to.regress)){

    cur_seurat <- ScaleData(cur_seurat, features=rownames(cur_seurat), model.use=scale.model.use)

  } else if(all(vars.to.regress %in% colnames(seurat_obj@meta.data))){

    cur_seurat <- ScaleData(cur_seurat, features=rownames(cur_seurat), model.use=scale.model.use, vars.to.regress=vars.to.regress)

  } else{

    stop(paste0("Some variables specified in vars.to.regress are not found in seurat_obj@meta.data"))

  }



  # compute average expression of each gene

  cur_expr <- GetAssayData(cur_seurat, slot='data')

  expr <- t(as.matrix(cur_expr))

  averExpr <- rowSums(expr) / ncol(expr)



  # run PCA with Seurat function

  cur_pca <- Seurat::RunPCA(

    seurat_obj,

    cur_seurat,

    features = cur_genes,

    reduction.key=paste0('pca', cur_mod),

    verbose=verbose, ...

  )@reductions$pca@cell.embeddings

  )@reductions$pca

  pc <- cur_pca@cell.embeddings[,1]

  pc_loadings <- cur_pca@feature.loadings[,1]



  # correlate average expression with eigengene

  pca_cor <- cor(averExpr, pc)



  # run harmony

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



    # add this PCA as its own reduction in the seurat object

    seurat_obj@reductions$ME <- Seurat::CreateDimReducObject(

    embeddings = cur_pca,

    assay = Seurat::DefaultAssay(seurat_obj)

      embeddings = cur_pca@cell.embeddings,

      assay = Seurat::DefaultAssay(cur_seurat)

    )



  # run harmony

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

    cur_harmony <- harmony::RunHarmony(

      seurat_obj,

      group.by.vars=group.by.vars,

      reduction="ME", verbose=verbose, ...

    )@reductions$harmony@cell.embeddings

    )@reductions$harmony

    ha <- cur_harmony@cell.embeddings[,1]

    ha_loadings <- cur_pca@feature.loadings[,1]



    if(pca_cor < 0){

      cur_harmony@cell.embeddings[,1] <- -ha

      ha_loadings <- -ha_loadings

    }



    # add harmonized PCA as its own reduction in the seurat object

    seurat_obj@reductions$ME_harmony <- Seurat::CreateDimReducObject(

      embeddings = cur_harmony,

      embeddings = cur_harmony@cell.embeddings,

      assay = Seurat::DefaultAssay(seurat_obj)

    )



    seurat_obj <- SetMELoadings(

      seurat_obj,

      loadings=ha_loadings,

      harmonized=TRUE,

      wgcna_name=wgcna_name

    )



  }



  if(pca_cor < 0){

    cur_pca@cell.embeddings[,1] <- -pc

    pc_loadings <- -pc_loadings

  }



  # add this PCA as its own reduction in the seurat object

  seurat_obj@reductions$ME <- Seurat::CreateDimReducObject(

    embeddings = cur_pca@cell.embeddings,

    assay = Seurat::DefaultAssay(cur_seurat)

  )



  seurat_obj <- SetMELoadings(

    seurat_obj,

    loadings=pc_loadings,

    harmonized=FALSE,

    wgcna_name=wgcna_name

  )



  # return seurat object

  seurat_obj
@@ -601,13 +546,23 @@ ComputeModuleEigengene <- function( 
#' Computes module eigengenes for all WGCNA co-expression modules

#'

#' @param seurat_obj A Seurat object

#' @param modules table containing module / gene assignments, as in GetModules(seurat_obj).

#' @param group.by.vars groups to harmonize by

#' @param verbose logical indicating whether to print messages

#' @param vars.to.regress character vector of variables in seurat_obj@meta.data to regress when running ScaleData

#' @param scale.model.use model to scale data when running ScaleData choices are "linear", "poisson", or "negbinom"

#' @param wgcna_name name of the WGCNA experiment

#' @keywords scRNA-seq

#' @export

#' @examples

#'  ModuleEigengenes(pbmc)

ModuleEigengenes <- function(

  seurat_obj, group.by.vars=NULL,

  modules=NULL, verbose=TRUE,

  seurat_obj,

  group.by.vars=NULL,

  modules=NULL,

  vars.to.regress = NULL,

  scale.model.use = 'linear',

  verbose=TRUE,

  wgcna_name=NULL, ...

){
@@ -620,6 +575,22 @@ ModuleEigengenes <- function( 
  me_list <- list()

  harmonized_me_list <- list()



  # re-set feature loadings:

  seurat_obj <- SetMELoadings(

    seurat_obj,

    loadings=c(""),

    harmonized=FALSE,

    wgcna_name=wgcna_name

  )

  if(harmonized){

    seurat_obj <- SetMELoadings(

      seurat_obj,

      loadings=c(""),

      harmonized=TRUE,

      wgcna_name=wgcna_name

    )

  }



  # get modules from Seurat object, else use provided modules

  if(is.null(modules)){

    modules <- GetModules(seurat_obj, wgcna_name)
@@ -638,9 +609,15 @@ ModuleEigengenes <- function( 


    # compute module eigengenes for this module

    seurat_obj <- ComputeModuleEigengene(

      seurat_obj = seurat_obj, cur_mod = cur_mod,

      modules=modules, group.by.vars=group.by.vars, verbose=verbose,

      wgcna_name=wgcna_name, ...

      seurat_obj = seurat_obj,

      cur_mod = cur_mod,

      modules=modules,

      group.by.vars=group.by.vars,

      vars.to.regress = vars.to.regress,

      scale.model.use = scale.model.use,

      verbose=verbose,

      wgcna_name=wgcna_name,

      ...

    )



    # add module eigengene to ongoing list
@@ -669,7 +646,6 @@ ModuleEigengenes <- function( 


  # set module factor levels based on order

  MEs <- GetMEs(seurat_obj, harmonized, wgcna_name)

  print(colnames(MEs))

  modules$module <- factor(

    as.character(modules$module),

    levels=colnames(MEs)
@@ -821,7 +797,7 @@ ModuleConnectivity <- function( 
  seurat_obj,

  harmonized=TRUE,

  assay = NULL,

  slot = NULL,

  slot = 'data',

  group.by = NULL,

  group_name = NULL,

  wgcna_name = NULL,
@@ -837,11 +813,10 @@ ModuleConnectivity <- function( 
  genes_use <- GetWGCNAGenes(seurat_obj, wgcna_name)

  params <- GetWGCNAParams(seurat_obj, wgcna_name)



  if(is.null(assay)){assay <- params$metacell_assay}

  if(is.null(slot)){slot <- params$metacell_slot}

  if(is.null(assay)){assay <- seurat_obj@active.assay}



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

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

    cells.use <- seurat_obj@meta.data %>% subset(get(group.by) %in% group_name) %>% rownames

    MEs <- MEs[cells.use,]

  } else{

    cells.use <- colnames(seurat_obj)
@@ -856,8 +831,6 @@ ModuleConnectivity <- function( 


  datExpr <- t(as.matrix(exp_mat))



  print('running signedKME:')



  kMEs <- WGCNA::signedKME(

    datExpr,

    MEs,
@@ -880,12 +853,12 @@ ModuleConnectivity <- function( 


#' RunEnrichr

#'

#' Run Enrichr gene set enrichment tests on scWGCNA modules

#' Run Enrichr gene set enrichment tests on hdWGCNA modules

#'

#' @param seurat_obj A Seurat object

#' @param dbs character vector 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 wgcna_name The name of the hdWGCNA experiment in the seurat_obj@misc slot

#' @keywords scRNA-seq

#' @export

#' @examples
@@ -936,13 +909,13 @@ RunEnrichr <- function( 


#' OverlapModulesDEGs

#'

#' Performs Fisher's Exact Test for overlap between DEGs and scWGCNA modules.

#' Performs Fisher's Exact Test for overlap between DEGs and hdWGCNA modules.

#'

#' @param seurat_obj A Seurat object

#' @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

#' @param wgcna_name The name of the hdWGCNA experiment in the seurat_obj@misc slot

#' @keywords scRNA-seq

#' @export

#' @examples
@@ -983,10 +956,8 @@ 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,
@@ -1032,7 +1003,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 The name of the scWGCNA experiment in the seurat_obj@misc slot

#' @param wgcna_name The name of the hdWGCNA experiment in the seurat_obj@misc slot

#' @keywords scRNA-seq

#' @export

#' @examples
@@ -1077,15 +1048,16 @@ ProjectModules <- function( 
    )

  }



  # scale the dataset if needed:

  if(!scale_genes & sum(genes_use %in% rownames(GetAssayData(seurat_obj, slot='scale.data'))) == length(genes_use)){

    print("Scaling already done.")

  } else if(scale_genes){

    print("Scaling dataset...")

    seurat_obj <- Seurat::ScaleData(

      seurat_obj, features = genes_use, ...

    )

  }

  #

  # # scale the dataset if needed:

  # if(!scale_genes & sum(genes_use %in% rownames(GetAssayData(seurat_obj, slot='scale.data'))) == length(genes_use)){

  #   print("Scaling already done.")

  # } else if(scale_genes){

  #   print("Scaling dataset...")

  #   seurat_obj <- Seurat::ScaleData(

  #     seurat_obj, features = genes_use, ...

  #   )

  # }





  # project modules:
@@ -1093,6 +1065,8 @@ ProjectModules <- function( 
    seurat_obj,

    group.by.vars=group.by.vars,

    modules = modules,

    vars.to.regress = vars.to.regress,

    scale.model.use = scale.model.use,

    wgcna_name = wgcna_name_proj,

    ...

  )
@@ -1264,7 +1238,6 @@ ComputeROC <- function( 
    stop('Invalid feature selection. Valid choices: hMEs, MEs, scores.')

  )



  print('here')



  # add group column to MEs:

  MEs <- as.data.frame(MEs)
@@ -1272,7 +1245,6 @@ ComputeROC <- function( 
  MEs_p <- as.data.frame(MEs_p)

  MEs_p$group <- seurat_test@meta.data[[group.by]]



  print('here')



  # only keep common groups:

  MEs <- subset(MEs, group %in% groups_common)
@@ -1445,7 +1417,7 @@ MotifScan <- function( 
#' Overlap modules with TF target genes

#'

#' @param seurat_obj A Seurat object

#' @param wgcna_name The name of the scWGCNA experiment in the seurat_obj@misc slot

#' @param wgcna_name The name of the hdWGCNA experiment in the seurat_obj@misc slot

#' @keywords scRNA-seq

#' @export

#' @examples
@@ -1581,7 +1553,7 @@ MotifTargetScore <- function( 
#' @param seurat_obj A Seurat object

#' @param n_hubs number of hub genes to use in the UMAP computation

#' @param exclude_grey logical indicating whether to include grey module

#' @param wgcna_name The name of the scWGCNA experiment in the seurat_obj@misc slot

#' @param wgcna_name The name of the hdWGCNA experiment in the seurat_obj@misc slot

#' @param ... Additional parameters supplied to uwot::umap

#' @keywords scRNA-seq

#' @export
@@ -1696,7 +1668,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 The name of the scWGCNA experiment in the seurat_obj@misc slot

#' @param wgcna_name The name of the hdWGCNA experiment in the seurat_obj@misc slot

#' @keywords scRNA-seq

#' @export

#' @examples
@@ -1863,9 +1835,9 @@ ModuleTraitCorrelation <- function( 
#' @param n_permutations Number of permutations for the module preservation test.

#' @param parallel logical determining whether to run preservation analysis in parallel

#' @param seed random seed for the permutation analysis.

#' @param return_raw if TRUE, returns the module preservation statistics, else returns seurat_obj with the stats added to the scWGCNA experiment.

#' @param wgcna_name The name of the scWGCNA experiment in the seurat_obj@misc slot

#' @param wgcna_name_ref The name of the scWGCNA experiment in the seurat_ref@misc slot

#' @param return_raw if TRUE, returns the module preservation statistics, else returns seurat_obj with the stats added to the hdWGCNA experiment.

#' @param wgcna_name The name of the hdWGCNA experiment in the seurat_obj@misc slot

#' @param wgcna_name_ref The name of the hdWGCNA experiment in the seurat_ref@misc slot

#' @keywords scRNA-seq

#' @export

#' @examples

R/getters_and_setters.R

+94 −44

File changed.

Preview size limit exceeded, changes collapsed.

R/plotting.R

+93 −16
Original line number Diff line number Diff line
@@ -485,7 +485,7 @@ ModuleFeaturePlot<- function( 
  seurat_obj, module_names=NULL, wgcna_name = NULL,

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

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

  label_legend = FALSE, ucell = FALSE

  label_legend = FALSE, ucell = FALSE, raster=FALSE, raster_dpi=500

){



  if(is.null(wgcna_name)){wgcna_name <- seurat_obj@misc$active_wgcna}
@@ -554,11 +554,17 @@ ModuleFeaturePlot<- function( 


    # plot with ggplot

    p <- cur_plot_df %>%

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

      # ggplot(aes(x=umap1, y=umap2, color=val))

      geom_point(size=point_size, alpha=alpha) +

      ggtitle(cur_mod) + umap_theme() +

      labs(color="")

      ggplot(aes_string(x=x_name, y=y_name, color="val"))



    # rasterise?

    if(raster){

      p <- p + ggrastr::rasterise(geom_point(size=point_size, alpha=alpha), dpi=raster_dpi)

    } else{

      p <- p + geom_point(size=point_size, alpha=alpha)

    }



    # add title and theme:

    p <- p + ggtitle(cur_mod) + umap_theme() + labs(color="")



    # UCell?

    if(!ucell){
@@ -601,7 +607,7 @@ ModuleFeaturePlot<- function( 
#' @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

#' @param wgcna_name The name of the hdWGCNA experiment in the seurat_obj@misc slot

#' @keywords scRNA-seq

#' @export

#' @examples
@@ -700,7 +706,7 @@ EnrichrBarPlot <- function( 
#' @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

#' @param wgcna_name The name of the hdWGCNA experiment in the seurat_obj@misc slot

#' @keywords scRNA-seq

#' @export

#' @examples
@@ -796,7 +802,7 @@ EnrichrDotPlot <- function( 
#' @param mods Names of the modules to plot. If mods = "all", all modules are plotted.

#' @param outdir The directory where the plots will be stored.

#' @param plot_size A vector containing the width and height of the network plots.

#' @param wgcna_name The name of the scWGCNA experiment in the seurat_obj@misc slot

#' @param wgcna_name The name of the hdWGCNA experiment in the seurat_obj@misc slot

#' @keywords scRNA-seq

#' @export

#' @examples
@@ -925,7 +931,7 @@ ModuleNetworkPlot <- function( 
#' @param vertex.label.cex The font size of the gene labels

#' @param hub.vertex.size The size of the hub gene nodes

#' @param other.vertex.size The size of the other gene nodes

#' @param wgcna_name The name of the scWGCNA experiment in the seurat_obj@misc slot

#' @param wgcna_name The name of the hdWGCNA experiment in the seurat_obj@misc slot

#' @keywords scRNA-seq

#' @export

#' @examples
@@ -1085,7 +1091,8 @@ HubGeneNetworkPlot <- function( 
#' @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

#' @param keep_grey_edges logical determining whether to show edges between genes in different modules (grey edges)

#' @param wgcna_name The name of the hdWGCNA experiment in the seurat_obj@misc slot

#' @keywords scRNA-seq

#' @export

#' @examples
@@ -1099,6 +1106,7 @@ ModuleUMAPPlot <- function( 
  vertex.label.cex=0.5,

  label_genes = NULL,

  return_graph = FALSE, # this returns the igraph object instead of plotting

  keep_grey_edges = TRUE,

  wgcna_name=NULL,

  ...

){
@@ -1159,6 +1167,11 @@ ModuleUMAPPlot <- function( 
    col

  })



  # keep grey edges?

  if(!keep_grey_edges){

    edge_df <- edge_df %>% subset(color != 'grey90')

  }



  # subset edges:

  groups <- unique(edge_df$color)

  if(sample_edges){
@@ -1383,7 +1396,7 @@ OverlapBarPlot <- 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 The name of the scWGCNA experiment in the seurat_obj@misc slot

#' @param wgcna_name The name of the hdWGCNA experiment in the seurat_obj@misc slot

#' @keywords scRNA-seq

#' @export

#' @examples
@@ -1454,7 +1467,7 @@ ROCCurves <- function( 
#' Displays the top n TFs in a set of modules as a bar plot

#'

#' @param seurat_obj A Seurat object

#' @param wgcna_name The name of the scWGCNA experiment in the seurat_obj@misc slot

#' @param wgcna_name The name of the hdWGCNA experiment in the seurat_obj@misc slot

#' @keywords scRNA-seq

#' @export

#' @examples
@@ -1559,7 +1572,7 @@ MotifOverlapBarPlot <- function( 
#'

#'

#' @param seurat_obj A Seurat object

#' @param wgcna_name The name of the scWGCNA experiment in the seurat_obj@misc slot

#' @param wgcna_name The name of the hdWGCNA experiment in the seurat_obj@misc slot

#' @keywords scRNA-seq

#' @export

#' @examples
@@ -1717,7 +1730,7 @@ DoHubGeneHeatmap <- function( 
#' @param plot_labels logical determining whether to plot the module labels

#' @param label_size the size of the module labels

#' @param mod_point_size the size of the points in each plot

#' @param wgcna_name The name of the scWGCNA experiment in the seurat_obj@misc slot

#' @param wgcna_name The name of the hdWGCNA experiment in the seurat_obj@misc slot

#' @keywords scRNA-seq

#' @export

#' @examples
@@ -1835,7 +1848,7 @@ PlotModulePreservation <- function( 
#' @param

#' @param

#' @param

#' @param plot_labels logical determining whether to plot the module labels#' @param wgcna_name The name of the scWGCNA experiment in the seurat_obj@misc slot

#' @param plot_labels logical determining whether to plot the module labels#' @param wgcna_name The name of the hdWGCNA experiment in the seurat_obj@misc slot

#' @keywords scRNA-seq

#' @export

#' @examples
@@ -2203,3 +2216,67 @@ ModuleTFNetwork <- function( 




}



PlotKMEs <- function(

  seurat_obj,

  n_hubs=10,

  text_size=2,

  ncol = 5,

  plot_widths = c(3,2),

  wgcna_name = NULL

){



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



  modules <- GetModules(seurat_obj, wgcna_name) %>% subset(module != 'grey')

  mods <- levels(modules$module); mods <- mods[mods != 'grey']

  mod_colors <- modules %>% subset(module %in% mods) %>%

    select(c(module, color)) %>%

    distinct





  #get hub genes:

  hub_df <- do.call(rbind, lapply(mods, function(cur_mod){

    print(cur_mod)

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

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

    names(cur)[3] <- 'kME'

    cur <- dplyr::arrange(cur, kME)

    top_genes <- cur %>% dplyr::top_n(n_hubs, wt=kME) %>% .$gene_name

    cur$lab <- ifelse(cur$gene_name %in% top_genes, cur$gene_name, "")

    cur

  }))

  head(hub_df)



  plot_list <- lapply(mods, function(x){

    print(x)

    cur_color <- subset(mod_colors, module == x) %>% .$color

    cur_df <- subset(hub_df, module == x)

    top_genes <- cur_df %>% dplyr::top_n(n_hubs, wt=kME) %>% .$gene_name

    p <- cur_df %>% ggplot(aes(x = reorder(gene_name, kME), y = kME)) +

      geom_bar(stat='identity', width=1, color = cur_color, fill=cur_color) +

      ggtitle(x) +

      #xlab(paste0('kME_', x)) +

      theme(

        axis.ticks.x = element_blank(),

        axis.text.x = element_blank(),

        plot.title = element_text(hjust=0.5),

        axis.title.x = element_blank(),

        axis.line.x = element_blank()

      )

    p_anno <- ggplot() + annotate(

      "label",

      x = 0,

      y = 0,

      label = paste0(top_genes, collapse="\n"),

      size=text_size,

      fontface = 'italic',

      label.size=0

    ) + theme_void()

    patch <- p + p_anno + plot_layout(widths=plot_widths)

    patch

  })



  wrap_plots(plot_list, ncol=ncol)



}

CRA Git | Maintained and supported by SUSTech CRA and CCSE