Key regulatory genes associated with 20 PCD patterns were sourced from multiple repositories, including the Kyoto Encyclopedia of Genes and Genomes (KEGG) Database (22), the GeneCards Database (23), the Molecular Signatures Database (24), the Reactome Database, FerrDb database, the Human Autophagy Database, and published articles. The comprehensive gene list comprising 20 distinct PCD patterns (Supplementary Table 2), encompassing 7 alkalosis-related genes (25), 338 anoikis-related genes (26, 27), 579 apoptosis-related genes, 19 Cuproptosis-related genes (28–30), 10 disulfidptosis-related genes (31), 23 entosis-related genes, 15 entotic cell death-related genes, 34 immune cell death (32), 220 lysosome dependent cell death-related genes, 8 methuosis-related genes, 39 MPT driven necrosis, 101 necroptosis-related genes, 24 NETosis-related genes, 8 netotic cell death-related genes, 5 oxeiptosis-related genes, 9 parthanatos-related genes, 59 parapoptosis-related genes (33), 52 pyroptosis-related genes (34), 367 autophagy-related genes and 88 ferroptosis-related genes.

The R package “limma” was employed to extract differentially expressed genes (DEGs) between the IPF cohort and health control cohort in the IPF dataset (p value < 0.05).

Further, the study performed the gene ontology (GO), and Kyoto encyclopedia of genes and genomes (KEGG) enrichment analysis of all the DEGs using the “ClusterProfiler” R package.

The GSE122960 (35) dataset provided scRNA-seq information from the lung tissue of three IPF patients and three healthy donors. To ensure the correctness and reliability of the information, the scRNA-seq information underwent stringent screening and quality control. It excluded cells with mitochondrial gene expression levels beyond 20% from the dataset. Furthermore, cells containing fewer than 300 genes and genes represented by fewer than five cells are excluded. The selected cells were subsequently subjected to downstream analysis. The scRNA-seq information was normalized by the “NormalizeData” function, which was then converted to Seurat objects and the first 1,500 highly variable genes were identified using the “FindVariableFeatures” function. Afterward, the “RunPCA” function of the “Seurat” R package was applied to perform principal component analysis (PCA) to reduce the dimensionality of the scRNA-seq data based on the top 1,500 genes. The functions “FindNeighbors” and “FindClusters” were used for cell clustering analysis.

Cell clusters were annotated using reference data from the Human Cell Atlas and were subsequently refined based on specific cell biomarkers including Type I Alveolar Epithelial (AT1) Cells (AGER and RTKN2), Type II Alveolar Epithelial (AT2) Cells (LAMP3), Club Cells (SCGB3A2), Ciliated Cells (TPPP3), Basal Cells (KRT5 and TP63), Macro/Mono (CD68 and CD14), Dendritic Cells (CLEC10A), T/NKT Cells (CD3D), Plasma Cells (IGHG4), B Cells (MS4A1), Mast Cells (TPSB2), Endothelial Cells (VWF), Fibroblasts (DCN).

Furthermore, the “Scissor” R package was employed to integrate scRNA-seq information, bulk RNA-seq information, and phenotypic information to identify cell subpopulations that display significant correlations with the prognosis of IPF patients. The “IPF-related genesEA” R package was employed to conduct gene set enrichment analysis (GSEA) on the scRNAseq dataset based on the “AUCell,” “UCell,” “singscore,” and “ssgsea” algorithms.

The datasets enrolled in the study were obtained from the Gene Expression Omnibus (GEO) database.¹ The prognosis signature was constructed based on four datasets, including GSE27957 (36), GSE28042 (36), GSE70866 (37), and GSE93606 (38). GSE70866 was used as the training set with various clinical features and samples, while the other three data sets were used as the validation sets. Each gene expression was transformed into a z-score across patients in all cohorts. The screening criteria for samples were as follows: (1) The individual received a diagnosis of IPF; (2) The IPF patient possessed comprehensive bulk-RNAseq information; (3) The IPF patient had a survival information record (refer to Supplementary Table 3 for detailed clinical parameters).

To enhance the precision and consistency of the prognostic signature, the study incorporated ten machine-learning algorithms into the analysis. These algorithms encompass random survival forest (RSF) (39), elastic network (Enet), Lasso, Ridge, Stepwise Cox, CoxBoost (40), partial least squares regression for Cox (plsRcox) (41), supervised principal components (SuperPC) (42), generalized boosted regression modeling (GBM) (43), and survival support vector machine (survival-SVM) (44) (refer to Supplementary Table 4 for details). Several algorithms have demonstrated the ability to perform feature selection, which was employed to screen crucial genes, including Lasso, stepwise Cox, Coxboost, and RSF algorithms. Ultimately, we developed 101 algorithmic combinations to enhance the identification of prognostic signatures with commendable accuracy and stability. Subsequent survival analysis and signature comparison were employed to assess the superiority and generalizability of the prognostic signature.

