A total of 196 patients with IPF, including their clinical information, and 20 healthy donors were included in this study. The gene expression of BAL cells (GSE70866) was obtained from Gene Expression Omnibus (GEO) database (23). Among them, RNA microarray chips from 20 healthy donors and 112 patients with IPF (from Freiburg and Siena) were performed using Agilent-028004 SurePrint G3 Human GE 8x60K Microarray, while RNA microarray chips from 64 IPF patients (from Leuven) were evaluated by Agilent-039494 SurePrint G3 Human GE v2 8x60K Microarray. The basic information of patients was shown in Table 1 . In addition, the “limma” and “sva” functions were performed to background-correct, normalize and remove batch effects from raw expression data (24, 25). The cohort of 112 patients with IPF was employed as the training set, and the cohort of 64 patients with IPF was used as a validation set. Ethical approval was not required in this study since all data was available in the public domain.

The “limma” package was applied to screen differentially expressed genes (p value<0.05). The gene sets of hallmarks were obtained from the Molecular Signatures Database (MSigDB), and the Z-score of hallmarks was quantified using a single-sample gene set enrichment analysis (ssGSEA) algorithm (R package ‘gsva’) based on transcriptome profiling data. Univariate Cox proportional-hazards (Cox-PH) regression was applied to estimate the significance of different hallmarks in IPF through the “survival” function.

Weighted Gene Co-expression Network Analysis (WGCNA) was carried to identify the module that was most correlated with hypoxia based on transcriptome profiling data and ssGSEA scores. Firstly, the gene expression matrix was transformed into a similarity matrix by using the Pearson test between pairwise genes. Secondly, the similarity matrix was transformed to an adjacency matrix by Topological overlap measure (TOM). Gene co‐expression networks were generated by the WGCNA package. Module membership (MM) represented the correlation between module eigengenes and gene expression profiles. 369 candidate genes were identified with a threshold of correlation coefficient of MM >0.5 and the p-value <0.001.

The univariate Cox regression was performed to assess the prognostic value of candidate genes. Then, the importance of the survival-related differential genes was calculated and ranked using the random forest via the “randomForestSRC” package. A risk model was built by multivariate Cox regression, and the risk score was calculated as follows: Risk Score = β1x1 +β2x2+β3x3+…+ βnxn.

Single-Sample GSEA (ssGSEA) and MCP counter algorithm were applied to estimate the infiltration levels of immune cells based on specific gene expression signatures of immune cells. Subsequently, we explored the correlation between risk score and immune cells. Furthermore, we also analyzed the relationship between risk score and NLRP3.

