Pulmonary fibrosis (PF) is a chronic, progressive tissue repair response, which leading to irreversible scarring and lung remodeling (King et al., 2011). PF can occur secondary to certain predisposing factors or diseases, such as radiation (He et al., 2019), asbestos (Pira et al., 2018), silica (Cao et al., 2020), drugs (Della Latta et al., 2015), autoimmune diseases (Fischer and Distler, 2019), etc. However, some patients with PF without a clear cause, which is called idiopathic pulmonary fibrosis (IPF). IPF is a chronic, age-related interstitial lung disease (ILD) characterized by excessively deposition of extracellular matrix (ECM) protein and irreversible loss of lung function, causing progressive respiratory failure (Richeldi et al., 2017; Barratt et al., 2018). The pathogeny of IPF is still unknown, but it likely related to heredity and environment. There are large regional differences in the incidence of IPF, ranging from 0.35 to 1.30 per 100,000 individuals in Asia–Pacific countries, 0.09 to 0.49 per 100,000 individuals in Europe, and 0.75 to 0.93 per 100,000 individuals in North America (Maher et al., 2021). IPF tends to occur in men between 40 and 50 years of age and has a poor prognosis. The average life expectancy of untreated IPF patients is only 3–5 years, and most patients die of acute exacerbations of IPF or respiratory failure. Actually, acute exacerbations of IPF can occur at any time during the course of the disease and are associated with extremely high mortality (Spagnolo and Wuyts, 2017). Although two antifibrotic drugs, nintedanib and pirfenidone, have been shown to delay the progression of IPF, there is currently no drug that can cure IPF (Raghu et al., 2015).

Epigenetics usually refers to the heritable modification of genetic material without changing gene sequence, including DNA methylation, RNA methylation, histone modification, chromosome remodeling, etc., which plays an important role in various diseases and tumors (Berger et al., 2009). At present, more than 100 kinds of RNA (mRNA, lncRNA, snRNA, etc.) have been found post-transcriptional modifications, among which N6-methyladenosine (m6A) is the most common (Yue et al., 2015; Boccaletto et al., 2018). M6A RNA modification is a dynamic and reversible post-transcriptional modification process mediated by m6A WER proteins (methyltransferase “writers”, demethylase “erasers”, binding proteins “readers”), which plays a crucial regulatory role in RNA metabolism, splicing, translation and other processes (Wang et al., 2020). Previous studies have shown that m6A is widely involved in the development of various diseases, such as pneumonia, lung cancer, colorectal cancer, breast cancer, nasopharyngeal cancer, systemic lupus erythematosus, etc. (Li et al., 2018; Chang et al., 2020; Yue et al., 2020; Maher et al., 2021; Meng et al., 2021; Li et al., 2022). For example, Li et al. (2021) found that SNHG4 promoted LPS-induced inflammation by inhibiting METTL3-mediated m6A level of STAT2 mRNA. And research pointed out that overexpressed FTO enhanced the expression of MZF1 by reducing the m6A modification level and stability of MZF1 mRNA, thereby promoting the development of lung cancer (Liu et al., 2018). Similarly, enhanced activity of methyltransferase METTL3 increased the m6A modification level of JUNB mRNA and accelerated the progression of TGF-β-induced lung adenocarcinoma (LUAD) (Wanna-Udom et al., 2020). These studies indicated that RNA methylation regulators could affect the development of the above diseases by regulating the m6A modification of RNA. M6A-related genes can also be used as diagnostic and prognostic markers for lung diseases. For example, studies found that m6A-related genes (EGFR, RFXAP, KHDRBS2, ADAMTS6, etc.) were determined to be associated with overall survival (OS) in patients with LUAD, in which RFXAP and KHDRBS2 exhibited independent prognostic value (Sun et al., 2021). Additionally, Jia et al. (2022) showed that three m6A-related genes (FAM71F1, MT1E, and MYEOV) were identified as prognostic genes in Lung Squamous Carcinoma (LUSC). However, there are few reports on m6A methylation modification in the occurrence and development of IPF. Therefore, it is of great significance to explore m6A-related genes and construct IPF-related prognostic risk model to assist in judging the progression and prognosis of IPF.