GSE150910 (45), GSE24206 (46), GSE28042 (36), GSE53845 (47), and GSE70866 (37) from the GEO database were enrolled to develop the diagnostic signature. The GSE150910 dataset served as the training cohort, while four datasets were utilized as the validation cohorts.

An efficient diagnostic signature for accurate prediction of IPF was developed by merging diverse machine-learning techniques based on PCD-related genes. The machine learning algorithms encompass Lasso, Ridge, Stepglm, XGBoost, Random Forest (RF) (39), Elastic Net (Enet), Partial Least Squares Regression for Generalized Linear Models (plsRglm), Generalized Boosted Regression Modeling (GBM), Naive Bayes (48), Linear Discriminant Analysis (LDA) (49), Generalized Linear Model Boosting (glmBoost), and Support Vector Machine (SVM). Ultimately, a total of 101 signatures were developed. For each signature, the area under the receiver operating characteristic curve (AUC) values were calculated across all validation datasets, and the signature exhibiting the highest average AUC in the validation cohort was deemed the superior diagnostic signature, owing to overfitting in the training cohort. The diagnostic model’s superiority was established by comparing its AUC value with that of the clinical features.

The “pRRophetic” (50) R package was applied to predict the therapeutic response of IPF patients to common drugs, and the value of the PCDI.prog signature in guiding the selection of drugs for IPF patients was assessed based on the IC50 values in different PCDI.prog score groupings.

Statistical differences between groups were determined by Student’s t-test for normally distributed variables, and for non-normally distributed variables, statistical differences between groups were determined by the Wilcoxon test. The statistical studies were conducted using the R project (version 4.3.3).

In this study, we collected a total of 20 PCD patterns and 2013 key regulatory genes from the existing published articles and online databases. We removed 449 duplicate gene symbols, resulting in 1,564 PCD-related genes for subsequent analysis (Figure 2A).

A total of 524 genes with significant differential expression (p-value < 0.05) were identified in the GSE70866 cohort (Figure 2B). Furthermore, the DEGs are associated with many PCD pathways and signal transduction pathways related to IPF, as demonstrated by the KEGG and GO enrichment studies (Figures 2C,D).

In addition, the study investigated the profiling of pathway activity in normal and fibrotic lung tissues to analyze PCD patterns (Figure 2E). The results suggested that except NETosis pattern, the activity of other PCD patterns decreased in fibrotic lung tissue. The vast majority of PCD patterns showed significant differences in activity in normal and fibrotic lung tissues.

A total of 9,295 cells, obtained from lung tissue samples of three IPF individuals and three healthy individuals, were grouped into 23 unique clusters using the “Seurat” and “clustree” R packages (Supplementary Figures 1A,B). A total of 13 unique cell types were distinguished using particular cell marker genes (Supplementary Figure 1C; Figure 3A). The population of epithelial cells in the lungs of individuals with IPF was shown to have a higher percentage of airway epithelial cells (specifically basal cells, ciliated cells, and club cells) and a notable drop in alveolar epithelial cells (AT1 cells and AT2 cells). This pattern aligns with prior findings (51, 52) (Figure 3B).

Furthermore, the study provided an initial description of the variations in PCD-related pathway activity among different cell subtypes in fibrotic lung tissues based on the integration of four algorithms (Figures 3C,D). As we all know, fibroblasts play a crucial role in the development of IPF by promoting the production of collagen and cell proliferation, leading to the formation of fibrous tissue (53). The study suggested that there is a difference in the behavior of PCD patterns in fibroblasts in fibrotic lung tissues, especially heightened activity in the anoikis pathway in fibroblasts. Previous studies have demonstrated that anoikis takes place in lung fibroblasts during the development of fibrosis and has been identified as a significant factor in the progression of IPF (54), indicating possible treatment targets for IPF.

A univariate regression analysis was conducted on a set of 524 DEGs to find 137 prognostic genes. Furthermore, the “Scissor” R package was employed to integrate the scRNA-seq GSE122960 dataset and the bulk-RNAseq GSE70866 dataset. We identified 1,934 cells positively associated with OS in patients with IPF and 2,041 cells that were negatively associated with OS in patients with IPF (Figure 3E). Additionally, the “Findmarkers” function was used to identify 1,392 marker genes for cells that had prognostic significance. Ultimately, a total of 44 IPF-related genes (Supplementary Table 5) were determined by identifying the common genes that serve as both prognostic genes and marker genes (Figure 3F).

Subsequently, we employed a total of 10 distinct machine-learning algorithms to develop and construct 101 prognostic signatures. The robustness of these signatures was evaluated using a 10-fold cross-validation approach in four distinct cohorts (GSE70866 as the training cohort, and three external validation cohorts including GSE27957, GSE28042, and GSE93606 cohort).