R software (version 3.6.3, http://www.r-project.org) was used to analyze data and plot graphs. Kaplan-Meier (K-M) curve plotted the relationship between score and clinical outcome, and the log-rank test was used to evaluate differences using “survival” and “survminer” packages. Univariate and multivariate Cox regressions were performed with the “survival” package. Subsequently, the Wilcox test was used to calculate the differences between the groups, and p<0.05 was considered significant. WGCNA was used to describe the correlation patterns among gene by following the protocol of WGCNA package. Time-dependent concordance index (C-index) and the receiver operating characteristic (ROC) analysis were applied to assess the predictive capacity of survival among different variables using R packages’ survConcordance’ and ‘survivalROC’. Then, gene set enrichment analysis (GSEA) was performed using GSEA software (http://software.broadinstitute.org/gsea). Finally, the calibration of the nomogram was plotted by the “rms” package.

The z-scores of various hallmarks were calculated based on the transcriptome profiling of the training set and gene sets of MSigDB using ssGSEA. In addition, the hazard ratios (HR) of various hallmarks were calculated and ranked based on the survival information. The results indicated that the inflammatory response, angiogenesis, epithelial-mesenchymal transition, hypoxia, immune response, and other signaling pathways impacted the overall survival of patients with IPF from the training cohort ( Figure 1A ). Among these hallmarks, inflammatory response and hypoxia exhibited more powerful effects on survival compared to others. Then, a co-expression network of the survival-related hallmarks was generated based on their ssGSEA z-scores. The result revealed that compared to other hallmarks, inflammatory response and hypoxia showed a higher correlation ( Figure 1B ). Additionally, multivariate Cox analysis demonstrated that both inflammatory response and hypoxia were prognostic factors of IPF independent from other clinical features (including age and sex) ( Figure 1C ). Based on the median value of ssGSEA score of inflammatory response and hypoxia, the training cohort was divided into the high- and low-score groups, with inflammatory response and hypoxia presenting worse clinical outcomes as indicated by their high z-score ( Figures 1D, E ).

In total, 4072 genes were identified as differentially expressed genes between IPF and healthy donors. Then, WGCNA was used to construct a scale-free co-expression network based on the average linkage method, and 20 modules were generated ( Figure 2A ). Among them, both green-yellow and light green modules showed a strong positive correlation with inflammatory response and hypoxia (r>0.6, p<0.0001) ( Figure 2B ). Subsequently, univariate Cox regression was performed to examine the prognostic value of the 369 candidate genes from the green-yellow and light green modules, which revealed 42 and 174 genes to be negatively and positively associated with prognosis, respectively ( Figure 2C ). Furthermore, the random forest analysis was applied to rank and screen relative important genes based on their impact on overall survival. Consequently, TFPI, CTSE, TLR2, CDC42EP2, CCL2, CECR6, SPP1, WFDC2, HS3ST1, and HAMP were identified as the top 10 genes ( Figure 2D ). Approximately 1023 risk models were calculated for the different combinations of the 10 genes and K-M analysis was performed to screen for the best risk model. The Log-rank p-value of all risk models was calculated and ranked, and the top 20 was selected for further study ( Figure 2E ). The five-gene signature was identified as the final risk model, as it exhibited the smallest p-value and fewer risk genes, with the risk score = 0.44× HS3ST1+ 0.27× WFDC2 + 0.16 × SPP1 + 1.12 × TFPI + 1.06 × CDC42EP2.

In the training set, the 112 patients were divided into high and low-risk groups based on the median of risk scores. The Bar plots demonstrated that the proportion of deaths in the high-risk group was 87%, significantly higher than the low-risk group ( Figure 3A ). Meanwhile, the risk score was significantly higher in patients who died during follow-up compared to survivors ( Figure 3B ). Furthermore, the K-M analysis indicated that the high-risk group had worse clinical outcomes than the low-risk group (HR=4.56, P<0.001, Figure 3C ). Area under the ROC curve (AUC) analysis revealed the reliable predictive ability of the model for the overall survival of IPF, with the 0.5-, 1-, 2-, 3-, and 5-years AUC values of 0.83, 0.82, 0.80, 0.81, and 0.94, respectively ( Figure 3D ). Furthermore, univariate and multivariate Cox analyses indicated that the five-gene signature was an independent risk factor (HR=2.78, P<0.0001, Figure 3E ). In addition, the C-index of the five-gene signature was significantly higher than that of age and sex, and was approximately 0.8 for 1-5 years ( Figure 3F ).

Similar to the training cohort, the high-risk group in the validation cohort had a higher mortality rate than the low-risk group, with a significant elevation of the risk score in patients who died ( Figures 4A, B ). The K-M plotter showed that patients with IPF in the high-risk group exhibited poorer overall survival than the low-risk group (HR=3.11, P=0.004, Figure 4C ). The ROC analyses indicated that the five-gene signature had a good prediction accuracy, with the AUC for 1-, 2-, 3-, and 5-years being 0.78, 0.79, 0.80, and 0.91, respectively ( Figure 4D ). The multivariate Cox regression analysis demonstrated that the five-gene signature served as an independent prognostic factor (HR=2.38, P<0.0001, Figure 4E ), which was further verified by the C-index that showed that the five-gene signature to have excellent prediction value in the validation cohort ( Figure 4F ).

Previous results have suggested that the ssGSEA score of inflammatory response and hypoxia was associated with the clinical outcome of patients with IPF in the training set. Following this, the analysis of the relationship between inflammatory response and hypoxia with the five-gene signature, demonstrated that the high-risk group had a higher ssGSEA score of inflammatory response and hypoxia than that of the low-risk group ( Figure 5A ). Meanwhile, the expression of inflammatory response and hypoxia-related genes, including HIF1A, CCR2, CCR5, TLR6, TLR2, IL1B, and IL6, was higher in the high-risk group than low-risk group ( Figure 5B ). In order to further explore the potential relationship of inflammatory response, hypoxia and related genes with the risk score, the correlation of risk score associated inflammatory response and hypoxia and related genes with clinical performance was analyzed, which showed that the patients with IPF having a high-risk score and high ssGSEA scores of hypoxias or inflammatory response presented worse clinical outcomes ( Figures 5C, D ). In addition, the patients with IPF showing a high expression of related-genes, including HIF1A, CCR2, CCR5, TLR2, IL1B, and IL6, and high risk-scores had a poorer prognosis ( Figures 5E–K ). Subsequent investigation of the five-gene signature on clinical outcome in different subgroups revealed that patients with IPF in the high-risk group had a poor prognosis in both male and female patients ( Figures 6A, B ). In both subgroups of patients ages ≥ 70 years and those aged < 70 years, the risk model retained its prognostic capacity to discriminate high-risk subset ( Figures 6C, D ). The GSEA analyses, which investigated the significantly altered KEGG pathway between the high- and low-risk groups, demonstrated that the chemokine signalling pathway, cytokine receptor interaction, JAK-STAT signalling pathway, leukocyte transendothelial migration, NOD-like receptor signalling pathway and VEGF signalling pathway were enriched in the high-risk group ( Figures 7A–F ). This further confirmed that the high-risk group was association with inflammation response and hypoxia status.

Although previous results indicated that the hallmarks of hypoxia and inflammatory response and five-gene signature heavily influences the prognosis of IPF patients, its potential function warranted additional explored. As reported, the BAL contains a large number of white blood cells, and they play crucial roles in the development and deterioration of organ fibrosis (26, 27). It is unclear whether there are differences in immune cell infiltration and inflammatory factors between high- and low- risk groups. Here, we analyzed the correlation between the infiltration level of immune cell and risk score, based on ssGSEA and MCP counter algorithm. Then, the results from ssGSEA analysis demonstrated that the high-risk group exhibited higher immune cell infiltration, including CD4+ T cell, natural killer cell, eosinophil, MDSC, Macrophage, neutrophil, dendritic cell, regulatory T cell, and follicular helper cell, compared to low-risk group ( Figure 8A ). Meanwhile, the results from MCP counter indicated that monocytic lineage and neutrophils abundance in the high-risk group was significantly higher that the low-risk group ( Figure 8B ). Furthermore, we also found that high-risk group had higher NLRP3 expression, which is known NLRP3 inflammasome drive chronic inflammation, compared to low-risk group ( Figure 8C ).

In order to quantify the risk assessment of individual patients with IPF and enhance the clinical applicability, a nomogram with five genes was constructed to predict the probability of 1-, 3-, and 5-years overall survival ( Figure 9A ). Calibration plots assessing the performance of the nomogram demonstrated its excellent prediction accuracy for 1-, 3- and 5-year survival probability (45° line, Figure 9B ).