Revealing hub pathway cross-talk for premature newborns with bronchopulmonary dysplasia by the integration of pathway analysis and Monte Carlo Cross-Validation

Introduction

Bronchopulmonary dysplasia (BPD), one of the most common complications of prematurity, is a lung injury in preterm infants resulting from oxygen and mechanical ventilation (1,2). Although progress has been made in neonatal care, the incidence of BPD is still increasing, mainly because the survival rate of younger and more premature babies has increased (3). Hence, effective prevention and treatments for BPD are on urgent need, and their necessary precondition is a comprehensive understanding of molecular and pathological mechanisms underlying this disease.

Previous literature has demonstrated that prenatal factors have important roles in the progression of BPD, such as surfactant deficiency and maternal infection (4,5). Recently, genetic variants were discovered in preterm infants of BPD (6,7). For instance, mutations of the lipid transporter ATP-binding cassette subfamily A member 3 (ABCA3), which is involved in surfactant synthesis in alveolar type II cells, have been revealed in BPD premature newborns (7). Moreover, these studies have highlighted a potential role for molecular pathways involved in inflammation and lung development in affected infants. Cross-talk among gene pathways can be means of regulatory interaction among different pathways or can express the gene overlap among pathways (8,9).

In this study, we proposed to reveal hub pathway cross-talk for premature newborns with BPD based on the integration of pathway related analyses and Monte Carlo Cross-Validation (MCCV) algorithm. Ultimately, hub cross-talk was extracted from seed cross-talk dependent on the MCCV method. The study of cross-talk involved in BPD could help in defining the molecular mechanisms leading to the onset and progression of the pathology (10).

Materials and methods

Data collection

Gene expression data with accession number E-GEOD-8586 (11) was collected from the ArrayExpress database for premature newborns at risk of BPD related studies. In E-GEOD-8586, there were 54 newborns with a gestational age from 23 to 27 weeks. These 54 newborns included 20 infants who later developed BPD (BPD group) and 34 infants who did not develop BPD (normal controls). There were 23 females and 31 males. Subsequently, in order to control the quality of the data, background corrections and standard normalizations were conducted.

Identification of differentially expressed genes (DEGs)

For investigating the differentially expressed levels of genes across BPD samples and normal controls, the empirical bayes (EB) approach was implemented to identify DEGs. In particular, the EB approach is equivalent to shrinkage of the estimated sample variances towards a pooled estimate, resulting in far more stable inference when the number of arrays is small (12). During this process, P-value for each gene was adjusted by Benjamini-Hochberg (BH) test (13). Importantly, only genes that met the thresholds of false discovery rate (FDR) <0.01 and |log2FoldChange| >2 simultaneously were considered as DEGs in BPD and normal control groups.

Determination of differentially expressed pathways (DEPs)

In the present study, all human biological pathways were derived from the Ingenuity Pathways Analysis (IPA) tool. Subsequently, DEGs of BPD were mapped to these biological pathways with an attempt to make them more correlated with the disease, and we gained pathways that were responsible for coordinating DEG activities, named as background pathways. Then we evaluated the significance between DEGs and genes in background pathways using Fisher's exact test (14). P-values were adjusted by BH test (13). In addition, the cut-off for DEPs across BPD patients and normal controls was P<0.01.

DEP cross-talk investigation

Based on the DEPs, we further focused on the relationships among any two of them. If a correlation existed between two DEPs, the DEP pair would be defined as a cross-talk for DEPs (10). Hence, we selected the discriminating score (DS) to access the strength of cross-talk among DEPs in the present study (10). Supposing that × and y were two DEPs in cross-talk, the DS was computed as following formula:

Seed cross-talk extraction

With an attempt to evaluate the functional activities of DEP cross-talk more deeply, a random forest (RF) classification model (15) was utilized to access the classified performance of these pairs between BPD and normal control groups. The method was comprised of three sections, first, drawing Ntree bootstrap samples dependent on the DS for each DEP cross-talk (Ntree = 500); secondly, growing a regression tree for each of the bootstrap samples; and finally, aggregating the predictions of the Ntree trees. The top 10 DEP cross-talk in descending order of AUC were defined as seed cross-talk.

MCCV for seed cross-talks

