Effect of M6A regulators on diagnosis, subtype classification, prognosis and novel therapeutic target development of idiopathic pulmonary fibrosis

Abstract

Molecular biology studies show that RNA N6-methyladenosine (m6A) modifications may take part in the incidence and development of idiopathic pulmonary fibrosis (IPF). Nonetheless, the roles of m6A regulators in IPF are not fully demonstrated. In this study, 12 significant m6A regulators were filtered out between healthy controls and IPF patients using GSE33566 dataset. Random forest algorithm was used to identify 11 candidate m6A regulators to predict the incidence of IPF. The 11 candidate m6A regulators included leucine-rich PPR motif-containing protein (LRPPRC), methyltransferase-like protein 3, FTO alpha-ketoglutarate dependent dioxygenase (FTO), methyltransferase-like 14/16, zinc finger CCCH domain-containing protein 13, protein virilizer homolog, Cbl proto-oncogene like 1, fragile X messenger ribonucleoprotein 1 and YTH domain containing 1/2. A nomogram model was constructed based on 11 candidate m6A regulators and considered beneficial to IPF patients using decision curve analysis. Consensus clustering method was used to distinctly divide IPF patients into two m6A patterns (clusterA and clusterB) based on 12 significant m6A regulators. M6A scores of all IPF patients were obtained using principal component analysis to quantify the m6A patterns. Patients in clusterB had higher m6A scores than those in clusterA. Furthermore, patients in clusterB were correlated with Th17 and Treg cell infiltration, innate immunity and Th1 immunity, while those in clusterA were correlated with adaptive immunity and Th2 immunity. Patients in clusterB also had higher expressions of mesenchymal markers and regulatory factors of fibrosis but lower expressions of epithelial markers. Lastly and interestingly, two m6A regulators, LRPPRC (p = 0.011) and FTO (p = 0.042), were identified as novel prognostic genes in IPF patients for the first time using an external GSE93606 dataset. Both of them had a positive correlation with a better prognosis and may serve as therapy targets. Thus, we conducted virtual screening to discover potential drugs targeting LRPPRC and FTO in the treatment of IPF. In conclusion, m6A regulators are crucial to the onset, development and prognosis of IPF. Our study on m6A patterns may provide clues for clinical diagnosis, prognosis and targeted therapeutic drugs development for IPF.

Keywords: consensus clustering; diagnostic model; idiopathic pulmonary fibrosis; m6A regulators; novel therapeutic targets; prognostic markers.

FIGURE 1

The study workflow.

FIGURE 2

Landscape of m6A regulators in IPF. (A) Expression heatmap of 12 significant m6A regulators in healthy controls (Type: con) and IPF patients (Type: treat). (B) Chromosomal positions of the 26 m6A regulators. (C) Differential expression histogram of 12 significant m6A regulators identified between healthy controls (Type: con) and IPF patients (Type: treat). *p < 0.05, **p < 0.01, and ***p < 0.001.

FIGURE 3

Correlation between writers and erasers in IPF (A–E) Writer genes: CBLL1, METTL14, METTL16, VIRMA and ZC3H13; eraser genes: ALKBH5 and FTO.

FIGURE 4

Construction of random forest (RF) model and support vector machine (SVM) model. (A) Reverse cumulative distribution of residual of RF and SVM model. (B) Boxplots of residual of RF and SVM model. (C) Cross-validation curve of RF model. (D) The importance of the 12 significant m6A regulators. (E) ROC curves of RF and SVM model.

FIGURE 5

Construction of the nomogram model. (A) Nomogram model constructed by the 11 candidate m6A regulators. (B) Predictive value of nomogram model through a calibration curve. (C) Decisions curve analysis of nomogram model showing benefits to IPF patients. (D) Clinical impact curve of nomogram model.

FIGURE 6

Consensus clustering of 12 significant m6A regulators in IPF. (A–D) Consensus matrices with cluster count from 2 to 5 showing an optimal cluster (clusterA and cluster B) with k = 2. (E) Expression heatmap in clusterA and clusterB. (F) Differential expression histogram in clusterA and clusterB. (G) Principal component analysis based on 12 significant m6A regulators showing a notable distinction between clusterA and clusterB. (H) Gene ontology enrichment showing potential biological functions of 402 m6A-related differentially expressed genes (DEGs) on the etiopathogenesis of IPF. *p < 0.05, **p < 0.01, and ***p < 0.001.

FIGURE 7

