The gene expression matrix and clinical characteristics of IPF patients were retrieved from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/), and a total of 505 IPF patients and 147 healthy volunteers were included. The bronchoalveolar lavage fluid (BALF) came from the GSE70866 series, including IPF patients (n = 176) and healthy volunteers (n = 20), which were detected using the GPL14550 and GPL17077 platforms, respectively (Prasse et al., 2019). The lung tissue came from the GSE47460 series, including interstitial lung disease patients (n = 254) and healthy volunteers (n = 108), which were tested using the GPL6480 and GPL14550 platforms, respectively (Kim et al., 2015; Peng et al., 2016; Anathy et al., 2018). The source of peripheral blood mononuclear cells (PBMCs) was the GSE28221 series, including IPF patients (n = 75) and healthy volunteers (n = 19), which were detected using the GPL6480 platform (Herazo-Maya et al., 2013; Huang et al., 2015). Due to the different detecting platforms of samples, the batch effect detected by the platforms was eliminated using the ComBat method in the “SVA” R package for data from the same histologic origin. Principal component analysis (PCA) was used to determine the goodness of batch correction.

Based on previous studies, 36 m6A regulators were considered, including 19 readers, 14 writers, and three erasers. The m6A regulators described above were extracted using the “limma” package in R software. The Wilcoxon test was used to determine whether there were significant changes between IPF patients and healthy volunteers. The “ggpubr” package was employed to draw expression boxplots, and the “pheatmap” package was used to draw a visual m6A expression heatmap.

We constructed machine learning classifiers of RF and SVM to predict the reliability of IPF occurrence (Noble, 2006; Altmann et al., 2010). A hyperplane was identified to distinguish IPF patients from healthy volunteers. Then, we established a nomogram to predict the occurrence of IPF in different GEO series based on five selected m6A regulators. Statistical analysis was performed using the R packages “rms” and “rmda.” To assess model accuracy and validity, plot calibration curves, decision curve analysis (DCA), and receiver operating characteristic (ROC) curves. Finally, the nomogram formula of the GSE47460 series was brought into GSE70866 and GSE28221 series for verification. ROC curves were drawn to calculate the area under the curve (AUC) and 95% confidence interval (CI).

The 31 m6A-related genes of the BALF-derived GSE70866 series were analyzed using univariate Cox regression analysis. Multi-variables faced the risk of over-fitting. LASSO regression was applied to analyze the influence of the survival status of these three variables, and the parameter λ was used to adjust the model complexity to obtain a less variable and more representative combination. Prognostic high-risk and low-risk clusters were obtained according to the risk score using LASSO regression. The 33 m6A-related genes of the PBMC-derived GSE28221 series were analyzed using univariate Cox regression analysis.

The risk score of IPF patients can be used as an independent prognostic factor. Among the clinical data, age, gender, and global alignment and proportion (GAP) score were included as independent variables in univariate and multivariate Cox regression analyses to clarify independent prognostic factors for IPF patients. To balance the statistical data, we used 1000× risk scores. To perform these analyses, we used the “forestplot,” “survival,” and “glmnet” packages.

Spearman correlation analysis was performed on the lung function using the diffusing capacity of the lungs for carbon monoxide (DLCO) and m6A regulators (|R| > 0.2, p < 0.001) in the GSE47460 series using the “limma” package. Spearman correlation analysis was performed on METTL14, G3BP2, and ZC3H13 with other m6A regulators (|R| > 0.4, p < 0.001) in the GSE70866 series. The “ggplot2,” “ggpubr,” and “ggExtra” packages were used for visual scatter plot drawing.

To determine whether the outcomes-related m6A regulators METTL14, G3BP2, and ZC3H13 are associated with IPF, series from GSE70866 were divided into groups according to the consensus levels of m6A regulators using the “ConsensusClusterPlus” package. Outputs included consensus cumulative distribution function (CDF) plots, the relative change in area under the CDF curve, and the consensus matrix. PCA was used to determine the fitness of classification. Kaplan-Meier survival analysis was used to evaluate the difference in overall survival between different m6A types using the “survival” package.

