multiWGCNA: the full workflow
2024-05-06
Contents
- 1 Introduction
- 2 Install the multiWGCNA R package
- 3 Load microarray data from human post-mortem brains
- 4 Perform network construction, module eigengene calculation, module-trait correlation
- 5 Compare modules by overlap across conditions
- 6 Perform differential module expression analysis
- 7 Perform the module preservation analysis
- 8 Summarize interesting results from the analyses
- 9 Print the session info
1 Introduction
The multiWGCNA R package is a WGCNA-based procedure designed for RNA-seq datasets with two biologically meaningful variables. We developed this package because we’ve found that the experimental design for network analysis can be ambiguous. For example, should I analyze my treatment and wildtype samples separately or together? We find that it is informative to perform all the possible design combinations, and subsequently integrate these results. For more information, please see Tommasini and Fogel. BMC Bioinformatics. 2023.
In this vignette, we will be showing how users can perform a full start-to-finish multiWGCNA analysis. We will be using the regional autism data from Voineagu et al. 2011 (https://www.nature.com/articles/nature10110). This dataset has two sample traits: 1) autism versus control, and 2) temporal cortex versus frontal cortex.
2 Install the multiWGCNA R package
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("multiWGCNA")The development version of multiWGCNA can be installed from GitHub like this:
if (!require("devtools", quietly = TRUE))
install.packages("devtools")
devtools::install_github("fogellab/multiWGCNA")Load multiWGCNA library
library(multiWGCNA)3 Load microarray data from human post-mortem brains
Since this is purely an example of how to run the entire analysis from start to finish, we are going to limit our analysis to 1500 randomly selected probes.
Input data:
- The expression data can be either a SummarizedExperiment object or a data.frame with genes as rows and samples as columns
- The sampleTable should be a data.frame with the first column named “Sample”. The second and third column names can be arbitrary, i.e. “disease” or “Genotype”
# Download data from the ExperimentHub
library(ExperimentHub)
eh = ExperimentHub()
# Note: this requires the SummarizedExperiment package to be installed
eh_query = query(eh, c("multiWGCNAdata"))
autism_se = eh_query[["EH8219"]]
#> see ?multiWGCNAdata and browseVignettes('multiWGCNAdata') for documentation
#> loading from cache
# Collect the metadata in the sampleTable
sampleTable = colData(autism_se)
# Randomly sample 1500 genes from the expression matrix
set.seed(1)
autism_se = autism_se[sample(rownames(autism_se), 1500),]
# Check the data
assays(autism_se)[[1]][1:5, 1:5]
#> GSM706412 GSM706413 GSM706414 GSM706415 GSM706416
#> ILMN_1672121 11.034264 10.446682 11.473705 11.732849 11.43105
#> ILMN_2151368 10.379812 9.969130 9.990030 9.542288 10.26247
#> ILMN_1757569 9.426955 9.050024 9.347505 9.235251 9.38837
#> ILMN_2400219 12.604047 12.886037 12.890658 12.446960 12.98925
#> ILMN_2222101 12.385019 12.748229 12.418027 11.690253 13.10915
# sampleTable$Status = paste0("test_", sampleTable$Status, "_test")
# sampleTable$Sample = paste0(sampleTable$Sample, "_test")
sampleTable
#> DataFrame with 58 rows and 3 columns
#> Sample Status Tissue
#> <character> <character> <character>
#> GSM706412 GSM706412 autism FC
#> GSM706413 GSM706413 autism FC
#> GSM706414 GSM706414 autism FC
#> GSM706415 GSM706415 autism FC
#> GSM706416 GSM706416 autism FC
#> ... ... ... ...
#> GSM706465 GSM706465 controls TC
#> GSM706466 GSM706466 controls TC
#> GSM706467 GSM706467 controls TC
#> GSM706468 GSM706468 controls TC
#> GSM706469 GSM706469 controls TC
# Define our conditions for trait 1 (disease) and 2 (brain region)
conditions1 = unique(sampleTable[,2])
conditions2 = unique(sampleTable[,3])4 Perform network construction, module eigengene calculation, module-trait correlation
Also, let’s find best trait for each module and identify outlier modules (ie modules driven by a single sample).
Let’s use power = 10 since Voineagu et al. 2011 used that for all their networks. They also constructed “unsigned” networks and a mergeCutHeight of 0.10 so that modules with similar expression are merged. This step may take a minute.
# Construct the combined networks and all the sub-networks (autism only, controls only, FC only, and TC only)
autism_networks = constructNetworks(autism_se, sampleTable, conditions1, conditions2,
networkType = "unsigned", power = 10,
minModuleSize = 40, maxBlockSize = 25000,
reassignThreshold = 0, minKMEtoStay = 0.7,
mergeCutHeight = 0.10, numericLabels = TRUE,
pamRespectsDendro = FALSE, verbose=3,
saveTOMs = TRUE)
#> Flagging genes and samples with too many missing values...
#> ..step 1
#> Calculating module eigengenes block-wise from all genes
#> Flagging genes and samples with too many missing values...
#> ..step 1
#> ..Working on block 1 .
#> TOM calculation: adjacency..
#> ..will not use multithreading.
#> Fraction of slow calculations: 0.000000
#> ..connectivity..
#> ..matrix multiplication (system BLAS)..
#> ..normalization..
#> ..done.
#> ..saving TOM for block 1 into file combined-block.1.RData
#> ....clustering..
#> ....detecting modules..
#> ....calculating module eigengenes..
#> ....checking kME in modules..
#> ..removing 163 genes from module 1 because their KME is too low.
#> ..removing 126 genes from module 2 because their KME is too low.
#> ..removing 59 genes from module 3 because their KME is too low.
#> ..removing 18 genes from module 4 because their KME is too low.
#> ..removing 28 genes from module 5 because their KME is too low.
#> ..removing 24 genes from module 6 because their KME is too low.
#> ..removing 9 genes from module 7 because their KME is too low.
#> ..merging modules that are too close..
#> mergeCloseModules: Merging modules whose distance is less than 0.1
#> Calculating new MEs...
#> softConnectivity: FYI: connecitivty of genes with less than 20 valid samples will be returned as NA.
#> ..calculating connectivities..
#> Flagging genes and samples with too many missing values...
#> ..step 1
#> Calculating module eigengenes block-wise from all genes
#> Flagging genes and samples with too many missing values...
#> ..step 1
#> ..Working on block 1 .
#> TOM calculation: adjacency..
#> ..will not use multithreading.
#> Fraction of slow calculations: 0.000000
#> ..connectivity..
#> ..matrix multiplication (system BLAS)..
#> ..normalization..
#> ..done.
#> ..saving TOM for block 1 into file autism-block.1.RData
#> ....clustering..
#> ....detecting modules..
#> ....calculating module eigengenes..
#> ....checking kME in modules..
#> ..removing 139 genes from module 1 because their KME is too low.
#> ..removing 82 genes from module 2 because their KME is too low.
#> ..removing 95 genes from module 3 because their KME is too low.
#> ..removing 37 genes from module 4 because their KME is too low.
#> ..removing 35 genes from module 5 because their KME is too low.
#> ..removing 36 genes from module 6 because their KME is too low.
#> ..removing 27 genes from module 7 because their KME is too low.
#> ..merging modules that are too close..
#> mergeCloseModules: Merging modules whose distance is less than 0.1
#> Calculating new MEs...
#> softConnectivity: FYI: connecitivty of genes with less than 10 valid samples will be returned as NA.
#> ..calculating connectivities..
#> Flagging genes and samples with too many missing values...
#> ..step 1
#> Calculating module eigengenes block-wise from all genes
#> Flagging genes and samples with too many missing values...
#> ..step 1
#> ..Working on block 1 .
#> TOM calculation: adjacency..
#> ..will not use multithreading.
#> Fraction of slow calculations: 0.000000
#> ..connectivity..
#> ..matrix multiplication (system BLAS)..
#> ..normalization..
#> ..done.
#> ..saving TOM for block 1 into file controls-block.1.RData
#> ....clustering..
#> ....detecting modules..
#> ....calculating module eigengenes..
#> ....checking kME in modules..
#> ..removing 212 genes from module 1 because their KME is too low.
#> ..removing 107 genes from module 2 because their KME is too low.
#> ..removing 65 genes from module 3 because their KME is too low.
#> ..removing 51 genes from module 4 because their KME is too low.
#> ..removing 50 genes from module 5 because their KME is too low.
#> ..removing 40 genes from module 6 because their KME is too low.
#> ..merging modules that are too close..
#> mergeCloseModules: Merging modules whose distance is less than 0.1
#> Calculating new MEs...
#> softConnectivity: FYI: connecitivty of genes with less than 10 valid samples will be returned as NA.
#> ..calculating connectivities..
#> Flagging genes and samples with too many missing values...
#> ..step 1
#> Calculating module eigengenes block-wise from all genes
#> Flagging genes and samples with too many missing values...
#> ..step 1
#> ..Working on block 1 .
#> TOM calculation: adjacency..
#> ..will not use multithreading.
#> Fraction of slow calculations: 0.000000
#> ..connectivity..
#> ..matrix multiplication (system BLAS)..
#> ..normalization..
#> ..done.
#> ..saving TOM for block 1 into file FC-block.1.RData
#> ....clustering..
#> ....detecting modules..
#> ....calculating module eigengenes..
#> ....checking kME in modules..
#> ..removing 164 genes from module 1 because their KME is too low.
#> ..removing 138 genes from module 2 because their KME is too low.
#> ..removing 86 genes from module 3 because their KME is too low.
#> ..removing 33 genes from module 4 because their KME is too low.
#> ..removing 35 genes from module 5 because their KME is too low.
#> ..removing 34 genes from module 6 because their KME is too low.
#> ..removing 19 genes from module 7 because their KME is too low.
#> ..merging modules that are too close..
#> mergeCloseModules: Merging modules whose distance is less than 0.1
#> Calculating new MEs...
#> softConnectivity: FYI: connecitivty of genes with less than 11 valid samples will be returned as NA.
#> ..calculating connectivities..
#> Flagging genes and samples with too many missing values...
#> ..step 1
#> Calculating module eigengenes block-wise from all genes
#> Flagging genes and samples with too many missing values...
#> ..step 1
#> ..Working on block 1 .
#> TOM calculation: adjacency..
#> ..will not use multithreading.
#> Fraction of slow calculations: 0.000000
#> ..connectivity..
#> ..matrix multiplication (system BLAS)..
#> ..normalization..
#> ..done.
#> ..saving TOM for block 1 into file TC-block.1.RData
#> ....clustering..
#> ....detecting modules..
#> ....calculating module eigengenes..
#> ....checking kME in modules..
#> ..removing 146 genes from module 1 because their KME is too low.
#> ..removing 119 genes from module 2 because their KME is too low.
#> ..removing 52 genes from module 3 because their KME is too low.
#> ..removing 31 genes from module 4 because their KME is too low.
#> ..removing 34 genes from module 5 because their KME is too low.
#> ..removing 12 genes from module 6 because their KME is too low.
#> ..removing 13 genes from module 7 because their KME is too low.
#> ..removing 10 genes from module 8 because their KME is too low.
#> ..merging modules that are too close..
#> mergeCloseModules: Merging modules whose distance is less than 0.1
#> Calculating new MEs...
#> softConnectivity: FYI: connecitivty of genes with less than 9 valid samples will be returned as NA.
#> ..calculating connectivities..5 Compare modules by overlap across conditions
For calculating significance of overlap, we use the hypergeometric test (also known as the one-sided Fisher exact test). We’ll save the results in a list. Then, let’s take a look at the mutual best matches between autism modules and control modules.
# Save results to a list
results=list()
results$overlaps=iterate(autism_networks, overlapComparisons, plot=TRUE)
#>
#> #### comparing combined and autism ####
#>
#> #### comparing combined and controls ####
#>
#> #### comparing combined and FC ####
#>
#> #### comparing combined and TC ####