Unraveling neutrophil dynamics in pulmonary tuberculosis with insights from transcriptome data and neutrophil extracellular traps-related genes

Data source and preprocessing

In the present study, a single-cell dataset (GEO number: GSE167232) [21] from mice infected with Mtb was obtained, comprising 17,101 cells and spanning 14 cell types: alveolar macrophages 1 (AM_1), AM_2, AM_3, AM_4, interstitial macrophages 1 (IM_1), IM_2, IM_3, IM_4, Monocytes, Neutrophils, Dendritic cells (DC), DC.103+11B-, Macrophages/T-cells, Macrophages/ B-cells. Additionally, the single-cell dataset encompasses three distinct disease states: Infected, Bystander, and Uninfected. Furthermore, 49 transcriptomic studies related to patients with TB and other clinical conditions were obtained from the “curatedTBData” package (https://github.com/wejlab/curatedTBData). To streamline integration, datasets with fewer than 15,000 genes or those not sourced from human peripheral blood were excluded. This refinement resulted in 31 high-quality datasets, including GSE101705 (n = 44) [22], GSE107104 (n = 33) [23], GSE107731 (n = 6), GSE112104 (n = 51) [24], GSE19435 (n = 33) [25], GSE19439 (n = 42) [25], GSE19442 (n = 51) [25], GSE19443 (n = 44) [25], GSE19444 (n = 54) [25], GSE22098 (n = 274) [25], GSE25534 (n = 102) [26], GSE28623 (n = 108) [27], GSE29536 (n = 15) [28], GSE34608 (n = 44) [29], GSE37250 (n = 537) [30], GSE39939 (n = 157) [31], GSE39940 (n = 334) [31], GSE40553 (n = 204) [32], GSE41055 (n = 27) [33], GSE42825 (n = 42) [34], GSE42826 (n = 102) [34], GSE42830 (n = 95) [34], GSE42832 (n = 90) [34], GSE50834 (n = 44) [35], GSE54992 (n = 39) [36], GSE56153 (n = 71) [37], GSE62147 (n = 52) [38], GSE62525 (n = 42) [39], GSE69581 (n = 50) [40], GSE73408 (n = 109) [41], GSE81746 (n = 8) [42], GSE83456 (n = 202) [43], GSE83892 (n = 116) [44], GSE84076 (n = 36) [45], GSE89403 (n = 453) [46], GSE94438 (n = 418) [47], GSETornheim (n = 100). The expression values for these datasets have been standardized in “curatedTBData” package. To mitigate potential batch effects across disparate datasets, the expression profiles corresponding to each dataset were amalgamated into a unified matrix. Subsequently, the Quantile Normalization strategy was implemented to systematically rectify batch effects. Additionally, 69 gene sets related to neutrophil extracellular traps were obtained from published studies [16, 48], with detailed information available in Supplementary Tables 1 and 2.

Significantly, our cohort comprised 31 human TB patients, stratified into four primary categories: healthy controls (Control), individuals with subclinical infection (Subclinical), those diagnosed with extrapulmonary tuberculosis (EPTB), and individuals with pulmonary tuberculosis (PTB), as well as those with latent tuberculosis infection (LTBI). Subclinical denotes a condition where individuals are infected but do not manifest overt symptoms, while EPTB refers to tuberculosis affecting organs or tissues outside the lungs. PTB is the predominant form, instigated by Mtb and typically characterized by pulmonary symptoms. LTBI describes individuals harboring the bacteria without active disease, signifying potential for disease progression.

Calculation of signature score

The Seurat tool’s AddModuleScore function [49] was employed to calculate the signature score for a list of genes, evaluating and integrating module scores for samples based on the expression levels within predefined gene sets.

Cell-cell interactions analysis

CellChat [50] (version 1.6.1) was employed to model intra- and intercellular communication networks based on the abundance of known ligand-receptor (L-R) transcript pairs. To quantify interactions between cell types, a series of functions, including identifyOver ExpressedGenes, identifyOverExpressedInteractions, projectData, computeCommunProb, filterCommunication, aggregateNet, and netAnalysis_computeCentrality, were sequentially applied with default parameters. The mergeCellChat function was utilized to compare the intensity differences between Infected and Bystander, as well as Infected and Uninfected states.

Random forest modeling for TB diagnostics

To assess the predictive capability of IL1B, G0S2, PTAFR, and CSF3R for different diagnostic types among TB patients, including Control, Subclinical, EPTB, PTB, and LTBI, a random forest model was employed with the following steps: (1) The 31 TB patient queues were split into training and testing sets in a 3:1 ratio. It is important to note that patient samples lacking the aforementioned five labels were excluded. (2) IL1B, G0S2, PTAFR, and CSF3R were used as features, and the corresponding diagnostic types of patient samples served as labels. (3) The randomForest function from the “randomForest” package (version 4.7-1.1) [51] was utilized to train a random forest multi-classification model, with the formula: TBStatus ~ IL1B + G0S2 + PTAFR + CSF3R. All parameters were set to their default values. (4) The performance of the random forest model was evaluated using the testing set, and an Area Under the Receiver Operating Characteristic Curve (AUC) curve was plotted.

The strategy for constructing a random forest multi-classification model for the 11-feature gene set, composed of IL1B, G0S2, PTAFR, CSF3R, and their highly interacting genes, remained consistent with the aforementioned approach. It is worth noting that to construct a model effectively predicting the progression stages of TB patients over time, one only needs to replace the labels of different diagnostic types with labels representing different time stages since the onset of TB.

Construction TBSig for predicting the progress of TB patients

To construct a signature related to NETs for predicting disease progression in TB patients, 11 key genes (including IL1R2, MYD88, PTAFR, CSF3R, IL1RAP, IL1R1, CASP1, IL1B, G0S2, IFNG, IL4) inferred by this study are crucial for this purpose. The specific methodology is outlined as follows: (1) Utilizing the “ssgsea” method in GSVA [52] tool to estimate the activity of the gene set comprising these 11 genes in samples from TB patients. (2) Considering that the activity scores in TB patients approximately follow a Gaussian distribution, applying the Maximum Likelihood Estimation (MLE) strategy to estimate the corresponding mean and standard deviation of this distribution. (3) Given the fitted normal distribution, the TBSig score was considered as a continuous score variable X with PDF f_X(x), and the cumulative distribution function (CDF) can be obtained from:

$$TBSig=P(X\le x)=F_x(x)={\displaystyle {\int }_0^x{f}_x(u)du}$$

According to the above formula, the TBSig score is ranged between 0 and 1.

Heterogeneity analysis in TB patients

The variability of TBSig scores across different diagnostic types in TB patients was assessed using the Gini coefficient and the interquartile range (IQR). The Gini coefficient, a statistical measure of economic inequality within a population, was computed using the Gini function from the “DescTools” R package (version 0.99.52) [53]. This coefficient offered insights into the distribution of TBSig scores among various diagnostic categories. Additionally, the IQR, calculated as the difference between the third quartile (Q3) and the first quartile (Q1), served as another measure to evaluate the spread of TBSig scores within the TB patient population.

Inference of PTB patient subtypes

Given the higher heterogeneity observed in PTB patients compared to other types of patients, including Control, Subclinical, EPTB, and LTBI, it suggests the presence of potential subtypes within this patient category. To address this, the Seurat [49] standard protocol was employed, utilizing the FindClusters function to achieve clustering of PTB patients. The parameters used were “algorithm = original Louvain algorithm; resolution = 0.5.” This resulted in the identification of three distinct subgroups, defined as C1 (n = 478), C2 (n = 408), and C3 (n = 501).

To validate the stability of unsupervised clustering, consensus clustering was performed using the “ConsensusClusterPlus” package (version 1.62.0) [54] with the parameters “maxK = 6, clusterAlg = pam, distance = maximum, pItem = 0.90”. The optimal value of k was determined using the highest and the elbow point of the Delta area plot.

Differentially expressed genes (DEGs) associated with the PTB subtypes

The identification of DEGs among subtypes involved the application of the FindAllMarkers function from the “Seurat” package (version 4.3.0) [55]. Notably, genes meeting the criteria of adjusted p-values (adjusted by Benjamini-Hochberg method) below 0.01 were considered as significantly differentially expressed.

Functional enrichment analysis

To scrutinize distinct biological processes and pathways exhibiting significant differential expression among PTB subtypes, the “clusterProfiler” package (version 4.0.5) [56] was employed. Gene symbols were converted into Entrez IDs using the bitr function for subsequent analysis. Following this, the enrichGO function with the “ont = BP” parameter and the enrichKEGG function were utilized for gene ontology (GO) and KEGG pathway enrichment analysis, respectively. Significantly enriched GO and KEGG terms were identified based on adjusted p-values below 0.05 (adjusted by Benjamini-Hochberg method).

Estimation of cellular compositions in PTB subtypes

Using CIBERSORT [57], the infiltration of immune cell types in different subtypes of PTB was estimated, with the LM22 expression profile as a reference containing 22 cell types and default parameters. Additionally, xCell (version 1.1.0) [58] was employed to assess the abundance of 67 cell types, also with default parameters.

Statistical analysis

Pearson correlations were employed to assess the concordance between groups. Standard statistical tests, including Student’s t-test, Wilcoxon rank-sum test, and Kruskal-Wallis test, were applied to analyze both clinical and expression data in the study. These statistical analyses were conducted using R version 4.3.2.

Availability of data and materials

Transcriptomic data and relevant clinical information for the TB cohorts were extracted from publicly accessible databases, as elaborated in Supplementary Table 1. LL authored and reviewed all the codes, which can be obtained upon request from the corresponding author.

Uncovering the intricacies of cellular changes in TB through the analysis of single-cell data and NET-related genes

We initially acquired single-cell data pertaining to Mtb infection in mice from the research conducted by Pisu et al. [21]. This dataset comprises 17,101 cells and includes comprehensive single-cell clustering and annotation information provided by the research team (Figure 1A; Supplementary Table 1). The dataset delineates three distinct states in mice: Uninfected, Infected, and Bystander (Figure 1B). To gain a nuanced understanding of the alterations in diverse cell types within the murine lung microenvironment prompted by Mtb infection, we used Bystander as a reference and conducted a comparative analysis of cell-type abundance between the Infected and Uninfected states. The findings underscore a pronounced transformation in the Neutrophils cell type during the Infected state (Figure 1C). To further substantiate the dynamic changes in Neutrophils across various states of TB, we conducted a thorough analysis of cell-cell communication (see Materials and Methods). The observed differences in cell-type interactions between states revealed an augmented interaction between Neutrophils and Macrophages/T-cells during the Infected state compared to the Uninfected and Bystander states (Figure 1D). This heightened interaction implies a more robust immune response when addressing Mtb infection, reinforcing the pivotal role of Neutrophils cells in the progression of TB.

Figure 1. Interplay of cell types and neutrophil extracellular traps (NETs) in the single-cell landscape of Mtb infected mouse lung tissue. (A) Uniform Manifold Approximation and Projection (UMAP) showing the single-cell atlas of Mtb-infected mouse lung tissue [21]. Each point represents an individual cell, with distinct cell types marked by color codes. (B) UMAP showing the cell distribution based on uninfected, infected, and bystander statuses in mice. (C) Bar plot showing the change ratio of cell type fractions related to bystander mice. Notably, neutrophils exhibit the highest changes compared to uninfected mice and other cell types. (D) Heatmap showcasing ligand-receptor interactions within the mouse lung microenvironment. Greater inclination towards yellow indicates higher interactions among different cell types in the Infected group compared to bystander or uninfected status. (E) Violin plot coupled with a boxplot, revealing the signature scores of a gene set composed of neutrophil extracellular trap-related genes (NETs) in various cell types. (F, G) Violin plots showing the signature scores of NETs in Neutrophils (F) and the entire single-cell dataset (G) across different Mtb-infected statuses. (H) Violin plot showing the expression levels of Il1b, G0s2, Ptafr, and Csf3r between infected and Others (i.e., bystander and uninfected) statuses. P-values were obtained by t-test. (I) Violin plots showing the expression levels of Il1b, G0s2, Ptafr, and Csf3r across different cell types.

Several recent studies [16, 17, 59, 60] posit that Neutrophil Extracellular Traps (NETs), extracellular structures released by neutrophils, wield a significant influence in immune defense and inflammation regulation. We curated a set of 69 genes [16] associated with NETs and evaluated their activity across different cell types within the lung tissue microenvironment during Mtb infection in mice, highlighting the most robust activity scores within Neutrophils (Figure 1E; Supplementary Table 2). To delve into the dynamic fluctuations of NETs in distinct states of TB, we conducted a comparative analysis of NETs signature scores among the Uninfected, Infected, and Bystander states. The discerned pattern revealed the highest scores in the Infected state and the lowest in the Uninfected state, a consistent observation across the entire dataset (Figure 1F, 1G). The escalated NETs scores in the Infected state may be ascribed to the activation of neutrophils, releasing NETs to impede bacterial dissemination during Mtb infection.

Particularly, in the analysis discerning which NETs-related genes are closely associated with the progression of TB during the Infected state in comparison to the Uninfected and Bystander states, we conducted a differential expression analysis in Neutrophils cells. The results highlighted significant alterations in genes such as IL1B, G0S2, PTAFR, and CSF3R, manifesting specific high expression levels in Neutrophils cells (Figure 1H, 1I; t-test; p < 0.05; Supplementary Table 2). These genes, encompassing IL1B, G0S2, PTAFR, and CSF3R, exert pivotal roles in TB, contributing to immune responses and inflammation regulation [61, 62]. Their intricate interplay underscores a complex yet critical impact on the pathogenesis of TB. A comprehensive exploration of these genes holds promise in revealing novel therapeutic avenues and immune modulation strategies.

Significance of IL1B, G0S2, PTAFR, and CSF3R in the progression of TB patients

Although the single-cell data from murine Mtb provides insights into the putative roles of IL1B, G0S2, PTAFR, and CSF3R in the progression of TB, a more extensive investigation is imperative to elucidate their roles in TB patients. Consequently, we acquired data from 31 human TB datasets (see Supplementary Table 1; see Materials and Methods), comprising a total of 4229 samples, and applied IL1B, G0S2, PTAFR, and CSF3R to score these samples. The findings reveal consistent high expression scores for these four genes in the PTB state across diverse stages of patients (23/31) (Figure 2A). Notably, at the expression level, these four NETs-related genes exhibit significantly elevated expression levels in PTB compared to other stages, aligning with observations from murine single-cell data (Figures 2B and 1F).

Figure 2. Multi-faceted analysis of NETs signature scores and gene expression patterns in TB patients. (A) Boxplot showing the distribution of signature scores for NETs, composed of IL1B, G0S2, PTAFR, and CSF3R, across different statuses in 31 publicly available TB datasets (with a total of 4209 samples). Red points above the box indicate that the scores in pulmonary tuberculosis (PTB) are higher than in other statuses, while green points below the box suggest lower scores in PTB compared to other statuses. (B) Heatmap depicting the expression levels of IL1B, G0S2, PTAFR, and CSF3R in patients from 31 TB datasets. Column annotations under “Source” represent the dataset to which the sample belongs, while column annotations under “Treatment” indicate the clinical status of the patients. (C) Violin plots showing the distribution of signature scores based on demographic features, including age status (i), gender (ii), smoking status (iii), and ethnicity status (iv). P-values were determined using the Wilcoxon or Kruskal-Wallis tests. (D) Violin plots showing the distribution of signature scores in TB patients with complications, including diabetes (i), and lung cancer (ii). P-values were obtained through the Wilcoxon test. (E) Violin plots showing the distribution of signature scores in TB patients with bacterial infection (i, ii) or autoimmune disease (ii). P-values were obtained through the Wilcoxon test or Kruskal-Wallis test. (F) Violin plots showing the distribution of signature scores in different health states of patients with pneumonia and sarcoidosis. P-values were determined using the Wilcoxon test or Kruskal-Wallis test. (G) Violin plots showing the distribution of signature scores at different stages of TB progression. P-values were obtained through the Wilcoxon test or Kruskal-Wallis test. (H) Heatmap presenting the co-expression relationships between the gene expression of IL1B, G0S2, PTAFR, and CSF3R, and inflammation-related genes. (I) Heatmap showing the co-expression relationships between the gene expression of IL1B, G0S2, PTAFR, and CSF3R, and genes associated with T cell exhaustion and reactive oxygen species (ROS) pathway.

Further investigation was conducted to examine the correlation of these four genes with clinical features, encompassing demographic characteristics, complications, and distinct stages of TB progression. In terms of demographic features, the signature scores associated with IL1B, G0S2, PTAFR, and CSF3R exhibited a notable increase in older patients compared to their younger counterparts, a trend also observed in female patients compared to males, and in smokers relative to non-smokers (Figure 2C). Among different ethnic cohorts, African Americans demonstrated the lowest scores, while individuals from the Indian subcontinent exhibited the highest scores (Figure 2C). Within the spectrum of complications, patients with diabetes and lung cancer displayed notably higher scores compared to their counterparts without these conditions (Figure 2D). Patients with Still’s complications exhibited the highest signature scores, whereas Strep was associated with the lowest scores (Figure 2E). In the context of pneumonia and sarcoidosis complications, positive patients consistently manifested higher signature scores than their negative counterparts (Figure 2F). Regarding the diverse stages of TB progression, those with positive progression exhibited higher scores than those with negative progression, and patients further removed from the onset of tuberculosis displayed correspondingly lower signature scores derived from these four genes (Figure 2G).

Considering the pivotal role of the inflammatory response pathway in the development of TB, the activation of inflammatory pathways such as NF-κB, IL-1, TNF-α, IL-6 upon Mtb infection in lung tissue leads to the activation and recruitment of immune cells, enhancing their antimycobacterial capabilities [1, 63]. The analysis of genes associated with the inflammatory response, including IL1B, G0S2, PTAFR, and CSF3R, underscored a significant overall correlation, signifying their synergistic role in regulating immune system inflammation (Figure 2H). During TB progression, the ROS pathway assumes a pivotal role [64, 65]. IL1B, G0S2, PTAFR, and CSF3R exhibited a robust positive correlation with genes linked to reactive oxygen species pathway (ROS) (average r > 0.5) (Figure 2I). Notably, these four genes also demonstrated a positive association with genes associated with exhaustion, mirroring the intricate integration of immune and inflammatory regulation within the organism (Figure 2I). Collectively, these findings suggest their joint participation in modulating immune response and depletion processes, ensuring an effective response to pathogens, thereby unveiling the intricate nature of the inflammatory regulation network.

Development of TBSig as a diagnostic biomarker for predicting the progressions of TB patients

We observed marked heterogeneity in the distribution of signature scores derived from IL1B, G0S2, PTAFR, and CSF3 among distinct diagnostic categories of TB patients (Figure 3A). Specifically, PTB demonstrated the highest average signature scores, exhibiting significant distinctions in comparison to both the Control and LTBI groups. However, when compared to the subclinical and EPTB groups, no statistically significant differences were observed. At the gene expression level, these four genes consistently demonstrated higher expression in PTB compared to other TB types (Figure 3B). This led us to question the potential of IL1B, G0S2, PTAFR, and CSF3 as effective discriminators for various TB patient types.

Figure 3. Unveiling TB progression insights with emphasis on TBSig score through NETs-related genes. (A) Violin plots showing the distribution of signature scores (IL1B, G0S2, PTAFR, and CSF3R) across various stages of progression in pulmonary tuberculosis (PTB). The stages include Control, Latent Tuberculosis Infection (LTBI), Subclinical, PTB, and Extrapulmonary Tuberculosis (EPTB). (B) Violin plots showing the expression distribution of IL1B, G0S2, PTAFR, and CSF3R in different stages of progression in PTB patients. (C) The Area Under the Curve (AUC) curves showcase the predictive performance of a random forest classifier built with IL1B, G0S2, PTAFR, and CSF3R as features, across different diagnostic type of TB patients. (D) Network diagram showing gene interactions from the STRING database associated with IL1B, G0S2, PTAFR, and CSF3R. Different lines represent the sources of gene interaction relationships. (E) Heatmap showing the expression correlations of 11 highly interacting genes. Colors trending towards yellow indicate stronger correlations. (F, G) AUC curves showing the predictive performance of a random forest classifier constructed with 11 genes as features for different diagnostic types of TB (F) and various time stages leading up to the occurrence of TB (G). (H) Density curve plots showing the distribution of predicted scores constructed from 11 genes (left side) and the fitted error distribution (right side). (I) Violin combined with box plots showing the distribution of TBSig scores across different stages of progression in TB patients. (J) Violin combined with box plots showing the distribution of TBSig scores at different time stages leading up to the occurrence of TB patients. (K) Violin combined with box plots showing the distribution of TBSig scores under different states after treatment. (L) Violin combined with box plots showing the variability in the average scores of 4-genes, 11-genes, and TBSig across diverse datasets. For each iteration, 500 samples were randomly subsampled from the dataset, and this process was repeated 100 times. TBSig consistently demonstrated the least perturbation variability across all datasets.

To explore this, we constructed a random forest multi-classification diagnostic model, utilizing these four genes as features and clinical diagnostic types as labels. In the test set, the model exhibited robust diagnostic accuracy, reaching 0.983 for EPTB, 0.938 for Subclinical, and maintaining a commendable level above 0.8 for PTB (Figure 3C; see Materials and Methods). This underscores the significance of IL1B, G0S2, PTAFR, and CSF3 as pivotal biomarkers for diverse TB classifications. In an effort to further refine the model’s predictive performance for PTB, we endeavored to incorporate additional genes closely interacting with IL1B, G0S2, PTAFR, and CSF3. Retrieving an additional 7 genes highly interactive with these four from the STRING database, the resulting interaction network graph revealed credible relationships among them (Figure 3D; see Materials and Methods). Co-expression pattern analysis across 31 TB datasets demonstrated predominantly strong positive interactions between most of these genes (Figure 3E). Reintegrating these 11 genes into the random forest model notably improved diagnostic performance for various TB types, notably achieving an AUC of 0.871 for PTB and an impressive 0.995 for EPTB in the test set (Figure 3F). Moreover, our exploration of time since PTB onset revealed the highest diagnostic AUC at 3 months (0.973), followed by baseline (0.925), and extending beyond a year (0.916). The diagnostic AUC values for other time points consistently remained above 0.72 (Figure 3G).

Subsequently, in the pursuit of a robust biomarker for TB diagnosis, we introduced a cumulative Gaussian probability distribution strategy based on the enrichment score inferred by the “ssgsea” method, resulting in the creation of a diagnostic biomarker named TBSig (Figure 3H; see Materials and Methods). TBSig scores exhibited a discernible trend of decreasing severity across different TB diagnoses and exposure times (Figure 3I, 3J). Notably, TBSig displayed sensitivity in distinguishing patients with varied treatment outcomes, with the highest scores in Not Cured patients and the lowest in Definite Cure patients, aligning seamlessly with clinical expectations (Figure 3K). Crucially, TBSig demonstrated remarkable comparability across diverse patient cohorts, a pivotal characteristic for its practical application in clinical settings. Subsampling analysis of 31 datasets for 4-gene, 11-gene, and TBSig scores consistently revealed TBSig’s minimal variability in average scores among different datasets, underscoring its suitability for cross-dataset comparisons (Figure 3L; see Materials and Methods).

In summary, IL1B, G0S2, PTAFR, CSF3, and their interacting genes showcased significant differences among TB patients, offering robust discriminative features for diverse TB clinical types. The integration of interacting genes and the creation of the TBSig biomarker not only enhanced the accuracy and predictive performance of TB patient diagnosis but also presented a potential valuable tool for practical implementation in clinical settings.

Distinct clinical subtypes identification in pulmonary tuberculosis patients

We conducted an analysis using the Gini index and Interquartile range (IQR) to explore the inherent heterogeneity of TBSig among different types of TB patients (see Materials and Methods). The findings revealed that the Gini index and IQR of PTB patients were consistently higher than those of other types, suggesting the presence of distinct clinical subtypes within the PTB cohort (Figure 4A). Employing the Louvain clustering algorithm, three subtypes of PTB patients, named C1, C2, and C3, were successfully identified (Figure 4B; see Materials and Methods; Supplementary Table 3). Additionally, consensus clustering was applied to subgroup the PTB patients, confirming that the optimal number of clusters is 3 (Supplementary Figures 1 and 2; see Materials and Methods). The overlap analysis revealed a high degree of agreement between the consensus clustering and Louvain clustering results, indicating that the clusters identified by the Louvain method are stable (Supplementary Figure 3). Moreover, the projection of TBSig scores onto UMAP revealed a gradual increase from C3 to C1, a pattern corroborated by Figure 4D, where C3 exhibited the lowest, and C1 had the highest (Figure 4C, 4D). When assessing different subtypes of PTB patients, we utilized an inflammation-related gene set for scoring, revealing that the C1 subtype exhibited the highest score, while C3 had the lowest (Figure 4E; t-test). This suggests an inclination towards early inflammatory responses in the C1 subtype. Examination of the distribution of the three PTB subtypes at various time points in the progression of the disease illustrated that C1 had a higher distribution at Baseline and Month1, gradually declining over time, whereas C3 demonstrated the opposite trend (Figure 4F). The Sankey diagram further confirmed the predominance of C1 at 0 months and 6 months, with minimal distribution at 18 months (Figure 4G). In contrast, C2 and C3 exhibited a notable number of patients at 18 months, reinforcing the association of C1 with early-stage PTB.

Figure 4. Exploration of inferred PTB subtypes based on NETs-related genes. (A) Bar plots showing the heterogeneity indexes of the diagnostic types of TB patients based on Gini index and Interquartile range. (B) UMAP showing the distribution of different PTB subtypes, with each point representing an individual sample, and distinct colors indicating PTB categories. (C) Projection of TBSig scores for PTB patients onto UMAP. Each dot represents an individual PTB patient. (D) Violin plots combined with box plots showing the distribution of TBSig scores within different PTB subtypes. (E) Violin plots showing the activity scores of inflammation-related gene sets across different subtypes, with p-values obtained from t-tests. (F) Proportions of different PTB subtypes at various time stages relative to the occurrence of TB. (G) Sankey diagram showing the mapping relationships between different PTB subtypes and stages of PTB outbreaks. (H) Box plots showing the proportions of distinct cell types infiltrating various PTB subtypes, as estimated by CIBERSORT. ^nsp > 0.05; ^*p < 0.01; ^**p < 0.001; ^***p < 0.0001; ^****p < 0.00001. (I) Bar plot showing the Pearson correlations between Neutrophils with other immune cell types. (J) Scatter plot with a fitted line showing the association between Neutrophils and CD8+ T cells. Each dot represents an individual PTB patient. “R” indicates the Pearson correlation, and the p-value is obtained by t-test. (K) Bar plot showing the Pearson correlations between Neutrophils and other cell types, as estimated by xCell. (L) Box plots showing the distribution of Neutrophils and CD8+ T cells at different time stages relative to the onset of PTB outbreak. (M) Box plot showing the difference of CD8+ T cell fractions between low and high Neutrophil infiltration. P-value was obtained by t-test.

To assess differences in the immune microenvironment among the three PTB subtypes, we utilized CIBERSORT [57] to estimate immune cell composition. The results indicated significant differences in the majority of immune cell types among the three subtypes (18/22; p < 0.05) (Figure 4H; see Materials and Methods; Supplementary Table 4). Monocytes, Neutrophils, and T cells CD8+ dominated the proportions, with Neutrophils gradually decreasing from C1 to C3, while T cells CD8+ increased gradually (Figure 4H; t-test). Further analysis revealed a robust negative correlation between Neutrophils and CD8+ T cells (r = −0.61; t-test, p < 0.00001; Figure 4I, 4J). Validation through the xCell tool [58] confirmed a substantial negative correlation between CD8+ T cells and Neutrophils (Figure 4K; see Materials and Methods; Supplementary Table 4). Based on these findings, we postulated that early PTB patients might manifest a higher abundance of Neutrophils, potentially inhibiting the activation of CD8+ T cells. Subsequent analysis demonstrated that Neutrophil abundance was elevated in the early stages, diminishing over time, whereas CD8+ T cell abundance followed the opposite trajectory, commencing at lower levels in the early stages and gradually increasing (Figure 4L). Stratifying PTB patients into high Neutrophils and low Neutrophils groups based on average abundance, the results revealed that Low Neutrophils patients had higher CD8+ T cell levels, while High Neutrophils patients had lower CD8+ T cell levels, with a significant distinction between the two groups (Figure 4M).

To summarize, our analysis unveiled internal heterogeneity and subtype differences within the TB patients. Specifically, the C1 subtype displayed a heightened inflammatory response in the early stages, linked with increased Neutrophil abundance and suppression of CD8+ T cells.

CLEC4D’s potential significance in suppressing CD8+ T cell activity by early neutrophils during TB

Through differential expression analysis, we observed significant variations in gene expression across different subtypes of PTB patients (Figure 5A; see Materials and Methods). Specifically, C1 exhibited a higher number of significantly differentially expressed genes compared to C2 and C3, while C2 displayed the fewest differentially expressed genes (Figure 5A; Supplementary Table 5). Additionally, there was an overlap in differentially expressed genes between C2 and both C1 and C3, indicating that C2 might be in an intermediate state during the progression of PTB from C3 to C1, providing insights into crucial nodes in the development of the disease. Further exploration of the biological functions and pathways involving the top 100 differentially expressed genes in different subtypes revealed that C1 was associated with biological functions related to inflammation, including pathways involving NK-kappaB and cytokines (Figure 5B, 5C; see Materials and Methods; Supplementary Tables 6 and 7). This suggests that C1 may play a vital role in immune responses. Conversely, C2 was predominantly involved in the electron transport chain and oxidative phosphorylation pathways, emphasizing its critical role in cellular energy production (Figure 5B, 5C; Supplementary Tables 6 and 7). On the other hand, C3 exhibited pathways related to the regulation and activation of T cells, highlighting its significance in modulating and activating T cells within the immune system (Figure 5B, 5C; Supplementary Tables 6 and 7). These findings unveil the involvement of specific functional pathways in each subtype, offering a new perspective for a deeper understanding of the progression of PTB.

Figure 5. Functional characteristics and dynamic gene expression patterns in PTB subtypes. (A) Venn diagram showing the overlap of differentially expressed genes among different PTB subtypes. (B, C) The most significant 100 genes expressed in different PTB subtypes enrich in biological processes (B) and KEGG pathways (C). (D) Violin plot showing the distribution of differentially expressed genes in C1 that have a correlation with TBSig greater than 0.4 across different PTB subtypes. (E) Bar plot showing the correlation between gene expression and the abundance of CD8+ T cells and Neutrophils (estimated by CIBERSORT). The right side corresponds to Neutrophils, while the left side corresponds to CD8+ T cells. (F) Scatter plot showing the association between the expression of GK and CLEC4D genes and TBSig score. Each point represents an individual PTB patient, and the blue line is the fitted linear regression line, with “R” indicating the Pearson correlation coefficient. (G) Violin plot showing the expression of GK and CLEC4D in different cell types in a mouse single-cell dataset (GSE167232) [21]. (H) Box plot shows the dynamic changes in CLEC4D expression over different time points of PTB occurrence. (I) Box plot displays the difference in CD8+ T cells abundance between the high CLEC4D expression group and the control group in all PTB datasets. The p-value is obtained from a t-test. (J) Box plot illustrates the difference in CD8+ T cells abundance between the high CLEC4D expression group and the low group in different TB datasets. The p-values were obtained from a t-test.

In the investigation of key genes within the C1 subtype during the progression of PTB, a correlation analysis with TBSig identified 13 genes with correlation coefficients greater than 0.4. These genes displayed a clear trend from high to low expression levels from C1 to C3, indicating their potential crucial roles in the disease progression of TB patients (Figure 5D). Given the potential suppressive effect of early Neutrophils on CD8+ T cell activity, we analyzed the correlation between these genes and the abundance of Neutrophils and CD8+ T cells. Results revealed a strong positive correlation between GK and CLEC4D with Neutrophil abundance (r > 0.5), and a strong negative correlation with CD8+ T cell abundance (Figure 5E, 5F; r < −0.5). Single-cell data analysis further demonstrated the specific expression of CLEC4D in Neutrophils (Figure 5G). Combining these findings with existing literature, it is suggested that CLEC4D may participate in the early immune response to PTB, with its expression trend aligning with that of Neutrophils.

To validate the role of CLEC4D, we analyzed its dynamic changes during different stages of PTB progression. CLEC4D exhibited higher expression levels in the early stages, gradually decreasing over time, consistent with the trend observed in Neutrophils (Figure 5H). Stratifying PTB patients into high and low CLEC4D expression groups revealed significant differences in CD8+ T cell abundance between the two groups (Figure 5I). Further analysis across multiple datasets confirmed this trend, reinforcing the potential role of CLEC4D in regulating CD8+ T cell abundance (Figure 5J; t-test). Overall, these findings provide crucial clues for a comprehensive understanding of the molecular mechanisms underlying the progression of PTB.