2.4. Weighted Gene Correlation Network Analysis (WGCNA)

The genes whose expression was not detected in most of the samples or with variance zero were filtered out from the normalized gene count matrix. To identify strong co-expression networks of genes related to SM FA composition, we employed a co-expression analysis approach using the “WGCNA” package [33] in the R environment [32]. The gene count matrix was analyzed together with the IMF contents of palmitic, palmitoleic, stearic, oleic, n-6, and n-3 FAs. These FAs were chosen for their relative abundance in IMF (palmitic, palmitoleic, stearic, and oleic are the saturated FAs (SFAs) and monounsaturated FAs (MUFAs) most represented in pork fat tissues [6,11]) and for their nutritional value (n-6 and n-3 PUFAs).

To obtain the scale-free undirected co-expression networks between the genes in the gene count matrix, an adjacency matrix was built using Pearson’s correlations between each gene couple. The soft threshold power values were estimated through the function pickSoftThreshold() in the WGCNA package, and a soft threshold power (β) of 6 was used to raise the adjacency matrix. This value was chosen as the scale-free topology index (R²) reached the peak (R² > 0.70) for the first time when β = 6, and the minimum module size was 30 genes (Supplementary Figure S1). After β was determined, the adjacency matrix was calculated using the topological overlap measure (TOM) and the corresponding dissimilarity (dissTOM = 1 − TOM). The latter was used as a distance for gene hierarchical cluster, and then DynamicTree Cut algorithm [33] was used to identify the modules of genes. The default minimum cluster merge height of 0.25 was retained. After the construction of the gene co-expression network, WGCNA restituted a list of gene modules, which were named using color labels and clustered highly interconnected genes. The principal component of each module was defined as the module eigengene (ME); MEs represented the expression value of each module and were used to detect modules that may comprise genes having a biologically relevant role in the variations of the traits. The module-trait relationship (module membership, MM) was calculated as the Pearson’s correlation between the ME and the traits of interest. For each gene, the gene MM and the relative p-values were obtained and carefully evaluated, as they indicated the importance of that gene in a module. Due to the high number of gene modules showing a significant correlation with the analyzed traits, we decided to choose the gene modules for further analysis based on the threshold reported by Pampouille et al. [34]. Based on this study, only the gene modules with an absolute value of module-trait correlation higher than 0.7 were further analyzed. After selecting the gene modules exceeding the set threshold, the genes participating in those modules with a MM p-value < 0.001 were further considered to perform the functional enrichment analysis. The intramodular connectivity and gene significance (i.e., a coefficient and relative significance value indicating the importance of that gene for the trait, based on the correlation between the gene expression profile and the trait) were used to identify key genes in the networks.