Introduction

Idiopathic pulmonary fibrosis (IPF) is a fatal lung disease of unknown etiology characterized by chronic, progressive fibrosis.1,2 IPF patients often present with non-specific symptoms such as dyspnea and dry cough, leading to delayed or missed diagnosis.3 The prognosis of patients with IPF is generally poor, with significant variability in survival rates among individuals.4 Many patients eventually succumb to progressive chronic hypoxic respiratory failure,5 and the median survival following diagnosis ranges from 3 to 5 years,6 which is lower than many cancer survival rates.7 Despite evidence for immune mechanisms in lung fibrosis, immunotherapies have been unsuccessful for major types of IPF.8 Currently, IPF treatment options are limited to Food and Drug Administration (FDA)-approved drugs with first-line drugs approved by the, such as Pirfenidone and Nintedanib.9 These drugs can delay the deterioration of lung function, and patients exhibit varied responses and side effects, including gastrointestinal irritation.10 Consequently, specific diagnostic and prognostic markers for IPF are urgently needed to enable early diagnosis and identify patients with poor prognosis, thereby enhancing early screening and treatment strategies for IPF.

The pathological manifestations of IPF include excessive extracellular matrix deposition, fibroblast foci, and inflammatory cell infiltration.11 Repetitive epithelial injury leads to excessive wound repair and remodeling, ultimately leading to fibrosis.12 The innate and adaptive immune systems play crucial roles in fibrosis development.13 Additionally, programmed cell death is closely associated with the immune system.14 The innate immune system possesses the ability to activate programmed cell death mechanisms, enabling it to swiftly respond to cellular stressors and foreign microbial pathogens.15

Previously, it was believed that different modes of cell death had their fixed and unique pathways, and were mutually independent. However, Kanneganti’s team discovered that the internal proteins NP and PB1 of influenza A virus (IAV) can bind to Z-DNA binding protein 1 (ZBP1) to facilitate the activation of NLRP3 inflammasome through the receptor-interacting protein kinase (RIPK)-CASP8 pathway, triggering apoptosis, necroptosis, and pyroptosis in mouse bone marrow-derived macrophages.16 Subsequently, Malireddi et al officially named this complex cell death crosstalk pattern as PANoptosis.17 PANoptosis is a unique inflammatory programmed cell death pathway that plays a crucial role in the immune inflammatory response.18 This process is governed by the PANoptosome complex regulation, which incorporates key characteristics of pyroptosis, apoptosis, and necroptosis. However, it cannot be expressed exclusively through any individual pathway of cell death, such as pyroptosis, apoptosis, or necroptosis.19 PANoptosis has been implicated in various respiratory diseases, including asthma, acute lung injury, and silicosis.20,21 However, research on the role of PANoptosis in the development of IPF is currently limited.

This study utilized machine learning and single-cell analysis to explore the role of PANoptosis in IPF. We developed and experimentally validated a diagnostic and prognostic models based on PANoptosis related genes. The diagnostic model included AKT1, PDCD4, PSMA2, and PPP3CB. Conversely, the prognostic model included TNFRSF12A, DAPK2, UACA, and DSP. Additionally, potential therapeutic drugs, including Metergoline, Candesartan, and Selumetinib, were identified based on four prognostic genes associated with PANoptosis and validated through molecular docking. A flow chart of the analysis is shown in Figure 1.

Figure 1 Study flowchart.

Materials and Methods Data Source and Preprocessing

PANoptosis genes included genes involved in apoptosis, pyroptosis, and necroptosis. We collected five pathways downloaded from the MSigDB (V7.4) (GSEA | MSigDB:gsea-msigdb.org) database for analysis, including HALLMARK_APOPTOSIS, KEGG_APOPTOSIS, REACTOME_APOPTOSIS, REACTOME_PYROPTOSIS, and KEGG_NECROPTOSIS. The union of all gene sets for the five pathways was identified as the PANoptosis gene set.