Single sample GSEA for immune infiltration. (A) Heatmap revealing relationship between immune cells and the 12 significant m6A regulators. (B) Distinction of immune cells infiltration between high and low LRPPRC expression subgroups. (C) Distinction of immune cells infiltration between high and low FTO expression subgroups. (D) Distinction of immune cells infiltration between clusterA and clusterB. *p < 0.05, **p < 0.01, and ***p < 0.001.

FIGURE 8

Consensus clustering of the 402 m6A-related DEGs in IPF. (A–D) Consensus matrices with cluster count from 2 to 5 showing an optimal cluster (gene clusterA and gene cluster B) with k = 2. (E) Differential expression histogram in gene clusterA and gene clusterB. (F) Distinction of immune cells infiltration between gene clusterA and gene clusterB. (G) Distinction of m6A score between clusterA and clusterB. (H) Distinction of m6A score between gene clusterA and gene clusterB. *p < 0.05, **p < 0.01, and ***p < 0.001.

FIGURE 9

Effect of m6A patterns on differentiation of IPF. (A) Sankey diagram revealing relevance between m6A patterns, m6A gene patterns, and m6A scores. (B) Different expression level of epithelial markers and mesenchymal markers between clusterA and clusterB. (C) Different expression level of epithelial markers and mesenchymal markers between gene clusterA and gene clusterB. (D) Different expression level of regulatory factors of fibrosis between clusterA and clusterB. (E) Different expression level of regulatory factors of fibrosis between gene clusterA and gene clusterB. *p < 0.05, **p < 0.01, and ***p < 0.001.

FIGURE 10

Clinical prognostic value of LRPPRC and FTO in IPF and their binding pockets for virtual screening. (A) Kaplan-Meier survival curve of the overall survival in high and low LRPPRC expressions subgroups. (B) Kaplan-Meier survival curve of the overall survival in high and low FTO expression subgroups. (C) Binding pocket of LRPPRC: core of pocket (−15.83, 4.206, 9.538) Å, size of pocket (40, 40, 40) Å. (D) Binding pocket of FTO: core of pocket (29.199, −7.339, −23.371) Å, size of pocket (26, 26, 26) Å.

FIGURE 11

Virtual Screening for Potential Drugs Targeting LRPPRC. (A) The top five docked compounds (Z109823102, Z79383944, Z18792881, Z31753778, Z16009222) from Enamine HTS potentially targeting LRPPRC. Interaction type between each compound and LRPPRC: Z109823102 (hydrophobic interaction, hydrogen bond, π-cation interaction), Z79383944 (hydrophobic interaction, hydrogen bond, π-cation interaction), Z18792881 (hydrophobic interaction,π-cation interaction), Z31753778 (hydrophobic interaction, hydrogen bond, π-π stacking), Z16009222 (hydrophobic interaction, hydrogen bond). (B) The top five docked natural products (ZINC68568380, ZINC68563949, ZINC70706523, ZINC85907291, ZINC70706097) from ZINC potentially targeting LRPPRC. Interaction type between each product and LRPPRC: ZINC68568380 (hydrophobic interaction), ZINC68563949 (hydrophobic interaction), ZINC70706523 (hydrophobic interaction, hydrogen bond), ZINC85907291 (hydrophobic interaction), ZINC70706097 (hydrophobic interaction, hydrogen bond).

FIGURE 12

Virtual Screening for Potential Drugs Targeting FTO. (A) The top five docked compounds (Z28140847, Z316147040, Z31323863, Z335602852, Z45588056) from Enamine HTS potentially targeting FTO. Interaction type between each compound and FTO: Z28140847 (hydrophobic interaction, hydrogen bond, π-π stacking), Z316147040 (hydrophobic interaction, π-π stacking), Z31323863 (hydrophobic interaction, π-cation interaction), Z335602852 (hydrophobic interaction, π-π stacking), Z31323863 (hydrophobic interaction, π-cation interaction, π-π stacking), Z45588056 (hydrophobic interaction, hydrogen bond, π-π stacking). (B) The top five docked natural products (ZINC03875800, ZINC70665164, ZINC04404594, ZINC68569433, ZINC05220992) from ZINC potentially targeting FTO. Interaction type between each product and FTO: ZINC03875800 (hydrophobic interaction, π-cation interaction, π-π stacking), ZINC70665164 (hydrophobic interaction, π-π stacking), ZINC04404594 (hydrophobic interaction, π-cation interaction), ZINC68569433 (hydrophobic interaction), ZINC05220992 (hydrophobic interaction, π-π stacking).