Three IPF-related cohorts and the corresponding clinical data were collected from the Gene Expression Omnibus database, containing the Freiburg, Siena, and Leuven cohorts, all of which belong to the series GSE70866. Given that data of the Freiburg cohort and Siena cohort were detected with the same platform GPL14550, they were combined as the Freiburg-Siena cohort for further analysis. The gene expression levels of Freiburg-Siena and Leuven cohorts were normalized with log₂ (x + 1). The two cohorts included multiple bronchoalveolar lavage samples obtained from IPF patients (n = 176) and normal healthy individuals (n = 20) (Figure 1). Considering the apparent difference in the number of samples between IPF (n of Freiburg-Siena = 112) and healthy (n of Freiburg-Siena = 20) groups, the SMOTE (synthetic minority oversampling technique) algorithm was used to balance the sample size. This process was completed by the “DMwR” package for R (v4.1.0), and the target samples of both the IPF group and the healthy group (for the new Freiburg-Siena cohort) were set to 120. The new data produced based on SMOTE were used to explore the expression difference and identification effect of m7GRGs between the IPF and healthy groups.

Three m7GRGs, BUD23, METTL1, and RNMT, were collected from the literature survey (Ouyang et al., 2021). METTL1 and RNMT were included in both the Freiburg-Siena and Leuven cohorts and thus were further analyzed in this study. Based on the Freiburg-Siena cohort (the samples were more than the Leuven cohort), the average method was used to identify various m7G subtypes related to m7GRGs, with clustering by the k-means algorithm and describing the distance metric by euclidean. The k values were set as 2–20, and the optimal k value was determined by the cumulative distribution function plot. The identification of m7G subtypes was performed by utilizing the “ConsensuClusterPlus” package (Wilkerson and Hayes, 2010) in R (v4.1.0).

The gene expression data of the Freiburg-Siena cohort was used to select differentially expressed genes (DEGs) in various m7G subtypes by the threshold of the absolute value of log₂ (fold change) > 1 and adjusted p-value < 0.05, and the step was finished with the help of “limma” package (Ritchie et al., 2015). Based on data from both the Freiburg-Siena and Leuven datasets, prognosis-related genes (PRGs) were chosen by the p values of univariate Cox regression analysis less than 0.05. Candidate genes for constructing the m7GPS were identified from the intersection of DEGs and PRGs. With the “sample” algorithm in R (v4.1.0), IPF patients from the Freiburg-Siena cohort were randomly classified into the discovery cohort (70%) and internal validation cohort (30%). The Leuven cohort was identified as the external validation cohort.

Multiple methods have been used in medical prognostic signatures, such as the least absolute shrinkage and selection operator (LASSO) Cox regression and Random Forest. Considering that LASSO Cox regression was more common in this field (Supplementary Figure S1), it was used to construct the m7GPS in our study. The expression levels of candidate genes, survival status, and survival time data of the discovery cohort were applied to establish the m7GPS (consisting of specific genes) using the LASSO Cox regression. In addition to two parameters (family = “cox” and maxit = 1,000), all parameters of the “glmnet” function (for fitting the signature) and the “cv.glmnet” function (for validating the signature by using k-fold cross) were set by default (Friedman et al., 2010; Simon et al., 2011).

Using time-dependent receiver operating characteristic (ROC) curves (Robin et al., 2011) produced by the “timeROC” package, the accuracy of m7GPS in predicting the prognosis was assessed through the area under the curve (AUC) values of the discovery, internal validation, and external validation cohorts. A high AUC value (ranging from 0 to 1) represented the conspicuous accuracy of the m7GPS in predicting the survival of IPF patients. Moreover, confusion matrixes (calculated by the “caret” package) and F1 score (ranging from 0 to 1) were also utilized to assess the precision of the m7GPS. For the F1 score, the increasing value indicated a more conspicuous classification performance of m7GPS. The relationship of m7GPS with the prognosis of IPF individuals was evaluated by both the Kaplan-Meier curves (with log-rank test) and Cox regression analyses. For Cox regression analyses, an indicator (e.g., a gene) with a hazard ratio (HR) of more than 1 indicated a poor prognosis for IPF patients. AUC values were also applied to explore the ability of m7GPS and its constitutive genes in screening IPF patients.

