The research strategy is presented in Figure 1 . We conducted a search of the GEO database (http://www.ncbi.nlm.nih.gov/geo/) (24) for expression microarrays related to sepsis. We included datasets from clinical studies of sepsis in adults using peripheral blood samples. Retrieved gene expression profile data for sepsis patients (GSE57065, GSE65682 and GSE95233) using the R package GEOquery, merged the gene expression groups, and split them into sepsis and healthy control groups. Details of the datasets are shown in Table S1 .

We first compiled a list of 40 RNA methylation regulators in 5 categories (including m5C, m6Am, m1A, m7G and Ψ) from published research (25–27), which included writers, readers, and erasers, showed in Table S2 . We then obtained a total of 34 RMGs by screening gene expression with detectable expression in the GEO datasets.

All datasets were downloaded and outputs from mRNA array were normal-exponential background corrected and quantile normalized between arrays using the limma R package (28). Expression was normalized using a weighted linear regression, and precision weights were multiplied with corresponding log2 values to yield final gene expression values. Genes with an absolute log fold change (log₂FC) > 1 and adjusted p-value < 0.05 were considered up-regulated, and those with log₂FC < - 1 and adjusted p-value < 0.05 were considered down-regulated. We also performed Spearman correlation analysis to investigate the relationship between each gene. Functional interactions of the RMGs were explored using the Search Tool for Retrieval of Interacting Genes (STRING 11.0, http://string-db.org/cgi/input.pl), and a Protein-protein interaction (PPI) network was constructed. We visualized the chromosomal localization of RMGs using the R circos package (29).

The consensus clustering algorithm works by performing clustering on the same dataset multiple times, using various clustering methods or different parameter settings. The next step is to create a consensus matrix or a consensus function that measures the degree of agreement among the individual clustering solutions. This consensus matrix represents the overall consensus for each sample across different clustering methods, and is used to identify robust and stable clusters. We used consensus k means clustering to identify highly heterogeneous RNA methylation subtypes in sepsis based on RNA methylation regulator expression profiles. Clustering was performed using 100 iterations, with each iteration containing 80% of samples. We determined the optimal number of clusters using CDF curves of the consensus score, clear separation of the consensus matrix heatmaps, characteristics of the consensus cumulative distribution function plots, and adequate pairwise-consensus values between cluster members. To assess the expression distribution of RNA methylation regulators in different subtypes, we utilized moderated t-tests.

To better understand the situation of infiltrating immune cells in C1-C3 subtypes, we used the CIBERSORT algorithm. CIBERSORT is a bioinformatics tool used to deconvolute relative cell type proportions and gene expression profiles from bulk RNA sequencing datasets. The Gene Set Variation Analysis (GSVA) R software package was downloaded from http://www.bioconductor.org, and the pathway with differential enrichment in the three groups was obtained. GSVA is used to evaluate the activity or enrichment of biological pathways, gene sets, or functional signatures in individual samples from high-throughput gene expression data, such as RNA-Seq. The analysis quantifies the degree to which a specific gene set is upregulated or downregulated in each individual sample, generating a continuous enrichment score for each sample-gene set pair.

To find gene-sets significantly correlated to RNA methylation subtypes, we used weighted correlation network analysis (WGCNA) with the WGCNA R package (30). WGCNA is a powerful bioinformatics method used to analyze gene expression data and identify co-expression patterns among genes. Genes from modules highly associated with RNA methylation subtypes (the maximum correlation coefficient and P < 0.05) were selected for further Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis, to investigate the biological functions and signaling pathways involved in sepsis.

We created Random forest model (RF), support vector machine model (SVM), eXtreme Gradient Boosting (XGB) and generalized linear model (GLM) using the integrated dataset based on RMGs. The RF, GLM, SVM and XGB models all belong to machine learning model which is used to assess the importance of variables. Here, we applied them in screening of diagnostic genes for sepsis. The different classes of sepsis were utilized as the response variable, and the RMGs were used as explanatory variables. We then used the explain feature of the “DALEX” package in R to analyze the aforementioned four models and residual distribution, and the receiver operating characteristics (ROC) was plotted to select the best model using the integrated dataset. Finally, to evaluate the diagnostic efficacy of the five most important explanatory variables in the best model, we calculated ROC curve using the ‘survivalROC’ package in the three databases respectively. The AUC was used to evaluate the diagnostic genes: 0.5-0.7 (moderate), 0.7-0.8 (good), and >0.9 (excellent).

Nomogram model is a visual aid that combines multiple variables and assigns them a weighted score to estimate the likelihood of a specific outcome. Using the “rms” package, we established a nomogram model to predict the occurrence of sepsis. “Points” indicate the score of the corresponding factor below, and “Total Points” indicate the summation of all the scores of factors above. We then used calibration curves to assess the predictive power of the nomogram model. Finally, we evaluated the clinical value of the model using decision curve analysis.

The co-expressed gene network of hub RMGs was constructed using COEXPEDIA (https://www.coexpedia.org). The potential miRNAs targeting key RMGs were downloaded from miRanda (https://www.cs.kent.ac.uk/people/staff/dat/miranda/) and miRDB (http://mirdb.org/) databases, and overlapping miRNAs were selected to depict network. The final gene-miRNA and the co-expression network associated with sepsis were visualized using Cytoscape software.

This study included 39 adult patients from our hospital, consisting of 26 septic patients diagnosed with sepsis according to the Third International Consensus Definitions for Sepsis and Septic Shock (Sepsis-3) (1), and 13 healthy controls. Healthy volunteers were those who came to East Hospital for routine physical examination. Peripheral blood samples and corresponding clinical data were collected upon admission. The characteristics of the patients are shown in Table 1 . PBMCs were isolated within 4 hours after collection.

This study was approved by the Research Ethics Board of East Hospital, Tongji University (Shanghai, China). Written informed consent was obtained from all recruited patients or their authorized family members.

Blood samples in EDTA-containing tubes were centrifuged at 3000 rpm for 5 min at room temperature. Serum was removed, and PBMCs were separated from the remaining blood using Ficoll-Paque density gradient centrifugation, following the instruction manual (Solarbio, Cat. NO. P8610). The PBMCs were then stored at -80°C until further testing.

RNA was extracted from cells using TRIzol (Gene Copoeia, MD, USA), and 1 µg total RNA from each sample was reverse transcribed to cDNA using PrimeScript^® RT Master Mix (Takara, Cat. NO. RR036A). SYBR^® Green (Applied Biosystems™, Cat. NO. 4309155) was used for qRT-PCR analysis. The ΔCt method was used to analyze mRNA levels relative to glyceraldehyde-3-phosphate dehydrogenase (GAPDH). The primer sequences were listed in Table S3 .

GraphPad Prism 9.4 (GraphPad Inc, San Diego, USA) was used for statistical analyses. Normally distributed data were compared using Student’s t-test or one-way analysis of variance (one-way ANOVA), and the results were shown as the mean ± SD. Nonnormally distributed data were analyzed using the non-parametric Mann-Whitney U test, and the results were expressed as the median and interquartile range (IQR). Categorical variables were compared with the chi-square or Fisher’s exact test, and the results were shown as numbers and percentages. All statistical tests were two-sided. P < 0.05 was regarded as statistically significant.

The overall strategy of this study was presented in Figure 1 . To analyze the expression of RMGs in PBMCs from patients with sepsis, we first performed differential expressed genes (DEGs) analysis of the integrated gene expression matrix using the limma package. DEGs in PBMCs from healthy donors and septic patients were shown in heatmap ( Figure 2A ). Most of the genes that mediate one certain type of RNA methylation modifications displayed a similar expression pattern. For instance, genes mediating RNA Ψ modification (TRUB1, TRUB2, PUS1, RPUSD3, RPUSD4, PUS7, RPUSD2) were mainly suppressed in sepsis ( Figure 2A , Table S4 ).

We then investigated the correlation among RMGs by using Spearman’s correlation analysis and found that there was a significantly close correlation among these differentially expressed RMGs. Of note, YTHDF3 expression was positively correlated with the expression of TET2, whereas PCIF1 exhibited a negative correlation with most other regulators ( Figure 2B ). The correlated expression among RMGs suggested their functional synergy or antagonism in the pathogenesis of sepsis. Moreover, we performed a PPI network analysis of RMGs. Proteins with an interaction score ≥0.4 were selected and visualized ( Figure 2C ). The network comprised 40 nodes and 257 edges, representing genes and interactions between genes, respectively ( Table S5 ), also indicating the functional association in sepsis. The chromosomal locations of RMGs were depicted in Figure 2D .

Given that the differential expression levels of RNA modification genes and clustering stability, we employed consensus clustering, an unsupervised clustering method to obtain a robust ranking for subsequent analysis, and classified sepsis into three subtypes (cluster 1 to 3). The cumulative distribution function (CDF) plot displayed the consensus distributions for each cluster ( Figure 3A ). The delta area plot showed the relative change in the area under the CDF curve ( Figure 3B ). The CDF distribution was smoother when k = 3 or 4, while the increase in area under the CDF curve was relatively less at k = 4 than that at k = 3, suggesting that k = 3 was the optimal number of clusters. As shown in the Consensus matrix heatmap ( Figure 3C , S1 ), 3 clusters showed clear boundaries, indicating good cluster stability over repeated iterations. The principal component analysis (PCA) results also indicated the justification of the three clusters ( Figure 3D ).

In addition, the expression of RMGs in each subtype displayed distinct patterns (P < 0.05) ( Figures 3E, F ). In cluster 1, FTO was relatively downregulated, and of all the highly expressed genes, NSUN7 and TRMT6 were most upregulated. These genes were involved in the regulation of m6Am, m5C and m1A RNA methylation modification. NSUN3, encoding a key methyltransferase for RNA m5C modification showed a significantly increased expression in cluster 2. In cluster 3, among the relatively highly expressed regulators, NSUN5 and NSUN6 significantly increased, which were also correlated with m5C modification. Those results indicated that m5C RNA methylation may play a role in sepsis classification.

Besides, we analyzed the basic clinical characteristics of sepsis patients in GSE65682 in 3 clusters ( Table S6 ) and found significant differences in Molecular Diagnosis and Risk Stratification of Sepsis (MARS) endotypes among the three clusters (P < 0.001). In cluster 1, MARS 2 endotypes accounted for the highest proportion (42.44%). The proportion of MARS 1 (47.12%) was significantly higher in cluster 2 than that in the other two groups (P < 0.05). Cluster 3 was dominated by MARS 3 (35.56%). Previous studies have revealed that MARS 1 was consistently correlated to a poor outcome (31). But we did not observe significant differences in mortality among the three clusters (P = 0.074) ( Table S6 ).

By using CIBERSORT and GSVA analyses, we further identified a significant heterogeneity in terms of immune cell infiltration and molecular pathways among the three subtypes. The abundance and proportions of infiltrating immune cells differed among the three types, of which neutrophils and monocytes accounted for higher proportions as compared to the other immune cells ( Figure 4A ). In cluster 1, the frequencies of neutrophils and gamma delta T (γδT) cells were relatively increased compared to other clusters, while CD8 T cells, regulatory T cells (Tregs), activated NK cells were down-regulated (P < 0.05). The distribution of immune cells in cluster 2 and 3 showed relatively opposed patterns as compared to that in cluster 1 (P < 0.05). In cluster 3, CD8 T cells and monocytes were relatively upregulated, while γδT cells and neutrophils were down-regulated (P < 0.05). Immune cells that were up- or down-regulated in cluster 3 showed the same trend in cluster 2 with a less degree than those in cluster 3. Collectively, these results suggested that the three clusters have different immune infiltration in septic patients.

We then evaluated the different signaling pathways in these clusters. Figures 4B–D showed relatively up-regulated (11) and down-regulated (10) signaling pathways in the three sepsis subgroups. Compared to other subtypes, cluster 1 showed enhanced inflammatory signaling pathways activation, such as the Leukocyte transendothelial migration and MAPK signaling pathway. The Down-regulated pathways in cluster 1 were involved in amino acid metabolism ( Figure 4B ). In cluster 2, we found that up-regulated signaling pathways were associated with amino acids, glycosphingolipids and unsaturated fatty acids metabolism. Down-regulated signaling pathways were mainly related to maintaining the stability of the genome and amino acid metabolism ( Figure 4C ). In cluster 3, up-regulated signaling pathways were mainly involved in carbohydrate metabolism pathways. Down-regulated signaling pathways were related to inflammatory diseases ( Figure 4D ). These functional analyses indicated that septic patients with different RMG expression patterns have distinct molecular pathways.

Next, we performed WGCNA to analyze the gene co-expression networks and identify biologically meaningful modules that corresponded to designated phenotype-related genes. A scale free co-expression network was established with the soft threshold power as 12 (scale-free R ^{2 =} 0.90) ( Figures 5A, B ) and cut height as 0.25 ( Figure 5C ). The cluster dendrogram was displayed in Figure 5D . WGCNA identified 18 modules in the sepsis cohort ( Figure 5E ). Cluster 1 was positively correlated with MEblack module (containing 390 genes) ( Table S7 ), with a correlation coefficient of 0.64 and P value of 6.00E-112 ( Figure S2A ). Cluster 2 was positively associated with MEblue module (containing 1510 genes) ( Table S7 ), with a correlation coefficient of 0.58 and P value of 1.00E-84 ( Figure S2B ). Cluster 3 was significantly positively correlated to MEyellow (containing 547 genes) ( Table S7 ), with a correlation coefficient of 0.62 and P value of 4.00E-100 ( Figure S2C ). To identify genes highly correlated with the key module in each cluster, we set module membership (MM) > 0.8 and gene significance (GS) > 0.3, and results were shown in Table S8 .

To further explore the biological function of each module, we conducted GO and KEGG enrichment analyses ( Figures 6A–F ). GO enrichment analysis demonstrated that the biological functions of MEblack module (correlated with cluster 1) were mainly enriched in Neutrophil activation ( Figure 6A ); the MEblue module (correlated with cluster 2) was mainly related to cellular catabolic processes and mRNA stability ( Figure 6B ); and the MEyellow module (correlated with cluster 3) was mainly associated with Ribonucleoprotein complex biogenesis and carbohydrate metabolism pathways ( Figure 6C ). KEGG enrichment analysis showed that MEblack module (correlated with cluster 1) was mainly enriched in MAPK signaling pathway and cytokine-cytokine receptor interaction signaling pathway ( Figure 6D ); the MEblue module (correlated with cluster 2) was mainly related to Cushing syndrome and autophagy ( Figure 6E ); and the MEyellow module (correlated with cluster 3) was mainly associated with ribosome and Coronavirus disease-COVID-19 ( Figure 6F ).

To identify hub RMGs with diagnostic value for sepsis, firstly, we refined the selection of RMGs by four algorithms, namely the RF, SVM, XGB and GLM, and used differentially expressed RMGs in Figure 2A to construct the four diagnostic models. We further evaluated the efficiency of the four models by usage of “DALEX” R package. As shown in Figures 7A–C , SVM model was the most suitable model for its least sample residual. In addition, the performance of the four models was evaluated by using ROC curves, SVM model displayed the highest AUC value (AUC=0.998) ( Figure 7D ). To further assess the accuracy and reliability of the SVM model, we conducted ROC analysis of three GEO databases. The AUC amount to 0.808 in GSE57065 ( Figure 8A ), 0.976 in GSE65682 ( Figure 8B ) and 0.983 in GSE95233 ( Figure 8C ). Therefore, the above results indicated that the SVM diagnostic model exhibited excellent predictive efficiency.

Subsequently, according to the feature importance of models in Figure 6B , the top 5 genes (NSUN7, FTO, NOP2, PUS1, and PUS3) selected by SVM were used for the analysis of diagnostic performance. AUC values of the ROC of these signature genes were 0.912 for PUS1, 0.962 for NSUN7, 0.903 for NOP2, 0.688 for PUS3, and 0.966 for FTO in GSE57065 ( Figure 8D ). In GSE65682, the AUC values of ROC were 0.855 for PUS1, 0.983 for NSUN7, 0.950 for NOP2, 0.730 for PUS3, and 0.997 for FTO, respectively ( Figure 8E ). Finally, the AUC values of ROC were 0.861 for PUS1, 0.975 for NSUN7, 0.974 for NOP2, 0.484 for PUS3, and 0.998 for FTO in GSE95233 ( Figure 8F ). Table S9 reports the specificity, sensitivity and optimal cutoff point for discriminating between sepsis and healthy controls in GSE57065, GSE65682, and GSE95233, respectively. These results indicated that the screened signature genes particularly NSUN7 and FTO exhibited remarkable diagnostic efficiency in sepsis.

To predict the occurrence of sepsis, we constructed a nomogram model for diagnosing sepsis based on the expression levels of five hub RMGs (NSUN7, FTO, NOP2, PUS1 and PUS3) ( Figure 8G ). Decision curve analysis (DCA) revealed that the “nomogram model” curve was superior to the gray line, suggesting that the patients could benefit from the nomogram model at high-risk threshold from 0 to 1 and the clinical advantage of the nomogram model was higher ( Figure 8H ). We further assessed the predictive power of nomogram model via a calibration curve. The calibration curve showed no significant difference between the actual sepsis risk and the predicted risk, indicating that this model was accurate in predicting sepsis ( Figure 8I ). These findings suggested that the five hub genes are efficient for the diagnosis of sepsis.

We further investigated whether the hub RMGs were associated with immune cell infiltration in sepsis ( Figure 9 ). We observed that PUS3 was negatively correlated with resting mast cells and neutrophils. PUS1 was negatively correlated with M1 macrophages, M2 macrophages, neutrophils, plasma cells and γδT cells, but positively correlated with monocytes, activated NK cells and CD8 T cells. Interestingly, similar to PUS1, NOP2 and FTO were also positively correlated with monocytes, activated NK cells and CD8 T cells. Conversely, NSUN7 displayed a negative correlation with monocytes, activated NK cells and CD8 T cells, but positive correlation with neutrophils, γδT cells and M1 macrophages. These results suggested that PUS1, NOP2, and FTO may act synergistically in sepsis, while NSUN7 may counteract the function of aforementioned three genes.

Furthermore, by determining the expression of hub genes from single-cell transcriptome in Human Protein Atlas database (https://www.proteinatlas.org ), we found that NOP2, FTO, PUS1, and PUS3 were all enriched in T cells, however, the NSUN7 gene is not expressed in any immune cells in PBMCs from healthy subjects ( Figure S3 ). In addition, to better understand the genes and microRNA associated with key RMGs, we constructed a gene co-expression analysis and mRNA-miRNAs networks ( Figure S4 , Table S10 - S11 ), suggesting that NOP2 and PUS1 may be functionally synergistic with each other. Genes associated with hub RMGs may linked to the progression of sepsis, NOP2 and NAT10 were co-expressed, and previous research reported that NAT10 improves the survival and ameliorates lung injury in septic mice by inhibiting neutrophil pyroptosis (32). We found miRNAs such as miR-21-3p (33) and miR-126-5p (34) were also crucial components in regulating the sepsis pathogenesis by interacting with key RMGs.

To validate the expression and diagnostic potential of the five hub RMGs in clinical samples, we prospectively recruited patients with sepsis and healthy controls. The general clinical characteristics of the patients was shown in Table 1 . The cohort was comprised of 11 patients with septic shock, 15 sepsis patients without shock (sepsis), and 13 healthy volunteers. There were no significant differences in sex or age between the two (septic shock and sepsis) (P = 0.0965, P = 0.0740) or three groups (septic shock, sepsis and health) (P = 0.0501, P = 0.0772). Additionally, there were no significant differences in comorbidities or sites of infection between the two groups (septic shock and sepsis). Compared to the sepsis group, the patients in the shock groups exhibited higher Sequential Organ Failure Assessment (SOFA) and Acute Physiology and Chronic Health Evaluation (APACH) II scores (P = 0.0020, P = 0.0001), longer Mechanical ventilation time and higher mortality rate at day 7, 14 and 28 (P = 0.0177, P = 0.0020, P = 0.0038, P = 0.0008).

Firstly, we determined the genes expression in PBMCs by qRT-PCR analysis, and found the expression of NSUN7, NOP2, PUS1 and PUS3 in septic patients (septic shock and sepsis) at day 1 were significantly higher than in healthy volunteers, and FTO expression was lower in septic patients ( Figure 10A ). However, the expression levels of the five genes displayed no significant difference in survivors and non-survivors in septic patients at day 28 ( Figure S5A ). In addition, their expression did not show any significant difference between septic shock and sepsis patients ( Figure S5B ). Notably, the expression of NOP2 and PUS1 was concurrently increased in patients with sepsis of group C3 in the integrated gene expression matrix ( Figure 2F ), which was consistent with the results observed in our clinic samples.

Then, we investigated the diagnostic value of hub RMGs in sepsis using ROC analysis. The AUC of NSUN7, NOP2, PUS1, PUS3 and FTO at day 1 was 0.828, 0.707, 0.846, 0.834 and 0.976 respectively ( Figure 10B ). Notably, FTO exhibited the highest diagnostic value among these genes, with sensitivity of 0.846 and specificity of 1 ( Table S12 ). Furthermore, we evaluated the combined diagnostic value of NSUN7, NOP2, PUS1, PUS3 and FTO (the integrated RMGs model) among non-septic shock and septic shock patients. We compared the AUC of integrated RMGs model with Lactate (Lac) and procalcitonin (PCT) at day 1 and found that the AUC of the integrated RMGs model was not inferior to Lac and PCT (0.72 vs. 0.68, P = 0.78; 0.72 vs. 0.52, P = 0.24) ( Figure 10C ), with a sensitivity of 0.467 and specificity of 1 ( Table S13 ).