As mentioned above, seed cross-talk was identified, and then the MCCV method was adopted to validate the biological functions of seed cross-talk in BPD samples (16). The 54 samples were randomly split into two sets, the calibration set (Sc), contained nc samples for fitting the models; and the validation set (Sv), included nv samples for validating the model.

Of which E represents the squared prediction error, R represents the procedure repeated times (R=50). Theoretically, the fewer samples used in model calibration, the more repeat times were needed. For each bootstrap, the DEGs, DEPs, cross-talk and their DS values were carried out. The repeat times of DEP cross-talk were statistically calculated. More repeat times might imply more significance of this cross-talk.

Results

DEGs

In the present study, there were 20,514 genes in the pretreated E-GEOD-8586. After performing EB algorithm, a total of 132 DEGs which satisfied the cut-off (P<0.01 and |log2FoldChange| >2) between BPD premature newborns and normal controls were obtained. All DEGs were ranked in descending order of their P-values (Table I). The smaller the P-value, the more significant the expression of the DEG was. The most significant five DEGs across BPD group and normal controls were EMILIN1 (P=4.29Ex10−6), MFAP51 (P=8.77Ex10−6), SNCA1 (P=9.23Ex10−6), interleukin-62 (IL-62) (P=3.14Ex10−5) and STRADB1 (P=7.56Ex10−5).

Table I. Differentially expressed genes (DEGs).

Table I.

Differentially expressed genes (DEGs).

Rank	Genes	Rank	Genes	Rank	Genes	Rank	Genes
1	EMILIN1	34	LINC00702	67	TRIM58	100	CXCL3
2	MFAP51	35	SLC25A391	68	BPGM	101	MORC4
3	SNCA1	36	STRADB2	69	SLC19A2	102	C3
4	IL-62	37	XK1	70	ZBTB16	103	SERPINE11
5	STRADB1	38	EPB422	71	SNCA	104	NPTN-IT1
6	CXCL12	39	AHSP1	72	DMTN	105	PTGS2
7	BPGM1	40	KRT13	73	EPB42	106	CXCL6
8	CXCL32	41	S100A91	74	PTX31	107	SGK1
9	TNFAIP32	42	SOD22	75	IL-61	108	CXCL8
10	SOD21	43	BOK	76	HEMGN	109	S100A8
11	EPB421	44	GYS1	77	KRT1	110	TNFAIP3
12	PTX32	45	CCL22	78	BLVRB	111	LIF
13	CCL21	46	S100A81	79	STRADB	112	CXCL1
14	HBM1	47	KRT141	80	FAM46C	113	FOSB1
15	CXCL62	48	MFAP52	81	UCP2	114	CCL20
16	KRT14	49	ZBTB163	82	SERPINE12	115	IL-6
17	PLEKHO1	50	CCL202	83	XK	116	SOD2
18	UCP21	51	CXCL33	84	AHSP	117	TNFAIP6
19	KAT2B1	52	HBM4	85	GCH1	118	PTX3
20	POC1B	53	CXCL63	86	CHRDL2	119	S100A9
21	ADIPOR1	54	PTX33	87	LIF1	120	CCL2
22	ID1	55	CXCL82	88	CXCL31	121	S100A12
23	HBQ1	56	CXCL13	89	TNFAIP61	122	MFAP5
24	FAM46C1	57	SQRDL	90	CCL201	123	EGR1
25	DCAF12	58	MAP3K8	91	SLC25A39	124	FOS
26	BPGM2	59	DCK	92	HBM	125	FOSB
27	C15orf48	60	SUSD6	93	CSF3	126	SERPINE1
28	PLA2G2A1	61	KAT2B	94	TNFAIP31	127	FOSB2
29	SERPINB2	62	GYPA	95	CXCL11	128	EGR11
30	ZBTB162	63	SGK11	96	CXCL61	129	FOS1
31	HBM3	64	KRT5	97	CXCL81	130	ZBTB161
32	HEMGN1	65	YOD1	98	PLA2G2A	131	HBM2
33	SNCA2	66	PIP4K2A	99	SELE	132	RPS4Y1

DEPs