Weighted gene co-expression network analysis (WGCNA) is a comprehensive analysis technique based on biological network, which can identify a class of genes (or proteins) that are co-expressed, and cluster genes with similar expression patterns through algorithms into different modules, analyze the association between modules and characteristic traits or phenotypes, use clustering modules to associate with phenotypes to build a co-expression network, and explore the core genes (or proteins) in the modules, so as to provide ideas for exploring the molecular mechanism of diseases (Presson et al., 2008; Yin et al., 2018). Compared with microarray and high-throughput sequencing analysis, WGCNA is suitable for multiple statistical tests to analyze the correlation between genes and avoid losing the trend information of genes according to a fixed threshold screening.

The Cox proportional hazards model is essentially a regression model commonly used in medical research statistics to study the association between a patient’s survival time and one or more predictor variables (Cox, 1972). It is applicable to quantitative predictor variables and categorical variables. It mainly includes univariate and multivariate Cox regression analysis. Univariate Cox analysis is usually used to remove collinearity, but may lead to synergistic effects caused by other variables, so multivariate Cox regression is performed to correct other factors, which is often used for data modeling in survival analysis (Huang and Liu, 2006; Li et al., 2016).

In this paper, the microarray data GSE28042 was downloaded from the Gene Expression Omnibus (GEO) database, and the gene expression profiles of peripheral blood mononuclear cell (PBMC) and the corresponding clinical data of 75 IPF samples and 19 normal samples were obtained. Through the analysis of differentially expressed genes, WGCNA and m6A correlation analysis method, a group of m6A-related candidate genes with high IPF correlation and differential expression were screened. The genes associated with the prognosis of IPF were screened out by univariate Cox regression analysis, LASSO analysis, and multivariate Cox regression analysis, and the multivariate Cox model of prognostic risk of related genes was constructed. Finally, the prognostic predictive ability and independent prognostic characteristics of the risk model were evaluated by survival analysis and independent prognostic analysis, and verified by GSE93606 dataset, which is intended to provide a basis for prognostic prediction of IPF patients (Figure 1).

The etiology of IPF is still not fully understood, but some studies have shown that its pathogenesis may be related to the abnormal damage and repair of alveolar epithelial cells, epithelial-to-mesenchymal transition (EMT), fibroblast-to-myofibroblast transformation (FMT), and inflammatory response (King et al., 2011). Worldwide, the incidence and mortality of IPF are on the rise. Lung transplantation is the only treatment for IPF that can prolong life expectancy (Kumar et al., 2018). Unfortunately, IPF patients without lung transplantation have a short median survival time. M6A is the most abundant post-transcriptional modification in mRNA and long non-coding RNA (lncRNA) in most eukaryotes. In addition, studies have reported that m6A is involved in post-transcriptional modification, cell differentiation, cell recoding, cell stress and other processes, and plays an important role in lung diseases such as lung cancer, pulmonary hypertension and chronic obstructive pulmonary disease through various mechanisms. However, there are few studies on m6A in IPF. Therefore, it is necessary to explore the prognostic value of m6A-related genes in IPF and establish a set of prediction models for evaluating the survival time of IPF and improving the prognosis of patients.

In this study, we downloaded GSE28042 dataset from GEO database, which included peripheral blood monocyte cell gene expression profiles and their corresponding clinical information of 75 IPF samples and 19 normal samples, and analyzed the obtained data. The gene expression matrix was used for differential gene analysis, and 606 up-regulated genes and 686 down-regulated genes were screened. The correlation between each module and the trait was obtained by WGCNA analysis combined with correlation heat map. The black and pink modules with the highest positive and negative correlations were selected, and 405 intersection genes were obtained by intersection of the DEG and the module genes with the highest correlation in the selected WGCNA. Then, 218 m6A-related candidate genes were screened out from the 405 IPF highly correlated differentially expressed genes by Pearson correlation analysis, and the enrichment analysis of these genes showed that the above genes were mainly enriched in herpes simplex virus type Ⅰ(HSV-1) infection, neutrophil degranulation, ciliary assembly and other pathways. Studies have shown that chronic viral infections, mainly herpes virus infections, may contribute to the development of IPF. And HSV-1 is a kind of herpes virus, it can enter the alveoli through the respiratory tract and spread with the blood, resulting in focal necrotizing pneumonia, followed by diffuse pulmonary fibrosis (Luyt, 2020). Neutrophil degranulation is one of the important links that neutrophils participate in the inflammatory response. As inflammatory cells, neutrophils participate in the progression of PF by promoting the proliferation of fibroblasts and enhancing the differentiation of myofibroblasts (Gregory et al., 2015; Klopf et al., 2021). Cilia is an organelle protruding from the cell surface. The abnormal structure and function of cilia can cause various diseases, such as bronchiectasis and reproductive infertility (Jain et al., 2012; Girardet et al., 2019). Moreover, studies have shown that pulmonary fibrosis is associated with bronchiectasis (Fitzgerald et al., 2017). The above relevant findings suggest that the m6A-related candidate genes screened were closely related to the occurrence and development of PF. Therefore, we hypothesized that the m6A-related candidate genes were associated with IPF.

