DME

Package: hdWGCNA

Package: hdWGCNA

Title: hdWGCNA

Title: hdWGCNA

Version: 0.1.1.9015

Version: 0.1.2.0000

Authors@R: c(

Authors@R: c(

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

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

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

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

export(DoHubGeneHeatmap)

export(DoHubGeneHeatmap)

export(EnrichrBarPlot)

export(EnrichrBarPlot)

export(EnrichrDotPlot)

export(EnrichrDotPlot)

export(FindAllDMEs)

export(FindDMEs)

export(GetActiveWGCNA)

export(GetActiveWGCNA)

export(GetAvgModuleExpr)

export(GetAvgModuleExpr)

export(GetDatExpr)

export(GetDatExpr)

export(OverlapDotPlot)

export(OverlapDotPlot)

export(OverlapModulesDEGs)

export(OverlapModulesDEGs)

export(OverlapModulesMotifs)

export(OverlapModulesMotifs)

export(PlotDMEsVolcano)

export(PlotDendrogram)

export(PlotDendrogram)

export(PlotKMEs)

export(PlotKMEs)

export(PlotModulePreservation)

export(PlotModulePreservation)

# hdWGCNA 0.1.2.0000 (2022-09-08)

## Added

- Differential Module Eigengene (DME) tutorial

- FindDMEs function

- FindAllDMEs function

- PlotDMEsVolcano function

## Changes

- None

# hdWGCNA 0.1.1.9015 (2022-09-06)

# hdWGCNA 0.1.1.9015 (2022-09-06)

## Added

## Added

- None

- None

## Changes

## Changes

- Bugfix in `ModuleConnectivity` that caused kMEs to be out of order.

- Bugfix in `ModuleConnectivity` that caused kMEs to be out of order.

# hdWGCNA 0.1.1.9013 (2022-08-25)

# hdWGCNA 0.1.1.9013 (2022-08-25)

## Added

## Added

- Tutorial for using SCTransform normalized data.

- Tutorial for using SCTransform normalized data.

  seurat_obj

  seurat_obj

}

}

#' FindAllDMEs

#'

#' Modification of the Seurat function FindMarkers used to perform iterative one-versus-all differential module eigengene testing given a group of cells (ie clusters, cell types, etc).

#'

#' @param seurat_obj A Seurat object

#' @param group.by column in seurat_obj@meta.data containing cell grouping information

#' @param harmonized logical determining whether or not to use harmonized MEs

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

#' @param add_missing logical determining whether or not to add missing modules back into the resulting dataframe with NA values.

#' @param ... Additional parameters for the FindMarkers function

#' @keywords scRNA-seq

#' @export

#' @return A dataframe contaning differential ME results

#' @examples

#' FindAllDMEs

FindAllDMEs <- function(

  seurat_obj,

  group.by,

  harmonized=TRUE,

  wgcna_name=NULL,

  add_missing=FALSE,

  test.use='wilcox',

  only.pos=FALSE,

  logfc.threshold = 0,

  min.pct=0,

  verbose=FALSE,

  pseudocount.use=0,

  ...

){

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

  # get list of groups

  groups <- as.character(unique(seurat_obj@meta.data[[group.by]]))

  # get module eigengenes, remove grey module

  MEs <- GetMEs(seurat_obj, harmonized, wgcna_name)

  MEs <- MEs[, colnames(MEs) != 'grey']

  # set all negative values to zero

  MEs[MEs < 0] <- 0

  # transpose MEs so cells are columns and rows are MEs

  MEs <- t(MEs)

  # create new assay for MEs:

  ME_assay <- Seurat::CreateAssayObject(MEs)

  DMEs_list <- list()

  for(cur_group in groups){

    print(cur_group)

    barcodes1 <- colnames(seurat_obj)[seurat_obj@meta.data[[group.by]] == cur_group]

    barcodes2 <- colnames(seurat_obj)[seurat_obj@meta.data[[group.by]] != cur_group]

    # run differential test on the ME assay

    DMEs <- FindMarkers(

      ME_assay,

      cells.1 = barcodes1,

      cells.2 = barcodes2,

      slot='counts',

      test.use=test.use,

      only.pos=only.pos,

      logfc.threshold=logfc.threshold,

      min.pct=min.pct,

      verbose=verbose,

      pseudocount.use=pseudocount.use,

      ...

    )

    DMEs$module <- rownames(DMEs)

    DMEs$group <- cur_group

    # add missing modules if specified

    if(add_missing){

      missing_mods <- rownames(MEs)[!(rownames(MEs) %in% DMEs$module)]

      for(cur_mod in missing_mods){

        DMEs[cur_mod,] <- NA

        DMEs[cur_mod,'module'] <- cur_mod

      }

    }

    # add to ongoing list:

    DMEs_list[[cur_group]] <- DMEs

  }

  DMEs <- do.call(rbind, DMEs_list)

  DMEs

}

