The experimental design is illustrated in Figure 1 .

Two microarray datasets (GSE83456 and GSE152532) were obtained from the Gene Expression Omnibus (GEO) database (www.ncbi.nlm.nih.gov/geo). Perl programming language was employed to preprocess the two microarray datasets. Microarray datasets were merged with probe platform.The GSE83456 dataset (GPL10058 platform, PMID: 27706152), including 61 blood samples from healthy people (control group) and 92 blood samples from TB patients (TB group), was selected for preliminary analysis. Diagnosis of TB was based on the following: positive mycobacterial culture result from the site of disease; or caseating granuloma on biopsy and/or clinical/radiological features consistent with TB and a good response to therapy (19). The GSE152532 dataset (GPL10058 platform, PMID: 34555657), including 11 blood samples from healthy people and 136 blood samples from TB patients, was selected for validation. The blood samples from TB patients in GSE152532 dataset contained 111 latent TB samples and 25 active TB samples. Latent TB was confirmed in participants with a positive interferon gamma release assay(IGRA) and the absence of clinical and radiographic signs of ATB or any other significant co-morbidity. Active TB was defined as either a suspected (based on clinical suspicion and radiological and/or histological evidence) or microbiologically confirmed (by Mtb culture) new diagnosis of pulmonary or extrapulmonary TB disease, in the absence of any other significant co-morbidity (20).The cuproptosis-related genes in the Molecular Signature Database (MsigDB) v7.0 database (http://www.gsea-msigdb.org/gsea/msigdb/) were combined with gene sets relevant to cuproptosis from a prior study (21). R programming language (version 4.1.3) was employed to conduct this study.

The R package “limma” was utilized to identify differentially expressed CRGs (threshold P<0.05) from the GSE83456 dataset. Subsequently, the R packages “ggpubr” and “pheatmap” were utilized to generate box plots and heat maps, respectively. The differential CRGs were correlated by the R package “corrplot “ to explore correlations between genes.

The common method for studying cell heterogeneity, such as flow cytometry, rely on a limited repertoire of phenotypic markers. However, tissue disaggregation before flow cytometry can lead to lost or damaged cells, altering results (22). Aaron M Newman et al. presented cell-type identification by estimating relative subsets of RNA transcripts (CIBERSORT), a computational approach that accurately resolves relative fractions of diverse cell subsets in GEPs from complex tissues (23). CIBERSORT (https://CIBERSORT.stanford.edu//) is an analytical tool from the Alizadeh Lab and Newman Lab to impute gene expression profiles and provide an estimation of the abundances of member cell types in a mixed cell population, using gene expression data. Based on the gene expression data of the GSE83456 dataset, the CIBERSORT algorithm and Leukocyte signature matrix (LM22) were used to estimate the relative abundance of 22 immune cells in each sample. LM22 is a gene matrix that contains 547 white blood cells characteristic genes to differentiate 22 types of immune cells, including myeloid subgroup, Natural killer (NK) cells, naive and memory B-cells, and seven types of T-cell (24). Monte Carlo sampling was used to obtain the inverse folded product p-value of each sample. Only the samples with p-value < 0.05 were considered accurate immune cell components. The sum of the 22 immune cells in each sample was 1.

The correlation coefficient between CRGs and the characteristics of relevant immune cells was analyzed to further demonstrate the correlation between CRGs expression and the relative percentage of immune cells. Spearman correlation analysis was conducted, and a p-value less than 0.05 indicated a significant correlation. Finally, the results were visualized using the R package “corrplot”.

Based on the expression profile of CRGs, the R package “ConsensusClusterPlus” was utilized to apply the unsupervised clustering analysis. The 92 TB samples were classified into clusters using the k-means algorithm with 1,000 iterations. k values are defined from 1 to 9 to generate different subtypes. The optimal number of clusters was selected according to the cluster consensus score. The principal component analysis (PCA) was used to visualize the distribution of the clusters.

GSVA contributes to the current need of Gene set enrichment methods for RNA-seq data (25). GSVA is an open source software package for R which forms part of the Bioconductor project and can be downloaded at http://www.bioconductor.org. The R package “GSVA” was used for GSVA enrichment analysis to elucidate the differences of enriched gene sets among different CRGs clusters. The gene matrix transposed (gmt) file,including “ c2.cp.kegg.v2022.1.Hs.symbols.gmt” and “ c5.all.v2022.1.Hs.symbols.gmt”, were downloaded from the MSigDB website database for further GSVA analysis. The R package “limma” was utilized to identify differentially expressed pathways and biological functions by comparing GSVA scores between different CRGs clusters. A p-value < 0.05 was statistically significant.

Weighted correlation network analysis (WGCNA) can be used for finding clusters (modules) of highly correlated genes, for summarizing such clusters using the module eigengene or an intramodular hub gene, for relating modules to one another and to external sample traits (using eigengene network methodology), and for calculating module membership measures (26).The R packages “WGCNA” was utilized to identify the co-expression modules. The top 25% genes with the highest variance were used for WGCNA analysis to guarantee the accuracy of quality results. The optimal soft power was selected to construct the weighted adjacency matrix, which was further transformed into a topological overlap matrix (TOM). When the minimum module size was set to 100, the TOM dissimilarity measure (1-TOM) based on the hierarchical clustering tree algorithm was used to obtain the module. Each module was assigned a random color. Modular characteristic genes represent the overall gene expression profile of each module. Gene significance (GS) represented the correlation between genes and clinical phenotypes. The relationship between modules and disease status was reflected by module significance (MS). MS, defined as the mean of the gene significance values of all genes within a module, represents the correlation between module genes and traits

Based on two different CRGs clusters, the R package “caret” was applied for establishing machine learning models, including random forest (RF), support vector machine (SVM), generalized linear model (GLM), and eXtreme Gradient Boosting (XGB). The R package “pROC” was utilized to visualize the area below the ROC curve. The five genes with the lowest dropout-loss values were considered as top five variables. Therefore, the optimal machine learning model was determined, and the top five variables were identified as the key predictive genes associated with TB. ROC curve analysis was utilized to verify the diagnostic value of the diagnostic model in the GSE152532 dataset. Finally, based on the key predictive genes of TB, Spearman correlation analysis was performed in the GSE152532 dataset to explore the associations between prediction model-related genes with latent TB and active TB. A p-value < 0.05 was statistically significant.

The the top five variables were considered as the predictors. Based on the predictors, a nomogram was constructed to identify TB clusters by the R package “rms”. Each predictor was attributed a score, and the “total score” was obtained by summing the scores of all predictors. Calibration curve analysis and decision curve analysis (DCA) were applied to evaluate the predictive performance of the nomogram.

R Programming Language (version 4.1.3) was used for data analysis and statistical analyses. Bilateral wilcoxon tests were utilized to evaluate statistical differences between the two groups. Spearman correlation, calculated by “cor.test” function,was applied to analyze the relationship between the expression level of cuproptosis-related genes and immune cells. A p-value < 0.05 was statistically significant.

To clarify the biological function of CRGs in the occurrence and development of TB, the GSE83456 dataset was utilized to assess the differential expression of CRGs between TB and healthy people. A total of 11 CRGs were identified as differentially expressed cuproptosis genes ( Figure 2A ). The expressions of NFE2L2, NLRP3, ATP7B, SLC31A1, MTF1, and DLD were upregulated, while LIAS, LIPT1, DLAT, GLS, and DBT were downregulated in TB patients ( Figure 2B ). Subsequently, correlation analysis between these differentially expressed CRGs was performed to explore whether cuproptosis regulators participated in the progression of TB ( Figure 2C ). Among these genes, significant positive correlations were found between NFE2L2 and MTF1 (correlation coefficient, r= 0.56), LIPT1 and DLAT (correlation coefficient,r= 0.59), DLD and DLAT (correlation coefficient,r = 0.51), as well as GLS and DBT (correlation coefficient,r = 0.55) ( Figure 2D ).

An immune infiltration analysis was conducted to clarify immune differences between the TB and control groups. CIBERSORT analysis revealed differences in the abundance of 22 infiltrating immune cell types between the TB and control groups ( Figure 3A ). CD8+ T cells, resting memory CD4+ T cells, and follicular helper T cells were significantly decreased, while monocytes, M0, M1, and M2 macrophages, activated dendritic cells, eosinophils, and neutrophils were increased in the TB group ( Figure 3B ). Meanwhile, resting dendritic cells, eosinophils, M0, M1, and M2 macrophages, activated mast cells, resting mast cells, monocytes, neutrophils, activated NK cells, resting NK cells, plasma cells, activated memory CD4+ T cells, resting memory CD4+ T cells, naïve CD4+ T cells, CD8+ T cells, follicular helper T cells, gamma delta T cells, and T cells regulatory (Tregs) were associated with CRGs ( Figure 3C ). These results suggest that CRGs are key factors in regulating molecular and immune-invasive states in TB patients.

Based on the expression of 11 CRGs, 92 TB samples were grouped using a consensus clustering algorithm. The optimal number of clusters was observed when the k value was set to 2 (k = 2), and the CDF curve fluctuated within the minimum range of the consensus index of 0.2 to 0.6 ( Figures 4A, B ). At k = 2 ~ 9, the area under the CDF curve exhibited differences between two CDF curves (k and k-1) ( Figure 4C ). At k = 2, the concordance score of each subtype was the highest ( Figure 4D ). PCA analysis revealed that the 92 TB patients could be divided into Cluster 1 (n = 42) and Cluster 2 (n = 50), which were significantly different ( Figure 4E ).

The expression differences of 11 CRGs between Cluster 1 and Cluster 2 were comprehensively evaluated to explore the molecular characteristics between clusters. Distinct CRGs expression profiles were observed between Cluster 1 and Cluster 2 ( Figure 5A ). The expression of NFE2L2, NLRP3, SLC31A1, LIPT1, DLD, DLAT, MTF1, GLS, and DBT was significantly upregulated in Cluster 2( Figure 5B ). In addition, the immune cell infiltration analysis showed significant differences in the immune microenvironment between Cluster 1 and Cluster 2 ( Figure 5C ). The abundance of M0 macrophages, eosinophils, and neutrophils was significantly increased ( Figure 5D ).

The pathway activities and biological functions associated with each group were identified by GSVA. The functional enrichment results showed that lipopolysaccharide-mediated signaling, regulation of myeloid cell differentiation, and protoporphyrinogen IX metabolic process were significantly enriched in Cluster 1 ( Figure 6A ). In contrast, cerebellar Purkinje cell granule cell precursor cell signaling involved in the regulation of granule cell precursor cell proliferation, neuropeptide hormone activity, and positive regulation of gastrulation were significantly enriched in Cluster 2 ( Figure 6A ). KEGG pathway analysis suggested that Chronic myeloid leukemia, apoptosis, and neurotrophin signaling pathways were significantly enriched in Cluster 1 ( Figure 6B ), while olfactory transduction, taste transduction, maturity-onset diabetes of the young, retinol metabolism, and neuroactive ligand-receptor interaction were significantly enriched in Cluster 2 ( Figure 6B ).

The WGCNA algorithm was applied to establish co-expression networks and modules in the normal population and TB patients to identify key gene modules associated with TB. The variance of each gene expression in the GSE83456 dataset was calculated to select the top 25% genes with the highest variance for further analysis. Ten co-expression modules with different colors were obtained by dynamic cutting algorithm, and a heat map of the topological overlap matrix (TOM) was generated. Among them, 1,225 genes in the blue module exhibited the most significant relationship with TB ( Figure 7A ). In addition, the WGCNA algorithm was conducted to analyze the key gene modules closely related to CRGs clustering. Module-clinical features (Cluster 1 and Cluster 2) relationship analysis demonstrated a high correlation between the blue module (929 genes) and TB clusters ( Figure 7B ). Finally, the intersection of genes in the two modules yielded 154 genes by the R package “Venn” ( Figure 7C ).

Four machine learning models, including RF, SVM, GLM, and XGB, were constructed based on cluster-specific DEGs in the TB training cohort. The R package “DALEX” was applied to interpret the four models. The residual distribution of each model was plotted in the validation set. The residuals of XGB and SVM machine learning models were lower ( Figures 8A, B ). Subsequently, the genes of the top 15 features of each model were sequenced according to root mean square error (RMSE) ( Figure 8C ). Moreover, the discriminative performance of the four machine learning algorithms was evaluated by calculating receiver operating characteristic (ROC) curves based on 5-fold cross-validation in the training set (GSE83456 dataset) ( Figure 8D ). The areas under the ROC curve (AUC) were obtained for the four models (RF, AUC=0.975; SVM, AUC=0.979; XGB, AUC=0.957; GLM, AUC=0.716). Based on the residual and AUC, the XGB machine learning model demonstrated the best performance in distinguishing TB patients with different clusters. Finally, the top five genes (C13orf18, PACS1, MGC18216, RNASE2, and PLAUR) from the XGB model were selected as predictor genes for further analysis.

A nomogram was constructed to estimate the risk of cuproptosis clusters in 92 TB patients ( Figure 9A ). The prediction performance of the nomogram was evaluated by calibration curves and DCA. The predicted outcomes were consistent with the actual outcomes ( Figure 9B ). DCA indicated that the accuracy of the nomogram was relatively high, providing a net clinical benefit ( Figure 9C ).

The GSE152532 dataset was utilized to validate the accuracy of the machine-learning model. In the GSE152532 dataset, the ROC curve of the five genes of the XGB model (C13orf18, PACS1, MGC18216, RNASE2, and PLAUR) exhibited good performance (AUC= 0.825) ( Figure 10A ). The TB patients in GSE152532 dataset could be divided into latent TB and active TB groups. Based on the clinical characteristics, the five genes were utilized to predict latent and active TB ( Figures 10B-F ), C13orf18 (R=-0.18) and MGC18216 (R=-0.23) were negatively correlated with active TB ( Figures 10B, C ). RNASE2 (R=0.25) was positively correlated with active TB ( Figure 10F ).