Two Homo sapiens gene array expression series matrix files of tuberculosis peripheral blood samples, including GSE54992 (14) and GSE83456 (15), were collected from GEO (Table 1). GSE54992, in which tests were performed on the Affymetrix Human Genome U133 Plus 2.0 Array (GPL570 platform, Affymetrix, Inc.), contains the expression profiles of 39 samples, which comprise 27 active pulmonary TB cases, six healthy controls, and six latent TB cases. A total of 33 TB cases and healthy controls were enrolled, and 6 samples of latent TB were excluded from our investigation. GSE83456(GPL10558) contains 45 humans with pulmonary TB, 47 humans with extra-pulmonary TB, 49 cases of pulmonary sarcoidosis, and 61 healthy human controls. A total of 106 patients with PTB and healthy controls were recruited for this study. The flowchart of the study is displayed in Figure 1. The m⁶A differentially expressed genes (DEGs) analysis between pulmonary TB and healthy controls were conducted by using “limma”R packages (version 3.64.1) (16). A volcano plot was used to visualize the expression of DEGs through “pheatmap”R packages (version 1.0.13) (17) and a heatmap diagram was performed to exhibit the expression of m⁶A regulators in TB and normal control by “ggplot” R packages (version 3.5.2) (18).

Twenty-two widely studied m⁶A methylation regulators containing Eraser: FTO, ALKBH5; Writer: VIRMA, WTAP, ZC3H13, METTL14, METTL3, CBLL1, RBM15, RBM15B; Reader: YTHDF1, YTHDF2, YTHDF3, YTHDC1, YTHDC2, HNRNPC, HNRNPA2BP1, IGF2BP1, IGF2BP2, IGF2BP3, ELAVL1, LRPPRC were confirmed from published reports. Spearman correlation analysis was conducted to evaluate the association between m⁶A RNA methylation regulators and TB by using “corrr”R packages (version 0.4.5), and the results were visualized and plotted through “ggplot2” and “pheatmap” R packages (version 3.5.2, version 1.0.13respecively). Random forest (RF) algorithm was used to screen and identify TB-related m⁶A RNA methylation regulators Least absolute shrinkage and selection operator (LASSO) regression was performed to construct the prognostic model (19). Receiver Operating Characteristic (ROC) was used to assess the diagnostic performance of the established model.

Single-sample Gene Set Enrichment Analysis (ssGSEA) (20) was employed to assess the relative abundance of specific infiltrating immune cells and the activity of immune responses. This analysis was conducted through the “GSVA”package (version 2.2.0) of R software (version 4.0.2), which employing its built-in ssgsea method for scoring. The list of genes of infiltrated immune cells, including activated CD8+ T cells, natural killer T cells, Regulatory T cells (Tregs), activated dendritic cells, and macrophages were obtained from previous studies (21). The list of immune response genes was obtained from the ImmPort database (22). The correlation between m⁶A RNA methylation regulators and the proportion of immune cells and immune reaction activity was evaluated using Spearman correlation analysis.

Based on the expression profiles of 22 m⁶A RNA methylation regulators, unsupervised consensus clustering analysis was performed on TB samples by using the R package ConsensusClusterPlus (Wilkerson & Hayes, 2010) (23). The clustering analysis employed Euclidean distance as the similarity measure, combined with the k-means algorithm for clustering, and generated a consensus matrix through 1,000 resampling iterations with approximately 80% of samples participating in each iteration to evaluate clustering stability. By comparing the cumulative distribution function (CDF) curves, delta area curves, and consensus heatmaps for different cluster numbers (k = 2–9), the optimal number of clusters was determined. Subsequently, principal component analysis (PCA) was performed based on the expression matrix of the 22 m6A regulators for dimensionality reduction, aiming to visualize the distribution differences between the two modification patterns and validate the clustering effect. Among different m6A modification patterns, the expression levels of m6A regulators, immune cell infiltration abundance, immune response scores, and HLA gene expression were compared. For normally distributed variables, student’s t-test was used for comparisons between the two groups. For non-normally distributed variables, the Mann–Whitney U test (Wilcoxon rank-sum test) was applied. All statistical tests were two-sided, with the significance level set at P < 0.05.