A nomogram was constructed based on prognosis-related factors identified from univariate Cox analysis, where the survival probabilities of IPF patients can be predicted. Calibration curves were utilized to detect the accuracy of the nomogram in predicting the survival condition of patients. Decision curve analysis (DCA) (Vickers and Elkin, 2006) was applied to evaluate the net benefits for IPF patients in predicting survival rates based on the nomogram.

Differential expression genes (DENGs) between the high-risk group and low-risk group were identified for exploring the potential molecular functions and KEGG (Kyoto Encyclopedia of Genes and Genomes) (Kanehisa et al., 2017) signaling pathways of m7GPS in IPF. The calculation processes were carried out by the “limma,” “clusterProfiler,” and “GOplot” packages (Yu et al., 2012; Walter et al., 2015). The CIBERSORT (Newman et al., 2019) algorithm was utilized to predict the composition of 22 types of immune cells for the discovery, internal validation, and external validation cohorts. Wilcoxon rank-sum tests were exerted to disclose the difference in infiltration levels of immune cells between the high-risk group and low-risk group.

The weighted gene co-expression network analysis (WGCNA) determined candidate hub genes of m7GPS. In detail, the median absolute deviation (MAD) of all genes in the discovery cohort was calculated, and genes with a MAD value of more than one were screened for performing WGCNA. The steps of WGCNA were as follows: 1) the top 5,000 genes with the highest MAD values were divided into various modules related to the m7GPS risk score using the “WGCNA” package (Langfelder and Horvath, 2008); 2) The scale-free topology fitting index R² was set as 0.85 and identified the minimum soft threshold; 3) the correlations between genes in each module and m7GPS risk score were evaluated by gene significance and module membership; 4) the genes of the essential module were output for selecting hub genes. Subsequently, hub genes of m7GPS were identified in Cytoscape (v3.9.0) based on the “degree” algorithm. The same steps were also utilized for internal and external validation cohorts.

The “oncoPredict” package (Maeser et al., 2021) was applied to predict potential drugs (from the GDSC2 database) sensitive to IPF patients with high risk scores (more than the median value of all risk scores in the corresponding cohort). The drug sensitivity was assessed by half-limiting dose values.

The expression difference of m7GRGs between IPF and control groups was identified by Wilcoxon rank-sum test. The test was also utilized to evaluate clinical characteristics—ages, genders, GAP (gender-age-physiologic) stages, and subtypes—between high and low risk score IPF groups. All the analyses in this study were carried out in R (v4.1.0) and Cytoscape (v3.9.0). The design of this study is shown in Figure 1.

According to the new data based on SMOTE (n = 240), the expression level of METTL1 in the IPF group was higher than that in the normal group (p < 0.05); in contrast, the expression of RNMT between the two groups was not statistically different (p = 0.42) (Figure 2A). IPF patients with high METTL1 expression had an unfavorable prognosis—shorter survival (p < 0.05, Figure 2B). Based on the Freiburg-Siena cohort, despite no statistical significance, patients with upregulated RNMT expression tended to have an optimistic prognosis—more prolonged survival (p = 0.086, Figure 2B). Based on the expression of the two m7GRGs—METTL1 and RNMT, IPF patients in the Freiburg-Siena cohort were classified into two subtypes (the optimal k = 2) by the unsupervised clustering algorithm (Figures 2C,D). Subtype 1 and subtype 2 contained 42 and 70 IPF patients, respectively (Supplementary Table S1). The analysis of the Kaplan-Meier curve showed that the clinical prognosis of these two m7G subtypes was different, and IPF patients in subtype 2 had more unfavorable survival (p < 0.05, Figure 2E).

