Clinical data from 74 cases diagnosed with NEC and discharged from the Children’s Hospital of Soochow University between 1 June 2017, and 1 June 2022, were collected. These cases were matched 1:1 with non-NEC cases from the same period based on gestational age (±3 days) and weight (±100 g), resulting in the formation of NEC and non-NEC groups. The retrospective collection of case data included laboratory indicators, maternal pregnancy indicators, basic information of the infants upon admission, comorbidities, and treatment status. Exclusion criteria included: (1) clear presence of genetic metabolic diseases and chromosomal abnormalities; (2) severe congenital structural malformations; (3) refusal by the infant’s guardian to participate in this study. Diagnostic criteria and definitions: NEC diagnosis and staging were based on the modified Bell’s staging criteria (Patel et al., 2020). Diagnoses of SGA, RDS, PDA, and sepsis were referenced from Avery’s Diseases of the Newborn (Mu and Wang, 2022). Early-onset sepsis (EOS) was defined as sepsis occurring within 72 h of birth, while late-onset sepsis (LOS) was defined as sepsis occurring after 72 h of birth (Glaser et al., 2024). To establish the temporal sequence and minimize reverse causality, we restricted the analysis to sepsis episodes that occurred before the onset of NEC; episodes occurring after NEC onset were excluded. Maternal underlying diseases and comorbidities were all clearly diagnosed in the hospital. This study was approved by the Ethics Committee of the Children’s Hospital of Soochow University (Ethics No. 2023CS130) and conforms to the ethical standards of the 1964 Declaration of Helsinki and its later amendments or comparable ethical standards. Informed consent was obtained from the guardians of the infants, who agreed to the use and disclosure of their clinical data.

For NS, we selected the GSE25504 dataset, which contains 44 NS samples and 44 healthy control samples derived from infant blood. For NEC, we selected the GSE46619 dataset, including 5 NEC and 5 healthy control samples from intestinal tissues. Both datasets were obtained from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) of the National Center for Biotechnology Information (NCBI).

The selection of public datasets was guided by the following criteria: (1) human neonatal samples; (2) clearly diagnosed NEC or sepsis cases with corresponding healthy controls; (3) samples derived from clinically accessible sources (intestinal tissue or blood); (4) datasets with good data quality and standardized preprocessing; and (5) sufficient sample size to allow for differential expression analysis. Following a systematic search of the GEO database, GSE297483, GSE69686, GSE25504 and GSE46619 met all of these inclusion criteria.

For the selected transcriptome data, gene symbols were mapped according to their respective platforms. When multiple probes matched a single gene, the median expression value was used. The expression matrix was normalized using the log2 (X + 1) transformation. After quality control, quantile normalization was performed with the normalizeBetweenArrays function in the limma package to ensure comparable distributions across samples and reduce technical variability.

DEGs analysis was performed in the GSE46619 and GSE25504 datasets using the limma package (Robinson et al., 2009; Ritchie et al., 2015) (for differential expression analysis in RNA sequencing and microarray studies), with a cutoff criterion of P. adj.value <0.05 and |LogFC| > 1. The Benjamini–Hochberg procedure was used to adjust p-values for multiple testing. Overlapping DEGs with the same direction in NS and NEC diseases were identified using a Venn diagram tool. After intersecting the differential genes from both disease datasets, 70 overlapping genes were obtained, of which 69 were upregulated and 1 was downregulated.

GO (Gene Ontology) and KEGG (Kyoto Encyclopedia of Genes and Genomes) enrichment analyses were performed for common driving genes using the clusterProfiler package (Yu et al., 2012) (an R package for comparing biological themes among gene clusters). GO was used to annotate gene biological processes, molecular functions, and cellular components. Gene pathways were annotated through KEGG. Enrichment was considered statistically significant when the P-value was less than 0.05.