To analyze the effect of m⁶A modification patterns on TB, we used the “limma”R packages (version 3.64.1) to analyze the DEGs of the two m⁶A modification patterns. The significance criteria for the determination of DEGs were as |log2FC (fold change) | >0.5 and P.adj< 0.01. Meanwhile, logFC>0.5 and P.ajj <0.01 was defined as upregulated genes. LogFC <-0.5 and P.adj<0.01 was considered as downregulated genes.

To explore the potential biological functions and signaling pathways of differentially expressed genes (DEGs) under different m⁶A modification patterns, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis were performed using the R package clusterProfiler (version 4.16.0). The GO analysis included three categories: Biological Process (BP), Molecular Function (MF), and Cellular Component (CC). For the enrichment analysis, the org.Hs.eg.db database was used as the human gene annotation reference, and all genes detected in the combined GEO dataset were set as the background gene set. Enrichment calculations were conducted using the enrich GO and enrich KEGG functions. Multiple testing correction was applied using the Benjamini–Hochberg method with a significance threshold at adjusted p-value (Padj) < 0.05. Significantly enriched functional categories and signaling pathways were visualized using dot plots (dotplot) and bar plots (barplot).

The “c2.cp.kegg.v7.5.1.entrez.gmt” (25) and “h.all.v 7.5.1. entrez.gmt” (26) data were acquired from Molecular Signatures Database (MLgDB), which was analyzed through GSEA by using R “clusterProfiler”package (24).

