Weighted correlation network analysis (WGCNA) for metabolome and transcriptome

In order to group the biochemicals that were highly correlated, we built the co-expression network using WGCNA [102]. The WGCNA is an efficient and robust method in grouping metabolomic and transcriptomic data [103, 104] and allowed us to summarize each module by its module eigenvalue. A one-sided Fisher test was used to determine if a pathway was enriched within the turquoise and blue modules in metabolomic data. P values were then adjusted using Benjamini-Hochberg method, and a cut-off of P < 0.05 and FDR adjusted-P < 0.05 were chosen to determine if a pathway was significantly enriched. We used Pearson’s correlation between expression profile of each gene and module eigenvalue to identify module membership. Using the module eigenvalue, the module-traits relationships were estimated by calculating Pearson’s correlations between the module eigenvalue and the traits of interest. We considered 0.90 as a correlation cut-off to choose soft-thresholding power and set the minimal module size as 20. For metabolome, the metabolites were clustered into 8 modules plus 43 unclustered metabolites. The transformed values of the unclustered metabolites were combined with standardized module eigenvalues in the following analysis. For the transcriptome data, 14 modules (defined as clusters of highly interconnected genes) were identified by using DynamicTree Cut algorithm. WGCNA led to 14 different modules by using DynamicTree Cut algorithm. Over-representation of genes in the blue module was characterized based on gene ontology biological process and KEGG pathways using clusterProfiler [105].