A flow chart illustrating the data preparation, processing, and analysis is displayed in Figure 1. We used ‘bronchopulmonary dysplasia’ as the key word to search the Gene Expression Omnibus (GEO) database and to select datasets containing samples from different pathological stages and normal controls. Finally, the dataset GSE32472¹, was found to meet our requirements and was therefore downloaded. To identify the molecular networks and hub genes related to the pathological progress of BPD, WGCNA was conducted. GSE32472 provided microarray profiles of blood samples of newborns with BPD, including microarray assessment of gene expression at approximately the 5th, 14th, and 28th days of life. To ensure the stability of the selection, 100 blood samples at about the 28th day were selected, when a more definite diagnosis of BPD had been made. These samples were obtained from 38 controls, and 38 mild, 10 moderate, and 14 severe BPD cases. The expression data were normalized using quantile normalization function in limma package of R software (Ritchie et al., 2015). The genes with the highest variance in expression values (top 25%) were selected for co-expression network construction. Cluster analysis using the Pearson’s correlation matrices and the average linkage method were conducted to detect whether outlier samples existed for the purpose of ensuring the reliability of the network construction. A brief design of the study is shown in Figure 1A.

A co-expression network was constructed using the WGCNA algorithm package in R (Presson et al., 2008). First, the Pearson’s correlation matrices were constructed for all pair-wise genes. Next, a weighted adjacency matrix was constructed by using the power function amn = | Cmn| ^β (Cmn = Pearson correlation between gene m and gene n; a_mn = adjacency between gene m and gene n) (Zhang and Horvath, 2005). The parameter β served as a soft threshold parameter to expand strong correlations and penalize weak correlations between genes. To ensure a scale-free topology of the network, β was selected when the scale independence value was equal to 0.9. The adjacency was transformed into a topological overlap matrix (TOM) to measure the network connectivity of a gene, which is defined as the sum of its adjacency with all other genes. Hierarchical clustering was performed according to TOM-based dissimilarity to distribute genes with similar expression patterns into modules with a minimum cluster size of 50 (Ravasz et al., 2002). Highly similar modules were merged with a cut-off of 0.25.

Module eigengenes (MEs) were considered the major component in the principal component analysis for each gene module, and the expression patterns of all genes could be summarized into a single characteristic expression profile within a given module. To identify modules significantly associated with BPD severity, the correlation between MEs and BPD stage was evaluated by the Pearson’s correlation test with p < 0.05 as the cut-off. The modules most significantly related to BPD severity were considered as key modules and subjected to further analysis.

A module hub gene is a highly connected in-module gene that has the highest module member (MM) score of its corresponding module (Horvath and Dong, 2008). The MM score for every gene was calculated by the WGCNA function KME, which correlates the expression profile of a gene with the ME of a module to quantify the relationship between a gene and a given module. The absolute value of gene significance (GS) represents the Pearson’s correlation between a given gene and clinical features. We removed hub genes based on the cut-off criteria (| MM| ≥ 0.85, | GS| ≥ 0.45). Further, all genes in key modules were uploaded to STRING² to acquire information about the interaction between genes. Protein–protein interaction (PPI) networks were constructed with the species limited to ‘Homo sapiens’ and a confidence > 0.9. In the PPI network, genes with a degree ≥ 10 were defined as hub nodes. Hub genes common in both co-expression network and PPI network were selected for candidate hub genes identification.

Hub genes common in both co-expression network and PPI network were analyzed by ROC curve, and area under the curve (AUC) was calculated to distinguish the control group from the BPD group. In addition, one-way ANOVA and Pearson’s correlation tests were conducted to explore the relevance of the hub genes common in both co-expression network and PPI network in terms of disease severity. Candidate hub genes were identified using the following criteria: (1) a significant P value in the one-way ANOVA test and the Pearson’s correlation and (2) an AUC > 0.8.

To validate the candidate hub genes obtained by WGCNA, pre-term infants with or without BPD blood samples were collected from the Department of Neonatology of the Sixth Affiliated Hospital of Sun Yat-sen University. This research was approved by the ethics review board of Sixth Affiliated Hospital of Sun Yat-sen University (2019ZSLYEC-80), and written informed consent was provided by the participants’ legal guardians. From each sample, 100 ng of cDNA was obtained for RT-PCR amplification reaction, and the expression of an endogenous control (housekeeping gene: GAPDH) was used for the determination of the relative expression levels of the hub genes. Primer sequences for related hub genes are listed in Supplementary Table S1. Real hub genes were identified if the results of RT-PCR have significant difference.

To gain further insights into the functions of hub genes in the module most related to BPD, we performed biological process analysis and KEGG pathway enrichment analysis with ‘c2.cp.kegg.v7.1.symbols’ as background³.

Mapping to KEGG (Kyoto Encyclopedia of Genes and Genomes) database⁴, GSEA⁵ (Subramanian et al., 2007) was performed between control and BPD groups to confirm the expression pattern of hub biological pathways.

Non-parametric tests or t-tests based on data distribution characteristics were used to analyze the statistical significance of the difference in hub gene expression levels between the two groups. Analyses were conducted in GraphPad Prism 8.0.2. P < 0.05 was considered statistically significant.