#' FindDMEs

#'

#' Modification of the Seurat function FindMarkers used to perform differential module eigengene testing between two sets of cell barcodes.

#'

#' @param seurat_obj A Seurat object

#' @param barcodes1 character vector containing cell barcodes for the first group to test. Positive fold-change means up-regulated in this group.

#' @param barcodes2 character vector containing cell barcodes for the second group to test. Negative fold-change means up-regulated in this group.

#' @param harmonized logical determining whether or not to use harmonized MEs

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

#' @param add_missing logical determining whether or not to add missing modules back into the resulting dataframe with NA values.

#' @param ... Additional parameters for the FindMarkers function

#' @keywords scRNA-seq

#' @export

#' @return A dataframe contaning differential ME results

#' @examples

#' FindDMEs

FindDMEs <- function(

  seurat_obj,

  barcodes1,

  barcodes2,

  harmonized=TRUE,

  wgcna_name=NULL,

  add_missing=FALSE,

  test.use='wilcox',

  only.pos=FALSE,

  logfc.threshold = 0,

  min.pct=0,

  verbose=FALSE,

  pseudocount.use=0,

  ...

){

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

  # ensure that selected barcodes are in the seurat obj

  if(!(all(barcodes1 %in% colnames(seurat_obj)))){

    stop('Some barcodes in barcodes1 not found in colnames(seurat_obj).')

  }

  if(!(all(barcodes2 %in% colnames(seurat_obj)))){

    stop('Some barcodes in barcodes2 not found in colnames(seurat_obj).')

  }

  # check for overlap in the two groups of bacrodes:

  if(length(intersect(barcodes1, barcodes2)) > 0){

    stop('Some barcodes overlap in barcodes1 and barcodes2')

  }

  # get module eigengenes, remove grey module

  MEs <- GetMEs(seurat_obj, harmonized, wgcna_name)

  MEs <- MEs[, colnames(MEs) != 'grey']

  # set all negative values to zero

  MEs[MEs < 0] <- 0

  # transpose MEs so cells are columns and rows are MEs

  MEs <- t(MEs)

  # create new assay for MEs:

  ME_assay <- Seurat::CreateAssayObject(MEs)

  # run differential test on the ME assay

  DMEs <- FindMarkers(

    ME_assay,

    cells.1 = barcodes1,

    cells.2 = barcodes2,

    slot='counts',

    test.use=test.use,

    only.pos=only.pos,

    logfc.threshold=logfc.threshold,

    min.pct=min.pct,

    verbose=verbose,

    pseudocount.use=pseudocount.use,

    ...

  )

  DMEs$module <- rownames(DMEs)

  # add missing modules if specified

  if(add_missing){

    missing_mods <- rownames(MEs)[!(rownames(MEs) %in% DMEs$module)]

    for(cur_mod in missing_mods){

      DMEs[cur_mod,] <- NA

      DMEs[cur_mod,'module'] <- cur_mod

    }

  }

  DMEs

}

}

}

#' PlotModulePreservatgion

#' PlotModulePreservation

#'

#'

#' Plotting function for Module Preservation statistics

#' Plotting function for Module Preservation statistics

#'

