multiWGCNA: visualizing condition-specific networks
2024-05-06
Contents
- 1 Introduction
- 2 Load multiWGCNA library
- 3 Load astrocyte Ribotag RNA-seq data
- 4 Network construction
- 5 Compare modules by overlap
- 6 Identify a module of interest
- 7 Draw the multiWGCNA network
- 8 Observe differential co-expression of top module genes
- 9 Follow up with a preservation analysis
- 10 Determining if preservation value is significant
- 11 Conclusion
1 Introduction
In this vignette, we will be showing various ways users can analyze their condition-specific network across the conditions of their experiment. We will be using the astrocyte Ribotag data from Itoh et al. PNAS. 2018 (https://doi.org/10.1073/pnas.1716032115). This data comes with the multiWGCNAdata package hosted on ExperimentHub and can be accessed as shown below.
2 Load multiWGCNA library
library(multiWGCNA)
#> Loading required package: ggalluvial
#> Loading required package: ggplot2
#> 3 Load astrocyte Ribotag RNA-seq data
# Download data from the ExperimentHub
library(ExperimentHub)
#> Loading required package: BiocGenerics
#>
#> Attaching package: 'BiocGenerics'
#> The following objects are masked from 'package:stats':
#>
#> IQR, mad, sd, var, xtabs
#> The following objects are masked from 'package:base':
#>
#> Filter, Find, Map, Position, Reduce, anyDuplicated, aperm, append,
#> as.data.frame, basename, cbind, colnames, dirname, do.call,
#> duplicated, eval, evalq, get, grep, grepl, intersect, is.unsorted,
#> lapply, mapply, match, mget, order, paste, pmax, pmax.int, pmin,
#> pmin.int, rank, rbind, rownames, sapply, setdiff, table, tapply,
#> union, unique, unsplit, which.max, which.min
#> Loading required package: AnnotationHub
#> Loading required package: BiocFileCache
#> Loading required package: dbplyr
eh = ExperimentHub()
# Note: this requires the SummarizedExperiment package to be installed
eh_query = query(eh, c("multiWGCNAdata"))
astrocyte_se = eh_query[["EH8223"]]
#> see ?multiWGCNAdata and browseVignettes('multiWGCNAdata') for documentation
#> loading from cache
#> require("SummarizedExperiment")
# Collect the metadata in the sampleTable; the first column must be named "Sample"
sampleTable = colData(astrocyte_se)
# Check the data
assays(astrocyte_se)[[1]][1:5, 1:5]
#> EAE_1 EAE_2 EAE_3 EAE_4 EAE_5
#> 0610007P14Rik 6.8295061 6.855911 6.254815 6.018663 7.082242
#> 0610009B22Rik 6.7316000 7.131389 5.512180 5.350221 6.612084
#> 0610009O20Rik 5.2374819 4.613763 5.163960 4.499471 5.376452
#> 0610010B08Rik 0.6621082 1.237812 2.602756 1.404309 2.646006
#> 0610010F05Rik 2.7476431 3.870872 3.311788 3.329422 3.424676
sampleTable
#> DataFrame with 36 rows and 3 columns
#> Sample Disease Region
#> <character> <character> <character>
#> EAE_1 EAE_1 EAE Cbl
#> EAE_2 EAE_2 EAE Cbl
#> EAE_3 EAE_3 EAE Cbl
#> EAE_4 EAE_4 EAE Cbl
#> EAE_5 EAE_5 EAE Cbl
#> ... ... ... ...
#> healthy_12 healthy_12 WT Hippo
#> healthy_13 healthy_13 WT Sc
#> healthy_14 healthy_14 WT Sc
#> healthy_15 healthy_15 WT Sc
#> healthy_16 healthy_16 WT Sc
# Define our conditions for trait 1 (disease) and 2 (brain region)
conditions1 = unique(sampleTable[,2])
conditions2 = unique(sampleTable[,3])4 Network construction
We now perform network construction, module eigengene calculation, module-trait correlation. Let’s use power = 12 since we used this in our paper (Tommasini and Fogel. BMC Bioinformatics. 2023.) for all the networks.
# Construct the combined networks and all the sub-networks (EAE, WT, and each region)
# Same parameters as Tommasini and Fogel. BMC Bioinformatics
astrocyte_networks = constructNetworks(astrocyte_se, sampleTable, conditions1, conditions2,
networkType = "signed", TOMType = "unsigned",
power = 12, minModuleSize = 100, maxBlockSize = 25000,
reassignThreshold = 0, minKMEtoStay = 0, mergeCutHeight = 0,
numericLabels = TRUE, pamRespectsDendro = FALSE,
deepSplit = 4, verbose = 3)
This step takes a while since it performs seven network constructions, so we also provide the WGCNA object list (astrocyte_networks) in a loadable format. These were generated from the function above.
# Load pre-computed astrocyte networks
astrocyte_networks = eh_query[["EH8222"]]
#> see ?multiWGCNAdata and browseVignettes('multiWGCNAdata') for documentation
#> loading from cache
# Check one of the WGCNA objects
astrocyte_networks[["combined"]]
#> ##### datExpr #####
#> X EAE_1 EAE_2 EAE_3 EAE_4
#> 0610007P14Rik 0610007P14Rik 6.8295061 6.855911 6.254815 6.018663
#> 0610009B22Rik 0610009B22Rik 6.7316000 7.131389 5.512180 5.350221
#> 0610009O20Rik 0610009O20Rik 5.2374819 4.613763 5.163960 4.499471
#> 0610010B08Rik 0610010B08Rik 0.6621082 1.237812 2.602756 1.404309
#> 0610010F05Rik 0610010F05Rik 2.7476431 3.870872 3.311788 3.329422
#> 0610010K14Rik 0610010K14Rik 6.7580660 5.741782 7.068028 6.066397
#>
#> ##### conditions #####
#> Sample Cbl Ctx Hippo Sc EAE WT
#> 1 EAE_1 1 2 2 2 1 2
#> 2 EAE_2 1 2 2 2 1 2
#> 3 EAE_3 1 2 2 2 1 2
#> 4 EAE_4 1 2 2 2 1 2
#> 5 EAE_5 1 2 2 2 1 2
#> 6 EAE_6 2 1 2 2 1 2
#>
#> ##### module-trait correlation #####
#> Module Cbl Ctx Hippo
#> 1 combined_000 -0.0364476003508565 -0.255256505885536 0.197531123492674
#> 2 combined_001 0.288897873627807 0.224331519793051 -0.951891696612219
#> 3 combined_002 -0.173307792846142 0.192824350579657 -0.737858151950135
#> 4 combined_003 -0.602638473395467 0.358865557684427 -0.443014688492424
#> 5 combined_004 -0.298759692763236 0.237517325231285 0.389483421290234
#> 6 combined_005 0.33717807439664 0.217125827452453 0.404748752737107
#> Sc EAE WT
#> 1 0.0941729827437187 -0.129401075697594 0.129401075697594
#> 2 0.438662303191361 0.0564118630816024 -0.0564118630816024
#> 3 0.718341594216621 0.124624344929484 -0.124624344929484
#> 4 0.686787604203464 0.118489907579839 -0.118489907579839
#> 5 -0.328241053758283 -0.27623235802516 0.27623235802516
#> 6 -0.959052654586201 0.0993302279700684 -0.0993302279700684
#> p.value.Cbl p.value.Ctx p.value.Hippo
#> 1 0.832857799237777 0.132966929718593 0.248179967153461
#> 2 0.0874757049142457 0.188400881435516 4.88276900785685e-19
#> 3 0.312100684420068 0.259856986766255 2.82359321810442e-07
#> 4 0.000100501116087253 0.0316046790137704 0.00681175845919204
#> 5 0.0767322030921174 0.16305699722782 0.0188747122157827
#> 6 0.0443263856584159 0.203366848448585 0.014341658701339
#> p.value.Sc p.value.EAE p.value.WT trait log10Pvalue
#> 1 0.58485059401596 0.451944819774715 0.451944819774715 None 0.8762564
#> 2 0.00744524874705859 0.743829515380077 0.743829515380077 Hippo 18.3113338
#> 3 8.09340902303281e-07 0.468950669314544 0.468950669314544 Hippo 6.5491979
#> 4 3.73966334787861e-06 0.491276518086413 0.491276518086413 Sc 5.4271675
#> 5 0.0506450923633608 0.102930451571833 0.102930451571833 Hippo 1.7241197
#> 6 3.33227384463383e-20 0.564360643230857 0.564360643230857 Sc 19.4772593
#>
#> ##### module eigengenes #####
#> EAE_1 EAE_2 EAE_3 EAE_4 EAE_5
#> combined_000 0.02299076 0.04566290 -0.016901593 -0.04957158 0.05120053
#> combined_001 -0.12249369 -0.04708563 -0.103509765 -0.08154738 -0.10047042
#> combined_002 -0.09416519 0.11637829 0.001182349 0.02840979 0.04835513
#> combined_003 0.01507925 0.20555343 0.141335547 0.13334496 0.21012775
#> combined_004 0.21553261 -0.19353170 0.303492247 0.17313968 0.25254140
#> combined_005 -0.07089446 -0.01972901 -0.159155387 -0.09926386 -0.15912974
#>
#> ##### outlier modules #####
#> [1] "combined_000" "combined_019" "combined_025" "combined_028" "combined_030"
#> [6] "combined_035" "combined_036" "combined_037" "combined_038"5 Compare modules by overlap
Next, we compare modules (by hypergeometric overlap) across conditions. We’ll save the results in a list.
# Save results to a list
results = list()
results$overlaps = iterate(astrocyte_networks, overlapComparisons, plot=FALSE)
#>
#> #### comparing combined and EAE ####
#>
#> #### comparing combined and WT ####
#>
#> #### comparing combined and Cbl ####
#>
#> #### comparing combined and Ctx ####
#>
#> #### comparing combined and Hippo ####
#>
#> #### comparing combined and Sc ####
#>
#> #### comparing EAE and WT ####
#>
#> #### comparing EAE and Cbl ####
#>
#> #### comparing EAE and Ctx ####
#>
#> #### comparing EAE and Hippo ####
#>
#> #### comparing EAE and Sc ####
#>
#> #### comparing WT and Cbl ####
#>
#> #### comparing WT and Ctx ####
#>
#> #### comparing WT and Hippo ####
#>
#> #### comparing WT and Sc ####
#>
#> #### comparing Cbl and Ctx ####
#>
#> #### comparing Cbl and Hippo ####
#>
#> #### comparing Cbl and Sc ####
#>
#> #### comparing Ctx and Hippo ####
#>
#> #### comparing Ctx and Sc ####
#>
#> #### comparing Hippo and Sc ####
# Check the overlaps, ie between the EAE and wildtype networks
head(results$overlaps$EAE_vs_WT$overlap)
#> mod1 mod2 mod1.size mod2.size overlap p.value p.adj
#> 1 EAE_000 WT_000 1 2 0 1 1
#> 2 EAE_000 WT_001 1 1001 0 1 1
#> 3 EAE_000 WT_002 1 885 0 1 1
#> 4 EAE_000 WT_003 1 883 0 1 1
#> 5 EAE_000 WT_004 1 767 0 1 1
#> 6 EAE_000 WT_005 1 727 0 1 16 Identify a module of interest
Then, we perform differential module expression analysis to detect modules with disease-associated expression patterns. This incorporates the linear model described in the paper and tests for significance using ANOVA.
# Run differential module expression analysis (DME) on combined networks
results$diffModExp = runDME(astrocyte_networks[["combined"]],
sampleTable,
p.adjust = "fdr",
refCondition = "Region",
testCondition = "Disease")
# plot=TRUE,
# out="ANOVA_DME.pdf")
# Check results sorted by disease association FDR
results$diffModExp[order(results$diffModExp$Disease),]
#> Disease Disease*Region Region
#> combined_013 0.004565611 0.98775812 5.707380e-04
#> combined_016 0.056027353 0.02147806 2.496359e-09
#> combined_021 0.056027353 0.98775812 1.591773e-01
#> combined_033 0.056027353 0.98775812 6.349427e-06
#> combined_032 0.075974114 0.98775812 1.228102e-01
#> combined_020 0.081916454 0.98775812 1.545306e-05
#> combined_026 0.125548873 0.98775812 2.274020e-02
#> combined_005 0.129290571 0.98775812 1.985491e-18
#> combined_017 0.129290571 0.98775812 3.390485e-03
#> combined_031 0.129290571 0.98775812 8.937602e-06
#> combined_003 0.201313016 0.98775812 5.082058e-14
#> combined_004 0.201313016 0.98775812 8.237774e-03
#> combined_002 0.204080192 0.98775812 7.594394e-13
#> combined_009 0.204080192 0.98775812 1.817818e-24
#> combined_022 0.235673317 0.98775812 1.522679e-17
#> combined_010 0.260509013 0.98775812 1.278032e-05
#> combined_015 0.260509013 0.98775812 1.286140e-03
#> combined_008 0.305563496 0.98775812 9.832001e-05
#> combined_011 0.305563496 0.98775812 1.488179e-19
#> combined_012 0.364451132 0.98775812 5.490254e-17
#> combined_014 0.364451132 0.98775812 7.311012e-22
#> combined_018 0.364451132 0.98775812 5.347187e-19
#> combined_023 0.364451132 0.98775812 3.631151e-09
#> combined_025 0.364451132 0.98775812 1.769155e-01
#> combined_001 0.387743702 0.98775812 4.956367e-17
#> combined_036 0.387743702 0.98775812 4.677576e-01
#> combined_029 0.433469339 0.98775812 5.875322e-15
#> combined_037 0.435725138 0.98775812 2.151891e-01
#> combined_006 0.442761880 0.98775812 9.999069e-17
#> combined_028 0.442761880 0.98775812 3.015503e-01
#> combined_035 0.442761880 0.98775812 4.374292e-01
#> combined_007 0.500999238 0.98775812 9.767417e-05
#> combined_030 0.500999238 0.98775812 3.015503e-01
#> combined_000 0.504429894 0.98775812 4.374292e-01
#> combined_034 0.504429894 0.98775812 5.707380e-04
#> combined_038 0.515358345 0.98775812 4.037062e-01
#> combined_024 0.692301058 0.98775812 3.840878e-18
#> combined_019 0.717692761 0.98775812 4.677576e-01
#> combined_027 0.717692761 0.98775812 1.378638e-15
# You can check the expression of module M13 from Tommasini and Fogel. BMC Bioinformatics. 2023 like this. Note that the values reported in the bottom panel title are p-values and not adjusted for multiple comparisons like in results$diffModExp
diffModuleExpression(astrocyte_networks[["combined"]],
geneList = topNGenes(astrocyte_networks[[1]], "combined_013"),
design = sampleTable,
test = "ANOVA",
plotTitle = "combined_013",
plot = TRUE)
#> #### plotting combined_013 ####
#> Factors p.value
#> 1 Disease 0.0001170670
#> 2 Region 0.0003394679
#> 3 Disease*Region 0.13864880057 Draw the multiWGCNA network
We can now check to see if M13 is present in any of the sub-networks. An easy way to do this is using the network-network correspondences from hypergeometric overlap. These are stored in results$overlaps. We can plot these in a convenient visualization scheme that also organizes the three levels of the multiWGCNA analysis: 1) combined network, 2) EAE and wildtype networks, and 3) the four regional networks.
drawMultiWGCNAnetwork(astrocyte_networks,
results$overlaps,
"combined_013",
design = sampleTable,
overlapCutoff = 0,
padjCutoff = 1,
removeOutliers = TRUE,
alpha = 1e-50,
layout = NULL,
hjust = 0.4,
vjust = 0.3,
width = 0.5)
#> NULLThis corresponds to Figure 2C from Tommasini and Fogel. BMC Bioinformatics. 2023. We see that M13 is really only present in the EAE network, but not any of the other sub-networks. Most importantly, it cannot be resolved in the wildtype network. This makes M13 a biologically interesting network, both in terms of differential expression and differential co-expression.
We can identify the EAE module that corresponds to M13 using the overlap analysis:
bidirectionalBestMatches(results$overlaps$combined_vs_EAE)