Commit f1f6b1af authored by smorabit's avatar smorabit
Browse files

fixed some bugs

parent d2766d6e

R/WGCNA_functions.R

+31 −17
Original line number Diff line number Diff line
@@ -819,12 +819,13 @@ ProjectModules <- function( 
    modules <- TransferModuleGenome(modules, gene_mapping, genome1_col, genome2_col)

  }



  print('here')



  # get genes that overlap between WGCNA genes & seurat_obj genes:

  gene_names <- modules$gene_name

  genes_use <- intersect(gene_names, rownames(seurat_obj))



  print('n genes:')

  print(length(genes_use))



  # subset modules by genes in this seurat object:

  modules <- modules %>% subset(gene_name %in% genes_use)
@@ -924,7 +925,8 @@ TransferModuleGenome <- function( 
#' @examples

#' ComputeROC

ComputeROC <- function(

  seurat_obj, group.by=NULL,

  seurat_obj,

  group.by=NULL, # this needs to be a factor!!!

  split_col=NULL, # needs to be a logical!!!

  features = 'hMEs',

  seurat_test=NULL,
@@ -953,7 +955,7 @@ ComputeROC <- function( 
  }



  # get names of different cell groupings

  groups <- as.character(unique(Idents(seurat_obj)))

  groups <- levels(seurat_obj@meta.data[[group.by]])

  groups <- groups[order(groups)]



  # split seurat object into training & testing if the testing is not provided
@@ -990,35 +992,43 @@ ComputeROC <- function( 


  }





  # get names of different cell groupings in test

  groups_test <- as.character(unique(Idents(seurat_test)))

  groups_test <- levels(seurat_test@meta.data[[group.by]])

  groups_test <- groups_test[order(groups_test)]



  # groups that are in common (some might not have been predicted):

  groups_common <- intersect(

    groups,

    as.character(unique(seurat_test@meta.data[[group.by]]))

  )



  # check if groups are equal:

  if(sum(groups == groups_test) != length(groups)){

    stop("Different groups present in train & test data. Idents likely do not match.")

  }



  # project modules for test

  # project modules for test if they haven't already been projected:

  if(is.null(GetModules(seurat_test, wgcna_name=wgcna_name_test))){

    wgcna_name_test <- "ROC"

    seurat_test <- ProjectModules(

      seurat_test,

      seurat_ref = seurat_obj,

      group.by.vars=harmony_group_vars,

    wgcna_name_proj="ROC",

      wgcna_name_proj=wgcna_name_test,

      scale_genes=scale_genes, verbose=verbose

    )

  }



  # get MEs from seurat object

  if(features == 'hMEs'){

    MEs <- GetMEs(seurat_train, TRUE, wgcna_name_train)

    MEs_p <- GetMEs(seurat_test, TRUE, "ROC")

    MEs_p <- GetMEs(seurat_test, TRUE, wgcna_name_test)

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

    MEs <- GetMEs(seurat_train, FALSE, wgcna_name_train)

    MEs_p <- GetMEs(seurat_test, FALSE, "ROC")

    MEs_p <- GetMEs(seurat_test, FALSE, wgcna_name_test)

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

    MEs <- GetModuleScores(seurat_train, wgcna_name_train)

    MEs_p <- GetModuleScores(seurat_test, "ROC")

    MEs_p <- GetModuleScores(seurat_test, wgcna_name_test)

    stop("Haven't implemented this one yet >.<")

  } else(

    stop('Invalid feature selection. Valid choices: hMEs, MEs, scores.')
@@ -1028,6 +1038,10 @@ ComputeROC <- function( 
  MEs <- as.data.frame(MEs) %>% mutate(group = seurat_train@meta.data[[group.by]])

  MEs_p <- as.data.frame(MEs_p) %>% mutate(group = seurat_test@meta.data[[group.by]])



  # only keep common groups:

  MEs <- subset(MEs, group %in% groups_common)

  MEs_p <- subset(MEs_p, group %in% groups_common)



  # compute average MEs in each group:

  avg_MEs <- MEs %>% group_by(group) %>% summarise(across(!!mods, mean))

  avg_MEs_p <- MEs_p %>% group_by(group) %>% summarise(across(!!mods, mean))

R/plotting.R

+25 −15
Original line number Diff line number Diff line
@@ -304,7 +304,7 @@ ModuleCorrNetwork <- function( 
ModuleFeaturePlot<- function(

  seurat_obj, module_names=NULL, wgcna_name = NULL,

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

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

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

  label_legend = FALSE, ucell = FALSE

){
@@ -352,7 +352,6 @@ ModuleFeaturePlot<- function( 


    # reset the range of the plot:

    plot_range <- plot_df[,cur_mod] %>% range

    print(plot_range)

    if(restrict_range){

      if(abs(plot_range[1]) > abs(plot_range[2])){

        plot_range[1] <- -1*plot_range[2]
@@ -363,16 +362,16 @@ ModuleFeaturePlot<- function( 
      plot_df[,cur_mod] <- ifelse(plot_df[,cur_mod] < plot_range[1], plot_range[1], plot_df[,cur_mod])

    }



    # order points:

    if(order == TRUE){

      plot_df <- plot_df %>% dplyr::arrange(!!cur_mod)

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

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

    }



    cur_plot_df <- plot_df[,c(colnames(umap), cur_mod)]

    colnames(cur_plot_df)[3] <- "val"



    # order points:

    if(order_points == TRUE){

      cur_plot_df <- cur_plot_df %>% dplyr::arrange(val)

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

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

    }



    # plot with ggplot

    p <- cur_plot_df %>%

      ggplot(aes_string(x=x_name, y=y_name, color="val")) +
@@ -618,6 +617,7 @@ EnrichrDotPlot <- function( 
ModuleNetworkPlot <- function(

  seurat_obj, mods="all", outdir="ModuleNetworks",

  plot_size = c(6,6), wgcna_name=NULL,

  label_center = FALSE, # only label the genes in the middle?

  edge.alpha=0.25, vertex.label.cex=1, vertex.size=6, ...

){
@@ -678,6 +678,9 @@ ModuleNetworkPlot <- function( 
    reducedTOM <- matrix(0,nrow(reducedTOM),ncol(reducedTOM));

    reducedTOM[connections2keep] <- 1;



    # only label the top 10 genes?

    if(label_center){cur$gene_name[11:25] <- ''}



    # top 10 as center

    gA <- graph.adjacency(as.matrix(reducedTOM[1:10,1:10]),mode="undirected",weighted=TRUE,diag=FALSE)

    gB <- graph.adjacency(as.matrix(reducedTOM[11:n_genes,11:n_genes]),mode="undirected",weighted=TRUE,diag=FALSE)
@@ -986,7 +989,9 @@ OverlapBarPlot <- function( 
#' ROCCurves

ROCCurves <- function(

  seurat_obj,

  roc_df=NULL, conf_df=NULL, wgcna_name=NULL

  roc_df=NULL,

  conf_df=NULL,

  wgcna_name=NULL

){



  if(is.null(wgcna_name)){wgcna_name <- seurat_obj@misc$active_wgcna}
@@ -1176,8 +1181,11 @@ DoHubGeneHeatmap <- function( 
    seurat_obj$temp_ident <- Idents(seurat_obj)

  }



  # drop if there are missing levels:

  seurat_plot@meta.data[[group.by]] <- droplevels(seurat_plot@meta.data[[group.by]])



  # get modules

  modules <- GetModules(seurat_obj)

  modules <- GetModules(seurat_obj, wgcna_name)

  modules <- modules %>% subset(module != 'grey') %>% mutate(module = droplevels(module))

  mods <- levels(modules$module)
@@ -1256,6 +1264,11 @@ DoHubGeneHeatmap <- function( 


  }



  # useful for ratios on colorbars

  n_total_cells <- ncol(seurat_plot)

  width_cbar <- n_total_cells / 50





  # module colorbar

  mod_colors$value <- n_hubs

  mod_colors$dummy <- 'colorbar'
@@ -1269,9 +1282,6 @@ DoHubGeneHeatmap <- function( 
  }

  p_cbar <- wrap_plots(cbar_list, ncol=1)



  n_total_cells <- ncol(seurat_plot)

  width_cbar <- n_total_cells / 50



  if(combine){

    out <- wrap_plots(plot_list, ncol=1) +plot_layout(guides='collect')

    out <- (p_cbar | out) + plot_layout(widths=c(width_cbar, n_total_cells))

scWGCNA_testing.Rmd

+2 −1
Original line number Diff line number Diff line
@@ -899,6 +899,7 @@ pdf(paste0(fig_dir, 'test_hubgene_heatmap.pdf'), width=10, height=10, useDingbat 
p

dev.off()





```

CRA Git | Maintained and supported by SUSTech CRA and CCSE