Based on the Freiburg-Siena cohort, there were 1,579 DEGs between the various subtypes (subtype 2 versus subtype 1), including 1,490 upregulated expression genes and 89 downregulated expression genes (Figure 3A). Using univariate Cox regression analyses, 3,995 and 1,563 PRGs were identified from the Freiburg-Siena and Leuven queues, respectively. Two hundred and sixty-three m7GPS candidate genes were obtained by the intersection of DEGs and PRGs (Figure 3B), and they would be utilized to establish the m7GPS subsequently. Seventy-nine (70%) patients in the Freiburg-Siena cohort were randomly assigned to the discovery cohort, and the remaining 33 patients (30%) were assigned to the internal validation cohort (Supplementary Tables S2, S3). The ratio of survival status (alive or dead) between the two cohorts was balanced, and such a condition can also be seen for GAP stages (p > 0.05; Supplementary Table S4), suggesting the random grouping of the discovery and internal validation cohorts was unbiased for IPF severity. Another queue, Leuven, was defined as the external validation cohort (Supplementary Table S5). Although there was a statistically significant difference in GAP stages between the discovery and external validation cohorts, the distribution of both 1-year and 2-year survival/death ratios was balanced (Supplementary Table S4). Moreover, considering that the discovery cohort and the external validation cohort came from different countries and cities (Prasse et al., 2019) and that it is meaningful to verify the effectiveness of the m7GPS by cohorts of different countries and cities, the external validation cohort was conserved in this study.

Based on the discovery cohort and the LASSO Cox regression algorithm, the m7GPS containing six candidate genes—CCL2, CCL7, HS3ST1, MRVI1, TM4SF1, and TPST1—was constructed (Figures 3C,D). The coefficients of the six genes for calculating m7GPS risk scores can be seen in Figure 3E. The expression level of the six genes in subtype 2 was higher than that in subtype 1 (Figure 3F), and the expression level of the six genes in subtypes 1 and 2 were significantly positively correlated, respectively (Figure 3G), suggesting the rationality of m7GPS subtypes to some extent. Moreover, the positive or negative relationship between m7GPS and the two m7GRGs (METTL1 and RNMT) can be seen in Supplementary Figure S2, implying the association of m7GPS with m7G.

By analyzing the discovery cohort, m7GPS can highly accurately predict the 1–5 years survival status of IPF patients (AUC = 0.852–0.971; Figure 4A). This conclusion was confirmed by the internal and external validation cohorts (AUC = 0.730–0.864; Figure 4A). Notably, the sample number of the external validation cohort was not enough to calculate the 5-year AUC. The confusion matrixes verified the reliability of the m7GPS, and F1 scores (0.611–0.839) also supported the conspicuous precision of the m7GPS in predicting the survival status of IPF patients (Supplementary Table S6). In addition, m7GPS can predict the 1-year survival status of IPF more accurately than almost all single genes of m7GPS (Figure 4B).

The discovery, internal validation, and external validation cohorts consistently indicated that IPF deaths were more common in patients with m7GPS high risk scores; IPF individuals with a high m7GPS risk score had elevated expression levels of six genes (Figure 4C). Compared with patients with low risk scores, the data of discovery, internal validation, and external validation cohorts supported that IPF patients with high risk scores had a shorter survival time (Figure 4D). Furthermore, the univariate Cox regression analysis of the three cohorts showed that GAP, risk score, and the six genes (CCL2, CCL7, HS3ST1, MRVI1, TM4SF1, and TPST1) were related to the poor prognosis of IPF individuals (HR > 1, p < 0.05; Figure 5A). Multivariate Cox regression analyses of GAP and risk score demonstrated that the m7GPS risk score was an independent risk factor for the prognosis of IPF patients (HR > 1, p < 0.05; Figure 5B). Additionally, the concordance indexes of GAP stages combined with m7GPS were 0.76–0.80 for the discovery, internal validation, and external validation cohorts (Figure 5B), verifying the accuracy of m7GPS in predicting the prognosis of IPF patients. Although Random Forest was not utilized to establish the m7GPS, the concordance index (0.74) of the prognostic signature based on this method (using the “randomForestSRC” package) was not conspicuously less than that of the m7GPS based on LASSO Cox regression (Supplementary Figures S1B,C); this verified the feasibility of constructing prognostic signature from the perspective of m7G.

Further, via ROC curves of the 240 samples based on SMOTE method, the constitutive genes (CCL2, CCL7, HS3ST1, MRVI1, TM4SF1, and TPST1) and risk score of m7GPS demonstrated conspicuous accuracy in distinguishing IPF patients and healthy individuals (AUC > 0.843, Figure 6). Compared to its constitutive genes, m7GPS had the highest accuracy in differentiating IPF patients (AUC = 0.960; Figure 6), presenting as an excellent marker in screening IPF based on bronchoalveolar lavage.