In order to explore the role of m6A-related candidate genes in the prognosis of IPF, we screened out 30 genes associated with patient prognosis by univariate Cox analysis, and then screened out 5 key genes (RBM11, RBM47, RIC3, TRAF5, ZNF14) by LASSO analysis and multivariate Cox analysis. The above studies indicate that the five key genes and 23 m6A regulators are significantly correlated and modified by their regulation. This regulation can be direct or indirect, but its specific mechanism is still unknown. The results of protein-protein interaction analysis also showed that RBM11 and RBM47 had protein-protein interactions with m6A regulators, and the m6A-Atlas analysis showed that all five key genes had m6A sites (Tang et al., 2021), which added confidence to our results. We construct a riskscore model as an indicator to predict the prognosis of IPF [riskscore = (−0.44084*RBM11) + (0.631579*RBM47) + (−0.01935*RIC3) + (−0.58291*TRAF5) + (−0.00528*ZNF14)], and then survival analysis was performed to assess the effect of the above genes on the prognosis of IPF patients. The results of single-gene survival analysis showed that high expression of RBM11, RIC3, TRAF5, ZNF14 was associated with good prognosis of IPF, while high expression of RBM47 was associated with poor prognosis; overall survival analysis of the risk prognostic model showed that high-risk patients had poor survival, while low-risk patients had higher survival, which preliminarily indicated the correctness of the model. Simultaneous univariate and multivariate independent prognostic analyses indicated that the risk model for these five key genes was independent of other clinicopathological parameters (gender, age). TRAF5 is an important signal transducer for a wide range of TNF receptor superfamily members, and it mainly mediates the activation of NF-κB pathway (Au and Yeh, 2007). Indeed, study has shown that overactivation of NF-κB pathway is associated with apoptosis of alveolar epithelial type II cells (AEC2) and the development of PF (Yang et al., 2018). Besides, Ben-David et al. (2016) demonstrated that inflammatory signals regulate the expression and splicing of RIC3, thereby influencing the α7 nA-ChR mediated cholinergic anti-inflammatory pathway. Although the role of inflammation in fibrosis is controversial, it is still considered to be an important component of IPF. Recently, Kim et al. (2019) pointed out that RBM47 promotes the EMT of cells by promoting TJP1-mediated alternative splicing. Globally, EMT is considered to be one of the key mechanisms of PF. When tissues are subjected to various injuries, a series of immune signals are generated, leading to inflammation and promoting EMT. In this process, macrophages, neutrophils and other immune cells are recruited and release proinflammatory cytokines to maintain inflammation and pulmonary fibrosis (Salton et al., 2019). In conclusion, we speculate that the above three genes are closely related to the progression of pulmonary fibrosis. However, studies on RBM11 and ZNF14 in lung diseases are rare.

These results indicated that the key genes screened by bioinformatics methods were highly correlated with the occurrence and development of IPF, and had a significant correlation with the prognosis of IPF patients. Therefore, the above five key genes can provide reference for the diagnosis and treatment of IPF. We also analyzed the risk model. By drawing the time-dependent ROC curve of the prognostic model, we found that the AUC value gradually increased with the increase of time, indicating that the accuracy of our prognostic model was good. Finally, the prognostic model was verified by the GSE93606 dataset. It can be seen that the prognostic model is also applicable to this dataset, which further confirms that the prognostic risk model has a strong and stable prediction efficiency.

However, the study also has certain limitations. First, our results are based on data from existing public databases. Therefore, a large-scale, prospective, multicenter study is needed to further validate our results. Secondly, our study population is mainly from European and American populations. Therefore, our findings may not be optimal for patients from other countries and ethnicities. Finally, the correlation between some key genes and the development and progression of IPF has not been confirmed by biological experiments. In follow-up studies, experimental validation will be performed to reveal the relationship between key genes and IPF. In this way, we can determine their suitability as new diagnostic and therapeutic targets to provide a rationale for the clinical diagnosis and treatment of IPF.