Overlapping genes were imported into the STRING database (http://string-db.org) (Franceschini et al., 2013) to construct a PPI network with complex interaction relationships (combined score >0.4), which was visualized in Cytoscape (version 3.8.1). Hub genes were identified by scoring differential genes using the degree method and the cytoHubba plugin (Chin et al., 2014). Enrichment analysis and co-expression network analysis of hub genes were performed using the “clusterProfiler” package and GeneMANIA (http://www.genemania.org/) (Warde-Farley et al., 2010), respectively. Molecular Complex Detection (MCODE) plugin in Cytoscape was used to deconstruct functional modules, with selection criteria: degree cutoff of 2, K-core of 2, node score cutoff of 0.2, and maximum depth of 100. Datasets GSE297483 and GSE69686 were used for external validation.

Based on the results of external validation, these core genes were successfully identified. To evaluate the diagnostic value of each core gene and multiple genes in NEC and NS, receiver operating characteristic (ROC) curve analyses were performed separately. The area under the ROC curve (AUC) in the GSE297483 and GSE69686 datasets was used to quantify the diagnostic ability of the core genes. The “pROC” R package was used to generate ROC curves (Robin et al., 2011) larger AUC value indicates stronger discriminative ability of the model. In our study, all validated core genes achieved AUC values above 0.7, which suggests at least moderate discriminative ability.

Transcription regulatory relationships unraveled by sentence-based text mining (TRRUST) were used to obtain candidate transcription factors (TFs) regulating core genes (Han et al., 2017). This database contains rich information about TFs associated with target genes and their regulatory relationships with TFs. We constructed a TF-mRNA regulatory network and visualized it using Cytoscape. For these TFs, we performed internal validation (datasets GSE46619 and GSE25504).

Continuous variables were expressed as mean ± standard deviation or median (interquartile range) based on data distribution and compared using Student’s t-test or Mann-Whitney U test. Categorical variables were presented as frequencies (percentages) and analyzed using Chi-square test or Fisher’s exact test. Univariate logistic regression analysis was performed to identify risk factors associated with NEC development, with results presented as odds ratios (OR) with 95% confidence intervals (CI). Variables with P < 0.05 in univariate analysis were considered statistically significant. All statistical analyses were conducted using R software (version 4.4.3).

This study included 74 infants in the NEC group and 74 infants in the non-NEC group for comparative analysis (Table 1). The two groups showed no significant differences in baseline characteristics or perinatal factors (all P > 0.05). In terms of comorbidities, the incidence of feeding intolerance (FI) was higher in the NEC group (82.43% vs. 50.00%, P < 0.001). Laboratory indicators revealed that the platelet count (Plt) in the NEC group was significantly lower than that in the non-NEC group (184.63 vs. 243.00, P = 0.002), whereas the white blood cell count (WBC) (13.41 vs. 10.52, P = 0.040) and procalcitonin (PCT) levels (0.59 vs. 0.29, P < 0.001) were significantly higher (Table 1).

Among the 70 infants diagnosed with sepsis, 45 (64.3%) developed NEC while 25 (35.7%) did not. Notably, early-onset sepsis was significantly more prevalent in the NEC group compared to the non-NEC group (57.78% vs. 4.00%, P < 0.001), while late-onset sepsis showed the opposite pattern (42.22% vs. 96.00%) (Table 2). Additionally, positive blood cultures were more frequently observed in sepsis patients who developed NEC (46.67% vs. 4.00%, P < 0.001).

In summary, compared to non-NEC infants, those with NEC exhibited more severe systemic inflammatory responses (elevated WBC, PCT, and C-reactive protein (CRP), and decreased Plt), a higher incidence of FI, anemia, and sepsis, increased transfusion requirements, and longer hospital stays. We performed univariate logistic regression analysis for NEC based on statistically significant variables, and the forest plot results indicated that FI, neonatal anemia, blood product use, and NS were risk factors, while Plt levels approached being a protective factor (Figure 1A).

The study flowchart is illustrated in Figure 1B. A total of 204 and 1420 DEGs were identified from the GSE25504 and GSE46619 datasets, respectively (Figures 2A,B). In the intersection of these datasets, 70 overlapping DEGs with concordant directions of change were identified (69 upregulated and 1 downregulated) (Figures 2C,D). The complete list of DEGs is available as supplementary data, specifically in Supplementary Tables S1-S3. To investigate potential biological functions, GO enrichment analysis and KEGG pathway analysis were conducted on the overlapping DEGs using R software (Figures 2E,F). The results of the GO analysis revealed significant enrichment in the biological process (BP) terms, including cytokine-mediated immune response, leukocyte migration, and positive regulation of immune effector processes. In the cellular component (CC) terms, significant enrichment was observed in extracellular exosomes, specific granule membranes, and granule lumens. Furthermore, in the molecular function (MF) terms, significant enrichment was noted in cytokine receptor binding, immune receptor activity, and integrin binding (Supplementary Tables S4, S5). Additionally, pathways related to pertussis, the tumor necrosis factor (TNF) signaling pathway, cytokine-cytokine receptor interaction, rheumatoid arthritis, acute myeloid leukemia, inflammatory bowel disease, glucose metabolism, and tuberculosis were significantly enriched in the KEGG analysis. These findings suggest that pathways associated with immune response, cytokine signaling, and infection may play crucial roles in the pathogenesis of NS and NEC.

The PPI network contained 43 nodes and 263 interacting pairs (Supplementary Tables S6). An important subnetwork within the hub gene network was screened using the MCODE module in Cytoscape, which included 14 nodes and 80 interacting pairs (Supplementary Tables S7). These 14 nodes were IL1RN, FPR1, FPR2, HCK, S100A12, CSF3R, IL1B, FCGR1A, FCER1G, CD163, ITGAM, SPI1, MMP9, and C5AR1, with specific scores detailed in Supplementary Tables S8.

Hub genes were identified using the degree algorithm in the CytoHubba plugin, as detailed in Supplementary Tables S9. Figure 3A presents the decomposition diagram of the PPI network alongside the network constructed from the hub genes. A total of 10 hub genes were identified, including MMP9, FPR1, FCER1G, CD163, S100A12, ITGAM, SPI1, IL1RN, IL1B, and CSF3R. The interaction network of these genes was constructed using the GeneMANIA database, revealing significant enrichment in functions related to innate immunity and inflammation, such as tertiary granule formation, humoral immune response, secretory granule membrane, reactive oxygen species metabolic process, positive regulation of defense response, ficolin-1-rich granule, and NAD(P)H oxidoreductase activity (Figure 3B). Subsequent enrichment analysis indicated significant enrichment in immune cell chemotaxis and activation, as well as receptor signaling, particularly in biological processes (e.g., neutrophil/granulocyte chemotaxis, migration, activation, and heterotypic cell adhesion), cellular components (e.g., ficolin-1-rich granules and their membranes, tertiary granules and their membranes, endocytic vesicles and their membranes, lysosomal membranes, integrin complexes), molecular functions (e.g., immune receptor activity, receptor for advanced glycation end products (RAGE) binding, interleukin-1 receptor binding, cytokine binding, integrin binding, complement component C3b binding, IgG binding, complement receptor activity), and KEGG pathways (e.g., cytokine-cytokine receptor interaction, leukocyte transendothelial migration, IL-17 signaling pathway, Staphylococcus aureus infection, tuberculosis, pertussis, leishmaniasis, type 1 diabetes) (Figures 3C,D; Supplementary Tables S10, S11).

To identify reliable core genes, we conducted external validation on the selected genes. In the NEC dataset (GSE297483), only the expression levels of FPR1, S100A12, and CSF3R were significantly elevated compared to the control group (Figure 4A), while no statistical differences were observed for the other genes. In NS samples, nine genes had significantly higher expression levels than the control group, while IL1B showed no significant difference (GSE69686, Figure 4B). Consequently, we identified three core genes: FPR1, S100A12, and CSF3R. Subsequently, we generated ROC curves using external validation datasets to further assess the diagnostic value of these validated core genes. In the NEC-related validation dataset, the AUC values for FPR1, S100A12, and CSF3R were 0.762, 0.762, and 0.746, respectively (Figures 4C–E). Furthermore, in the NS-related validation dataset, the AUC values for all validated core genes exceeded 0.7, measuring 0.723, 0.813, and 0.807, respectively (Figures 4G–I).

The TRRUST database provides insights into the regulatory relationships between transcription factors and their target genes, elucidating their interactions. In this study, we utilized transcription factor binding site information from TRRUST to identify key transcription factors and target genes associated with NEC and NS. We identified a total of 35 associations involving 14 transcription factors (SPI1, ETS2, CEBPA, NFKBIA, JUN, SIRT1, FOS, CEBPB, RELA, NFKB1, ETS1, STAT1, STAT3, and SP1) and 6 core genes (IL1B, ITGAM, CSF3R, CD163, MMP9, and IL1RN), with specific regulatory relationships detailed in Supplementary Tables S12. Based on these findings, we constructed a transcription factor-messenger RNA regulatory network using Cytoscape software (Figure 5A) and validated it against the NEC and NS datasets. By analyzing gene expression data from the GSE46619 (NEC) and GSE25504 (NS) datasets, we confirmed the expression changes of these transcription factors in both conditions. Validation results indicated that several transcription factors exhibited significant expression changes in NEC and NS, further supporting their putative regulatory roles in these diseases (Figures 5B,C). We distinguished mRNA and TFs through different shapes and colors, effectively illustrating this regulatory network and suggesting that these TFs may play crucial roles in the pathological processes of NEC and NS.