added correlation functions (83d24a9f) · Commits · github_fork / HdWGCNA

R/WGCNA_functions.R

+16 −3

Original line number	Diff line number	Diff line
		@@ -348,6 +348,9 @@ ModuleEigengenes <- function(seurat_obj, group.by.vars=NULL, verbose=TRUE, wgcna
		# set as active assay if wgcna_name is not given
		if(is.null(wgcna_name)){wgcna_name <- seurat_obj@misc$active_wgcna}

		# are we going to run Harmony?
		harmonized = !is.null(group.by.vars)

		me_list <- list()
		harmonized_me_list <- list()
		modules <- GetModules(seurat_obj, wgcna_name)
		@@ -369,28 +372,38 @@ ModuleEigengenes <- function(seurat_obj, group.by.vars=NULL, verbose=TRUE, wgcna
		me_list[[cur_mod]] <- cur_me

		# run harmony
		if(!is.null(group.by.vars)){
		if(harmonized){
		# add module eigengene to ongoing list
		cur_harmonized_me <- seurat_obj@reductions$ME_harmony@cell.embeddings[,1]
		harmonized_me_list[[cur_mod]] <- cur_harmonized_me
		}
		}

		# merge module eigengene lists into a dataframe, add to Seurat obj
		# merge module eigengene lists into a dataframe, order modules, add to Seurat obj
		me_df <- do.call(cbind, me_list)
		me_df <- WGCNA::orderMEs(me_df)
		seurat_obj <- SetMEs(seurat_obj, me_df, harmonized=FALSE, wgcna_name)

		# merge harmonized module eigengene lists into a dataframe, add to Seurat obj
		if(!is.null(group.by.vars)){
		hme_df <- do.call(cbind, harmonized_me_list)
		hme_df <- WGCNA::orderMEs(hme_df)
		seurat_obj <- SetMEs(seurat_obj, hme_df, harmonized=TRUE, wgcna_name)
		}

		# set module factor levels based on order
		MEs <- GetMEs(cur_seurat, harmonized, wgcna_name)
		modules$module <- factor(
		as.character(modules$module),
		levels=colnames(MEs)
		)
		seurat_obj <- SetModules(seurat_obj, modules, wgcna_name)

		# remove temp dim reductions by setting to NULL
		seurat_obj@reductions$ME <- NULL
		seurat_obj@reductions$ME_harmony <- NULL

		# return seurat object f
		# return seurat object
		seurat_obj
		}

R/plotting.R

+216 −11

Original line number	Diff line number	Diff line
		@@ -9,14 +9,21 @@
		#' @examples
		#' PlotDendrogram
		PlotDendrogram <- function(
		seurat_obj, groupLabels="Module colors",
		seurat_obj, groupLabels="Module colors", wgcna_name=NULL,
		dendroLabels = FALSE, hang = 0.03, addGuide = TRUE, guideHang = 0.05,
		main = "", ...
		){

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

		# get WGCNA network and module data
		net <- GetNetworkData(seurat_obj, wgcna_name)
		modules <- GetModules(seurat_obj, wgcna_name)

		# plot dendrogram
		WGCNA::plotDendroAndColors(
		seurat_obj@misc[[seurat_obj@misc$active_wgcna]]$wgcna_net$dendrograms[[1]],
		as.character(seurat_obj@misc[[seurat_obj@misc$active_wgcna]]$wgcna_net$colors),
		net$dendrograms[[1]],
		as.character(modules$color),
		groupLabels=groupLabels,
		dendroLabels = dendroLabels,
		hang = hang,
		@@ -27,6 +34,196 @@ PlotDendrogram <- function(
		)
		}

		#' MECorrelogram
		#'
		#' Plot Module Eigengene correlogram
		#'
		#' @param seurat_obj A Seurat object
		#' @keywords scRNA-seq
		#' @export
		#' @examples
		#' MECorrelogram
		MECorrelogram <- function(
		seurat_obj, wgcna_name = NULL, exclude_grey=TRUE,
		type = 'upper', order='original', method='ellipse',
		tl.col = 'black', tl.srt=45,
		sig.level = c(0.0001, 0.001, 0.01, 0.05),
		insig='label_sig', pch.cex=0.7, col=NULL, ncolors=200, ...){

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

		# get MEs and convert to matrix
		MEs <- as.matrix(GetMEs(seurat_obj, wgcna_name))

		# exclude grey
		if(exclude_grey){
		MEs <- MEs[,colnames(MEs) != 'grey']
		}

		# setup color scheme:
		if(is.null(col)){
		colfunc <- grDevices::colorRampPalette(c('seagreen', 'white', 'darkorchid1'))
		col = colfunc(ncolors)
		}

		# perform correlation
		res <- Hmisc::rcorr(MEs)
		resP <- corrplot::cor.mtest(MEs, conf.level=0.95)$p

		# plot correlogram
		corrplot::corrplot(
		res$r, p.mat = resP,
		type=type, order=order,
		method=method, tl.col=tl.col,
		tl.srt=tl.srt, sig.level=sig.level,
		insig=insig, pch.cex=pch.cex,
		col = col, ...
		)

		}


		#' ModuleCorrNetworks
		#'
		#' Plot Module Eigengene correlogram
		#'
		#' @param seurat_obj A Seurat object
		#' @keywords scRNA-seq
		#' @export
		#' @examples
		#' ModuleCorrNetwork
		ModuleCorrNetwork <- function(
		seurat_obj, wgcna_name=NULL, cluster_col=NULL, exclude_grey=TRUE,
		reduction='umap', cor_cutoff=0.2, label_vertices=FALSE, edge_scale=5,
		vertex_size=15, niter=100, vertex_frame=FALSE
		){

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

		# Get module eigengenes
		MEs <- GetMEs(seurat_obj, wgcna_name)
		modules <- GetModules(seurat_obj, wgcna_name)

		# exclude grey
		if(exclude_grey){
		MEs <- MEs[,colnames(MEs) != 'grey']
		modules <- modules %>% subset(color != 'grey')
		}

		# get list of modules
		mods <- colnames(MEs)

		# what clusters are we using?
		if(is.null(clusters)){
		cluster_col <- Idents(seurat_obj)
		} else{
		}
		clusters <- droplevels(seurat_obj$annotation)
		MEs$cluster <- clusters

		# compute average MEs for each cluster
		cluster_ME_av <- do.call(
		rbind, lapply(
		split(MEs, MEs$cluster),
		function(x){colMeans(x[,mods])
		})) %>% as.data.frame

		# which cluster mas highest expression of each module?
		top_clusters <- sapply(mods, function(x){
		rownames(cluster_ME_av[cluster_ME_av[,x] == max(cluster_ME_av[,x]),])
		})

		# get UMAP / tSNE centroids for these clusters to use as starting coordinates
		red_df <- as.data.frame(seurat_obj@reductions[[reduction]]@cell.embeddings)
		red_df$cluster <- clusters

		# compute average coords for each cluster
		red_av <- do.call(
		rbind, lapply(
		split(red_df, red_df$cluster),
		function(x){colMeans(x[,1:2])
		})) %>% as.data.frame



		# get correlation matrix
		cor_mat <- Hmisc::rcorr(as.matrix(MEs[,1:ncol(MEs)-1]))$r
		cor_mat[lower.tri(cor_mat)] <- NA

		# melt matrix and remove NAs
		cor_df <- reshape2::melt(cor_mat) %>% na.omit

		# remove self-edges
		cor_df <- cor_df %>% subset(!(Var1 == Var2))

		# remove weak edges:
		cor_df <- cor_df %>% subset(abs(value) >= cor_cutoff)

		# vertex df:
		v_df <- data.frame(
		name = mods,
		cluster = as.character(top_clusters)
		)

		# add module colors:
		unique_mods <- distinct(modules[,c('module', 'color')])
		rownames(unique_mods) <- unique_mods$module
		v_df$color <- unique_mods[v_df$name,'color']

		# add reduction coords:
		v_df$x <- red_av[v_df$cluster, 1]
		v_df$y <- red_av[v_df$cluster, 2]

		# make graph:
		g <- igraph::graph_from_data_frame(
		cor_df,
		directed=FALSE,
		vertices=v_df
		)

		# igraph layout
		e <- get.edgelist(g, name=FALSE)
		l <- igraph::layout_with_fr(
		g, weights=E(g)$value,
		coords = as.matrix(data.frame(x=V(g)$x, y=V(g)$y)),
		niter=niter
		)

		# ggplot just to get the colors
		plot_df <- rbind(cor_df, data.frame(Var1=c('x', 'y'), Var2=c('y', 'x'), value=c(-1,1)))
		temp <- ggplot(plot_df, aes(x=value, y=value, color=value)) +
		geom_point() + scale_color_gradient2(high='darkorchid1', mid='white', low='seagreen', midpoint=0)
		temp <- ggplot_build(temp)
		E(g)$color <- temp$data[[1]]$colour[1:nrow(cor_df)]

		# label the vertices?
		if(label_vertices){labels <- V(g)$name} else{labels <- NA}

		# vertex_frame
		if(vertex_frame){frame_color <- 'black'} else{frame_color <- V(g)$color}

		# plot the graph
		plot(
		g, layout=l,
		edge.color=E(g)$color,
		edge.curved=0,
		edge.width=abs(E(g)$value) * edge_scale,
		vertex.color=V(g)$color,
		vertex.frame.color=frame_color,
		vertex.label=labels,
		vertex.label.family='Helvetica',
		vertex.label.color = 'black',
		vertex.label.cex=0.5,
		vertex.size=vertex_size
		)

		# plot colorbar:
		colfunc <- colorRampPalette(c( "seagreen", "white", "darkorchid1"))
		image.plot(legend.only=T, zlim=c(-1,1), col=colfunc(256) )


		}


		#' MEFeaturePlot
		#'
		@@ -38,14 +235,23 @@ PlotDendrogram <- function(
		#' @examples
		#' MEFeaturePlot
		MEFeaturePlot<- function(
		seurat_obj, modules,
		seurat_obj, module_names=NULL, wgcna_name = NULL,
		harmonized=TRUE, reduction='umap',
		order=TRUE, restrict_range=TRUE, point_size = 0.5,
		label_legend = FALSE
		){

		# get MEs from seurat object
		MEs <- GetMEs(seurat_obj, harmonized)
		if(is.null(wgcna_name)){wgcna_name <- seurat_obj@misc$active_wgcna}

		# get MEs, module data from seurat object
		MEs <- GetMEs(seurat_obj, harmonized, wgcna_name)
		modules <- GetModules(seurat_obj, wgcna_name)

		# use all modules except gray if not specified by the user
		if(is.null(module_names)){
		module_names <- colnames(MEs)
		module_names <- module_names[module_names != 'grey']
		}

		# get reduction from seurat obj
		umap <- seurat_obj@reductions[[reduction]]@cell.embeddings
		@@ -56,11 +262,10 @@ MEFeaturePlot<- function(
		plot_df <- cbind(umap, MEs) %>% as.data.frame()

		plot_list <- list()
		for(cur_mod in modules){
		for(cur_mod in module_names){

		cur_color <- cur_mod
		print(cur_color)
		print(head(plot_df[,cur_mod]))
		# get the color for this module:
		cur_color <- modules %>% subset(module == cur_mod) %>% .$color %>% unique

		# reset the range of the plot:
		plot_range <- plot_df[,cur_mod] %>% range
		@@ -89,7 +294,7 @@ MEFeaturePlot<- function(
		# plot with ggplot
		plot_list[[cur_mod]] <- plot_df %>%
		ggplot(aes_string(x=x_name, y=y_name, color=cur_mod)) +
		geom_point(size=0.5) +
		geom_point(size=point_size) +
		scale_color_gradient2(
		low='grey75', mid='grey95', high=cur_color,
		breaks = c(plot_range[1], 0, plot_range[2]),

Admin message