Network comparison: consensus modules and module preservation
2 October 2024
Contents
1 Installation
if(!requireNamespace('BiocManager', quietly = TRUE))
install.packages('BiocManager')
BiocManager::install("BioNERO")# Load package after installation
library(BioNERO)
set.seed(12) # for reproducibility2 Introduction
Comparing different coexpression networks can reveal relevant biological
patterns. For instance, seeking for consensus modules can identify
coexpression modules that occur in all data sets regardless of natural
variation and, hence, are core players of the studied phenotype.
Additionally, module preservation within and across species can reveal
patterns of conservation and divergence between transcriptomes.
In this vignette, we will explore consensus modules and module preservation
analyses with BioNERO. Although they seem similar, their goals are opposite:
while consensus modules identification focuses on the commonalities, module
preservation focuses on the divergences.
3 Data loading and description
We will use RNA-seq data of maize (Zea mays) and rice (Oryza sativa) obtained from Shin et al. (2020).
data(zma.se)
zma.se
## class: SummarizedExperiment
## dim: 10802 28
## metadata(0):
## assays(1): ''
## rownames(10802): ZeamMp030 ZeamMp044 ... Zm00001d054106 Zm00001d054107
## rowData names(0):
## colnames(28): SRX339756 SRX339757 ... SRX2792103 SRX2792104
## colData names(1): Tissue
data(osa.se)
osa.se
## class: SummarizedExperiment
## dim: 7647 27
## metadata(0):
## assays(1): ''
## rownames(7647): Os01g0100700 Os01g0100900 ... Os12g0641400 Os12g0641500
## rowData names(0):
## colnames(27): SRX831140 SRX831141 ... SRX263041 SRX1544234
## colData names(1): TissueAll BioNERO’s functions for consensus modules and module preservation
analyses require the expression data to be in a list. Each element of the
list can be a SummarizedExperiment object (recommended) or an expression data
frame with genes in row names and samples in column names.
4 Consensus modules
The most common objective in consensus modules identification is to find core modules across different tissues or treatments for the same species. For instance, one can infer GCNs for different types of cancer in human tissues (say prostate and liver) and identify modules that occur in all sets, which are likely core components of cancer biology. Likewise, one can also identify consensus modules across samples from different geographical origins to find modules that are not affected by population structure or kinship.
4.1 Data preprocessing
Here, we will subset 22 random samples from the maize data twice and find consensus modules between the two sets.
# Preprocess data and keep top 2000 genes with highest variances
filt_zma <- exp_preprocess(zma.se, variance_filter = TRUE, n = 2000)
## Number of removed samples: 1
# Create different subsets by resampling data
zma_set1 <- filt_zma[, sample(colnames(filt_zma), size=22, replace=FALSE)]
zma_set2 <- filt_zma[, sample(colnames(filt_zma), size=22, replace=FALSE)]
colnames(zma_set1)
## [1] "SRX3804716" "SRX2792102" "SRX2527287" "SRX2792108" "SRX3804723"
## [6] "SRX2792107" "SRX2792103" "SRX2792104" "SRX3804715" "SRX339808"
## [11] "SRX339758" "SRX339756" "SRX339809" "SRX2792105" "SRX339757"
## [16] "SRX2641029" "SRX339762" "SRX2792110" "SRX339807" "SRX3804718"
## [21] "ERX2154032" "SRX339764"
colnames(zma_set2)
## [1] "SRX2792111" "SRX339756" "SRX3804715" "SRX2792108" "SRX2792103"
## [6] "SRX339762" "SRX339807" "ERX2154030" "SRX2792105" "SRX2792107"
## [11] "SRX2527287" "ERX2154032" "SRX2792104" "SRX2527288" "SRX339809"
## [16] "SRX3804718" "SRX3804716" "SRX2792109" "SRX3804723" "SRX2792102"
## [21] "SRX339808" "SRX339758"
# Create list
zma_list <- list(set1 = zma_set1, set2 = zma_set2)
length(zma_list)
## [1] 24.2 Identification of consensus modules
As described in the first vignette, before inferring the GCNs, we need to
identify the optimal \(\beta\) power that makes the network closer to a scale-free
topology. We can do that with consensus_SFT_fit().
cons_sft <- consensus_SFT_fit(zma_list, setLabels = c("Maize 1", "Maize 2"),
cor_method = "pearson")
## Power SFT.R.sq slope truncated.R.sq mean.k. median.k. max.k.
## 1 5 0.658 -0.650 0.656 107.00 90.50 304.0
## 2 6 0.761 -0.706 0.774 79.70 62.90 244.0
## 3 7 0.833 -0.757 0.863 61.30 45.10 200.0
## 4 8 0.825 -0.830 0.868 48.30 32.80 169.0
## 5 9 0.822 -0.916 0.887 38.70 24.50 145.0
## 6 10 0.838 -0.975 0.910 31.50 18.90 126.0
## 7 11 0.824 -1.040 0.914 26.00 15.00 111.0
## 8 12 0.837 -1.090 0.932 21.70 12.20 97.6
## 9 13 0.851 -1.130 0.944 18.30 10.10 86.6
## 10 14 0.842 -1.190 0.941 15.50 8.42 77.3
## 11 15 0.844 -1.230 0.947 13.30 7.17 69.4
## 12 16 0.842 -1.270 0.951 11.50 6.13 62.5
## 13 17 0.855 -1.300 0.962 9.99 5.24 56.5
## 14 18 0.864 -1.320 0.967 8.73 4.54 51.3
## 15 19 0.869 -1.330 0.970 7.67 3.92 46.7
## 16 20 0.875 -1.340 0.977 6.78 3.38 42.7
## Power SFT.R.sq slope truncated.R.sq mean.k. median.k. max.k.
## 1 5 0.592 -0.621 0.574 94.70 82.40 265.0
## 2 6 0.702 -0.666 0.710 69.90 56.80 209.0
## 3 7 0.757 -0.713 0.771 53.20 40.30 167.0
## 4 8 0.829 -0.729 0.853 41.50 29.50 136.0
## 5 9 0.855 -0.782 0.884 32.90 22.10 113.0
## 6 10 0.868 -0.826 0.911 26.60 17.20 96.6
## 7 11 0.873 -0.868 0.931 21.80 13.70 83.2
## 8 12 0.837 -0.943 0.914 18.10 11.30 73.0
## 9 13 0.826 -1.000 0.917 15.20 9.30 64.7
## 10 14 0.815 -1.070 0.917 12.80 7.73 57.6
## 11 15 0.810 -1.120 0.926 11.00 6.53 51.5
## 12 16 0.823 -1.150 0.943 9.42 5.45 46.3
## 13 17 0.815 -1.210 0.942 8.16 4.57 41.8
## 14 18 0.827 -1.240 0.954 7.11 3.92 37.8
## 15 19 0.833 -1.260 0.964 6.23 3.42 34.3
## 16 20 0.831 -1.280 0.965 5.49 2.99 31.2This function returns a list with the optimal powers and a summary plot,
exactly as SFT_fit() does.
powers <- cons_sft$power
powers
## set1 set2
## 7 8
cons_sft$plot
Now, we can infer GCNs and identify consensus modules across data sets.
consensus <- consensus_modules(zma_list, power = powers, cor_method = "pearson")
## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.
## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.
## ..done.
## multiSetMEs: Calculating module MEs.
## Working on set 1 ...
## Working on set 2 ...
names(consensus)
## [1] "consMEs" "exprSize" "sampleInfo"
## [4] "genes_cmodules" "dendro_plot_objects"
head(consensus$genes_cmodules)
## Genes Cons_modules
## 1 ZeamMp030 lightcyan
## 2 ZeamMp044 grey60
## 3 ZeamMp092 grey60
## 4 ZeamMp108 blue
## 5 ZeamMp116 grey60
## 6 ZeamMp158 blueFinally, we can correlate consensus module eigengenes to sample metadata (here, plant tissues).1 NOTE: Blank grey cells in the heatmap represent correlation values that have opposite sign in the expression sets. For each correlation pair, consensus correlations are calculated by selecting the minimum value across matrices of consensus module-trait correlations.
consensus_trait <- consensus_trait_cor(consensus, cor_method = "pearson")
head(consensus_trait)
## trait ME cor pvalue group
## 1 endosperm MEblack NA NA Tissue
## 2 endosperm MEgrey60 0.05143433 0.82244990 Tissue
## 3 endosperm MEdarkred 0.03438255 0.88082052 Tissue
## 4 endosperm MEblue -0.09885816 0.66550539 Tissue
## 5 endosperm MElightcyan 0.09297754 0.68440932 Tissue
## 6 endosperm MEdarkturquoise -0.45285214 0.03330611 TissueAs with the output of module_trait_cor(), users can visualize the output
of consensus_trait_cor() with the function plot_module_trait_cor().
plot_module_trait_cor(consensus_trait)
5 Module preservation
Module preservation is often used to study patterns of evolutionary conservation and divergence across transcriptomes, an approach named phylotranscriptomics. This way, one can investigate how evolution shaped the expression profiles for particular gene families across taxa.
5.1 Data preprocessing
To calculate module preservation statistics, gene IDs must be shared by
the expression sets. For intraspecies comparisons, this is an easy task,
as gene IDs are the same. However, for interspecies comparisons, users need
to identify orthogroups between the different species and collapse the
gene-level expression values to orthogroup-level expression values.
This way, all expression sets will have common row names. We recommend
identifying orthogroups with OrthoFinder (Emms and Kelly 2015), as it is simple
to use and widely used.2 PRO TIP: If you identify orthogroups with OrthoFinder,
BioNERO has a helper function named parse_orthofinder() that parses
the Orthogroups.tsv file generated by OrthoFinder into a suitable data frame
for module preservation analysis. See ?parse_orthofinder for more details. Here, we will compare maize and rice expression
profiles. The orthogroups between these species were downloaded from the
PLAZA 4.0 Monocots database (Van Bel et al. 2018).
data(og.zma.osa)
head(og.zma.osa)
## Family Species Gene
## 1548 ORTHO04M000001 osa Os01g0100700
## 1549 ORTHO04M000001 osa Os01g0100900
## 4824 ORTHO04M000001 zma Zm00001d009743
## 4854 ORTHO04M000001 zma Zm00001d020834
## 4874 ORTHO04M000001 zma Zm00001d026672
## 4921 ORTHO04M000001 zma Zm00001d039873As you can see, the orthogroup object for BioNERO must be a data frame with
orthogroups, species IDs and gene IDs, respectively. Let’s collapse gene-level
expression to orthogroup-level with exp_genes2orthogroups(). By default, if
there is more than one gene in the same orthogroup for a given species, their
expression levels are summarized to the median. Users can also summarize
to the mean.
# Store SummarizedExperiment objects in a list
zma_osa_list <- list(osa = osa.se, zma = zma.se)
# Collapse gene-level expression to orthogroup-level
ortho_exp <- exp_genes2orthogroups(zma_osa_list, og.zma.osa, summarize = "mean")
# Inspect new expression data
ortho_exp$osa[1:5, 1:5]
## SRX831140 SRX831141 SRX831137 SRX831138 SRX831134
## ORTHO04M000001 6.909420 7.258330 94.20870 92.85195 123.01060
## ORTHO04M000002 9.203498 8.709974 66.45512 44.97913 33.86936
## ORTHO04M000003 9.417930 9.444861 42.57513 66.02237 55.37741
## ORTHO04M000004 9.019436 8.920091 96.22074 62.56506 109.32262
## ORTHO04M000005 40.845040 41.844234 52.33474 31.31474 22.42236
ortho_exp$zma[1:5, 1:5]
## SRX339756 SRX339757 SRX339758 SRX339762 SRX339763
## ORTHO04M000001 26.02510 15.07917 14.91571 13.82989 8.080476
## ORTHO04M000002 19.28281 13.94254 13.57854 12.96032 13.522525
## ORTHO04M000003 45.17294 48.63796 54.22404 42.12135 10.779117
## ORTHO04M000004 28.05475 38.53734 39.48070 27.13272 2.978207
## ORTHO04M000005 67.58868 34.87009 21.46280 12.79565 7.452068Now, we will preprocess both expression sets and keep only the top 1000 orthogroups with the highest variances for demonstration purposes.
# Preprocess data and keep top 1000 genes with highest variances
ortho_exp <- lapply(ortho_exp, exp_preprocess, variance_filter=TRUE, n=1000)
## Number of removed samples: 2
## Number of removed samples: 2
# Check orthogroup number
sapply(ortho_exp, nrow)
## osa zma
## 1000 10005.2 Calculating module preservation statistics
Now that row names are comparable, we can infer GCNs for each set. We will do that iteratively with lapply.
# Calculate SFT power
power_ortho <- lapply(ortho_exp, SFT_fit, cor_method="pearson")
## Power SFT.R.sq slope truncated.R.sq mean.k. median.k. max.k.
## 1 3 0.15600 1.740 0.846 191.00 192.00 249.0
## 2 4 0.12400 1.080 0.844 129.00 130.00 178.0
## 3 5 0.10300 0.784 0.817 91.70 92.90 134.0
## 4 6 0.06140 0.528 0.818 67.80 68.90 106.0
## 5 7 0.04550 0.397