adding support for consensus WGCNA (60b7f293) · Commits · github_fork / HdWGCNA

R/WGCNA_functions.R

+8 −3

Original line number	Diff line number	Diff line
		@@ -153,7 +153,7 @@ TestSoftPowers <- function(

		# add datExpr if not already added:
		if(!("datExpr" %in% names(GetActiveWGCNA(seurat_obj))) \| setDatExpr == TRUE){
		seurat_obj <- SetDatExpr(seurat_obj, use_metacells, group.by=group.by, group_name=group_name)
		seurat_obj <- SetDatExpr(seurat_obj, use_metacells=use_metacells, group.by=group.by, group_name=group_name)
		}
		# get datExpr
		datExpr <- GetDatExpr(seurat_obj)
		@@ -238,7 +238,13 @@ ConstructNetwork <- function(

		# add datExpr if not already added:
		if(!("datExpr" %in% names(GetActiveWGCNA(seurat_obj))) \| setDatExpr == TRUE){
		seurat_obj <- SetDatExpr(seurat_obj, use_metacells, group.by=group.by, group_name=group_name)
		seurat_obj <- SetDatExpr(
		seurat_obj,
		group_name = group_name,
		group.by=group.by,
		use_metacells=use_metacells,
		return_seurat=TRUE
		)
		}

		# get datExpr
		@@ -626,7 +632,6 @@ ModuleConnectivity <- function(seurat_obj, harmonized=TRUE, wgcna_name=NULL, ...
		slot=params$metacell_slot
		))[,genes_use]


		# get MEs:
		MEs <- GetMEs(seurat_obj, harmonized, wgcna_name)

R/getters_and_setters.R

+154 −13

Original line number	Diff line number	Diff line
		@@ -103,11 +103,27 @@ GetWGCNAGenes <- function(seurat_obj, wgcna_name=NULL){
		#' This function sets up the expression matrix from the metacell object.
		#'
		#' @param seurat_obj A Seurat object
		#' @param group_name A string containing a group present in the provided group.by column or in the Seurat Idents.
		#' @param use_metacells A logical determining if we use the metacells (TRUE) or the full expression matrix (FALSE)
		#' @param group.by A string containing the name of a column in the Seurat object with cell groups (clusters, cell types, etc). If NULL (default), scWGCNA uses the Seurat Idents as the group.
		#' @param multi.group.by A string containing the name of a column in the Seurat object with groups for consensus WGCNA (dataset, sample, condition, etc)
		#' @param multi_group_name A string containing the name of a group present in the multi.group.by column.
		#' @param wgcna_name A string containing the name of the WGCNA slot in seurat_obj@misc. Default = NULL which retrieves the currently active WGCNA data
		#' @keywords scRNA-seq
		#' @export
		#' @examples
		#' SetDatExpr(pbmc)
		SetDatExpr <- function(seurat_obj, use_metacells=TRUE, wgcna_name=NULL, group.by=NULL, group_name=NULL, ...){
		SetDatExpr <- function(
		seurat_obj,
		group_name,
		use_metacells=TRUE,
		group.by=NULL,
		multi.group.by = NULL,
		multi_group_name = NULL,
		return_seurat = TRUE,
		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}
		@@ -117,6 +133,10 @@ SetDatExpr <- function(seurat_obj, use_metacells=TRUE, wgcna_name=NULL, group.by
		genes_use <- GetWGCNAGenes(seurat_obj, wgcna_name)
		assay <- params$metacell_assay

		print('n_genes:')
		print(length(genes_use))
		print(head(genes_use))

		# use metacells or whole seurat object?
		if(use_metacells){
		s_obj <- GetMetacellObject(seurat_obj, wgcna_name)
		@@ -124,12 +144,27 @@ SetDatExpr <- function(seurat_obj, use_metacells=TRUE, wgcna_name=NULL, group.by
		s_obj <- seurat_obj
		}

		# columns to group by
		# get the metadata from the seurat object:
		seurat_meta <- s_obj@meta.data
		print(dim(seurat_meta))

		# columns to group by for cluster/celltype
		if(!is.null(group.by)){
		cells <- s_obj@meta.data %>% subset(get(group.by) == group_name) %>% rownames
		} else{
		cells <- colnames(s_obj)
		seurat_meta <- seurat_meta %>% subset(get(group.by) == group_name)
		}
		print(dim(seurat_meta))

		# subset further if multiExpr:
		if(!is.null(multi.group.by)){
		seurat_meta <- seurat_meta %>% subset(get(multi.group.by) == multi_group_name)
		}
		print(dim(seurat_meta))

		# get list of cells to use
		cells <- rownames(seurat_meta)
		print('cells:')
		print(head(cells))
		print(length(cells))

		# get expression data from seurat obj
		datExpr <- as.data.frame(
		@@ -143,19 +178,28 @@ SetDatExpr <- function(seurat_obj, use_metacells=TRUE, wgcna_name=NULL, group.by
		# transpose data
		datExpr <- as.data.frame(t(datExpr))

		print(dim(datExpr))

		# only get good genes:
		gene_list = GetWGCNAGenes(seurat_obj)[WGCNA::goodGenes(datExpr, ...)]
		if(is.null(multi.group.by)){
		gene_list = genes_use[WGCNA::goodGenes(datExpr, ...)]
		datExpr <- datExpr[,gene_list]
		}

		print(dim(datExpr))

		if(return_seurat){

		# update the WGCNA gene list:
		seurat_obj <- SetWGCNAGenes(seurat_obj, gene_list, wgcna_name)

		datExpr <- datExpr[,gene_list]

		# set the datExpr in the Seurat object
		seurat_obj@misc[[wgcna_name]]$datExpr <- datExpr

		# return seurat obj
		seurat_obj
		out <- seurat_obj
		} else{
		out <- datExpr
		}
		out
		}


		@@ -176,6 +220,103 @@ GetDatExpr <- function(seurat_obj, wgcna_name=NULL){

		}


		#' SetMultiExpr
		#'
		#' This function sets up the expression matrix input for consensus WGCNA based on
		#' the metacell expression matrix, or from the full expression matrix.
		#'
		#' @param seurat_obj A Seurat object
		#' @param group_name A string containing a group present in the provided group.by column or in the Seurat Idents.
		#' @param use_metacells A logical determining if we use the metacells (TRUE) or the full expression matrix (FALSE)
		#' @param group.by A string containing the name of a column in the Seurat object with cell groups (clusters, cell types, etc). If NULL (default), scWGCNA uses the Seurat Idents as the group.
		#' @param multi.group.by A string containing the name of a column in the Seurat object with groups for consensus WGCNA (dataset, sample, condition, etc)
		#' @param multi_groups A character vecrtor containing the names of
		#' @param wgcna_name A string containing the name of the WGCNA slot in seurat_obj@misc. Default = NULL which retrieves the currently active WGCNA data
		#' @keywords scRNA-seq
		#' @export
		#' @examples
		#' SetDatExpr(pbmc)
		SetMultiExpr <- function(
		seurat_obj,
		group_name,
		use_metacells=TRUE,
		group.by=NULL,
		multi.group.by = NULL,
		multi_groups = NULL,
		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 the WGCNA genes:
		gene_names <- GetWGCNAGenes(seurat_obj, wgcna_name)

		# get the different groups present if not specified by the user:
		if(is.null(multi_groups)){
		multi_groups <- unique(seurat_obj@meta.data[[multi.group.by]])
		} else{
		seurat_groups <- unique(seurat_obj@meta.data[[multi.group.by]])
		if(sum(multi_groups %in% seurat_groups) != length(multi_groups)){
		stop('Some or all groups specified in multi_groups not found in seurat_obj@meta.data[,multi.group.by]')
		}
		}

		# get the datExpr for each group
		datExpr_list <- lapply(multi_groups, function(x){
		SetDatExpr(
		seurat_obj,
		group_name = group_name,
		group.by = group.by,
		multi.group.by = multi.group.by,
		multi_group_name = x,
		return_seurat = FALSE,
		wgcna_name = wgcna_name
		) %>% as.matrix
		})
		datExpr_list

		# convert to multiExpr, get good genes:
		multiExpr <- WGCNA::list2multiData(datExpr_list)
		genes_use <- WGCNA::goodGenesMS(multiExpr)
		gene_names <- gene_names[genes_use]

		# subset the multiExpr by the good genes::
		datExpr_list <- lapply(1:length(multiExpr), function(i){
		multiExpr[[i]]$data[,genes_use]
		})
		multiExpr <- WGCNA::list2multiData(datExpr_list)
		names(multiExpr) <- multi_groups

		# update the WGCNA gene list:
		seurat_obj <- SetWGCNAGenes(seurat_obj, gene_names, wgcna_name)

		# set the multiExpr in the Seurat object
		seurat_obj@misc[[wgcna_name]]$multiExpr <- multiExpr
		seurat_obj

		}


		#' GetMultiExpr
		#'
		#' This function gets the expression matrix from the metacell object.
		#'
		#' @param seurat_obj A Seurat object
		#' @keywords scRNA-seq
		#' @export
		#' @examples
		#' GetDatExpr(pbmc)
		GetMultiExpr <- function(seurat_obj, 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]]$multiExpr

		}

		############################
		# WGCNA params
		###########################

consensus_testing.Rmd

0 → 100644

+257 −0

Original line number	Diff line number	Diff line
		```{r eval=FALSE}


		TestSoftPowersConsensus <- function(
		seurat_obj,
		use_metacells = TRUE,
		group.by=NULL, group_name=NULL,
		multi.group.by = NULL,
		multi_groups = NULL,
		setDatExpr = TRUE,
		powers=c(seq(1,10,by=1), seq(12,30, by=2)),
		make_plot=TRUE, outfile="softpower", figsize=c(7,7)
		){

		# add multiExpr if not already added:
		if(!("multiExpr" %in% names(GetActiveWGCNA(seurat_obj))) \| setDatExpr == TRUE){

		# set
		seurat_obj <- SetMultiExpr(
		seurat_obj,
		group_name = group_name,
		group.by = group.by,
		multi.group.by = multi.group.by,
		multi_groups = multi_groups
		)
		}

		multiExpr <- GetMultiExpr(seurat_obj)

		# pick soft thresh for each consensus group:
		powerTables <- list()
		for(i in 1:length(multiExpr)){

		cur_group <- names(multiExpr)[i]
		print(cur_group)

		# Call the network topology analysis function for each set in turn
		powerTable = list(
		data = WGCNA::pickSoftThreshold(
		multiExpr[[cur_group]]$data,
		powerVector=powers,
		verbose = 100,
		networkType="signed",
		corFnc="bicor"
		)[[2]]
		);
		powerTables[[cur_group]] <- powerTable

		# Plot the results:
		if(make_plot){
		pdf(paste0(outfile, '_', cur_group, '.pdf'), height=figsize[1], width=figsize[2], useDingbats=FALSE)

		colors = c("blue", "red","black")
		# Will plot these columns of the returned scale free analysis tables
		plotCols = c(2,5,6,7)
		colNames = c("Scale Free Topology Model Fit", "Mean connectivity", "Mean connectivity",
		"Max connectivity");

		# Get the minima and maxima of the plotted points
		ylim = matrix(NA, nrow = 2, ncol = 4);
		for (col in 1:length(plotCols)){
		ylim[1, col] = min(ylim[1, col], powerTable$data[, plotCols[col]], na.rm = TRUE);
		ylim[2, col] = max(ylim[2, col], powerTable$data[, plotCols[col]], na.rm = TRUE);
		}

		# Plot the quantities in the chosen columns vs. the soft thresholding power
		par(mfcol = c(2,2));
		par(mar = c(4.2, 4.2 , 2.2, 0.5))
		cex1 = 0.7;

		for (col in 1:length(plotCols)){
		plot(powerTable$data[,1], -sign(powerTable$data[,3])*powerTable$data[,2],
		xlab="Soft Threshold (power)",ylab=colNames[col],type="n", ylim = ylim[, col],
		main = colNames[col]);
		addGrid();

		if (col==1){
		text(powerTable$data[,1], -sign(powerTable$data[,3])*powerTable$data[,2],
		labels=powers,cex=cex1,col=colors[1]);
		} else
		text(powerTable$data[,1], powerTable$data[,plotCols[col]],
		labels=powers,cex=cex1,col=colors[1]);
		}
		dev.off()
		}
		}

		# set the power table in Seurat object:
		seurat_obj <- SetPowerTable(seurat_obj, powerTable$data)
		seurat_obj

		}


		#' ConstructNetwork
		#'
		#' This function constructs a co-expression network from a Seurat object
		#'
		#' @param seurat_obj A Seurat object
		#' @param soft_power
		#' @param
		#' @param
		#' @param
		#' @param
		#' @param
		#' @param
		#' @param
		#' @param
		#' @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,
		consensus = FALSE,
		multi.group.by = NULL,
		multi_groups = 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, ...
		){

		# constructing network on multiple datasets (consensus WGCNA)
		if(consensus){

		# add multiExpr if not already added:
		if(!("multiExpr" %in% names(GetActiveWGCNA(seurat_obj))) \| setDatExpr == TRUE){

		# set
		seurat_obj <- SetMultiExpr(
		seurat_obj,
		group_name = group_name,
		group.by = group.by,
		multi.group.by = multi.group.by,
		multi_groups = multi_groups
		)
		}

		multiExpr <- GetMultiExpr(seurat_obj)
		checkSets(multiExpr) # check data size

		# constructing network on a single dataset
		} else{

		# add datExpr if not already added:
		if(!("datExpr" %in% names(GetActiveWGCNA(seurat_obj))) \| setDatExpr == TRUE){
		print('in here')
		seurat_obj <- SetDatExpr(
		seurat_obj,
		group_name = group_name,
		group.by=group.by,
		use_metacells=use_metacells,
		return_seurat=TRUE
		)
		print('out here')

		}

		# get datExpr from seurat object
		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
		}


		# make output dir for the TOM
		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

		}


		```

scWGCNA_testing.Rmd

+100 −130

Original line number	Diff line number	Diff line
		@@ -30,7 +30,7 @@ theme_set(theme_cowplot())
		setwd("/dfs3b/swaruplab/smorabit/collab/woodlab/cocaine_mouse_2021/Nurr2c_vs_GFP/scWGCNA")

		# source all of the scWGCNA scripts:
		scripts <- dir("bin/")
		scripts <- dir("/dfs3b/swaruplab/smorabit/analysis/scWGCNA/bin/")
		scripts <- scripts[scripts != 'scWGCNA.R']
		for(script in scripts){
		source(paste0("bin/", script))
		@@ -402,155 +402,125 @@ dev.off()

		```

		ROC Curve comparing original modules & projected modules in a shared cell space

		Testing consensus WGCNA support:

		Need to make a function "SetMultiExpr" for consensus

		Also will pretty much do a whole workflow to check if the consensus works for
		the downstream tasks.

		```{r eval=FALSE}

		group.by <- 'annotation'
		groups <- as.character(unique(cur_seurat@meta.data[[group.by]]))
		# re-load seurat obj
		cur_seurat <- readRDS(file='data/test_wgcna_seurat.rds')
		cur_seurat <- SetModules(cur_seurat, NULL)

		# get Modules
		modules <- GetModules(cur_seurat)
		mods <- levels(modules$module)
		mods <- mods[mods != 'grey']
		######################################################
		# MultiExpr
		######################################################

		# get MEs:
		MEs <- GetMEs(cur_seurat, harmonized=TRUE) %>% mutate(group = cur_seurat@meta.data[[group.by]])
		MEs_p <- as.data.frame(GetMEs(seurat_proj, harmonized=TRUE)) %>% mutate(group = seurat_proj@meta.data[[group.by]])

		# 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))
		groups <- avg_MEs$group

		# scale each column between 0 & 1
		avg_MEs <- avg_MEs %>% summarise(across(!!mods, scale01))
		avg_MEs_p <- avg_MEs_p %>% summarise(across(!!mods, scale01))

		# convert to binary labels:
		thresh <- 0.75
		labels <- avg_MEs %>% purrr::map(~ifelse(. >= thresh, TRUE, FALSE))
		labels <- as.data.frame(do.call(cbind, labels));
		rownames(labels) <- as.character(groups)

		predictions <- avg_MEs_p %>% purrr::map(~ifelse(. >= thresh, TRUE, FALSE))
		predictions <- as.data.frame(do.call(cbind, predictions));
		rownames(predictions) <- as.character(groups)

		# loop over modules to compute ROC curves:
		plot_df <- data.frame()
		conf_df <- data.frame()
		auc_list <- list()
		mod_colors <- list()
		for(cur_mod in mods){
		print(cur_mod)
		cur_color <- modules %>% subset(module == cur_mod) %>% .$color %>% unique
		mod_colors[[cur_mod]] <- cur_color

		# compute ROC
		rocobj <- pROC::roc(labels[,cur_mod], avg_MEs_p[[cur_mod]])
		auc_list[[cur_mod]] <- as.numeric(rocobj$auc)
		# confidence intervals:
		# sens_ci <- ci.se(rocobj)

		cur_df <- data.frame(
		specificity = 1-rocobj$specificities,
		sensitivity = rocobj$sensitivities,
		#auc = rocobj$auc,
		module = cur_mod,
		color = cur_color,
		auc = as.numeric(rocobj$auc)
		# ci_lo = sens_ci[,1],
		# ci_med = sens_ci[,2],
		# ci_hi = sens_ci[,3]
		# just run SetMultiExpr by itself
		test_seurat <- SetMultiExpr(
		cur_seurat,
		group_name = "MHb-Neuron",
		group.by = "cell_type",
		multi.group.by = "Sample",
		multi_groups = NULL
		)
		test_mExpr <- GetMultiExpr(test_seurat)


		######################################################
		# test consensus
		######################################################

		cur_seurat <- TestSoftPowersConsensus(
		cur_seurat,
		group.by='cell_type', group_name="MHb-Neuron",
		multi.group.by = 'Sample',
		multi_groups = NULL,
		setDatExpr=TRUE
		)
		plot_df <- rbind(plot_df, cur_df)

		cur_conf <- as.data.frame(ci.se(rocobj))
		cur_conf$sensitivity <- 1-as.numeric(rownames(cur_conf))
		cur_conf$module <- cur_mod
		cur_conf$color <- cur_color
		conf_df <- rbind(conf_df, cur_conf)
		cur_seurat <- ConstructNetwork(
		cur_seurat,
		soft_power=c(5,8), # soft power can be a single number of a vector with a value for each datExpr in multiExpr
		group.by='cell_type', group_name="MHb-Neuron",
		multi.group.by = 'Sample',
		multi_groups = NULL,
		setDatExpr=FALSE,
		consensus=TRUE
		)

		}
		colnames(conf_df)[1:3] <- c('lo', 'mid', 'hi')
		# plot consensus dendro:
		pdf("figures/test_dendro_consensus.pdf",height=5, width=8)
		PlotDendrogram(cur_seurat, main='Test Consensus')
		dev.off()

		# compute all MEs in the full single-cell dataset
		cur_seurat <- ModuleEigengenes(
		cur_seurat,
		group.by.vars="Sample"
		)

		plot_list <- ModuleFeaturePlot(cur_seurat, order=TRUE)
		pdf("figures/test_consensus_hMEs.pdf",height=9, width=15)
		wrap_plots(plot_list, ncol=5)
		dev.off()

		library(pROC)
		rocobj <- pROC::roc(labels$red, avg_MEs_p$red)
		rocobj$auc
		# compute module connectivity:
		cur_seurat <- ModuleConnectivity(cur_seurat)

		sens_ci <- ci.se(rocobj)


		p <- ggroc(rocobj)
		library(igraph)

		pdf(paste0(fig_dir, 'test_ROC_blue.pdf'), width=6, height=6)
		p
		dev.off()
		#
		#
		# # confidence intervals:
		# sens_ci <- ci.se(rocobj)
		#
		#
		# plot_df <- data.frame(
		# specificity = rocobj$specificities,
		# sensitivity = 1-rocobj$sensitivities,
		# #auc = rocobj$auc,
		# module = "blue",
		# color = "blue"
		# # ci_lo = sens_ci[,1],
		# # ci_med = sens_ci[,2],
		# # ci_hi = sens_ci[,3]
		#
		# )

		plot_df <- plot_df %>% group_by(module) %>% arrange(sensitivity)
		conf_df <- conf_df %>% group_by(module) %>% arrange(sensitivity)
		auc_df <- distinct(plot_df[,c('module', 'auc')])

		p <- plot_df %>% ggplot(
		aes(x=specificity, y=sensitivity, color=module, fill=module),
		) +
		geom_line() +
		geom_ribbon(
		data=conf_df,
		aes(x = sensitivity, ymin=lo, ymax=hi, fill=module),
		inherit.aes=FALSE, alpha=0.4
		) +
		scale_color_manual(values = unlist(mod_colors)) +
		scale_fill_manual(values = unlist(mod_colors)) +
		scale_x_continuous(breaks = c(0, 0.5, 1), labels=c("0", "0.5", "1")) +
		scale_y_continuous(breaks = c(0, 0.5, 1), labels=c("0", "0.5", "1")) +
		xlab("1 - Specificity (FPR)") + ylab("Sensitivity (TPR)") +
		geom_text(
		data = auc_df,
		aes(color=module),
		x=0.75, y=0.1, label=paste0("AUC: ", format(auc_df$auc, digits=2)),
		inherit.aes=FALSE, size=4, color='black'
		# network plots:
		ModuleNetworkPlot(
		cur_seurat,
		outdir = paste0(fig_dir, 'test_consensus_networks/')
		)
		#annotate("text", x=0.75, y=0.1, label=paste0("AUC: ", plot_df$auc))
		# theme(
		# axis.text = element_blank(),
		# axis.ticks = element_blank()
		# )



		pdf(paste0(fig_dir, 'test_ROC.pdf'), width=9, height=6)
		p + facet_wrap(~module, ncol=5) + NoLegend()

		######################################################
		# see if old function still works
		######################################################

		# re-load seurat obj
		cur_seurat <- readRDS(file='data/test_wgcna_seurat.rds')

		# see if old function still works
		cur_seurat <- ConstructNetwork(
		cur_seurat, soft_power=8,
		group.by='cell_type', group_name="MHb-Neuron",
		setDatExpr=TRUE
		)

		# plot dendro:
		pdf("figures/test_dendro.pdf",height=5, width=8)
		PlotDendrogram(cur_seurat, main='Test')
		dev.off()


		unlist(auc_list)
		# network plots:
		ModuleNetworkPlot(
		cur_seurat,
		outdir = paste0(fig_dir, 'test_networks/')
		)






		```


		Project onto human AD dataset:

		Project modules from mouse onto human AD dataset:

		```{r eval=FALSE}

test_ROC.Rmd

+0 −13

Original line number	Diff line number	Diff line
		@@ -237,19 +237,6 @@ dev.off()

		```

		Analysis for Fig 1 of the scWGCNA paper:

		1. Integrate Zhou et al AD dataset & Morabito/Miyoshi et al AD dataset
		2. Train/Test scWGCNA on Zhou with 80/20 split, use Morabito et al as validation set

		```{r eval=FALSE}

		# TODO


		```


		Project onto human AD dataset:

		```{r eval=FALSE}

Admin message