A total of five independent public datasets were downloaded from the GEO database (http://www.ncbi.nlm.nih.gov/geo/). Specifically, GSE110147 includes 22 lung tissue samples from patients with IPF and 11 samples from normal lung tissues, which are utilized to identify differential genes and develop diagnostic prediction models. GSE53845 comprises 40 lung tissue samples from IPF patients and 8 normal tissue samples, which are employed to validate the diagnostic prediction model. GSE70866 comprises bronchial alveolar lavage (BAL) cell sequencing data from 176 patients with IPF and is utilized to develop a prognostic prediction model. GSE93606 includes peripheral whole blood sequencing data from 60 IPF patients, which is employed to assess the effectiveness of the prognostic prediction model. GSE122960 is a single-cell RNA sequencing dataset that contains lung tissue samples from 8 transplant donors and 4 IPF patients, aimed at investigating the role of key PANoptosis genes in pulmonary fibrosis.

Identification of Differentially Expressed Genes Related to PANoptosis

We utilized R to annotate, normalize, and perform log2 transformations on various datasets. Specifically, the “FactoMineR” and “factoextra” R packages were employed for principal component analysis (PCA). The “limma” package was used to conduct differential gene analysis (DEGs) on the processed dataset, with thresholds set at |log2 Fold change| > 1 and adjusted p-value < 0.05. Visualization of DEGs was achieved through heatmaps and volcano plots, generated using the “ggplot2” and “pheatmap” packages. Additionally, we employed “ggvenn” to filter PANoptosis differentially expressed genes (PANDEGs). For datasets measured across different platforms, batch effects were corrected using the “sva” R package.

Functional and Pathway Enrichment Analysis

We utilized the “clusterProfiler” R package to conduct Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses of PANDEGs. Additionally, we employed the “enrichplot” R package to visualize the top 10 enriched terms.

Immune Infiltration Analysis and Correlation Between PANDEGs and Immune Cells

ImmuCellAI is a tool designed for the accurate estimation of the abundance of 24 immune cell types from gene expression data, utilizing a gene set signature immune cell abundance assessment method.22 The “ImmuCellAI” R package was employed to evaluate the abundance of immune cell infiltration between the IPF and normal groups in the GSE110147 dataset. The “corrplot” and “ggplotify” R packages were utilized to determine and visualize correlations between immune cells. The differences in immune cell abundance are illustrated using a box plot, while the immune cell content is represented by a histogram, and the correlations between immune cells are depicted in a correlation matrix diagram.

Construction and Evaluation of Diagnostic Prediction Model

Least absolute shrinkage and selection operator (LASSO) is a dimensionality reduction technique for generalized linear regression that facilitates variable selection and parameter estimation through the incorporation of an L1 penalty term.23 GSE110147 was utilized as the training set, and LASSO regression was employed to identify the most relevant features in PANDEGs. GSE53845 served as a validation set to assess the expression differences of the selected characteristic genes. Subsequently, these genes were further validated through qPCR, ultimately leading to the identification of the PANoptosis diagnostic gene. Multivariate logistic regression analysis was conducted for modeling, using GSE110147 as the training set to develop the model, while GSE53845 was used as the validation set. The ROC curve was generated, and the area under the curve (AUC) was calculated. Finally, the “rms package” was employed to construct a nomogram, and a calibration curve was utilized to further assess the stability and reliability of the model.

Construction and Evaluation of Prognostic Prediction Model

The expression matrix of PANDEGs in GSE70866 was extracted, and the dataset contains survival information. Initially, univariate Cox regression analysis was employed to identify prognostic differential genes. Subsequently, LASSO regression and support vector machine-recursive feature elimination (SVM-RFE) were utilized to screen for characteristic genes. SVM-RFE scores and ranks the feature genes, selecting the top 15 genes with the lowest error rates. The results from the two machine learning methods were then intersected. Based on the identified independent prognostic factors, the “rms” R package was used to establish a prognostic nomogram, draw the ROC curve, and calculate the area under the curve (AUC). In generating the nomogram, patients in both the training and validation cohorts were categorized into high and low mortality risk groups according to the optimal cutoff value of the risk score. To ensure the accuracy of risk typing based on PANoptosis prognostic genes, we also conducted random permutation tests. The significance of the performance of these PANoptosis prognostic genes was then evaluated by comparing mean difference in survival time of the prognostic prediction model data and permutated datasets which have labels that were randomly shuffled 1000 times.

Single-Cell RNA Statistical Processing

Single-cell sequencing dataset analysis has been used to explore IPF with promising results.24–26 We analyzed the single-cell RNA sequencing (scRNA-seq) dataset GSE122960 following the standard procedures outlined in the “Seurat” package.27 To ensure the quality of the dataset, we excluded cells with gene expression levels below 300 genes or above 7000 genes, as well as those with more than 20% mitochondrial genes. We employed the “harmony” R package to integrate functions and mitigate batch effects. Unsupervised cell clustering was performed using graph-based methods and visualized through uniform manifold approximation and projection (UMAP). A similar methodology was applied for subcluster analysis. Annotation of cell clusters was conducted based on prior studies.25,28 We utilized the “AddModuleScore” function to calculate the PANoptosis prognostic gene score for each cell, and compared the differences between the normal group and the IPF group.

Pseudotime Trajectory Analysis and Cell–Cell Communication Analysis

Dimensionality reduction and trajectory reconstruction were conducted using the “Monocle” R package with the DDR-Tree algorithm to investigate macrophage trajectories. To identify potential interactions both among different macrophage groupings and with other cell populations, we employed the “CellChat” R package for a comprehensive analysis of intercellular communication molecules.

Validation of Expression of PANoptosis Genes in Bleomycin Induced Pulmonary Fibrosis Mouse Model

C57BL/6N male mice, aged 6–8 weeks and weighing an average of 20–25 g, were utilized for this study. Following anesthesia, the experimental mice were intubated and administered bleomycin (5 mg/kg) to induce pulmonary fibrosis, while the control group received normal saline. Twenty-one days post-bleomycin instillation, lung tissues from the mice were harvested. A portion of these tissues was subjected to hematoxylin and eosin (HE) and Masson staining, with lung tissue paraffin sections prepared and stained according to the manufacturer’s instructions. The remaining tissue was utilized for quantitative polymerase chain reaction (qPCR) analysis. The β-actin gene was used to normalize the expression of various genes. The primers used to detect mRNA levels are listed in Supplemental Table S1.

Identification of Candidate Drugs

L1000 Fireworks Display (L1000FWD) offers an interactive visualization of over 16,000 gene expression feature sets induced by drugs and small molecules (https://maayanlab.cloud/l1000fwd/).29 We utilized the L1000 FWD tool to analyze four PANoptosis prognostic genes that exhibit either high or low expression levels in patients with IPF and are correlated with poorer prognosis. Through this analysis, we identified small molecule drugs that are inversely associated with these genes, thereby highlighting potential drug candidates for the treatment of IPF.

Preparation of Ligand-Receptor Structure and Molecular Docking

Ligand (candidate drug molecule) structure preparation: the 2D structure files of the primary active compounds were obtained from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/). Receptor (target protein) preparation: the PDB structure file of the target protein was sourced from the PDB database (http://www.rcsb.org/). For target proteins not available in the PDB database, Alphafold2 was employed for homology modeling.30,31 Schrodinger molecular docking software was utilized to assess ligand-receptor interactions. The ligand files were imported into Schrodinger Maestro software, where the LigPrep module was employed with the default configuration of the OPLS4 force field, allowing for potential ionization at physiological pH and consideration of optical isomers of the ligands to obtain optimized ligands. Subsequently, the PDB file of the receptor was imported into Schrodinger Maestro software, and the protein preparation wizard was utilized to prepare the receptor with default settings. Finally, molecular docking studies were conducted using the Glide module of Schrodinger.

Statistic Analysis

All statistical analyses in this study were completed using R v4.3.2 software. The Wilcoxon test was utilized for comparing two groups. Spearman correlation test is used for correlation analysis. Kaplan-Meier (KM) curves and forest plots were employed to visualize the results of both univariate and multivariate Cox regression analyses. The Log rank test was used to assess prognosis-related differences. A p-value < 0.05 was considered statistically significant.

Results Identification of PANDEGs Between IPF and Control

First, principal component analysis (PCA) was conducted on the IPF group and the normal group within the GSE110147 data set, revealing a clear contrast between the two groups (Figure 2A). Based on the established threshold, a total of 3175 differentially expressed genes (DEGs) were identified in both the control and IPF groups, comprising 1198 down-regulated genes and 1977 up-regulated genes. The expression of these DEGs was visualized using heat maps and volcano plots (Figure 2B and C). The PANoptosis gene set included genes associated with pyroptosis, necrosis, and apoptosis, totaling 485 genes (Supplemental Table S2). The intersection of the 3175 DEGs and the PANoptosis gene set yielded 104 PANDEGs (Figure 2D). The expression of these PANDEGs was visualized through heatmaps (Figure 2E).

Figure 2 Identification of PANDEGs in IPF. (A) Principal component analysis results of IPF and control group in GSE110147, PCA suggests indicating clear discrimination between IPF patients and control. (B) DEGs heatmap obtained from differential analysis of GSE110147 dataset. (C) The volcano plot of DEGs, contains 1198 upregulated genes and 1977 downregulated genes. (D) The Venn diagram shows the cross genes between GSE110147 DEGs and PANoptosis related genes, there are 104 intersecting genes. (E) Clustering heatmap of PANDEGs.

Functional Enrichment Analysis of PANDEGs

GO and KEGG enrichment analyses were conducted to identify the relevant signaling pathways and biological functions associated with the PANDEGs. The results of the GO enrichment analysis indicated that these biological processes primarily involved the regulation of apoptosis signaling pathways, cytokine production, viral responses, and protein catabolism, among others (Figure 3A). The analysis of cellular components revealed a significant focus on proteasome complexes, endopeptidase complexes, secretory lumen granules, and cytoplasmic vacuoles (Figure 3B). Furthermore, the enrichment analysis of molecular functions highlighted activities such as protein serine/threoninase activity, ubiquitin-like protein ligase activity, and cytokine receptor binding (Figure 3C). The KEGG analysis results demonstrated that PANDEGs were notably enriched in pathways related to Alzheimer’s disease, necroptosis, Epstein-Barr virus infection, prion disease, influenza A, and the NOD-like receptor signaling pathway, among others (Figure 3D).

Figure 3 Functional Enrichment analysis of PANDEGs. (A) The top 10 GO enrichment analyses of biological processes. (B) The top 10 GO enrichment analyses of cellular component analysis. (C) The top 10 GO enrichment analyses of molecular function analysis. (D) The top 10 pathways for KEGG enrichment analysis.

Correlation of Immune Infiltration to PANDEGs

Immune cells play a crucial role in the onset and progression of IPF.32,33 Enrichment analysis indicated that PANDEGs are associated with various immune-related signaling pathways. Consequently, we assessed the correlation between immune cell infiltration and PANDEGs. The results revealed that among the 24 immune cell types examined, 17 exhibited significant differences between the IPF group and the normal group (Figure 4A and B). Additionally, the immune cells demonstrated notable correlations (Figure 4C). Specifically, macrophages, CD4 T cells, neutrophils, and nTreg cells et al were predominantly enriched in the IPF group, while the control group showed a higher enrichment of DC cells, B cells, and monocytes et al. Due to the crucial role of macrophages in the occurrence and development of IPF,34 we focus on macrophages. In IPF patients, the median enrichment score of macrophages is 0.02255, with a distribution range of (Lower Whisker: 0.019, Upper Whisker: 0.282). In normal patients, the median enrichment score of macrophages is 0.064, with a distribution range of (Lower Whisker: 0, Upper Whisker: 0.08). Therefore, the enrichment degree of macrophages is significantly increased in IPF patients. The results of the correlation test indicated that PANDEGs were primarily positively correlated with macrophages, CD4 T cells, and neutrophils et al. In contrast, a negative correlation was observed with monocytes, NK cells, B cells, and DC cells et al (Figure 4D).

Figure 4 Analysis results of immune infiltration. (A) The box plot displays immune cell infiltration and the differences in immune cell infiltration between IPF and normal samples. The boxes represents the quartiles of the data, whiskers indicating the degree of dispersion of the data, and the dots represent outliers (The box plot’s legend in the following figures also applies). (B) Cluster heatmap of 24 immune cell proportions in GSE110147 dataset. (C) Matrix diagram of 24 immune cell correlations. Green is positively correlated, red is negatively correlated, and the depth of the color indicates the strength of the correlation, × means no statistically significant. (D) Heatmap of the correlation between PANDEGs and 24 types of immune cells.

Machine Learning-Based Selection of Diagnostic Genes

Machine learning methods are increasingly employed in the diagnosis and prognosis of IPF.35–37 Initially, ten predicted signature (PPP3CB, PDCD4, PSMD8, AKT1, PSMA2, SOD1, PSMA4, CHMP2B, PSMD7, PSMD14) genes were selected from PANDEGs by LASSO regression analysis (Figure 5A and B). Subsequently, these genes were validated in the GSE53845 dataset. The results indicated that a total of 6 genes exhibited significant differences, with 5 genes displaying consistent trends in expression change (Figure 5C). QPCR was employed to confirm the expression levels of the 5 previously identified genes, and the results indicated that the expression levels of 4 genes—AKT1, PDCD4, PSMA2, and PPP3CB—were statistically significant (Figure 5D–H). Consequently, these 4 genes were ultimately designated as PANoptosis diagnostic genes in IPF.

Figure 5 Screening of PANoptosis diagnostic genes. (A) Lasso regression coefficient versus Log (λ) variation curve. (B) Mean square error in Lasso regression varies with Log (λ). (C) Validation of differential expression of candidate PANoptosis diagnostic genes in external dataset. (D-H) Differential expression verification of candidate PANoptosis diagnostic genes was performed using qPCR.PPP3CB, PDCD4, AKT1 and PSMA2 have statistically significant.

Establishment and Evaluation of the IPF Diagnostic Prediction Model

Initially, a correlation heat map and a circular diagram illustrating the correlation between the 4 diagnostic genes associated with IPF (Figure 6A and B). Subsequently, the ROC curve for these 4 genes was generated using the validation set GSE53845, and the AUC was calculated. AKT1 (AUC=0.818), PDCD4 (AUC=0.876), PSMA2 (AUC=0.817) and PPP3CB (AUC=0.873) were obtained (Figure 6C–F). We conducted a multivariate logistic regression analysis on 4 genes and utilized ROC curves to assess the models. In the training set, the model achieved an AUC of 1 (Figure 6G), while in the validation set, the AUC was 0.981 (Figure 6H). These results demonstrate that the diagnostic model possesses an exceptionally high diagnostic value. Subsequently, we constructed a diagnostic nomogram based on the characteristic genes (Figure 6I), and evaluated the model’s predictive performance with a calibration curve (Figure 6J). The results from the calibration curve indicate that the predicted probabilities of the model closely align with the actual outcomes, suggesting that the model exhibits considerable accuracy.

Figure 6 Diagnostic model of IPF was constructed and evaluated. (A-B) Correlation analysis of 4 PANoptosis diagnostic genes. (D-F) ROC curves and AUC of 4 genes in the verification set. (G) ROC curve and AUC of diagnostic model of IPF. (H) ROC curve and AUC of validation set. (I) Prediction of IPF using nomogram. (J) Calibration curve of nomogram.

Machine Learning-Based Selection of Prognostic Prediction Model

A univariate Cox regression analysis was conducted on 104 PANDEGs to identify prognosis-related genes, resulting in a forest plot that revealed 26 genes with statistically significant differences (Figure 7A). Subsequently, LASSO regression was applied to select 16 characteristic genes (Figure 7B and C). The SVM-RFE algorithm identified the top 15 genes with the lowest error rate, which were then intersected with the genes selected by LASSO regression, yielding a final set of 8 genes (Figure 7D). These 8 genes underwent multivariate Cox regression analysis, revealing that TNFRSF12A, DAPK2, UACA, and DSP were independent prognostic factors (Figure 7E). Kaplan-Meier curves illustrated survival differences between high and low expression groups for each gene (Supplementary Figure 1A-D). Consequently, these four genes are classified as PANoptosis prognostic genes.

Figure 7 Screening of PANoptosis prognostic genes. (A) Forest plot of PANDEGs with P <0.05 by univariate Cox regression analysis. (B) Lasso regression coefficient versus Log (λ) variation curve. (C) Mean square error in Lasso regression varies with Log (λ). (D) Venn diagram shows the 16 genes selected by LASSO regression analysis and the top 15 genes with the highest accuracy in SVM-REF algorithm. (E) Multivariate Cox regression analysis of intersecting genes.

Establishment and Evaluation of the IPF Prognostic Prediction Model

Based on the results of multivariate Cox regression analysis, a prognostic nomogram was established (Figure 8A). According to the expression levels of 4 genes and the corresponding coefficients of the prognostic model, allowing for the calculation of a risk score for each patient. Patients were classified into high-risk and low-risk groups according to the median risk score. In the training set, the prognosis of the high-risk group was significantly worse than that of the low-risk group (Figure 8B). The ROC curve indicates that the prognostic model possesses good predictive value, with an AUC of 0.741 (Figure 8C). Furthermore, clinical information was included in the GSE700886 dataset for multivariate Cox regression analysis, we observed that the risk score remained an independent prognostic factor for IPF patients (Supplementary Table S3). The GSE93606 dataset was utilized as the validation set for the model, and the KM curve further demonstrated that the prognosis of the high-risk group was significantly poorer than that of the low-risk group, with an AUC of 0.66 (Figure 8D and E). To evaluate the validity of the result and ensure that the good predictive performance is not due to random chance, we used random permutation test to further evaluate. The results showed that the red dashed line representing the prognostic model fell on the far right side of the image, significantly larger than the dataset with random permutation (p<0.0001)(Figures 8F). Therefore, the prognostic model was significantly better than random prediction. To further validate the expression of prognostic genes, we constructed a bleomycin-induced pulmonary fibrosis mouse model. HE and Masson staining results revealed that the lung structure of the control group appeared normal, with no evident inflammatory cell infiltration or fibrosis. Conversely, in the bleomycin group, the alveolar walls were notably thickened, a substantial infiltration of inflammatory cells was observed, and significant collagen deposition was noted (Figure 8G). QPCR was employed to confirm the expression of these 4 PANoptosis prognostic genes, and the results indicated statistically significant differences (Figure 8H-K).

Figure 8 Prognostic model of IPF was constructed and evaluated. (A) The nomogram for IPF survival. (B) Kaplan-Meier curves of high-risk and low-risk groups in the training set. (C) ROC curve and AUC of prognostic model of IPF. (D) Kaplan-Meier curves of high-risk and low-risk groups in the validation set. (E) ROC curve and AUC of of validation set. (F) Distribution of mean difference in survival time of randomly permutated datasets (n=1000), the validation on the prognostic prediction model data is indicated by the red dotted line. (G) HE and Masson staining results of control group and bleomycin group. (H-K) Differences in prognostic gene expression between control group and ipf group.

Single-Cell Level Expression of PANoptosis

We utilized the GSE122960 scRNA-seq dataset to elucidate the inherent cellular heterogeneity present in lung tissue, which comprised samples from 8 transplant donors and 4 patients with IPF. Following quality control procedures outlined in the methods section, we obtained a total of 57,050 cells for subsequent analysis. Drawing on marker genes identified in prior studies,25 we annotated a total of ten distinct cell types, including AT II cells, AT I cells, Macrophages, Monocytes, B cells, T&NK, Ciliated cells, Endothelial cells, Club cells, and Fibroblasts (Figure 9A and B). We employed the “AddModuleScore” function to evaluate the prognostic scores of PANoptosis-related genes across individual cells. The results indicated significant differences in scores among AT II cells, AT I cells, Macrophages, Monocytes, and T&NK cells between donor and IPF patients (Figure 9C). Given that immune infiltration analysis reveals a predominant enrichment of macrophages in IPF, and recognizing the critical role that macrophages play in the onset and progression of IPF,38–40 we focused our analysis specifically on macrophages. We classified macrophages into 6 distinct types: lung resident Macrophages, clec4e high Macrophages, profibrotic Macrophages, proinflammatory1 Macrophages, proinflammatory2 Macrophages and MTs high Macrophages (Figure 10A and B). It should be noted that because proinflammatory macrophages are distributed far apart on the UMAP, they are subdivided into 2 subgroups based on their hypervariable genes. To explore the sequential development of different macrophage subpopulations in IPF, we conducted a pseudotime trajectory analysis of macrophages. At the beginning of the trajectory, proinflammatory2 macrophages and lung resident macrophages were predominantly observed. Conversely, profibrotic macrophages were identified at the end of the trajectory (Figure 10C-E). Each cell was evaluated for PANoptosis prognostic genes using the “AddModuleScore” function (Figure 10F). The results indicated that lung resident macrophages, clec4e-high macrophages, and MTs-high macrophages exhibited higher scores in IPF patients, whereas proinflammatory2 macrophages displayed lower scores in the IPF group. Specifically, the Median PANoptosis score for lung resident macrophages in IPF patients is −0.0247 (Lower Whisker: −0.0734, Upper Whisker: 0.0224), donor group is −0.0431 (Lower Whisker: −0.1, Upper Whisker: 0.0235). And, the Median PANoptosis score for proinflammatory2 macrophages in IPF patients is −0.05 (Lower Whisker: −0.127, Upper Whisker: 0.171), donor group is −0.036, (Lower Whisker: −0.124, Upper Whisker: 0.386).

Figure 9 The expression of PANoptosis at the IPF cellular level. (A) UMAP plot is colored by different cell types, annotated 10 different cell types. (B) The expression of marker genes in 10 different cell types. The color is determined by the average expression level, with yellow being higher and blue being lower. The size of the circle represents the percentage of expression. (C) The box plot display the scoring of PANoptosis prognostic genes in different cell types of IPF and Donor groups results.

Figure 10 The expression of PANoptosi and trajectory of macrophage population. (A) UMAP plot is colored by different macrophage subclusters, annotated 7 different cell types. (B) The expression of marker genes in 7 different macrophage subclusters. The color is determined by the average expression level, with yellow being higher and blue being lower. The size of the circle represents the percentage of expression. (C) Macrophage development trajectory color coded by status. (D) Macrophage development trajectory color coded by cell subclusters. (E) Macrophage development trajectory color coded by pseudotime. (F) The box plot display the scoring of PANoptosis prognostic genes in different macrophage subclusters of IPF and donor groups results.

Cell-Cell Communication

To elucidate potential interactions between macrophages and other cell populations, we conducted a cell-cell communication analysis, which is based on calculations of ligand-receptor gene expression. We categorized macrophages into 2 types according to the previously established macrophage PANoptosis prognostic genes score, utilizing the median value to distinguish between PANoptosis high macrophages and PANoptosis low macrophages. The construction of the cell-cell communication network is founded on the number of interactions and interaction weights (Figure 11A). The results indicate that PANoptosis high macrophages exhibit enhanced intercellular communication capabilities. The outgoing and incoming signaling patterns of macrophages and other cell populations reveal that the signal transmission intensity of PANoptosis high macrophages is significantly greater than that of PANoptosis low macrophages (Figure 11B and C). In comparison to PANoptosis low macrophages, PANoptosis high macrophages are capable of additional cellular communication via the MIF pathway and the GALECTIN pathway.

Figure 11 Cell communication analysis in macrophage subclusters. (A) Circle plot depicts the number and strength of ligand receptor interactions between paired cell populations. The thicker the line, the stronger the interaction. The same color of the line and the point indicates that the point is a source cell, and the other end is a target cell. (B) The outgoing and incoming signal patterns of different cell populations. The darker the color of the square, the stronger the outgoing and incoming signal patterns. (C) Signal pathways involved in the interaction between incoming and outgoing in two distinct subtypes of macrophages. The dots represent the communication probability, which gradually increases from blue to red.

L1000FWD Screening for Drug Candidate Molecules

We utilized the L1000 FWD online platform to perform a database search for 3 upregulated genes and 1 downregulated gene in patients with IPF who exhibit poor prognosis. Subsequently, candidate small molecule drugs with opposing correlations were identified, with a focus on the comprehensive score, similarity score, and p-value. The top five candidate drugs are listed in Table 1. Taking into account the ranking and availability of these drugs, we ultimately selected the top 3 candidates drugs (Metergoline, Candesartan, Selumetinib) for further molecular docking.

Table 1 Top Five Drugs of Opposite Relevance in the L1000FWD

Molecular Docking Validation

We obtained comprehensive information regarding the drug candidates from the PubChem database. Specifically, Metergoline (PubChem ID: 28693) has a molecular weight of 403.5 g/mol and a molecular formula of C25H29N3O2. Candesartan (PubChem ID: 2541) has a molecular weight of 440.5 g/mol and a molecular formula of C24H20N6O3. Selumetinib (PubChem ID: 10127622) has a molecular weight of 457.7 g/mol and a molecular formula of C17H15BrClFN4O3. The three-dimensional structures of the target proteins are as follows: DSP (PDB: 1LM5), DAPK2 (PDB: 1Z9X), TNFRSF12A (PDB: 2RPJ), and UACA (AlphaFoldDB: AF-D3ZGS5-F1). We conducted high-precision molecular docking utilizing Schrodinger’s Glide module. Generally, a lower binding energy between the ligand and receptor indicates a more stable structure.41,42 A binding energy of less than −1 kcal/mol indicates that the ligand and receptor can spontaneously bind. The binding energies between candidate drugs and target proteins are presented in Table 2. Specifically, the binding energy of Selumetinib to each target protein is lower than −4 kcal/mol, indicating excellent binding activity. The binding energy of Metergoline to DSP is −4.072 kcal/mol, while its binding energy to DAPK2 is −7.124 kcal/mol, demonstrating extremely strong binding activity with DAPK2. In contrast, the binding energy of Candesartan to each target protein is less than −2, suggesting a certain level of binding activity. We utilized Schrodinger software to visualize the interactions between each candidate small molecule drug and the target protein (Figure 12). Table 3 presents the drug candidates along with their respective residues that interact through hydrogen bonds. The results indicate that the candidate small molecule drugs and the target proteins are primarily connected by hydrogen bonds, thereby forming a stable complex.

Table 2 Binding Energy of Drugs and Targets

Table 3 Molecular Docking Hydrogen Bonds and Residues

Figure 12 Validation of molecular docking of candidate drugs with its target. (A-C). Molecular docking between DSP and Candesartan(A), Selumetinib(B), Metergoline(C). (D-F) Molecular docking between DAPK2 and Candesartan(D), Selumetinib(E), Metergoline(F). (G-I) Molecular docking between DAPK2 and Candesartan(G), Selumetinib(H), Metergoline(I). (J-L) Molecular docking between TNFRSF12A and Candesartan(J), Selumetinib(K), Metergoline(L).

Discussion

PANoptosis, a distinctive form of programmed cell death introduced by Malireddi et al in 2019,17 is triggered by innate immunity and regulated by the PANoptosome complex. This process incorporates pyroptosis, apoptosis, and necroptosis, leading to biological effects that arise from overlapping regulatory mechanisms among other programmed cell death pathways. However, these effects cannot be solely attributed to any one of these individual pathways.43 Previous studies have demonstrated that PANoptosis contributes to the development of various acute and chronic pulmonary diseases.20,21 Additionally, PANoptosis has been implicated in the progression of several fibrotic diseases, such as myocardial fibrosis and liver fibrosis,44,45 and is being explored as a potential therapeutic strategy. However, its specific role in IPF remains unclear. This study represents the first comprehensive analysis of PANoptosis in the development and progression of IPF. We employed machine learning techniques to construct diagnostic and prognostic models associated with PANoptosis and conducted a detailed analysis of its role. Potential small-molecule drugs (Metergoline, Candesartan, Selumetinib) were identified using the L1000FWD application based on prognostic genes related to PANoptosis and subsequently validated through molecular docking studies. To our knowledge, this is the first bioinformatics investigation to elucidate the relationship between PANoptosis and IPF in human samples.

In the present study, we initially employed differential analysis to identify 104 PANDEGs. We then conducted GO and KEGG enrichment analyses to elucidate the relevant signaling pathways and biological functions associated with PANDEGs, revealing predominant enrichment in the immune-inflammatory response pathway. Subsequently, we assessed the expression profiles of immune cells in IPF and examined their correlation with PANDEGs. The results showed that PANDEGs predominantly positively correlated with macrophages, CD4+T cells, and neutrophils. Previous research has demonstrated that macrophages undergo distinct phenotypic and functional transitions at various stages of IPF progression.46,47

Machine learning and bioinformatics are widely used in the diagnosis and prognostic model construction of IPF patients. In terms of diagnostic models for IPF, Shi et al constructed a cuprotosis related diagnostic model with AUC values of 0.729 and 0.700 in the training and validation sets.48 Similarly, Liao et al constructed a Lipid Related Biomarker diagnostic model for IPF, this model has excellent diagnostic ability, with AUC values greater than 0.9 in multiple datasets and an average AUC value of 0.966.49 In this study, we identified four PANoptosis diagnostic genes (AKT1, PDCD4, PSMA2, and PPP3CB) using machine learning and validated with qPCR. An IPF diagnostic nomogram was then constructed based on these genes. The AUC value of the ROC curve for this model in the training group was 1, while it reached 0.981 in the validation group, demonstrating the model’s exceptional diagnostic performance. Protein kinase B (AKT1) is a serine/threonine protein kinase that is a crucial component of the PI3K/AKT/mTOR signaling pathway.50 Regulation of AKT1 pathway has been shown to improve pulmonary fibrosis.51–53 Programmed Cell Death 4 (PDCD4) is a pro-apoptotic protein that inhibits tumor transformation,54 and promotes fibroblast differentiation in both liver and renal fibrosis through TGF-β and other regulatory pathways.55,56 PSMA2, is one of the seven alpha subunits of the 20S proteasome and facilitates intracellular proteolysis in an ATP- and ubiquitin-independent manner.57 It can also mitigate the oxidative stress response induced by influenza A virus by downregulating NRF2.58 PPP3CB, a member of the phosphoprotein phosphatases (PPPs) group,59 plays a critical role in regulating cell migration and epithelial-mesenchymal transition (EMT) process.60 In fibrotic diseases, PPP3CB causes cyclosporine A-induced secretion of MMP-9, which contributes to renal fibrosis.61 Although our diagnostic model has excellent diagnostic capabilities, it still needs to be noted that the sample size included is relatively small. In addition, the construction of diagnostic models is based on RNA sequencing of lung tissue, which requires invasive biopsy and may be a significant harm to patients.

The prognosis for patients with IPF is poor, with survival rates exhibiting significant variability among individuals.4 Therefore, identifying patients at high risk of developing IPF is crucial. Currently, some studies have screened prognostic biomarkers for IPF. Liao et al developed a prognostic model based on anoikis genes, which predicts AUC values of 0.784, 0.779, and 0.788 for the 1, 2, and 3 year survival rates of IPF patients.62 Fan et al constructed a prognostic model based on five autophagy related genes,63 Zhu et al developed a prognostic model based on endoplasmic reticulum stress-related genes,64 both have good predictive ability. Our prognostic prediction model for IPF demonstrated the capability to accurately identify high-risk patients, with an AUC value of 0.741 in the training set and 0.66 in the validation set. However, the model developed in this study still has certain limitations, and it is important to acknowledge the differences in sample sources between the training and validation datasets. Due to the fact that the sequencing results used in the prognostic model come from bronchoalveolar lavage fluid and blood, which are relatively easy to obtain, it is worth further promotion and research. Single-cell analysis enables the comparison of expression profile differences across samples at the cellular level as well as the investigation of differentiation trajectories and their interactions among key cell populations.65–67 Our single-cell analyses, revealed that the PANoptosis prognostic gene scores varied among AT II cells, AT I cells, and macrophages, and that the macrophage score in the IPF group was significantly higher than that in the normal group. Under healthy conditions, tissue-resident macrophages (TR-AMs) and interstitial macrophages (IMs) are predominantly found in the lungs. Under pathological conditions, such as during the acute inflammatory reaction induced by lung infection, circulating monocytes are recruited into lung tissue, where they migrate into the alveoli, and differentiate into macrophages.68 Using pseudotime trajectory analysis, we identified a differentiation trajectory primarily characterized by proinflammatory macrophages in the early stages of inflammation. As the disease progresses, lung-resident macrophages gradually differentiate into profibrotic macrophages. Scores for macrophage subclusters revealed that lung-resident macrophages exhibited higher scores in IPF patients, whereas proinflammatory macrophages demonstrated lower scores in the IPF group. These findings suggest that PANoptosis may influence the conversion of lung resident macrophages into profibrotic macrophages. Additionally, we compared the communication between these two macrophage types and other cell types. In comparison to low PANoptosis macrophages, high PANoptosis macrophages exhibit enhanced cell-cell communication via the MIF and GALECTIN pathways. MIF inhibition has been shown to mitigate bleomycin-induced pulmonary fibrosis.69,70 Additionally, elevated MIF levels have been identified as an independent risk factor through clinical investigations for 3-month mortality in IPF patients during acute exacerbations.71 GALECTIN, an S-type lectin, encompasses a family of molecules that includes GAL-1, −2, −3, −5, −7, −8, −9, and −10, among others.72 Notably, GAL-1 and GAL-3 have been implicated in the development of fibrosis based largely on studies in galectin-deficient mice.73,74 TD139, a Gal-3 inhibitor, has shown antifibrotic efficacy in preclinical models of pulmonary fibrosis and in patients with IPF.69,70,72 Consequently, variations in PANoptosis activity within macrophages may influence other cell types, thereby modulating IPF progression through these specific pathways.

Current understanding suggests that the pathogenesis of IPF involves multiple pathways, which explains the reason for single-target therapies being ineffective against this condition. Similar to cancer treatment, the future direction of IPF therapy is likely to involve multidrug and multitarget combination strategies.75 Consequently, we employed the L1000FWD online tool to search database for 3 up regulated genes and 1 down regulated gene associated with poor prognosis in IPF patients. Subsequently, drugs exhibiting negative correlations were screened, and based on their availability and ranking, Metergoline, Candesartan, and Selumetinib were selected for further investigation. Metergoline, an ergot derivative, is used to treat hyperprolactinemic amenorrhea in women.76 Numerous studies have demonstrated that metergoline irreversibly blocks 5-HT7 receptors.77–79 SB-269970, another 5-HT7 antagonist, has been shown to mitigates bleomycin-induced pulmonary fibrosis by down regulating oxidative stress and inflammation.80 Candesartan, an angiotensin II receptor antagonist, is primarily used to treat hypertension and chronic heart failure.81 The renin-angiotensin-aldosterone system (RAAS) is extensively implicated in the development of fibrosis in various organs including the heart, liver, kidneys, and lungs.82 Studies have shown that candesartan mitigates bleomycin- and silica-induced pulmonary fibrosis.83,84 Selumetinib, an inhibitor that targets the downstream effector protein MEK within the RAS-RAF-MEK-ERK (MAPK) signaling pathway, is currently used for treating various tumors.85 The MAPK pathway regulates cellular processes associated with fibrosis, including cell proliferation, apoptosis, and myofibroblast transformation.86,87 Studies have indicated that MEK inhibition can reduce bleomycin-induced pulmonary fibrosis.88,89 To further investigate the effects of these three small-molecule inhibitors on PANoptosis prognosis, we conducted molecular docking studies. These results indicated that these small-molecule inhibitors possessed the capacity to spontaneously bind to the four PANoptosis prognostic genes. However, it should be noted that the restrictive sampling of ligand and receptor conformations and the use of approximate scoring functions may produce results that are not correlated with the actual experimental binding affinity.

Our study provides valuable insights into the relationship between specific PANoptosis genes and the development of IPF; however, it has certain limitations. Firstly, the relatively small sample size, may impact the robustness of our findings. Although we performed remove batch effect, heterogeneity in the data may be present. Second, we lacked clinical samples from patients with IPF for validation purposes. Finally, the mechanisms by which these characteristic genes influence IPF through PANoptosis warrant further investigation.

Conclusion

In summary, we explored the potential role of PANoptosis in IPF and developed both diagnostic and prognostic models for the disease. These models were externally validated and corroborated in animal models using qPCR. Furthermore, our single-cell analysis revealed alterations in PANoptosis activity in IPF and its impact on cell-cell communication. Additionally, we identified potential therapeutic drugs for patients with IPF based on PANoptosis prognostic genes and assessed their effectiveness through molecular docking studies. However, the clinical efficacy of these drug candidates requires further investigation.

Data Sharing Statement

Publicly available datasets were analyzed in this study. This data can be found here: https://www.ncbi.nlm.nih.gov/geo/; GSE110147, GSE53845, GSE70866, GSE93606 and GSE122960. The original contributions presented in the study are included in the article/supplementary material. Further inquiries can be directed to the corresponding authors.

Ethics Approval

This study received approval from the Ethics Committee of Zunyi Medical University Hospital (KLL-2020-165), and all applicable institutional and governmental regulations concerning the ethical use of animals were followed. Our study is exempt from approval based on national legislation guidelines regarding human data. Specifically, it is item 1 and 2 of Article 32 of the Measures for Ethical Review of Life Science and Medical Research Involving Human Subjects dated February 18, 2023, China. The legislative content of these two items is that research conducted using legally obtained public data, data generated through observation without interfering with public behavior, or research conducted using anonymized information data can be exempted from ethical approval.

Acknowledgments

We thank the Chongqing University Central Hospital, for providing the experimental platform. We also thank Figdraw (www.figdraw.com) for the assistance in creating graphical abstract.

Author Contributions

All authors made a significant contribution to the work reported, whether that is in the conception, study design, execution, acquisition of data, analysis and interpretation, or in all these areas; took part in drafting, revising or critically reviewing the article; gave final approval of the version to be published; have agreed on the journal to which the article has been submitted; and agree to be accountable for all aspects of the work.

Funding

This study was supported by National Trauma Regional Medical Center Major Research Project (Co-built by the Municipal Commission) (jjzx2021-gjcsqyylzx01), research Project of Chongqing Talent Program (cstc2022ycjh-bgzxm0245), Eagle Talent Project of Chongqing Emergency Medical Center(SYRCCY20230312), key project co-organized by the Health Commission and the Science & Technology Bureau of Chongqing Province (2024ZDXM024), General Project of Chongqing Province Natural Science Foundation (CSTB2024NSCQ-MSX0873), Chongqing Key Laboratory of Emergency Medicine(2023-KFKT-03).

Disclosure

The authors report no conflicts of interest in this work.

References

1. Richeldi L, Collard HR, Jones MG. Idiopathic pulmonary fibrosis. Lancet. 2017;389(10082):1941–1952. doi:10.1016/S0140-6736(17)30866-8

2. Maher TM. Interstitial lung disease: a review. JAMA. 2024;331(19):1655–1665. doi:10.1001/jama.2024.3669

3. Lu W