Weighted gene co‐expression Network analysis (WGCNA) to identify networks of co‐expressed transcripts

To obtain networks of co-expressed transcripts that were categorized as modules, we performed a Weighted Gene Co-expression Network Analysis (WGCNA). The counts obtained with Kallisto (v0.43.0) [79] were transformed using variance stabilizing transformation (vst) as implemented in DESeq2 (v1.18.1) [81]. The vst transformed counts were used to perform a co-expression network analysis with the R package WGCNA (v1.63) [85]. For more details on the methodology, see [85–87]. In short, (Additional File 1, Figure S1, workflow, right side), a similarity matrix was built by calculating Pearson correlations between the expression values of all pairs of transcripts. Using the similarity matrix, a signed weighted adjacency matrix was obtained as described by the formula:

Where cor_ij is the Pearson correlation between the expression pattern of transcript ‘i’ and transcript ‘j’ (the similarity value). The value of β was chosen based on the soft-thresholding approach [85]. With this value of β, we obtained a weighted network with an approximate scale-free topology (β = 14, scale-free topology R² = 0.84). In a signed weighted adjacency matrix negative and small positive correlations get negligibly small adjacency values shifting the focus on strong positive correlations. Seeing the adjacency matrix as a network, the nodes correspond to the transcripts and the connections between nodes correspond to the adjacency values (transformed correlation coefficients). A topological overlap matrix (TOM), which in addition to the adjacency matrix considers topological similarity (shared neighbors reinforce the connection strength between two nodes), was constructed using the adjacency matrix [88]. To define transcript modules, a hierarchical clustering tree was constructed using the dissimilarity measure (1-TOM). Transcript modules were defined by cutting the branches of the tree using the Dynamic Hybrid Tree Cut algorithm [89] and the minimum module size was set to 30 transcripts. Transcript modules with similar expression profiles were merged by hierarchical clustering of the eigengene correlation values. Briefly, a hierarchical clustering tree was created with the eigengene dissimilarity measure (1-correlation coefficient of eigengenes) and a tree height cut of 0.20 was used (corresponds to an eigengene cor ≥ 0.80). Eigengenes were calculated with the function moduleEigengenes (default settings) [85]. A module eigengene corresponds to the module’s first principal component and can be seen as a weighted average expression profile [85]. To find significantly associated modules with age, correlations between age and eigengenes of the merged modules were calculated. Each module was named after a color by WGCNA.

The adjacency matrix of the WGCNA was visualized using Cytoscape (v3.7.1) [90], only including pairs of nodes with a cor_ij ≥ 0.90. Each module’s color corresponds with the respective module name (e.g., saddlebrown color for the saddlebrown module).

To identify co-expression modules containing age-related DETs, we looked for age-related DETs from the Iso-MaSigPro analysis in the WGCNA modules. Those modules that were significantly correlated with age and held the highest number of Iso-MaSigPro DETs were further inspected.