The “limma” package and the “VennDiagram” package were used to identify DEGs associated with different m6A modifications. Adjusted p-values <0.05 and |log2FC| > 1 were considered as DEGs, and volcano plots were drawn for all genes. The “clusterProfiler” package was used to analyze GO functional annotation, including the biological process (BP), celluar component (CC), molecular function (MF) terms, and KEGG pathway analysis. The significance criterion was set at p-value <0.05. The results were plotted using a bioinformatics tool (https://www.bioinformatics.com.cn) based on the R language (version 4.03) (Ashburner et al., 2000; Kanehisa and Goto, 2000; Gene Ontology, 2021).

Initially, we used single-sample gene set enrichment analysis (ssGSEA) to quantify the relative abundance of 23 immune cell types in the immune microenvironment across different m6A clusters. Panels of specific functional genes are used to mark each immune cell type. The relative abundance of each immune cell type was calculated from the enrichment fraction of gene set expression in transcriptome sequencing in ssGSEA analysis and normalized to a uniform distribution. Next, based on the results of the immune-related analysis, a boxplot of differentially expressed immune cells of different m6A clusters was drawn, and an association heatmap was drawn by using m6A regulators and immune cell expression abundance for correlation analysis.

The “pRRophetic” package was used to predict potential drugs sensitive to IPF patients with different m6A risk groups. Drug susceptibility was assessed by the half-maximal inhibitory concentration.

To investigate the differentially expressed m6A gene between IPF patients and healthy people, according to the detection ratio of the groups, we selected the GSE47460 series derived from lung tissue for differentially expressed m6A gene analysis. After batch correction, the two data sets with the same sample source were merged, and the m6A genes were extracted (Supplementary Figures S1A–H). Due to the small number of healthy volunteers, only differentially expressed m6A genes with opposite expression trends in the GSE70866 and GSE47460 series were excluded. We extracted 21 differentially expressed m6A genes from the GSE47460 series, eight DEGs with the opposite expression trend, and one undetected m6A gene of GSE47460 was removed (Figure 1A). We obtained 12 reliable and significantly differentially expressed m6A regulators, including six writers and six readers, and presented the heatmap of differentially expressed m6A regulators (Figure 1B).

Modeling of IPF disease-associated m6A regulators was performed using differentially expressed m6A regulators. Simultaneous RF and SVM models were established in the GSE47460 series to identify m6A regulators characterized by IPF. “ROC curve,” “residual box plot,” and “residual reverse cumulative distribution” suggested that the RF model was more suitable than the SVM model for constructing m6A genes related to pathogenesis (Figure 1C and Supplementary Figures S1I,J). The point with the smallest cross-validation error constructed by RF trees was used as the best model, and the importance scores of 12 m6A regulators were calculated (Figure 1D). The top five m6A regulators (PCIF1, RBM15B, CBLL1, SND1, and FMR1) were used to build a nomogram to predict prevalence in healthy volunteers and IPF patients (Figure 1E).

The calibration curve can be seen that the “apparent” dotted line corresponding to the entire cohort and the bootstrap corrected solid line are close to the “Ideal” dotted line, indicating that the model has predictive performance (Figure 2A). The DCA method was used to draw a “decision curve” to evaluate the relationship between the nomogram and the gene score to predict the benefits and risks of various cut points in the IPF prevalence model. Figure 2B shows that the selected m6A genes have better benefits at various threshold probabilities.

An m6A gene-associated nomogram of the GSE47460 derived from lung tissue series was then drawn (Figure 2C). The ROC curve indicated the sensitivity and specificity of the diagnosis. The AUC of the SE47460 series was 89.5% (95% CI: 86.19%–92.89%), suggesting a high prediction accuracy (Figure 2D). Then we brought the five m6A regulators into the GSE70866 and GSE28221 to establish the diagnosis-related nomogram construction. The AUC value of the ROC curve was 78.5% (Supplementary Figures S1K,L) and 80.8% (Supplementary Figures S1M,N).

To validate the nomogram constructed based on the lung tissue-derived GSE47460 series (Figure 2C), we used 196 patients from the BALF-derived GSE70866 series and 94 patients from the PBMC-derived GSE28221 series to perform IPF incidence prediction validation and draw ROC curves. The AUCs for the GSE70866 and GSE28221 series were 69.8% and 68.9%, respectively (Figures 2E,F). The source of BALF is lung tissue, and the predictive accuracy was better. There were few samples in the GSE70866 and GSE28221 series; a larger sample size is needed. These findings suggest that PCIF1, RBM15B, CBLL1, SND1, and FMR1 can be used as the characteristic m6A regulators of IPF, and they have diagnostic significance in lung tissue and PBMC.

To clarify the impact of m6A-related genes on the outcomes of IPF, the clinical information provided by the BALF-derived GSE70866 series was used to perform univariate Cox regression analysis on the 31 m6A regulators extracted from the data set (Figure 3A). Three protective m6A regulators (p < 0.05) (METTL14, G3BP2, and ZC3H13) were identified. LASSO regression was employed to identify these regulators, and risk scores were calculated (Figures 3B,C). We used univariate Cox regression models to compare the effects of age, gender, GAP, and m6A regulator risk score on patient survival. Figure 3D shows that GAP (HR 1.395, 95% CI: 1.237–1.574, p < 0.001) and the m6A regulator risk score (HR 1.430, 95% CI: 1.181–1.732, p < 0.001) predict IPF prognosis.

Multivariate Cox regression analysis was performed to control for the influence of confounding factors and to analyze the influencing factors with independent effects. Figure 3E shows that the GAP score (HR 1.474, 95% CI: 1.287–1.687, p < 0.001) and m6A regulator risk score (HR 1.420, 95% CI: 1.180–1.709, p < 0.001) remained IPF predictors. The BALF-derived GSE28221 series was used to perform univariate Cox regression analysis on the 33 m6A regulators extracted from the data set (Supplementary Figures S2A). YTHDC1 was a protective m6A regulator (p = 0.048). Considering the different sources of PBMC and BALF, PBMC may be affected by several organs and tissues through blood circulation, and the sample size of GSE28221 is small. The evidence for identifying IPF outcomes related m6A regulators in PBMC is insufficient, and more clinical data support is needed.

To identify the optimal IPF prognosis-related m6A modulators, we used lung function data from the GSE47460 series. Lower DLCO indicates worse lung function and is associated with poor outcomes. Spearman correlation analysis was used to analyze the correlation between DLCO and m6A regulators (|R| > 0.2, p < 0.01), and we found six genes in the GSE47460 series that were weakly correlated with DLCO. KIAA1429, FMR1, PCIF1, YTHDC1, and METTL14 were positively correlated with DLCO and considered potential protective factors for prognosis. YTHCF1 negatively correlated with DLCO and was a potential risk factor for poor prognosis (Supplementary Figures S3A–E). The expression of METTL14 was positively correlated with DLCO (R = 0.23) (Figure 3F) and was identified as a critical variable for outcomes in the GSE70866 and GSE47460 series. The expression of YTHDC1 was positively correlated with DLCO, consistent with the PMBC-derived GSE28221 series.

We determined whether METTL14, G3BP2, and ZC3H13 levels correlate with other m6A genes. METTL14, G3BP2, and ZC3H13 are protective factors for prognosis, and METTL14 is positive correlations with ZC3H13 and G3BP2 (Figures 3G,H). Besides, five m6A genes were moderately correlated with the METTL14 expression level, including three readers (IGF2BP3, LRPPRC, and YTHDC2) and the two writers (KIAA1429 and WTAP) (Supplementary Figures S3F–J). G3BP2 found correlated with LRPPRC, YTHDC2, YTHDF2, YTHDF3, KIAA1429, and WTAP by Spearman correlation analysis (|R| > 0.4, p < 0.001) (Supplementary Figures S3K–P). In summary, for the three outcomes-related m6A regulators, we believe that METTL14 has the strongest regulation on outcomes, including the regulation of lung function, and positively correlates with m6A writers and readers to promote RNA methylation, followed by G3BP2, which is also positively correlated to m6A writers and readers.

We applied the “consensus clustering method” to explore the impact of different m6A patterns on prognosis based on three m6A prognostic regulators. Compared with clusters 1 to 9, the growth rate in cluster two was flat in the CDF plot (Supplementary Figure S4A). Kaplan-Meier survival analysis of IPF patients with two different m6A clusters showed that the prognosis of clusters A and B were not significant (p = 0.123) (Supplementary Figure S4B). We believe that this m6A model does not correlate well with IPF outcomes. According to these conclusions, we selected METTL14 and G3BP2 with strong prognostic correlations as variables to distinguish m6A patterns. Compared with clusters 1 to 9, the growth rate in cluster three was flat in the CDF plot (Figure 4A). Figure 4B shows that the relative change in area under the CDF curve increased insignificantly for cluster 3. In the consistency matrix of cluster 3, the intra-group correlation was higher, and the inter-group correlation was low (Figure 4C and Supplementary Figures S4D–J). In summary, we identified three m6A patterns. Kaplan-Meier survival analysis of IPF patients showed significant differences in outcomes among the three clusters (p = 0.005) (Figure 4D). Clusters B (89 cases) and C (17 cases) had poor prognosis; therefore, we believe that cluster A (70 cases) is the m6A low-risk group (containing 70 cases); clusters B and C are the m6A high-risk group (containing 106 cases) for Kaplan-Meier survival analysis. The low-risk group had significantly better outcomes than the high-risk group (p = 0.008) (Figure 4E). The PCA scatter diagram was obtained for visualization, which distinguished between m6A low-risk and high-risk (Figure 4F). Boxplots and heatmaps of m6A-related genes were drawn based on different m6A risk groups to demonstrate the differential expression levels of 31 m6A regulators (Figures 4G,H); 18 m6A regulators were significantly differentially expressed.

A clinical heatmap was drawn for the five prognostic m6A genes, gender, age (older than 65 years old), survival status, and GAP score data of patients in the m6A high- and low-risk groups (Figure 5A). The chi-square test was performed on the clinical information, and there was no significant difference in gender and age among patients between risk groups, while the survival status (p < 0.001) and GAP scores (p < 0.05) were significantly different. Using the “limma” package to analyze the DEGs of the m6A high-risk versus low-risk group, 195 DEGs were obtained (Supplementary Table S1), including 27 downregulated genes and 168 upregulated genes. Supplementary Figure S4K is a volcano plot of various m6A clusters, and Figure 5B is an enhanced volcano plot of various m6A risk groups, showing the primary DEGs and marking the genes related to cell adhesion. The KEGG pathway analysis enriched adherens junction and cell adhesion molecules (Figure 5C). YES1, ITGB8, SORBS1, CLDN3, LRRC4, NRXN3, and PARD3 were DEGs related to adherens junction and cell adhesion molecules. Three downregulated genes and four upregulated genes in the high-risk group are shown in Figure 5B. Functions related to cell adhesion were observed in the GO annotation. Regulation of cell junction assembly was the affected BP (Figure 5E). These findings suggest that cell adhesion and junction functions might be regulated by m6A regulators, mediating IPF progression. Using TRRUST prediction, SMAD4 and SMAD3 were the main transcription factors of the TGF-β/SMAD pathway and are considered the essential pathway for IPF progression (Figure 5D). Myofibril and microtubule cytoskeleton were the affected CCs associated with fibroblasts.

To explore the influence of immune cells on IPF patients with different m6A types, we applied ssGSEA to visualize the abundance of immune cells in IPF samples. We found that neutrophils in the high-risk group were significantly upregulated compared with the low-risk group. Clinical study also showed that the upregulation of neutrophils signals poor IPF outcomes (Jegal, 2022). The high-risk group was markedly enriched in activated CD56dim natural killer cells, myeloid-derived suppressor cells, and type 2 T helper cells (Figure 6A). We applied Spearman correlation analysis to assess the relationships between m6A regulators and immune cell infiltration (Figure 6B). We evaluated the correlation between three m6A regulators regulating immune cells and found that G3BP2 and METTL14 were m6A regulators with negative correlations with neutrophil expression. We also explored differential immune cell infiltration between patients with high and low G3BP2 or METTL14 expression (Supplementary Figures S5A,B). WTAP is an m6A methyltransferase, and various immune cells positively correlate with WTAP. WTAP upregulates T cells, including activated CD4 T cell, regulatory T cell, T follicular helper cell, and type 2 T helper cell (Supplementary Figures S5C).

The “pRRophetic” package is a computational model that predicts chemotherapy responses based on tumor gene expression data (Altmann et al., 2010). IPF drugs are associated with the development of anti-tumor drugs such as nintedanib, an oral small molecule tyrosine kinase inhibitor initially developed for lung cancer and approved for the treatment of IPF. Pirfenidone has a sensitizing effect on the combination of paclitaxel and carboplatin used in clinical practice (Branco et al., 2022). We believe that the identification of related anti-tumor drugs has potential significance for IPF treatment. We obtained six drugs potentially associated with m6A based on a two-sided p-value <0.01. Drugs sensitive to the m6A high-risk group included erlotinib, XAV939, WZ-1-84, 681640, and A-770041 (Figures 7A–E). KIN001-135 was more sensitive to the m6A low-risk group (Figures 7F). The m6A low-risk group had better prognosis, while the m6A high-risk group was associated with more potentially sensitive drugs.