To investigate the biological functional differences between different m⁶A modification-based TB subtypes, Gene Set Variation Analysis (GSVA) was applied to quantitatively evaluate pathway activation levels. The gene sets used in the analysis were collected from the HALLMARK gene sets in the Molecular Signatures Database (MSigDB, https://www.gsea-msigdb.org/gsea/index.jsp) (25). GSVA analysis was performed using the R package “GSVA” (version 2.2.0) (20) and employing a non-parametric kernel method to calculate an enrichment score for each sample in each specific pathway. The main parameters were set as follows: method = “gsva”, kcdf = “Gaussian”, mx.diff = TRUE. Subsequently, the R package “limma” (version 3.64.1) was used to compare the GSVA pathway scores between different m⁶A modification-based TB subtypes. The activation score for each pathway was input as the dependent variable into a linear model for different test without other covariates. To control the false discovery rate (FDR) from multiple hypothesis test, the results were adjusted using the Benjamini–Hochberg method with a corrected FDR < 0.05 as the threshold for statistical significance. The significant differences were visualized with a volcano plot.

This case-control investigation included 34 patients with pulmonary tuberculosis who were diagnosed according to clinical laboratory tests including blood, sputum or bronchoalveolar lavage fluid, simple skin tests, and histological findings. 34 age- and sex-matched healthy individuals who underwent a physical examination at the First Affiliated Hospital of Chengdu Medical College were free of diabetes, hypertension, heart, liver, kidney, and other organ diseases or dysfunctions, and had a history of malignant tumors and severe organ dysfunction. None of all TB patients with TB received the standard antituberculosis treatment.

Peripheral venous blood samples were collected from all the participants in the early morning under fasting conditions. Two milliliters of peripheral blood was collected into EDTA anticoagulant tubes. The samples were either immediately analyzed or aliquoted for a single use and stored at -80°C until further use. The study was approved by the Ethics Committee of the First Affiliated Hospital of Chengdu Medical College (Sichuan, China) and conducted in accordance with the Declaration of Helsinki. Written informed consent was obtained from all participants.

Total RNA was isolated from PBMC of TB patients and healthy controls using a whole blood total RNA extraction kit (Simgen, China) according to the manufacturer’s protocols. The concentration and purity of each total RNA were detected (using the A260/A280 and A260/A230 ratios) through a NaoDrop ND-1000 spectrophotometer (Invitrogen). For PCR analysis, 1 mg of total RNA was used to synthesize cDNA by reverse transcription using a PrimeScript TM RT reagent kit (Tiangen Biotech, China) following the manufacturer’s protocol. The product was used as a template for PCR in a CFX-96 real-time PCR system that employed SYBR VRPremix Ex TaqTM II (Tiangen Biotech, China). The primer sequences used for amplification are listed in Table 2. Relative expression of each gene was determined using the 2^-△△Ct method.

All data analyses were performed in R software (version 4.0.2). For gene expression analysis, after normalization and log transformation, the data approximated a normal distribution under the large sample assumption. Therefore, differential analysis was conducted using the limma package (version 3.64.1). For comparisons of continuous variables between independent samples, student’s t-test was applied when the data met the assumptions of normal distribution and homogeneity of variance. If the variables did not follow a normal distribution (e.g., Shapiro–Wilk test p < 0.05) or the sample size was small with markedly skewed distribution, the Mann–Whitney U test was used. All statistical tests were two-sided, and a p-value < 0.05 was considered statistically significant.

The expression correlations among m6A regulators were calculated using Spearman’s rank correlation analysis. The analysis was based on the paired expression values of each regulator, and the results were expressed as correlation coefficients (ρ). The significance of correlations was determined based on p-values without adjustment for multiple tests. The correlation coefficient matrix and the significance results were visualized via heatmaps and network plots.

The receiver operating characteristic (ROC) curve and area under the ROC curve (AUC) were calculated to evaluate the feasibility of m⁶A regulators as potential markers for TB diagnosis. Meanwhile, the discriminative performance of the model was evaluated using the ROC curve. To compare the predictive ability of the model across different datasets, ROC curves were calculated based on the training set and the validation set, respectively. The area under the curve (AUC) and its 95% confidence interval were computed using the R package “pROC” (Version 1.18.5). To assess the robustness of the model, the study samples were randomly divided into a training set and a validation set at a ratio of 7:3 for model fitting and validation without additional resampling or cross-validation.

The workflow of this study is illustrated in Figure 1. After the integration of the GEO datasets GSE54992 and GSE83456, the raw expression data were first subjected to background correction and normalization using “limma”R packages (version 3.64.1), followed by batch-effect adjustment via the empirical Bayesian algorithm implemented in the “sva” package to ensure comparability across datasets. The data of boxplots of normalized expression values demonstrated that after normalization, the median levels of the boxplots were well aligned, indicating comparable distributions across samples (Figures 2A-D). Subsequently, principal component analysis (PCA) was conducted to evaluate overall clustering characteristics and group consistency. The PCA results (Figures 2E, F) revealed clear separation between TB group and healthy control group in the principal component space. PC1 and PC2 explained 13.52% and 9.74% of the total variance, respectively, indicating that the first two principal components can effectively capture the primary variation structure among the samples.

Then, the landscape of genetic expression of m⁶A RNA methylation genes between the TB group and control group was analyzed through “limma” R packages. Volcano plots of DEGs in the above two datasets showed that METTL3, VIRMA, RBM15, RBM15B, YTHDF1, YTHDC1, YTHDC2, ELAVL1, LRPPRC and ALKBH5 were downregulated in TB (Figure 3A). The data further demonstrated that the results of the heat map (Figure 3B), box plot (Figure 3C) and chromosome map (Figure 3D) were the same as those of the volcano plot.

Moreover, the validation and clinical relevance of different m⁶A regulators in TB were evaluated. We examined the mRNA expression of METTL3, VIRMA, RBM15, RBM15B, YTHDF1, YTHDC1, YTHDC2, ELAVL1, LRPPRC and ALKBH5 in peripheral blood by using qRT-PCR. METTL3, VIRMA, YTHDF1, YTHDC1, YTHDC2, ELAVL1 and LRPPRC mRNA levels were significantly downregulated in TB samples compared to those in the control group (Figures 4A-J). In addition, we evaluated the association between serum METTL3, VIRMA, YTHDF1, YTHDC1, YTHDC2, ELAVL1, LRPPRC mRNA expression, and clinical markers in TB patients. As shown in Figure 5, there was a close association between LRPPRC and Lymphocyte percentage (Lymph%) (r=0.358, P = 0.0002) and number (Lymph#) (r=0.415, P<0.0001) (Figures 5A-P).

The correlation between the expression levels of the 22 m⁶A genes in TB and normal samples was analyzed. The results were visualized using a heat map (Figure 6A) and network map (Figure 6B). The upper right section presents the expression correlation of m⁶A -related genes in all samples, whereas the lower left section shows the expression correlation of m⁶A -related genes in TB samples. These data revealed that YTHDF2 was strongly associated with YTHDF1 in the TB group (Figure 6C). YTHDF1 was highly correlated with VIRMA in all the samples ((Figure 6D).

To further explore the diagnostic ability of m⁶A -related regulators in TB, the random forest method was used (Figures 7A, B). The samples were randomly divided into the training (70%) and validation (30%) sets. Boxplots (Figures 7C, D) revealed significant differences in the model scores between the TB and healthy groups in both the training and validation sets. The ROC curve (Figure 7E) demonstrated that the constructed model exhibited excellent diagnostic performance for TB, indicating that m⁶A genes had strong predictive power.

LASSO regression analysis was used to screen variables and construct a prediction model as follows: risk scores = METTL3 × (-1.265) + METTL14 × 0.598 + WTAP × (-1.559) + RBM15 × (-0.926) + RBM15B × (-0.432) + CBLL1 × 1.623 + YTHDF2 × 0.564 + YTHDC1 × (-0.744) + YTHDC2 × (-1.844) + HNRNPC × (-1.119) + HNRNPA2B1 × 1.109 + IGF2BP1 × 0.559 + IGF2BP3 × (-2.372) + ELAVL1 × (-1.273) + LRPPRC × (-2.031) + ALKBH5 × 1.297. As showed in Figures 8A, B, the results were the same as those in the above conclusion. The boxplot shows a significant difference in the risk scores between the TB and healthy groups (Figure 8C). The ROC curve (Figure 8D) indicated that the risk model had strong diagnostic capability for TB.

In addition, the predictive values of METTL3, VIRMA, YTHDF1, YTHDC1, YTHDC2, ELAVL1and LRPPRC individually and in combination with Ziehl-Neelsen staining were evaluated by using a binary logistic regression model. Receiver operating characteristic (ROC) curves were plotted to analyze and compare their predictive values (Figure 9). The areas under the ROC for predicting TB using METTL3, VIRMA, YTHDF1, YTHDC1, YTHDC2, ELAVL1and LRPPRC were 0.543 (95% CI: 0.331–0.755), 0.564 (95% CI: 0.354–0.774), 0.462 (95% CI: 0.254–0.669), 0.521(CI:0.313-0.730), 0.556(CI:0.350-0.761), 0.564(CI:0.358-0.77) and 0.436(CI:0.222-0.649) respectively, indicating that none of the individually identified genes can possess sufficient diagnostic accuracy for clinical application. When compared to individual predictions, the combined ROC area under the curve for these seven markers and Ziehl-Neelsen staining was 0.953 (95% CI: 0.874–1.00) (P<0.001).

Immune cell infiltration is a vital component of the tumor microenvironment and is closely associated with the development of various diseases (Gajewski et al., 2013). Therefore, the relationship between m⁶A regulators, immune cell infiltration, and immune response was explored in TB (Figures 10A, B). FTO and METTL3 were positively correlated with most immune cells, whereas LRPPRC exhibited a negative correlation. Type 1 T helper cells were positively associated with most m⁶A genes, whereas gamma delta T cells were negatively correlated. In terms of immune processes, ALKBH5 and METTL3 levels were negatively correlated with immune processes. These findings suggested that m⁶A regulators could serve as predictors of the immune microenvironment and immune processes in TB, with METTL3 acting as a significant immunosuppressive regulator, and HNRNPA2B1 acting as an immune-activating gene.

Based on the expression of 22 m⁶A regulators, unsupervised cluster analysis was applied to classify the different m⁶A modification patterns in TB (Figures 11A-C). The data demonstrated that the optimal clustering stability was k = 2 in the consensus clustering. Patients with TB were divided into two subtypes (clusters 1 and 2). Furthermore, PCA revealed a statistically significant difference between the two m⁶A molecular subtypes (Figure 11D). Heat maps ((Figure 11E) and boxplots ((Figure 11F) demonstrated the expression specificity of m⁶A regulators between these two molecular subtypes.

m⁶A modification may affect the translation efficiency or stability of immune-related genes, leading to differences in immune response intensity between the two subtypes. To explore the differences in immune microenvironment features between different m⁶A modification patterns, we evaluated immune cell infiltration, immune response, and HLA expression under different m⁶A modification patterns. As shown in Figure 12A, the number of immune cells differed between the two groups. Pattern 1 exhibited a relatively high proportion of activated immune cells, including gamma delta T cells, neutrophils, and NK cells. Pattern 2 showed a higher proportion of type 1 T-helper cells. Chemokine and cytokine processes were relatively more active in pattern 1, whereas the TCR signaling pathway was highly active in pattern 2 (Figure 12B), which aligns with previous analysis of immune processes. To enhance the reliability of the results, CIBERSORT was used to calculate the immune cell content in the different modes (Figure 12C). Similar to ssGSEA, there was a high level of Gamma Delta T cells in Pattern 1. However, there was no significant difference in HLA family gene expression between the two patterns (Figure 12D). These findings suggested that pattern 1 mediates an active immune response, whereas pattern 2 regulates a mild immune response, revealing the important role of m⁶A methylation in regulating the formation of different immune microenvironments in TB.

The two subtypes may represent different biological states. We use “limma” package to investigate the biological responses under the two m⁶A modification patterns. Volcano and heat maps were used to display the results of the difference analysis (Figures 13A, B). We set the thresholds for differentially expressed genes (DEGs) as |log2 fold change (logFC)| > 0.5 and adjusted P-value (Padj) < 0.01. Enrichment analysis was conducted on the DEGs (Figures 13C, D). GSEA was used to conduct an enrichment analysis for the two patterns (Figures 13E, F). The KEGG and HALLMARK results showed that pattern 2 was more biologically active (Tables 3, 4).

To further elucidate the impact of m⁶A modification on the transcriptomic landscape of tuberculosis (TB), this study conducted a secondary subtype analysis based on differentially expressed genes (DEGs), which was built upon the previously identified m⁶A modification patterns derived from 22 m⁶A regulators. This analysis aimed to extend from the m⁶A regulatory level to the gene expression level and systematically uncovered the downstream effects of m⁶A modification on TB molecular heterogeneity. The clustering results (Figures 14A–C) showed that the m6A-related signature genes robustly distinguished two expression subtypes among TB samples. PCA analysis (Figure 14D) further validated clear separation between the two subtypes in transcriptomic space, suggesting that m6A modification may drive distinct downstream transcriptional responses. The heatmap (Figure 14E) displayed differential expression patterns of the signature genes across the two subtypes, reflecting potential functional divergences in metabolic regulation, immune responses, and signaling pathways. The Sankey diagram (Figure 14F) revealed the relationship between the m⁶A modification patterns and expression subtypes, indicating that different m⁶A modification patterns may shape distinct transcriptional states by regulating specific gene networks.

As shown in Figure 15A, no significant differences were observed in immune cell infiltration, immune processes, or expression of HLA family genes (Figures 15B-E), suggesting that the differences between the two subtypes might not be driven by immune infiltration. We then employed the GSVA package to convert the expression matrix into a pathway activation score matrix in the “h.all.v7.5. symbols” gene set. The R package limma was used to compare pathway activation scores between the two subtypes, and a volcano plot was generated (Figure 15F). Subtype A primarily activates the hallmark heme metabolism pathway, whereas subtype B mainly activates the hallmark TGF beta signaling pathway.