Although there was no statistical difference in ages and genders between m7G subtypes, GAP stages were statistically different among various m7G subtypes, implying the underlying clinical value (similar to GAP stages) of the m7GPS (based on m7G subtypes) (Figure 7A). To explore the potential clinical application value of m7GPS, this study constructed a nomogram based on GAP and risk score (both of which were related to the prognosis of patients in univariate Cox analysis) (Figure 7B). The 1-, 3-, and 5-year survival probabilities of IPF patients can be predicted according to the nomogram. For example, for the IPF patient (ID: GSM1820740) with GAP scores equaling to 5 and m7GPS risk scores equaling to 2.605, the predicted probability of his survival being less than 5 years, 3 years, and 1 year was 0.865, 0.613, and 0.116, respectively (Figure 7B). According to the calibration curve, the prediction accuracy of the nomogram was evaluated. The 1-year survival probability predicted by the nomogram was similar to the actual survival probability of IPF patients, although the probability of long-term survival (3- and 5-years) predicted by the nomogram was relatively inaccurate (Figure 7C). DCA revealed that compared with the non-implementation strategy and all strategies, the nomogram based on GAP and m7GPS provided net benefits for IPF patients in predicting 1-, 3-, and 5-year survival rates (Figure 7D).

The DENGs between the high-risk and low-risk groups (classified based on the median risk scores in the cohort) were determined by discovery cohort data, and they were used to perform the molecular functions and KEGG analyses. The same calculations were performed for the internal and external validation cohorts. The results showed that 13 functions were screened simultaneously from three cohorts (Figure 8A), such as receptor ligand activity, signaling receptor activator activity, cytokine activity, and chemokine activity (Figure 8B). The results of the KEGG signaling pathway also suggested similar results (Figure 8C). The signaling pathway, viral protein interaction with cytokine and cytokine receptor, was found in all discovery, internal validation, and external validation cohorts (Figure 8D), suggesting that the potential molecular mechanism of m7GPS may be related to this pathway.

The data analysis results of the three cohorts showed that IPF patients with high m7GPS risk scores had increasing levels of activated mast cells (Figure 8E). In addition, increasing NK cells and decreasing M0 macrophages cells and dendritic cells were detected in patients with high risk scores (Figure 8E).

For further understating the underlying mechanisms of m7GPS in IPF, the WGCNA determined candidate hub genes of m7GPS. For the discovery cohort, five samples with distinct outliers were excluded, and the remaining 74 samples were utilized for subsequent analysis (Figure 9A). With the optimal scale-free R² (equaling 0.86), the soft threshold was set as 12 to construct a scale-free network (Figures 9B,C). Eight modules were determined in hierarchical clustering (Figure 9D). The “brown” module had the strongest correlation with the m7GPS risk score (Pearson correlation coefficient = 0.86, p < 0.05; Figure 9E), which was also supported by the significant positive correlation between gene significance of risk score and module membership in the “brown” module (Pearson correlation coefficient = 0.62, p < 0.05; Figure 9F). CYR61, C1ORF116, MAL2, SFTA2, and SUSD2 were determined as the hub genes of the “brown” module in the Cytoscape (Figure 9G), indicating the five genes were essential for the m7GPS risk score in IPF based on the discovery cohort. Additionally, CYR61 was identified as the hub gene of m7GPS risk score in the internal cohort (Supplementary Figure S3), while hub genes in the external cohort were PPP1R14A and C11ORF9 (Supplementary Figure S4).

Using the “oncoPredict” package, this study predicted drugs potentially applicable to IPF patients with high risk scores. In this research, seven drugs that may be related to m7GPS were obtained from the discovery and validation cohorts based on the bilateral p-value < 0.01. Furthermore, the unilateral p values of these seven drugs were calculated to select drugs sensitive to high-risk IPF patients. As a result, in both the discovery and validation cohorts, ZM447439-1050, AZD1332-1463, and Ribociclib-1632 were identified as drugs that might apply to patients with high-risk scores (Figures 10A,B).