The best-performing predictive signature was determined as the signature with the greatest mean C-index in three external validation cohorts, due to overfitting in the training cohort (Figure 4A). The findings indicated that the CoxBoost + Enet[alpha = 0.1] algorithm combination demonstrated the highest average C-index (0.676), making it the optimal algorithm combination for developing the programmed cell death index prognostic (PCDI.prog) signature.

Afterward, the PCDI.prog score for each patient was calculated. The samples were then categorized into the low PCDI.prog score group and the high-PCDI.prog score group based on the optimal threshold value of the PCDI.prog score determined using the R package “survminer.” Subsequently, KM survival analysis and assessment of prognostic performance were conducted. Demonstrating a substantial difference in OS between the low-PCDI.prog score group and the high-PCDI.prog score group in all four cohorts (Figure 4B).

The accuracy and reliability of the PCDI.prog signature predicting 1-,2-, 3-, 4-, and 5-year survival of IPF patients was supported by empirical evidence that the area under the curve (AUC) values exceeded 0.65 in multiple distinct cohorts (Figure 4C).

Furthermore, due to the scarcity of prognostic models for non-tumor diseases compared to tumors, and the limited availability of datasets including comprehensive gene expression data for genes associated with signatures, a total of 10 prognostic models published in IPF were ultimately gathered from the existing literature. The features encompass a range of biological processes observed in the IPF cohort, such as hypoxia, autophagy, pyroptosis, epithelial-mesenchymal transition, epigenetic regulation, and inflammation (Supplementary Table 6). These features were then compared to the C-index of the PCDI.prog signature. The findings indicated that the PCDI.prog signature exhibited superior performance compared to the majority of the signatures within their respective categories (Figure 4D).

In addition, we compared the predictive value of the PCDI.prog signature with other clinical variables (Figure 5A). The C-index of the PCDI.prog signature was significantly higher than other clinical variables, covering GAP score, age, and gender. The univariate Cox regression analysis showed that compared with other features, PCDI.prog signature was regarded as a risk factor (Figure 5B). Finally, to facilitate clinical application, a nomogram was created, integrating the factors of age, gender, GAP score, and PCDI.prog score (Figure 5C).

The “pRRophetic” R package is a computational model that predicts chemotherapy responses based on gene expression data (55). Nintedanib, an orally administered small molecule tyrosine kinase inhibitor initially designed for lung cancer, has been approved for the treatment of IPF. It serves as an example of an anti-tumor medication associated with medications used for IPF. Furthermore, the addition of pirfenidone (56) can diminish the efficacy of the paclitaxel and carboplatin combination. Therefore, we believe that the identification of related anti-tumor drugs has potential significance for IPF treatment. Eight potential drugs for IPF patients were obtained through drug susceptibility analysis (Figure 5D). Individuals with low PCDI.prog scores exhibited a notable rise in sensitivity to staurosporine, linsitinib, sapitinib, taselisib, and alpelisib. Individuals with high PCDI.prog scores exhibited a notable rise in sensitivity to palbociclib, sabutoclax, and mitoxantrone. It is suggested that the PCDI.prog signature has a potential guiding effect on the treatment of IPF patients.

A total of 101 diagnostic signatures of GSE150910 were evaluated using 10-fold cross-validation (Figure 6A). The AUC values for each signature were computed across all validation datasets, including GSE110147, GSE24206, GSE28042, GSE53845, and GSE70866. The combination of Stepglm[both] and NaiveBayes was recognized as the optimal signature, yielding the greatest average AUC value of 0.856. The Stepglm[both] algorithm identified 9 crucial genes (SLC39A8, HIF1A, TIMP1, ACSL1, ALOX5, MET, IL1R1, HTRA1, TP53INP1) (Figure 6B), whereas the NaiveBayes algorithm developed the programmed cell death diagnostic (PCDI.diag) signature.

The diagnostic effectiveness of the PCDI.diag signature in predicting IPF was compared with other clinical factors in both the training and validation cohorts (Figure 6C). Regrettably, the GSE110147 and GSE24206 cohorts lack comprehensive clinical data. The findings indicated that the PCDI.diag signature exhibited markedly superior accuracy compared to other clinical features in GSE150910, GSE28042, GSE53845, and GSE70866, including age, rs35705950_genotype, and gender. The results suggest that the PCDI.diag signature, derived from nine crucial genes, may offer novel insights into the preliminary diagnosis of IPF. It is worth noting that the AUC values of gender in multiple datasets are weak, which means that gender cannot well distinguish the IPF patient group from the healthy control group. Previous studies have also shown that Although IPF seems to affect older men predominantly, the true prevalence of IPF in women is difficult to establish, and women may be underdiagnosed while men are overdiagnosed with IPF based on gender alone (57).