Figure 1 shows the flow chart of data preparation, processing, analysis, and validation for this study. The data were normalized using the limma package of R software (Figure 2A and R code in Supplementary Table S2). The co-expression analysis included 100 samples with clinical information, sample information, and expression matrix. Input files are provided in Supplementary Tables S3, S4. Sample clustering was performed based on Pearson’s correlation matrices and the average linkage method. No outliers were detected (Figure 2B). The genes showing the highest expression variance (top 25%) were selected for subsequent WGCNA using the WGCNA package in R software. Genes with similar expression patterns were then grouped by average linkage hierarchical clustering. In our study, β = 23 (scale-free R² > 0.901) was selected as the soft threshold to ensure a scale-free network (Figures 3A,B). Next, we constructed a systematic clustering tree using the WGCNA package. In Figure 3C, each short vertical line represents a gene, and each color represents one module composed of genes with similar expression patterns. The genes shown in gray were the genes that could not be merged. A total of 10 modules were identified (Figure 3C).

We tested the relevance of each module for BPD clinical information, focusing on different BPD stages. As displayed in Figure 3D, the yellow module (P = 5e-08, R² = 0.51) was most significantly and positively correlated with BPD severity, whereas the red module (P = 3e-11, R² = −0.60) showed the opposite result. The correlation between the yellow module and BPD severity gradually increased and finally became positive. The red module showed the opposite pattern. Based on the above findings, the red and yellow modules were identified as key modules correlated with BPD severity and were thus, further analyzed.

PPI networks of the red and yellow modules were constructed with a cutoff confidence > 0.9 (Figures 4A,B). A total of 31 genes in the red module and 41 genes in the yellow module were identified with a degree ≥ 10 as hub genes in each PPI network. Based on | MM| ≥ 0.85 and | GS| ≥ 0.45, a total of 76 genes in the red module and 77 genes in the yellow module were selected as hub genes in each co-expression network (Figures 4C,D). A total of 21 genes in the red module and 13 genes in the yellow module were identified in both the PPI and co-expression networks (Figures 4E,F). All GS, MM, and intramodule connectivity values of each identified module are listed in Supplementary Tables S5–S7.

All the 21 candidate hub genes in the red module showed significance in the one-way ANOVA. A total of nine genes had an AUC ≥ 0.80, and 19 showed a significant correlation with disease severity in the Pearson’s correlation analysis. Eventually, nine genes in the red module with an AUC ≥ 0.80 and significant P values in the Pearson’s correlation as well as one-way ANOVA were regarded as candidate hub genes (Figure 4G). Similarly, of the 13 candidate hub genes in the yellow module, 12 showed significance in the one-way ANOVA, 11 genes had an AUC ≥ 0.80, and 10 showed a significant correlation with disease severity in the Pearson’s correlation analysis. Nine genes in the yellow module had an AUC ≥ 0.80 and significant P values in the Pearson’s correlation as well as in the one-way ANOVA, and were thus, selected as candidate hub genes (Figure 4H). Detailed information about the red and yellow modules in relation to the Pearson’s correlation, ROC, and one-way ANOVA has been provided in Supplementary Tables S8–S11. The severity plot for the candidate hub genes is shown in Figure 5A. The expression levels of candidate hub genes in the yellow module increased with disease severity, and the expression levels of candidate hub genes in this module were significantly increased in different BPD severity conditions compared with those of normal controls. In contrast, candidate hub genes in the red module showed decreasing expression levels with greater disease severity, and markedly decreased expression levels in different BPD severity conditions (Figure 5B). To further clarify the clinical significance and identify real hub genes, we collected the BPD patient’s blood for qRT-PCR validation in vitro. The results showed that most of these genes had statistically significant differences and were considered as real hub genes, except for MAPK14, CEACAM3, CSF2RB, and CD3G (Figure 6).

To learn more about the function of the identified hub genes, they were subjected to perform the biological process and KEGG pathway enrichment analyses. Real hub genes in the red module, which exhibited a negative correlation with disease severity, were significantly enriched in 20 BPs: T cell activation, positive regulation of leukocyte cell-cell adhesion, positive regulation of cell-cell adhesion, regulation of leukocyte cell-cell adhesion, leukocyte cell-cell adhesion, positive regulation of cell adhesion, regulation of cell-cell adhesion, positive regulation of T cell activation, regulation of cell-cell adhesion, positive regulation of T cell activation, positive regulation of lymphocyte activation, positive regulation of leukocyte activation, regulation of T cell activation, positive regulation of cell activation, regulation of lymphocyte activation, regulation of leukocyte activation, T cell differentiation, lymphocyte differentiation, interleukin-4 production, T cell differentiation in thymus, positive regulation of T cell differentiation in thymus, and regulation of cell-cell adhesion mediated by integrin (Figure 7A). The real hub genes were also enriched in three KEGG pathways: T cell receptor signaling pathway, primary immunodeficiency, and hematopoietic cell lineage (Figure 7B). The real hub genes in the yellow module, which showed a positive correlation with disease severity, were enriched in 4 BPs, neutrophil degranulation, neutrophil activation involved in immune response, neutrophil mediated immunity, and neutrophil activation (Figure 7C). KEGG pathway enrichment analysis showed that real hub genes in the yellow module were enriched in the bladder cancer pathway, and leukocyte transendothelial migration (Figure 7D).

As bladder cancer pathway is not related to BPD, it was not considered in GSEA confirmation. According to the results of the GSEA, the terms ‘T cell receptor signaling pathway’ and ‘primary immunodeficiency’ were significantly enriched in the control group while the term ‘hematopoietic cell lineag’ was not (Figures 8A,B and Supplementary Table S12). On the contrary, the term ‘leukocyte transendothelial migration’ was significantly enriched in the BPD group (Figure 8C and Supplementary Table S13). These results successfully confirmed the expression pattern of hub biological pathways.