Integration of transcriptome and machine learning to identify the potential key genes and regulatory networks affecting drip loss in pork

Abstract

Low level of drip loss (DL) is an important quality characteristic of meat with high economic value. However, the key genes and regulatory networks contributing to DL in pork remain largely unknown. To accurately identify the key genes affecting DL in muscles postmortem, 12 Duroc × (Landrace × Yorkshire) pigs with extremely high (n = 6, H group) and low (n = 6, L group) DL at both 24 and 48 h postmortem were selected for transcriptome sequencing. The analysis of differentially expressed genes and weighted gene co-expression network analysis (WGCNA) were performed to find the overlapping genes using the transcriptome data, and functional enrichment and protein-protein interaction (PPI) network analysis were conducted using the overlapping genes. Moreover, we used machine learning to identify the key genes and regulatory networks related to DL based on the interactive genes of the PPI network. Finally, nine potential key genes (IRS1, ESR1, HSPA6, INSR, SPOP, MSTN, LGALS4, MYLK2, and FRMD4B) mainly associated with the MAPK signaling pathway, the insulin signaling pathway, and the calcium signaling pathway were identified, and a single-gene set enrichment analysis (GSEA) was performed to further annotate the functions of these potential key genes. The GSEA results showed that these genes are mainly related to ubiquitin-mediated proteolysis and oxidative reactions. Taken together, our results indicate that the potential key genes influencing DL are mainly related to insulin signaling mediated differences in glycolysis and ubiquitin-mediated changes in muscle structure and improve the understanding of gene expression and regulation related to DL and contribute to future molecular breeding for improving pork quality.

Keywords: RNA-seq; WGCNA; drip loss; machine learning; meat quality; single-gene GSEA.

Figure 1.

Identification and verification of DEGs. (A) PCA analysis of H and L groups. The samples in H group and L group were clearly clustered into two branches, indicating significant differences between the H and L groups and good reproducibility within the H or L group.(B) Volcano plot of DEGs. Red and blue represent significantly up- and downregulated genes, respectively (FDR < 0.05). (C) The heat maps of DEGs in H and groups. Each row represents a gene, and each column represents a sample. Identification of DEGs and Hierarchical cluster analysis. (D) Verification of RNA-seq data by qRT‒PCR. X-axis represents 10 selected genes, and Y-axis represents the expression levels of genes from RNA-seq and qRT-PCR. (E) Linear regression analysis of expression levels between RNA-seq and qRT‒PCR data. X-axis represents the log2foldchange of RNA-seq, and Y-axis indicates the log2foldchange of qRT-PCR.

Figure 2.

Weighted gene co-expression network analysis (WGCNA). (A) Analysis of network topology for various soft-thresholding powers. The left part of (A) shows the scale-free fit index (Y-axis) as a function of the soft-thresholding power (X-axis) and the right part of (A) represents the mean connectivity (degree, Y axis) as a function of the soft-thresholding power (X-axis). (B) Gene dendrogram of co-expression modules. It was obtained by clustering the dissimilarity based on consensus topological overlap, and the corresponding module colors indicated by the color row and the highly related gene modules are shown in different colors. (C) Module-trait associations. Each row corresponds to a module eigengene, and each column represents a trait; each cell contains the corresponding correlation coefficient (top) and P-value (bottom), and the red and blue represent positive and negative correlation, respectively.

Figure 3.

Function enrichment analysis of intersection genes. (A) Venn diagrams of DEGs and the genes derived from WGCNA. (B) GO enrichment analysis. The X-axis shows the number of DEGs, while the Y-axis represents the GO terms. The bar color corresponds to different GO categories, and the blue for biological process, the green for cellular component and the red for molecular function. (C) KEGG pathway enrichment analysis. X-axis shows the rich ratio of the intersection genes, while the Y-axis represents the KEGG pathways. The size of the bubble represents the number of intersection genes annotated in each pathway, and the color represents the P-value.

Figure 4.

Identification of potential key genes. (A) PPI network. The genes in the network are the intersection genes derived from the DEGs and the WGCNA. (B) Venn diagrams of intersection genes derived from the top 10 genes for 12 algorithms in cytoHubba plugin. (C) LASSO regression algorithm. The left part of (C) represents the confidence interval of each lambda, and the horizontal axis shows the logarithm of the lambda, and the vertical axis shows the mean-squared error. The right part of (C) represents the distribution of the LASSO coefficient. Each color line shows the changing tendency of each gene coefficient chosen by the LASSO algorithm, and the horizontal axis shows the log value of lambda, the vertical axis shows the coefficient corresponding to lambda, and the numeral on the uppermost axis shows several genes whose coefficient is not zero at different log lambda values. (D) SVM-RFE algorithm. In the left part of (D), the horizontal axis represents the number of critical genes, and the vertical axis represents the generalization error under 10-fold cross-validation. In the right part of (D), the trend of the line graph represents the relationship between the number of critical genes and the generalization error. (E) Venn diagram represents the number of intersection genes extracted from LASSO and SVM-RFE algorithms.

Figure 5.

GSEA of potential key genes. (A) The function annotation for ESR1 gene. (B) The function annotation for INSR gene, (C) the function annotation for IRS1 gene. (D) The function annotation for LGALS4 gene. (E) The function annotation for MSTN gene. (F) The function annotation for MYLK2 gene. (G) The function annotation for SPOPL gene. (H) The function annotation for FRDM4B gene. (I) The function annotation for HSPA6 gene.

Figure 6.

Correlations among the expression of the potential key genes and drip loss phenotypes.