m6A Regulator-Mediated Methylation Modification Patterns and Characteristics in COVID-19 Patients

Abstract

Background: RNA N6-methyladenosine (m6A) regulators may be necessary for diverse viral infectious diseases, and serve pivotal roles in various physiological functions. However, the potential roles of m6A regulators in coronavirus disease 2019 (COVID-19) remain unclear.

Methods: The gene expression profile of patients with or without COVID-19 was acquired from Gene Expression Omnibus (GEO) database, and bioinformatics analysis of differentially expressed genes was conducted. Random forest modal and nomogram were established to predict the occurrence of COVID-19. Afterward, the consensus clustering method was utilized to establish two different m6A subtypes, and associations between subtypes and immunity were explored.

Results: Based on the transcriptional data from GSE157103, we observed that the m6A modification level was markedly enriched in the COVID-19 patients than those in the non-COVID-19 patients. And 18 essential m6A regulators were identified with differential analysis between patients with or without COVID-19. The random forest model was utilized to determine 8 optimal m6A regulators for predicting the emergence of COVID-19. We then established a nomogram based on these regulators, and its predictive reliability was validated by decision curve analysis. The consensus clustering algorithm was conducted to categorize COVID-19 patients into two m6A subtypes from the identified m6A regulators. The patients in cluster A were correlated with activated T-cell functions and may have a superior prognosis.

Conclusions: Collectively, m6A regulators may be involved in the prevalence of COVID-19 patients. Our exploration of m6A subtypes may benefit the development of subsequent treatment modalities for COVID-19.

Keywords: COVID-19; consensus clustering; diagnostic biomarkers; m6A methylation modification; m6A regulators.

Figure 1

Landscape of the 26 m6A regulators in COVID-19. (A) Differential expression analysis of the 26 m6A regulators identified between samples with different COVID-19 status. (B) Expression heat map of the 26 m6A regulators in samples. (C) GSVA enrichment analysis between Non-COVID-19 and Non-ICU-COVID-19 samples. (D) GSVA enrichment analysis between Non-COVID-19 and ICU-COVID-19 samples. (E) GSVA enrichment analysis between Non-ICU-COVID-19 and ICU-COVID-19 samples. (F) The PPI network analysis among the differentially expressed genes. (G) Chromosomal positions of the 26 m6A regulators. *p < 0.05, **p < 0.01, and ***p < 0.001.

Figure 2

Correlation between m6A regulators in COVID-19. (A) Correlation plot of 26 m6A regulators. (B–J) Correlation between writers and erasers in COVID-19. Writer genes: METTL3, METTL14, METTL16, RBM15B, VIRMA, CBLL1, and ZC3H1; eraser genes: ALKBH5 and FTO. *p < 0.05, **p < 0.01, and ***p < 0.001.

Figure 3

Establishment of RF model and SVM model. (A) Reverse cumulative distribution of residual was displayed to demonstrate the residual distribution of RF and SVM model. (B) Boxplots of residual was displayed to demonstrate the residual distribution of RF and SVM model. (C) The influence of the number of decision trees on the error rate. (D) The importance of the 26 m6A regulators based on the RF model. (E) ROC curves revealed the accuracy of the RF and SVM model.

Figure 4

Establishment of the nomogram model. (A) Establishment of the nomogram model based on the 8 selected m6A regulators. (B) Predictive robustness of the nomogram model as disclosed by the calibration curve. (C) Decisions based on the nomogram model may benefit COVID-19 patients. (D) Clinical impact of the nomogram model as evaluated by the clinical impact curve.

Figure 5

Consensus clustering of the 18 significant m6A regulators in COVID-19. (A) Consensus matrices of the 18 significant m6A regulators for k = 2. (B) Differential expression analysis of the 18 significant m6A regulators in cluster A and cluster B. (C) Expression heatmap of the 18 significant m6A regulators in cluster A and cluster B. (D) PCA for the expression data of the 18 significant m6A regulators that indicates an obvious difference in transcriptomes between the two m6A subtypes. (E) GO analysis that investigates the potential mechanism underlying the effect of the 139 m6A-related DEGs on the occurrence and development of COVID-19. *p < 0.05, **p < 0.01, and ***p < 0.001.

Figure 6

Single sample gene set enrichment analysis. (A) Correlation between infiltrating immune cells and the 18 significant m6A regulators. (B) Difference in the abundance of infiltrating immune cells between high and low METTL3 expression groups. (C) Differential immune cell infiltration between cluster A and cluster B. *p < 0.05, **p < 0.01, and ***p < 0.001.

Figure 7

Consensus clustering of the 139 m6A-related DEGs in COVID-19. (A) Consensus matrices of the 139 m6A-related DEGs for k = 2. (B) Expression heat map of the 139 m6A-related DEGs in gene cluster A and gene cluster B. (C) Differential expression of the 18 significant m6A regulators in gene cluster A and gene cluster B. (D) Differential immune cell infiltration between gene cluster A and gene cluster B. (E) Differences in m6A score between cluster A and cluster B. (F) Differences in m6A score between gene cluster A and gene cluster B. *p < 0.05, **p < 0.01, and ***p < 0.001.

Figure 8

Role of m6A subtypes in distinguishing COVID-19. (A) Sankey diagram demonstrating the relationship between m6A subtypes, m6A gene subtypes, and m6A scores. (B) Differential expression levels of cytokines between cluster A and cluster B. (C) Differential expression levels of cytokines between gene cluster A and gene cluster B. *p < 0.05, **p < 0.01, and ***p < 0.001.