As described above, the IPA tool with 589 biological pathways (covering 5,169 genes) was implemented to execute pathway enrichment analysis for 132 DEGs of BPD infants, and the Fisher's exact test was used to evaluate significance of the pathways. Consequently, a total of 137 DEPs with P<0.01 were identified between BPD samples and normal controls, of which role of IL-17A in psoriasis (P=5.96Ex10−13) was the most significant one. Furthermore, DEPs with larger number of DEGs might be more important than that of the smaller ones. The DEPs with DEG number ≥4 are listed in Table II. Among these DEPs, flucocorticoid receptor signaling (P=7.11Ex10−7) contains 6 DEGs. Interestingly, the number of DEGs enriched in 6 DEPs was 5, including granulocyte adhesion and diapedesis (P=1.55Ex10−6), role of IL-17A in psoriasis (P=5.96Ex10−13), glucocorticoid receptor signaling (P=5.32Ex10−6), role of IL-17A in arthritis (P=2.24E-10−8), IL-17A signaling in airway cells (P=4.54Ex10−8) and flucocorticoid receptor signaling (P=4.70Ex10×10−5).

Table II. Differentially expressed pathways (DEPs) with DEG number ≥4.

Table II.

Differentially expressed pathways (DEPs) with DEG number ≥4.

DEPs	P-value	Number of genes in pathway	Number of DEGs
Glucocorticoid receptor signaling	7.11×10−7	255	6
Granulocyte adhesion and diapedesis	1.55×10−6	163	5
Role of IL-17A in psoriasis	5.96×10−13	13	5
Glucocorticoid receptor signaling	5.32×10−6	255	5
Role of IL-17A in arthritis	2.24×10−8	54	5
IL-17A signaling in airway cells	4.54×10−8	63	5
Glucocorticoid receptor signaling	4.70×10−5	255	5
Role of IL-17A in arthritis	6.74×10−7	54	4
IL-17A signaling in airway cells	1.18×10−6	63	4
Agranulocyte adhesion and diapedesis	5.89×10−5	173	4
Role of IL-17A in arthritis	2.72×10−7	54	4
IL-17A signaling in airway cells	4.78×10−7	63	4
Glucocorticoid receptor signaling	2.25×10−5	255	4
Atherosclerosis signaling	6.60×10−8	119	4
Role of IL-17A in psoriasis	2.63×10−9	13	4
Granulocyte adhesion and diapedesis	1.07×10−4	163	4
Agranulocyte adhesion and diapedesis	1.36×10−4	173	4

Seed cross-talk

In this section, first of all, DS was assigned to any pair of DEPs or cross-talk by comparing their gene expression levels. In consequence, 609 cross-talk among 137 DEPs were identified, and the DS distribution is illustrated in Fig. 1. Two of the cross-talk had higher DS than the others: Role of IL-17A in psoriasis; IL-17A signaling in fibroblasts (DS=6.328) and LXR/RXR activation; role of IL-17A in psoriasis (DS=6.134). Subsequently, we applied the RF classification of cross-talk for each sample based on the DS, and chose the AUC to evaluate the classification performance of each cross-talk. The cross-talk was ranked according to AUC in a descending order and the top 10 cross-talk were denoted as seed cross-talks (Table III). The best cross-talk was the pair of pathways acute phase response signaling and gepatic fibrosis/hepatic stellate cell activation (AUC=0.915). Besides, the AUC for the following two cross-talks were ranged from 0.8 to 0.9: The pair of pathways role of IL-17F in allergic inflammatory airway diseases and role of IL-17A in psoriasis (AUC=0.857), and the pair of pathways LXR/RXR activation and hepatic fibrosis/hepatic stellate cell activation (AUC=0.812).

Figure 1. DS was assigned to any pair of DEPs or cross-talk by comparing their gene expression levels. In consequence, 609 cross-talk among 137 DEPs were identified, and the DS distribution are illustrated. DEP, differentially expressed pathway; DS, computing discriminating score.

Table III. Seed cross-talk with AUC value for random forest (RF) classification.

Table III.

Seed cross-talk with AUC value for random forest (RF) classification.