#'

    # don't include grey & gold:

    # don't include grey & gold:

    plot_df <- plot_df %>% subset(!(module %in% c('grey', 'gold')))

    plot_df <- plot_df %>% subset(!(module %in% c('grey', 'gold')))

    if(grepl("Rank", statistic)){

    if(grepl("Rank", statistic)){

      cur_p <-  plot_df %>% ggplot(aes(x=size, y=value, fill=module, color=module)) +

      cur_p <-  plot_df %>% ggplot(aes(x=size, y=value, fill=module, color=module)) +

        geom_point(size=mod_point_size, pch=21, color='black') +

        geom_point(size=mod_point_size, pch=21, color='black') +

    if(plot_labels){

    if(plot_labels){

      cur_p <- cur_p + geom_text_repel(label = plot_df$module, size=label_size)

      cur_p <- cur_p + geom_text_repel(label = plot_df$module, size=label_size, max.overlaps=Inf, color='black')

    }

    }

    plot_list[[statistic]] <- cur_p

    plot_list[[statistic]] <- cur_p

  wrap_plots(plot_list, ncol=ncol)

  wrap_plots(plot_list, ncol=ncol)

}

}

#' PlotDMEsVolcano

#'

#' Plotting function for the results of FindDMEs and FindAllDMEs

#'

#' @param DMEs dataframe output from FindDMEs or FindAllDMEs

#' @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 show_cutoff logical determining whether to plot the significance cutoff. Should set this to FALSE if using facet_wrap.

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

#' @keywords scRNA-seq

#' @export

#' @return A ggplot object

#' @examples

#' PlotDMEsVolcano

PlotDMEsVolcano <- function(

  DMEs,

  plot_labels=TRUE,

  mod_point_size=4,

  label_size=4,

  show_cutoff = TRUE,

  wgcna_name=NULL

){

  # remove NAs:

  DMEs <- na.omit(DMEs)

  # lowest non-zero value

  lowest <- DMEs %>% subset(p_val_adj != 0) %>% top_n(-1, wt=p_val_adj) %>% .$p_val_adj

  DMEs$p_val_adj <- ifelse(DMEs$p_val_adj == 0, lowest, DMEs$p_val_adj)

  # fix infinite fold change

  max_fc <- max(abs(DMEs$avg_log2FC))

  max_fc <- DMEs %>% subset(abs(avg_log2FC) != Inf) %>% .$avg_log2FC %>% max

  DMEs$avg_log2FC <- ifelse(DMEs$avg_log2FC == -Inf, -1*max_fc, DMEs$avg_log2FC)

  DMEs$avg_log2FC <- ifelse(DMEs$avg_log2FC == Inf, max_fc, DMEs$avg_log2FC)

  # get modules and module colors

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

  module_colors <- modules %>% dplyr::select(c(module, color)) %>% distinct

  # module names

  mods <- levels(modules$module)

  mods <- mods[mods %in% DMEs$module]

  mod_colors <- module_colors$color; names(mod_colors) <- as.character(module_colors$module)

  # annotate modules with significant DME

  DMEs$anno <- ifelse(DMEs$p_val_adj < 0.05, DMEs$module, '')

  # set x-axis limit

  xmax <- max_fc

  # plot basics

  p <- DMEs %>%

    ggplot(aes(x=avg_log2FC, y=-log10(p_val_adj), fill=module, color=module))

  if(show_cutoff){

    p <- p +

      geom_vline(xintercept=0, linetype='dashed', color='grey75', alpha=0.8) +

      geom_rect(

        data=DMEs[1,],

        aes(xmin=-Inf, xmax=Inf, ymin=-Inf, ymax=-log10(0.05)), fill='grey75', alpha=0.8, color=NA)

  }

  # add points:

  p <- p + geom_point(size=mod_point_size, pch=21, color='black')

  # label points?

  if(plot_labels){

    p <- p + geom_text_repel(aes(label=anno), color='black', min.segment.length=0, max.overlaps=Inf, size=label_size)

  }

  p <- p +

    scale_fill_manual(values=mod_colors) +

    scale_color_manual(values=mod_colors) +

     xlim((-1*xmax)-0.1, xmax+0.1) +

    xlab(bquote("Average log"[2]~"(Fold Change)")) +

    ylab(bquote("-log"[10]~"(Adj. P-value)")) +

    theme(

     panel.border = element_rect(color='black', fill=NA, size=1),

     panel.grid.major = element_blank(),

     axis.line = element_blank(),

     plot.title = element_text(hjust = 0.5),

     legend.position='bottom'

   ) + NoLegend()

   p + NoLegend()

}