#' Pick power to fit network to a scale-free topology
#'
#' @param exp A gene expression data frame with genes in row names and
#' samples in column names or a `SummarizedExperiment` object.
#' @param net_type Network type. One of 'signed', 'signed hybrid' or
#' 'unsigned'. Default is signed.
#' @param rsquared R squared cutoff. Default is 0.8.
#' @param cor_method Correlation method. One of "pearson", "biweight"
#' or "spearman". Default is "spearman".
#'
#' @return A list containing: \itemize{
#' \item{power}{Optimal power based on scale-free topology fit}
#' \item{plot}{A ggplot object displaying main statistics of the SFT fit test}
#' }
#'
#' @author Fabricio Almeida-Silva
#' @seealso
#' \code{\link[WGCNA]{pickSoftThreshold}}
#' @rdname SFT_fit
#' @export
#' @importFrom WGCNA pickSoftThreshold
#' @importFrom ggplot2 theme element_text geom_hline
#' @examples
#' data(filt.se)
#' sft <- SFT_fit(filt.se, cor_method = "pearson")
SFT_fit <- function(exp, net_type = "signed", rsquared = 0.8,
cor_method = "spearman") {
exp <- handleSE(exp)
texp <- t(exp)
if(cor_method == "pearson") {
sft <- WGCNA::pickSoftThreshold(
texp, networkType = net_type, powerVector = 3:20,
RsquaredCut = rsquared
)
} else if(cor_method == "biweight") {
sft <- WGCNA::pickSoftThreshold(
texp, networkType = net_type, powerVector = 3:20,
RsquaredCut = rsquared, corFnc = bicor,
corOptions = list(use = 'p', maxPOutliers = 0.05)
)
} else if (cor_method == "spearman"){
sft <- WGCNA::pickSoftThreshold(
texp, networkType = net_type, powerVector=3:20,
RsquaredCut = rsquared,
corOptions = list(use = 'p', method = "spearman")
)
} else {
stop("Please, specify a correlation method (one of 'spearman', 'pearson' or 'biweight').")
}
wgcna_power <- sft$powerEstimate
if(is.na(wgcna_power)){
message("No power reached R-squared cut-off, now choosing max R-squared based power")
wgcna_power <- sft$fitIndices$Power[which(sft$fitIndices$SFT.R.sq == max(sft$fitIndices$SFT.R.sq))]
}
# Create data frame with indices to plot
sft_df <- data.frame(
power = sft$fitIndices$Power,
fit = -sign(sft$fitIndices$slope) * sft$fitIndices$SFT.R.sq,
meank = sft$fitIndices$mean.k.
)
sft_df <- sft_df[sft_df$fit > 0, ]
# Plot 1
p1 <- ggplot(sft_df, aes_(x = ~power, y = ~fit)) +
geom_point(color = "gray10") +
ggrepel::geom_text_repel(aes_(label = ~power)) +
labs(
x = "Soft threshold (power)",
y = expression(paste("Scale-free topology fit - ", R^{2})),
title = "Scale independence"
) +
theme_bw() +
ylim(c(0,1)) +
geom_hline(yintercept = rsquared, color="brown3") +
theme(legend.position = "none")
# Plot 2
p2 <- ggplot(sft_df, aes_(x = ~power, y = ~meank)) +
geom_point(color = "gray10") +
ggrepel::geom_text_repel(aes_(label = ~power)) +
labs(
x = "Soft threshold (power)",
y = "Mean connectivity (k)",
title = "Mean connectivity"
) +
theme_bw()
# Combined plot
sft_plot <- patchwork::wrap_plots(p1, p2, ncol = 2)
result <- list(power = as.numeric(wgcna_power), plot = sft_plot)
return(result)
}
#' Reconstruct gene coexpression network from gene expression
#'
#' @param exp A gene expression data frame with genes in row names and
#' samples in column names or a `SummarizedExperiment` object.
#' @param net_type Network type. One of 'signed', 'signed hybrid' or
#' 'unsigned'. Default: 'signed'.
#' @param module_merging_threshold Correlation threshold to merge similar
#' modules into a single one. Default: 0.8.
#' @param SFTpower SFT power generated by the function \code{SFT_fit}.
#' @param cor_method Correlation method. One of "pearson", "biweight" or
#' "spearman". Default is "spearman".
#' @param verbose Logical indicating whether to display progress
#' messages or not. Default: FALSE.
#'
#' @return List containing: \itemize{
#' \item Adjacency matrix
#' \item Data frame of module eigengenes
#' \item Data frame of genes and their corresponding modules
#' \item Data frame of intramodular connectivity
#' \item Correlation matrix
#' \item Parameters used for network reconstruction
#' \item Objects to plot the dendrogram in \code{plot_dendro_and_colors}.
#' }
#' @author Fabricio Almeida-Silva
#' @seealso
#' \code{\link[WGCNA]{adjacency.fromSimilarity}},\code{\link[WGCNA]{TOMsimilarity}},\code{\link[WGCNA]{standardColors}},\code{\link[WGCNA]{labels2colors}},\code{\link[WGCNA]{moduleEigengenes}},\code{\link[WGCNA]{plotEigengeneNetworks}},\code{\link[WGCNA]{mergeCloseModules}},\code{\link[WGCNA]{plotDendroAndColors}},\code{\link[WGCNA]{intramodularConnectivity}}
#' \code{\link[dynamicTreeCut]{cutreeDynamicTree}}
#' @rdname exp2gcn
#' @export
#' @importFrom WGCNA TOMsimilarity standardColors labels2colors moduleEigengenes plotEigengeneNetworks mergeCloseModules plotDendroAndColors intramodularConnectivity
#' @importFrom dynamicTreeCut cutreeDynamicTree
#' @importFrom stats as.dist median cor fisher.test hclust na.omit prcomp qnorm qqplot quantile sd var
#' @importFrom grDevices colorRampPalette
#' @examples
#' data(filt.se)
#' # The SFT fit was previously calculated and the optimal power was 16
#' gcn <- exp2gcn(filt.se, SFTpower = 18, cor_method = "pearson")
exp2gcn <- function(exp, net_type="signed",
module_merging_threshold = 0.8,
SFTpower = NULL, cor_method = "spearman",
verbose = FALSE) {
params <- list(
net_type = net_type,
module_merging_threshold = module_merging_threshold,
SFTpower = SFTpower,
cor_method = cor_method
)
norm.exp <- handleSE(exp)
if(is.null(SFTpower)) { stop("Please, specify the SFT power.") }
if(verbose) { message("Calculating adjacency matrix...") }
cor_matrix <- calculate_cor_adj(cor_method, norm.exp, SFTpower, net_type)$cor
adj_matrix <- calculate_cor_adj(cor_method, norm.exp, SFTpower, net_type)$adj
#Convert to matrix
gene_ids <- rownames(adj_matrix)
adj_matrix <- matrix(adj_matrix, nrow=nrow(adj_matrix))
rownames(adj_matrix) <- gene_ids
colnames(adj_matrix) <- gene_ids
#Calculate TOM from adjacency matrix
if(verbose) { message("Calculating topological overlap matrix (TOM)...") }
tomtype <- get_TOMtype(net_type)
TOM <- WGCNA::TOMsimilarity(adj_matrix, TOMType = tomtype)
#Hierarchically cluster genes
dissTOM <- 1-TOM #hclust takes a distance structure
geneTree <- hclust(as.dist(dissTOM), method="average")
#Detecting coexpression modules
if(verbose) { message("Detecting coexpression modules...") }
old.module_labels <- dynamicTreeCut::cutreeDynamicTree(
dendro = geneTree, minModuleSize = 30, deepSplit = TRUE
)
nmod <- length(unique(old.module_labels))
palette <- rev(WGCNA::standardColors(nmod))
old.module_colors <- WGCNA::labels2colors(
old.module_labels, colorSeq = palette
)
#Calculate eigengenes and merge the ones who are highly correlated
if(verbose) { message("Calculating module eigengenes (MEs)...") }
old.MElist <- WGCNA::moduleEigengenes(
t(norm.exp), colors = old.module_colors, softPower = SFTpower
)
old.MEs <- old.MElist$eigengenes
#Calculate dissimilarity of module eigengenes
MEDiss1 <- 1-cor(old.MEs)
#Hierarchically cluster module eigengenes to see how they're related
old.METree <- hclust(as.dist(MEDiss1), method="average")
#Then, choose a height cut
MEDissThreshold <- 1-module_merging_threshold
#Merge the modules.
if(verbose) { message("Merging similar modules...") }
if(cor_method == "pearson") {
merge1 <- WGCNA::mergeCloseModules(
t(norm.exp), old.module_colors, cutHeight = MEDissThreshold,
verbose = 0, colorSeq = palette
)
} else if(cor_method == "spearman") {
merge1 <- WGCNA::mergeCloseModules(
t(norm.exp), old.module_colors, cutHeight = MEDissThreshold,
verbose = 0, corOptions = list(use = "p", method = "spearman"),
colorSeq = palette
)
} else if(cor_method == "biweight") {
merge1 <- WGCNA::mergeCloseModules(
t(norm.exp), old.module_colors, cutHeight = MEDissThreshold,
verbose = 0, corFnc = bicor, colorSeq = palette
)
} else {
stop("Please, specify a correlation method. One of 'spearman', 'pearson' or 'biweight'.")
}
new.module_colors <- merge1$colors
new.MEs <- merge1$newMEs
#Plot dendrogram of merged modules eigengenes
new.METree <- hclust(as.dist(1-cor(new.MEs)), method="average")
#Get data frame with genes and modules they belong to
genes_and_modules <- as.data.frame(cbind(gene_ids, new.module_colors))
colnames(genes_and_modules) <- c("Genes", "Modules")
if(verbose) { message("Calculating intramodular connectivity...") }
kwithin <- WGCNA::intramodularConnectivity(adj_matrix, new.module_colors)
result.list <- list(
adjacency_matrix = adj_matrix,
MEs = new.MEs,
genes_and_modules = genes_and_modules,
kIN = kwithin,
correlation_matrix = cor_matrix,
params = params,
dendro_plot_objects = list(
tree = geneTree, unmerged = old.module_colors
)
)
return(result.list)
}
#' Plot eigengene network
#'
#' @param gcn List object returned by \code{exp2gcn}.
#'
#' @return A base plot with the eigengene network
#' @importFrom WGCNA plotEigengeneNetworks
#' @importFrom graphics par layout
#' @export
#' @rdname plot_eigengene_network
#' @examples
#' data(filt.se)
#' gcn <- exp2gcn(filt.se, SFTpower = 18, cor_method = "pearson")
#' plot_eigengene_network(gcn)
plot_eigengene_network <- function(gcn) {
on.exit(graphics::layout(1))
opar <- par(no.readonly = TRUE)
on.exit(par(opar), add = TRUE, after = FALSE)
p <- WGCNA::plotEigengeneNetworks(
gcn$MEs, "", marDendro = c(3,5,2,6), marHeatmap = c(3,4,2,2)
)
return(p)
}
#' Plot dendrogram of genes and modules
#'
#' @param gcn List object returned by \code{exp2gcn}.
#'
#' @return A base plot with the gene dendrogram and modules.
#' @importFrom WGCNA plotDendroAndColors
#' @importFrom graphics layout par
#' @export
#' @rdname plot_dendro_and_colors
#' @examples
#' data(filt.se)
#' gcn <- exp2gcn(filt.se, SFTpower = 18, cor_method = "pearson")
#' plot_dendro_and_colors(gcn)
plot_dendro_and_colors <- function(gcn) {
on.exit(graphics::layout(1))
opar <- par(no.readonly = TRUE)
on.exit(par(opar), add = TRUE, after = FALSE)
p <- WGCNA::plotDendroAndColors(
gcn$dendro_plot_objects$tree,
cbind(
gcn$dendro_plot_objects$unmerged,
gcn$moduleColors
),
c("Unmerged", "Merged"),
dendroLabels = FALSE, hang = 0.03, addGuide = TRUE, guideHang = 0.05
)
return(NULL)
}
#' Perform module stability analysis
#'
#' @param exp A gene expression data frame with genes in row names and
#' samples in column names or a `SummarizedExperiment` object.
#' @param net List object returned by \code{exp2gcn}.
#' @param nRuns Number of times to resample. Default is 20.
#'
#' @return A base plot with the module stability results.
#' @seealso
#' \code{\link[WGCNA]{sampledBlockwiseModules}},\code{\link[WGCNA]{matchLabels}},\code{\link[WGCNA]{plotDendroAndColors}}
#' @rdname module_stability
#' @export
#' @importFrom WGCNA sampledBlockwiseModules matchLabels plotDendroAndColors
#' @importFrom graphics par layout
#' @examples
#' data(filt.se)
#' filt <- filt.se[1:100, ] # reducing even further for testing purposes
#' # The SFT fit was previously calculated and the optimal power was 16
#' gcn <- exp2gcn(filt, SFTpower = 16, cor_method = "pearson")
#' # For simplicity, only 2 runs
#' module_stability(exp = filt, net = gcn, nRuns = 2)
module_stability <- function(exp, net, nRuns = 20) {
norm.exp <- handleSE(exp)
expr <- as.matrix(t(norm.exp))
net_type <- net$params$net_type
SFTpower <- net$params$SFTpower
cor_method <- net$params$cor_method
module_merging_threshold <- 1 - net$params$module_merging_threshold
TOMType <- get_TOMtype(net_type)
mods0 <- WGCNA::sampledBlockwiseModules(
nRuns = nRuns, replace = FALSE, datExpr = expr,
maxBlockSize = 5000, checkSoftPower = FALSE, corType = cor_method,
networkType = net_type, TOMType = TOMType, TOMDenom = "mean",
mergeCutHeight = module_merging_threshold,
reassignThreshold = 0, numericLabels = FALSE,
checkMissingData = FALSE, quickCor = 0, verbose = 2
)
nGenes <- ncol(expr)
# Define a matrix of labels for the original and all resampling runs
labels <- matrix(0, nGenes, nRuns + 1)
labels[, 1] <- mods0[[1]]$mods$colors
# Relabel modules in each of the resampling runs
labs <- Reduce(cbind, lapply(2:(nRuns+1), function(x) {
labels[, x] <- WGCNA::matchLabels(mods0[[x-1]]$mods$colors, labels[, 1])
}))
labels <- cbind(labels[,1], labs)
on.exit(graphics::layout(1))
opar <- par(no.readonly = TRUE)
on.exit(par(opar), add = TRUE, after = FALSE)
p <- WGCNA::plotDendroAndColors(
mods0[[1]]$mods$dendrograms[[1]], labels,
c("Full data set", paste("Resampling", seq_len(nRuns))),
main = "Dendrogram and modules: resampled data",
autoColorHeight = FALSE, colorHeight = 0.65,
dendroLabels = FALSE, hang = 0.03, guideHang = 0.05,
addGuide = TRUE, guideAll = FALSE, cex.main = 1, cex.lab = 0.8,
cex.colorLabels = 0.7, marAll = c(0, 5, 3, 0)
)
return(NULL)
}
#' Correlate module eigengenes to trait
#'
#' @param exp A gene expression data frame with genes in row names and
#' samples in column names or a `SummarizedExperiment` object.
#' @param metadata A data frame containing sample names in row names and
#' sample annotation in the first column. Ignored if `exp` is
#' a `SummarizedExperiment` object, since the function will extract colData.
#' @param MEs Module eigengenes. It is the 2nd element of the result list
#' generated by the function \code{exp2gcn}.
#' @param cor_method Method to calculate correlation. One of 'pearson',
#' 'spearman' or 'kendall'. Default is 'spearman'.
#' @param transpose Logical indicating whether to transpose the heatmap of not.
#' Default is FALSE.
#' @param palette RColorBrewer's color palette to use. Default is "RdYlBu",
#' a palette ranging from blue to red.
#' @param continuous_trait Logical indicating if trait is a continuous variable.
#' Default is FALSE.
#' @param cex.lab.x Font size for x axis labels. Default: 0.6.
#' @param cex.lab.y Font size for y axis labels. Default: 0.6.
#' @param cex.text Font size for numbers inside matrix. Default: 0.6.
#'
#' @return A data frame with correlation and correlation p-values for each pair
#' of ME and trait along with a heatmap.
#' @details Significance levels:
#' 1 asterisk: significant at alpha = 0.05.
#' 2 asterisks: significant at alpha = 0.01.
#' 3 asterisks: significant at alpha = 0.001.
#' no asterisk: not significant.
#'
#' @author Fabricio Almeida-Silva
#' @rdname module_trait_cor
#' @export
#' @importFrom WGCNA corPvalueStudent labeledHeatmap
#' @importFrom reshape2 melt
#' @importFrom RColorBrewer brewer.pal
#' @importFrom SummarizedExperiment colData
#' @importFrom graphics par
#' @examples
#' data(filt.se)
#' gcn <- exp2gcn(filt.se, SFTpower = 18, cor_method = "pearson")
#' module_trait_cor(filt.se, MEs=gcn$MEs)
module_trait_cor <- function(exp, metadata, MEs, cor_method = "spearman",
transpose = FALSE, palette = "RdYlBu",
continuous_trait = FALSE,
cex.lab.x = 0.6, cex.lab.y = 0.6,
cex.text = 0.6) {
if(is(exp, "SummarizedExperiment")) {
metadata <- as.data.frame(SummarizedExperiment::colData(exp))
}
exp <- handleSE(exp)
metadata <- metadata[colnames(exp), , drop=FALSE]
trait <- handle_trait_type(metadata, continuous_trait)
modtraitcor <- cor(as.matrix(MEs), trait, use = "p", method = cor_method)
nSamples <- ncol(exp)
modtraitpvalue <- WGCNA::corPvalueStudent(modtraitcor, nSamples)
modtraitsymbol <- pval2symbol(modtraitpvalue)
modtrait_long <- reshape2::melt(modtraitcor)
pvals_long <- reshape2::melt(modtraitpvalue)
combined_long <- merge(modtrait_long, pvals_long, by = c("Var1", "Var2"))
colnames(combined_long) <- c("ME", "trait", "cor", "pvalue")
combined_long$ME <- as.character(combined_long$ME)
combined_long$trait <- as.character(combined_long$trait)
textMatrix <- paste(signif(modtraitcor, 2), modtraitsymbol, sep = "")
dim(textMatrix) <- dim(modtraitcor)
if(transpose) {
modtraitcor <- t(modtraitcor)
textMatrix <- t(textMatrix)
yLabels <- colnames(trait)
xLabels <- names(MEs)
xSymbols <- names(MEs)
ySymbols <- NULL
xColorLabels <- TRUE
par(mar = c(5, 5, 1, 1))
} else {
par(mar = c(6, 8.5, 3, 3))
yLabels <- names(MEs)
xLabels <- colnames(trait)
xColorLabels <- FALSE
xSymbols <- NULL
ySymbols <- names(MEs)
}
cols <- colorRampPalette(rev(RColorBrewer::brewer.pal(10, palette)))(100)
hm <- WGCNA::labeledHeatmap(
Matrix = modtraitcor, yLabels = yLabels, xLabels = xLabels,
ySymbols = ySymbols, xSymbols = xSymbols, colorLabels = FALSE,
colors = cols, textMatrix = textMatrix, setStdMargins = FALSE,
cex.text = cex.text, cex.lab.x = cex.lab.x, cex.lab.y = cex.lab.y,
zlim = c(-1,1), cex.main = 1,
main = paste("Module-trait relationships")
)
return(combined_long)
}
#' Calculate gene significance for a given group of genes
#'
#' @param exp A gene expression data frame with genes in row names and
#' samples in column names or a `SummarizedExperiment` object.
#' @param metadata A data frame containing sample names in row names and
#' sample annotation in the first column. Ignored if `exp` is
#' a `SummarizedExperiment` object, since the function will extract colData.
#' @param genes Character vector of genes to be correlated with traits.
#' If not given, all genes in `exp` will be considered.
#' @param alpha Significance level. Default is 0.05.
#' @param min_cor Minimum correlation coefficient. Default is 0.2.
#' @param use_abs Logical indicating whether to filter by correlation using
#' absolute value or not. If TRUE, a \code{min_cor} of say 0.2 would keep all
#' correlations above 0.2 and below -0.2. Default is TRUE.
#' @param palette RColorBrewer's color palette to use. Default is "RdYlBu",
#' a palette ranging from blue to red.
#' @param show_rownames Logical indicating whether to show row names or not.
#' Default is FALSE.
#' @param continuous_trait Logical indicating if trait is a continuous variable.
#' Default is FALSE.
#'
#' @return A heatmap of gene significance (GS) and a list containing: \itemize{
#' \item{filtered_corandp}{Filtered matrix of correlation and p-values}
#' \item{raw_GS}{Raw matrix of gene significances}
#' }
#' @seealso
#' \code{\link[WGCNA]{corPvalueStudent}}
#' \code{\link[RColorBrewer]{RColorBrewer}}
#' @rdname gene_significance
#' @export
#' @importFrom reshape2 melt
#' @importFrom WGCNA corPvalueStudent
#' @importFrom ComplexHeatmap pheatmap
#' @importFrom RColorBrewer brewer.pal
#' @importFrom SummarizedExperiment colData
#' @examples
#' data(filt.se)
#' gs <- gene_significance(filt.se)
gene_significance <- function(exp, metadata, genes = NULL,
alpha = 0.05, min_cor = 0.2,
use_abs = TRUE,
palette = "RdYlBu", show_rownames = FALSE,
continuous_trait = FALSE) {
if(is(exp, "SummarizedExperiment")) {
metadata <- as.data.frame(SummarizedExperiment::colData(exp))
}
exp <- handleSE(exp)
final_exp <- exp[, rownames(metadata)]
if(!is.null(genes)) {
final_exp <- final_exp[genes, , drop = FALSE]
}
trait <- handle_trait_type(metadata, continuous_trait)
GS <- cor(as.matrix(t(final_exp)), trait, use = "p")
# Filter by correlation coefficient and p-value
melt.cor <- reshape2::melt(GS)
nSamples <- ncol(final_exp)
GS.pvalue <- WGCNA::corPvalueStudent(GS, nSamples)
melt.pvalue <- reshape2::melt(GS.pvalue)
corandp <- merge(melt.cor, melt.pvalue, by = c("Var1", "Var2"))
corandp$Var1 <- as.character(corandp$Var1)
corandp$Var2 <- as.character(corandp$Var2)
colnames(corandp) <- c("gene", "trait", "cor", "pvalue")
if(use_abs) {
corandp <- corandp[corandp$pvalue < alpha & abs(corandp$cor) > min_cor, ]
} else {
corandp <- corandp[corandp$pvalue < alpha & corandp$cor > min_cor, ]
}
cols <- colorRampPalette(rev(RColorBrewer::brewer.pal(10, palette)))(100)
p <- ComplexHeatmap::pheatmap(
as.matrix(GS), border_color = NA, color = cols,
show_rownames = show_rownames, main = "Gene-trait correlations"
)
resultlist <- list(filtered_corandp = corandp, raw_GS = GS)
return(resultlist)
}
#' Get GCN hubs
#'
#' @param exp A gene expression data frame with genes in row names and
#' samples in column names or a `SummarizedExperiment` object.
#' @param net List object returned by \code{exp2gcn}.
#'
#' @return Data frame containing gene IDs, modules and
#' intramodular connectivity of all hubs.
#' @author Fabricio Almeida-Silva
#' @seealso
#' \code{\link[WGCNA]{signedKME}}
#' @rdname get_hubs_gcn
#' @export
#' @importFrom WGCNA signedKME
#' @examples
#' data(filt.se)
#' gcn <- exp2gcn(filt.se, SFTpower = 18, cor_method = "pearson")
#' hubs <- get_hubs_gcn(filt.se, gcn)
get_hubs_gcn <- function(exp, net) {
exp <- handleSE(exp)
cor_method <- net$params$cor_method
genes_modules <- net$genes_and_modules
MEs <- net$MEs
kIN <- net$kIN[, 2, drop=FALSE]
# Add kIN info
genes_modulesk <- merge(
genes_modules, kIN, by.x = "Genes", by.y = "row.names"
)
# Calculate kME
if(cor_method == "spearman") {
MM <- WGCNA::signedKME(
t(exp), MEs, outputColumnName = "",
corOptions = "use = 'p', method = 'spearman'"
)
} else if(cor_method == "pearson") {
MM <- WGCNA::signedKME(t(exp), MEs, outputColumnName = "")
} else if(cor_method == "biweight") {
MM <- WGCNA::signedKME(
t(exp), MEs, corFnc = "bicor", outputColumnName = "",
corOptions = "maxPOutliers = 0.1"
)
}
# Add kME info
kme_info <- merge(genes_modulesk, MM, by.x="Genes", by.y="row.names")
# Keep top 10% genes with the highest degree for each module
genes_mod_list <- split(kme_info, kme_info[,2])
genes_mod_list$grey <- NULL
top_10 <- lapply(genes_mod_list, function(x) {
indices <- round(nrow(x) * 0.1)
return(x[order(x$kWithin, decreasing=TRUE), ][seq_len(indices), ])
})
# Pick genes from the top 10% degree with kME above 0.8
hubs <- Reduce(rbind, lapply(top_10, function(x) {
cols <- c("Genes", "Modules", "kWithin", unique(x$Modules))
h <- x[, cols]
h <- h[h[,4] > 0.8, c(1,2,3)]
}))
rownames(hubs) <- seq_len(nrow(hubs))
colnames(hubs) <- c("Gene", "Module", "kWithin")
return(hubs)
}
#' Helper function to perform Fisher's Exact Test with parallel computing
#'
#' @param genes Character vector containing genes on which enrichment will
#' be tested.
#' @param reference Character vector containing genes to be used as background.
#' @param genesets List of functional annotation categories.
#' (e.g., GO, pathway, etc.) with their associated genes.
#' @param adj Multiple testing correction method.
#' @param bp_param BiocParallel back-end to be used.
#' Default: BiocParallel::SerialParam()
#'
#' @return Results of Fisher's Exact Test in a data frame with TermID,
#' number of associated genes, number of genes in reference set,
#' P-value and adjusted P-value.
#'
#' @importFrom stats p.adjust fisher.test
#' @importFrom BiocParallel bplapply SerialParam
#' @noRd
par_enrich <- function(genes, reference, genesets, adj = "BH",
bp_param = BiocParallel::SerialParam()) {
reference <- reference[!reference %in% genes]
tab <- BiocParallel::bplapply(seq_along(genesets), function(i) {
RinSet <- sum(reference %in% genesets[[i]])
RninSet <- length(reference) - RinSet
GinSet <- sum(genes %in% genesets[[i]])
GninSet <- length(genes) - GinSet
fmat <- matrix(
c(GinSet, RinSet, GninSet, RninSet),
nrow = 2, ncol = 2, byrow = FALSE
)
colnames(fmat) <- c("inSet", "ninSet")
rownames(fmat) <- c("genes", "reference")
fish <- stats::fisher.test(fmat, alternative = "greater")
pval <- fish$p.value
inSet <- RinSet + GinSet
res <- c(GinSet, inSet, pval)
return(res)
}, BPPARAM = bp_param)
rtab <- do.call(rbind, tab)
rtab <- data.frame(as.vector(names(genesets)), rtab)
rtab <- rtab[order(rtab[, 4]), ]
colnames(rtab) <- c("TermID", "genes", "all", "pval")
padj <- stats::p.adjust(rtab$pval, method = adj)
tab.out <- data.frame(rtab, padj)
return(tab.out)
}
#' Perform enrichment analysis for a set of genes
#'
#' @param genes Character vector containing genes for overrepresentation
#' analysis.
#' @param background_genes Character vector of genes to be used as background
#' for the Fisher's Exact Test.
#' @param annotation Annotation data frame with genes in the first column and
#' functional annotation in the other columns. This data frame can be exported
#' from Biomart or similar databases.
#' @param column Column or columns of \code{annotation} to be used for
#' enrichment.
#' Both character or numeric values with column indices can be used.
#' If users want to supply more than one column, input a character or
#' numeric vector. Default: all columns from \code{annotation}.
#' @param correction Multiple testing correction method. One of "holm",
#' "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr" or "none".
#' Default is "BH".
#' @param p P-value threshold. P-values below this threshold will be
#' considered significant. Default is 0.05.
#' @param bp_param BiocParallel back-end to be used.
#' Default: BiocParallel::SerialParam()
#'
#' @return Data frame containing significant terms, p-values and
#' associated genes.
#'
#' @author Fabricio Almeida-Silva
#' @rdname enrichment_analysis
#' @export
#' @importFrom BiocParallel bplapply
#' @examples
#' \donttest{
#' data(filt.se)
#' data(zma.interpro)
#' genes <- rownames(filt.se)[1:50]
#' background_genes <- rownames(filt.se)
#' annotation <- zma.interpro
#' # Using p = 1 to show all results
#' enrich <- enrichment_analysis(genes, background_genes, annotation, p = 1)
#' }
enrichment_analysis <- function(genes, background_genes, annotation,
column = NULL,
correction = "BH", p = 0.05,
bp_param = BiocParallel::SerialParam()) {
## Get a dataframe of expressed genes and their annotations
gene_column <- colnames(annotation)[1]
# Handle missing values
annotation[is.na(annotation)] <- "Unannotated"
if(is.null(column)) { #all columns?
gene_annotation <- annotation
} else {
if(is.numeric(column)) { # Input is column indices
gene_annotation <- annotation[, c(1, column)]
} else { #Input is column names
gene_annotation <- annotation[, c(paste(gene_column), column)]
}
}
background <- gene_annotation[gene_annotation[,1] %in% background_genes, ]
if(ncol(background) == 2) { #only one annotation category
# Named list of expressed genes and their annotations
annotation_list <- split(background[,1], background[,2])
# Perform the enrichment analysis
enrich.table <- par_enrich(
genes, background_genes, annotation_list, adj = correction,
bp_param = bp_param
)
signif_enrich <- as.data.frame(enrich.table[enrich.table$padj <= p, ])
if(nrow(signif_enrich) > 0) {
signif_terms_genes <- annotation_list[as.character(signif_enrich[,1])]
signif_terms_genes <- lapply(signif_terms_genes, function(x) unique(x[x %in% genes]))
# Save final result with a column containing gene IDs of associated genes
signif_enrich$GeneID <- unlist(
lapply(signif_terms_genes, paste, collapse = ",")
)
df_signif_enrich <- signif_enrich
} else {
df_signif_enrich <- NULL
}
} else { # more than 1 annotation category
annotation_list <- lapply(seq_along(background)[-1], function(x) {
return(split(background[,1], background[,x]))
})
# Remove unannotated genes
annotation_list_final <- lapply(annotation_list, function(x) {
return(x[names(x) != "Unannotated"])
})
# Perform the enrichment analysis
signif_enrich <- bplapply(annotation_list_final, function(x) {
enrich.table <- par_enrich(
genes, background_genes, x, adj = correction,
bp_param = bp_param
)
sig_enrich <- as.data.frame(enrich.table[enrich.table$padj < p, ])
return(sig_enrich)
}, BPPARAM = bp_param)
# Remove empty data frames from list
signif_enrich <- signif_enrich[vapply(signif_enrich, nrow, numeric(1)) > 0]
if(length(signif_enrich) > 0) {
# Create a data frame containing annotations and the annotation class
annot_correspondence <- Reduce(rbind, lapply(seq_along(background)[-1], function(x) {
annot_correspondence <- data.frame(
TermID = background[, x],
Category = names(background)[x]
)
annot_correspondence <- annot_correspondence[!duplicated(annot_correspondence[,c(1,2)]),]
return(annot_correspondence)
}))
# Add column containing the annotation class
list_signif_enrich <- lapply(signif_enrich, function(x) {
merge(x, annot_correspondence)
})
# Reduce list of data frames to a single data frame
df_signif_enrich <- Reduce(rbind, list_signif_enrich)
# Get genes associated to each term
annotation_list_concat <- unlist(annotation_list, recursive = FALSE) # create list from list of lists
signif_terms_genes <- annotation_list_concat[as.character(df_signif_enrich[,1])]
signif_terms_genes <- lapply(signif_terms_genes, function(x) unique(x[x %in% genes]))
# Get final table with enriched terms and a column containing the associated genes
df_signif_enrich$GeneID <- unlist(
lapply(signif_terms_genes, paste, collapse = ",")
)
} else {
df_signif_enrich <- NULL
}
}
return(df_signif_enrich)
}
#' Perform enrichment analysis for coexpression network modules
#'
#' @param net List object returned by \code{exp2gcn}.
#' @param background_genes Character vector of genes to be used as background
#' for the Fisher's Exact Test.
#' @param annotation Annotation data frame with genes in the first column
#' and functional annotation in the other columns. This data frame can
#' be exported from Biomart or similar databases.
#' @param column Column or columns of \code{annotation} to be used for enrichment.
#' Both character or numeric values with column indices can be used.
#' If users want to supply more than one column, input a character or
#' numeric vector. Default: all columns from \code{annotation}.
#' @param correction Multiple testing correction method. One of "holm",
#' "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr" or "none".
#' Default is "BH".
#' @param p P-value threshold. P-values below this threshold will be considered
#' significant. Default is 0.05.
#' @param bp_param BiocParallel back-end to be used.
#' Default: BiocParallel::SerialParam()
#' @return A data frame containing enriched terms, p-values, gene IDs
#' and module names.
#' @author Fabricio Almeida-Silva
#' @rdname module_enrichment
#' @importFrom BiocParallel bplapply SerialParam
#' @export
#' @examples
#' \donttest{
#' data(filt.se)
#' data(zma.interpro)
#' background <- rownames(filt.se)
#' gcn <- exp2gcn(filt.se, SFTpower = 18, cor_method = "pearson")
#' mod_enrich <- module_enrichment(gcn, background, zma.interpro, p=1)
#' }
module_enrichment <- function(net = NULL, background_genes, annotation,
column = NULL,
correction = "BH", p = 0.05,
bp_param = BiocParallel::SerialParam()) {
# Divide modules in different data frames of a list
genes.modules <- net$genes_and_modules
list.gmodules <- split(genes.modules, genes.modules$Modules)
list.gmodules <- list.gmodules[names(list.gmodules) != "grey"]
enrichment_allmodules <- BiocParallel::bplapply(seq_along(list.gmodules), function(x) {
message("Enrichment analysis for module ",
names(list.gmodules)[x], "...")
l <- enrichment_analysis(
genes = as.character(list.gmodules[[x]][, 1]),
background_genes = background_genes,
annotation = annotation,
column = column, correction = correction, p = p
)
return(l)
}, BPPARAM = bp_param)
names(enrichment_allmodules) <- names(list.gmodules)
# Remove NULL elements from list
enrichment_filtered <- enrichment_allmodules[!vapply(enrichment_allmodules,
is.null, logical(1))]
if(length(enrichment_filtered) == 0) {
message("None of the modules had significant enrichment.")
enrichment_final <- NULL
} else {
# Add module name to each data frame
enrichment_modnames <- lapply(seq_along(enrichment_filtered), function(x) {
return(cbind(
enrichment_filtered[[x]], Module=names(enrichment_filtered)[x]
))
})
# Reduce list of data frames to a single data frame
enrichment_final <- Reduce(rbind, enrichment_modnames)
}
return(enrichment_final)
}
#' Get 1st-order neighbors of a given gene or group of genes
#'
#' @param genes Character vector containing genes from which
#' direct neighbors will be extracted.
#' @param net List object returned by \code{exp2gcn}.
#' @param cor_threshold Correlation threshold to filter connections.
#' As a weighted network is a fully connected graph, a cutoff must be selected.
#' Default is 0.7.
#'
#' @return List containing 1st-order neighbors for each input gene.
#' @seealso \code{exp2gcn} \code{SFT_fit}
#' @author Fabricio Almeida-Silva
#' @export
#' @rdname get_neighbors
#' @examples
#' data(filt.se)
#' genes <- rownames(filt.se)[1:10]
#' gcn <- exp2gcn(filt.se, SFTpower = 18, cor_method = "pearson")
#' neighbors <- get_neighbors(genes, gcn)
get_neighbors <- function(genes, net, cor_threshold = 0.7) {
net_type <- net$params$net_type
power <- net$params$SFTpower
edges <- net$correlation_matrix
edges <- cormat_to_edgelist(edges)
colnames(edges) <- c("Gene1", "Gene2", "Weight")
filt_edges <- edges[edges$Gene1 %in% genes | edges$Gene2 %in% genes, ]
if(net_type == "signed") {
filt_edges <- filt_edges[abs(filt_edges$Weight) >= cor_threshold, ]
} else {
filt_edges <- filt_edges[filt_edges$Weight >= cor_threshold, ]
}
result_list <- lapply(genes, function(x) {
y <- filt_edges[rowSums(filt_edges == x) > 0, ]
y <- c(as.character(y$Gene1), as.character(y$Gene2))
y <- unique(y[y != x])
return(y)
})
names(result_list) <- genes
return(result_list)
}
#' Get edge list from an adjacency matrix for a group of genes
#'
#' @param net List object returned by \code{exp2gcn}.
#' @param genes Character vector containing a subset of genes from which
#' edges will be extracted. It can be ignored if the user wants to extract
#' an edge list for a given module instead of individual genes.
#' @param module Character with module name from which edges will be extracted.
#' To include 2 or more modules, input the names in a character vector.
#' @param filter Logical indicating whether to filter the edge list or not.
#' @param method Method to filter spurious correlations. One of "Zscore",
#' "optimalSFT", "pvalue" or "min_cor". See details for more information
#' on the methods. Default: 'optimalSFT'
#' @param r_optimal_test Numeric vector with the correlation thresholds
#' to be tested for optimal scale-free topology fit. Only valid
#' if \code{method} equals "optimalSFT". Default: seq(0.4, 0.9, by = 0.1)
#' @param Zcutoff Minimum Z-score threshold. Only valid if
#' \code{method} equals "Zscore". Default: 1.96
#' @param pvalue_cutoff Maximum P-value threshold. Only valid
#' if \code{method} equals "pvalue". Default: 0.05
#' @param rcutoff Minimum correlation threshold. Only valid
#' if \code{method} equals "min_cor". Default: 0.7
#' @param nSamples Number of samples in the data set from which
#' the correlation matrix was calculated. Only required
#' if \code{method} equals "pvalue".
#' @param check_SFT Logical indicating whether to test if the resulting network
#' is close to a scale-free topology or not. Default: FALSE.
#' @param bp_param BiocParallel back-end to be used.
#' Default: BiocParallel::SerialParam()
#'
#' @return Data frame with edge list for the input genes.
#' @details The default method ("optimalSFT") will create several different
#' edge lists by filtering the original correlation matrix by the thresholds
#' specified in \code{r_optimal_test}. Then, it will calculate a scale-free
#' topology fit index for each of the possible networks and return the network
#' that best fits the scale-free topology.
#' The method "Zscore" will apply a Fisher Z-transformation for the correlation
#' coefficients and remove the Z-scores below the threshold specified
#' in \code{Zcutoff}.
#' The method "pvalue" will calculate Student asymptotic p-value for the
#' correlations and remove correlations whose p-values are above the threshold
#' specified in \code{pvalue_cutoff}.
#' The method "min_cor" will remove correlations below the minimum correlation
#' threshold specified in \code{rcutoff}.
#' @seealso
#' \code{\link[WGCNA]{scaleFreeFitIndex}},\code{\link[WGCNA]{corPvalueStudent}}
#' \code{\link[igraph]{fit_power_law}}
#' @seealso \code{SFT_fit}
#' @seealso \code{exp2gcn}.
#' @author Fabricio Almeida-Silva
#' @rdname get_edge_list
#' @export
#' @importFrom WGCNA scaleFreeFitIndex corPvalueStudent
#' @importFrom ggplot2 theme geom_point geom_line theme_bw
#' @importFrom BiocParallel bplapply SerialParam
#' @examples
#' data(filt.se)
#' gcn <- exp2gcn(filt.se, SFTpower = 18, cor_method = "pearson")
#' genes <- rownames(filt.se)[1:50]
#' edges <- get_edge_list(gcn, genes=genes, filter = FALSE)
get_edge_list <- function(net, genes = NULL, module = NULL,
filter = FALSE, method = "optimalSFT",
r_optimal_test = seq(0.4, 0.9, by=0.1),
Zcutoff = 1.96, pvalue_cutoff = 0.05, rcutoff = 0.7,
nSamples = NULL,
check_SFT = FALSE,
bp_param = BiocParallel::SerialParam()) {
# Define objects containing correlation matrix and data frame of genes and modules
cor_matrix <- net$correlation_matrix
# Should we extract genes in a module?
if(!is.null(module)) {
genes_modules <- net$genes_and_modules
keep <- genes_modules[genes_modules$Modules %in% module, 1]
cor_matrix <- cor_matrix[keep, keep]
}
# Should we extract a user-defined gene set?
if(!is.null(genes)) {
cor_matrix <- cor_matrix[genes, genes]
}
# Should we filter the matrix?
if(filter) {
# Create edge list from correlation matrix
edges <- cormat_to_edgelist(cor_matrix)
colnames(edges) <- c("Gene1", "Gene2", "Weight")
if(method == "Zscore") {
# Apply Fisher-Z transformation to correlation values
edgesZ <- edges
edgesZ$Weight <- 0.5 * log((1+edges$Weight) / (1-edges$Weight))
edgelist <- edgesZ[edgesZ$Weight >= Zcutoff, ]
} else if(method == "optimalSFT") {
cutoff <- r_optimal_test
# Create list of 2 elements: edge lists, each with a different correlation threshold, and degree
list_mat <- replicate(length(cutoff), cor_matrix, simplify = FALSE)
list_cormat_filtered <- BiocParallel::bplapply(seq_along(list_mat), function(x) {
# Convert r values below threshold to NA and set diagonals to 0
matrix <- list_mat[[x]]
matrix[matrix < cutoff[x] ] <- NA
diag(matrix) <- 0
# Calculate degree
degree <- rowSums(matrix, na.rm=TRUE)
# Convert lower triangle and diagonals to NA
matrix[lower.tri(matrix, diag=TRUE)] <- NA
# Convert symmetrical matrix to edge list (Gene1, Gene2, Weight)
matrix <- na.omit(data.frame(as.table(matrix)))
result <- list(matrix = matrix, degree = degree)
return(result)
}, BPPARAM = bp_param)
degree_list <- lapply(list_cormat_filtered, function(x) return(x[[2]]))
# Calculate scale-free topology
sft.rsquared <- unlist(lapply(degree_list, function(x) WGCNA::scaleFreeFitIndex(x)$Rsquared.SFT))
max.index <- which.max(sft.rsquared)
# Plot scale-free topology fit for r values
plot.data <- data.frame(x = cutoff, y = sft.rsquared)
plot <- ggplot(plot.data, aes_(x = ~x, y = ~y, group = 1)) +
geom_point(color = "firebrick", size = 4) +
geom_line(color = "firebrick", size = 2) +
labs(
x = "Correlation (r) values",
y = expression(paste("Scale-free topology fit - ", R^{2})),
title = "Scale-free topology fit for given r values"
) +
theme_bw()
print(plot)
optimalr <- cutoff[max.index]
message("The correlation threshold that best fits the scale-free topology is ", optimalr)
edgelist <- list_cormat_filtered[[max.index]][[1]]
} else if(method == "pvalue") {
if(is.null(nSamples)) {
stop("Please, specify the number of samples used to calculate the correlation matrix")
}
# Calculate Student asymptotic p-value for given correlations
cor.pvalue <- WGCNA::corPvalueStudent(edges$Weight, nSamples)
# Create a data frame of correlations and p-values
corandp <- edges
corandp$pvalue <- cor.pvalue
# Create a final edge list containing only significant correlations
edgelist <- corandp[corandp$pvalue < pvalue_cutoff, c(1,2,3)]
} else if(method == "min_cor") {
edgelist <- edges[edges$Weight >= rcutoff, ]
} else{
stop("Please, specify a valid filtering method.")
}
} else {
# Create edge list from correlation matrix without filtering
edgelist <- cormat_to_edgelist(cor_matrix)
colnames(edgelist) <- c("Gene1", "Gene2", "Weight")
}
# Check scale-free topology fit
if(check_SFT) {
test <- check_SFT(edgelist, net_type="gcn")
}
return(edgelist)
}