ID	Cross-talk	AUC
1	Acute phase response signaling; hepatic fibrosis/hepatic stellate cell activation	0.915
2	Role of IL-17F in allergic inflammatory airway diseases; role of IL-17A in psoriasis	0.857
3	LXR/RXR activation; hepatic fibrosis/hepatic stellate cell activation	0.812
4	LXR/RXR activation; acute phase response signaling	0.799
5	Role of IL-17A in psoriasis; IL-17A signaling in fibroblasts	0.790
6	IL-17A signaling in airway cells; role of hypercytokinemia/hyperchemokinemia in the pathogenesis of influenza	0.786
7	IL-17 signaling; differential regulation of cytokine production in intestinal epithelial cells by IL-17A and IL-17F	0.786
8	Role of IL-17F in allergic inflammatory airway diseases; atherosclerosis signaling	0.782
9	LXR/RXR activation; glucocorticoid receptor signaling	0.782
10	IL-17A signaling in gastric cells; role of hypercytokinemia/hyperchemokinemia in the pathogenesis of influenza	0.778

Hub cross-talk

Herein, we divided 20 BPD samples and 30 normal controls of total 54 samples in the gene expression data into two sets according to the ratio of 6:4. Among the 10 seed cross-talk, we uncovered that 3 cross-talk were repeated many times. The pair of pathways role of IL-17F in allergic inflammatory airway diseases and role of IL-17A in psoriasis was repeated 29 times. The pair of pathways role of IL-17A in psoriasis and IL-17A signaling in fibroblasts was repeated 28 times. The pair of pathways IL-17A signaling in airway cells and role of hypercytokinemia/hyperchemokinemia in the pathogenesis of influenza was repeated 28 times. Hence we defined these three cross-talk as hub cross-talk for BPD patients, which might play more crucial roles in the progression of BPD samples than the others.

Discussion

In the present study, we investigated hub cross-talk of DEPs for BPD samples by integrating a series of bioinformatics analyses, identification of DEGs, extraction of DEPs, exploration of seed cross-talk and investigation of hub cross-talk. Although the scenario might be intuitive, the biological functions and activities of DEPs have not been well studies (17,18). In order to increase the confidence and feasibility of our study, we implemented the MCCV to screen the hub cross-talk (10,19).

Consequently, a total of 3 hub cross-talk for premature newborns with BPD were gained. IL-17F is a cytokine that shares sequence similarity with IL-17 and is expressed by activated T cells (20). It inhibits angiogenesis of endothelial cells and induces endothelial cells to produce IL-2, TGFB1/TGFB, and signaling and function in inflammation (21). The involvement of IL-17A/IL-17F in inflammatory immune responses and host defense mechanisms, and the correlation between IL-17A/IL-17F with other cytokines has been reviewed previously (22,23). The review indicated that IL-17A and IL-17F possibly work together in human disease. Coincidently, in our results, the hub cross-talk (role of IL-17F in allergic inflammatory airway diseases and role of IL-17A in psoriasis) also revealed this phenomenon in BPD development.

IL-17A has been demonstrated to take vital part in allergic airway inflammation, and promote inflammation via inducing various proinflammatory cytokines and chemokines, recruiting neutrophils, increasing antibody production, and activating T cells (24). We found that IL-17A was correlated to all of the three hub cross-talk of BPD infants, including the pair of pathways role of IL-17F in allergic inflammatory airway diseases and role of IL-17A in psoriasis; the pair of pathways role of IL-17A in psoriasis and IL-17A signaling in fibroblasts; the pair of pathways IL-17A signaling in airway cells and role of hypercytokinemia/hyperchemokinemia in the pathogenesis of influenza. Importantly, IL-17A has been reported to exhibit a vital role in the development of pulmonary fibrosis, suggesting it is a potential treatment target for inflammation-induced fibrosis (25). Taken together, we inferred that the three hub cross-talk closely correlated to BPD.

In summary, 3 hub cross-talk among DEPs were identified in our study, which might give insight into revealing the molecular mechanism underlying the premature newborns with BPD.

Acknowledgements

Not applicable.

Funding

No funding was received.

Availability of data and materials

The datasets used in this study are available from the corresponding author on reasonable request.

Authors' contributions

CW, RP and BZ conceived the study. MC and GC acquired, and analyzed the data. CW wrote the manuscript. RP, GC and MX explained the results and revised the manuscript. All authors read and approved the final manuscript.

Ethics approval and consent to participate

Not applicable.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References