Identification and validation of ubiquitination-related signature and subgroups in immune microenvironment of tuberculosis

Data collection

The GEO database (http://www.ncbi.nlm.nih.gov/geo) is a repository for public genome datasets. For this study, gene expression profiles were downloaded from GEO, specifically GSE157657, GSE62525, GSE83456, GSE93272, GSE76925, GSE166253, and GSE31210. The platform contains an explanation document that aims to link the probes with their respective genes.

Identification of DEGs and gene ontology enrichment analysis

To identify differentially expressed genes (DEGs) between the TB and HC, the R package “limma” was utilized [12]. The criteria of filtering were established as follows: P < 0.05 and |log2 fold change| > 0.264 or |log2 fold change | > 0.585. To undertake Gene Ontology (GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis for the DEGs, we utilized the DAVID tool (https://david.ncifcrf.gov/). Enriched GO pathways and KEGG pathways were determined based on the P-value cut-off criterion of < 0.05.

ssGSEA algorithm calculates the function score

We retrieved ubiquitination terms from the iuucd database (http://iuucd.biocuckoo.org/) to specifically select the relevant functional items for assessing disparities in ubiquitination-related functions within the samples. Next, we applied the single sample gene set enrichment analysis (ssGSEA) algorithm using the R package “GSVA” to score each piece. A higher score indicates a greater relative gene expression level associated with the specific functional item. Based on this scoring metric, we evaluated the activation state in relation to the particular functional item of interest.

Machine learning

To improve the regularity, interpretability, and predictive accuracy of predictive models and to select relevant variables for model incorporation, we utilized a Least Absolute Shrinkage and Selection Operator (LASSO) regression approach [13]. Support Vector Machine (SVM) method, a technique that establishes a threshold between categories for prediction based on one or more feature vectors, was implemented for analyzing the data [14]. Random forest (RF) approaches were also employed for predicting continuous variables accurately, precisely without significant fluctuations [15]. eXtreme Gradient Boosting (Xgboost) ensemble learning algorithm was also used, which utilizes decision trees as base learners [16]. These machine learning analyses (LASSO, SVM-RFE, Xgboost and RF) were performed using the R packages “glmnet,” “kernlab,” “randomForest,” and “xgboost,” respectively. The genes identified through the intersection of these analyses were considered to be hub genes related to ubiquitination in the diagnosis of TB.

Constructing the nomogram model

The carefully selected prediction variables (clinical relevance and utility) were utilized to build and validate the nomogram model [17]. Following the final multivariable logistic regression analysis, the nomogram was created in order to provide a visual instrument for evaluating individual risk projections for constipation in TB patients. The effectiveness of the prediction model was evaluated by calculating the area under the receiver operating characteristic curve. Furthermore, we assessed the clinical practicality of the prediction model by utilizing the decision curve analysis (DCA) method to determine net benefits.

Subclusters analysis with 11 Ub-related hub genes

The TB samples were subject to an unsupervised hierarchical clustering analysis using the “ConsensusClusterPlus” R package [18] based on the expression of 11 Ub-related genes. Subsequently, we employed GSVA [19] to elucidate the function between the Ub-related subclusters identified through the clustering mentioned above analysis.

Gene set variation analysis (GSVA)

To perform GSVA analysis, we downloaded files “h.all.v7.5.1.symbols”, “c2.cp.kegg.v7.5.1.symbols” and “c2.cp.reactome.v7.5.1.symbols” from MSigDB.

GSVA R package was used to calculate the activities of two sets for each sample. DEGs between the two Ub-related subclusters were identified using a significance cut-off of | log2 fold change (FC)| > 0.585 and P < 0.05. These DEGs were illustrated using a volcano plot. Subsequently, GO and KEGG enrichment analyses were carried out using the DAVID tool to uncover the biological functions of the identified DEGs.

Evaluation of immune cell infiltration

The CIBERSORT algorithm was utilized by us to establish the linked cellular immune infiltration with the help of standard LM22 gene signature. Marker information for 22 immune cells was obtained for the purpose of computing the immune cell infiltration in each dataset [20]. Spearman correlation analysis was performed to investigate the correlations between immune cells and hub genes, enabling the exploration of the relationship between these hub genes and immune cells [21].

Quality control and processing of single cell data

We obtained single cell raw sequencing data (10X Genomics) from previous study, which can be accessed under the accession numbers SRR11038989 and SRR11038994 [22]. The number of expressed genes was calculated for each single cell. Following preliminary quality control procedures, we utilized the R package “Seurat” to perform a standardized analysis workflow. Firstly, we standardized the data using “NormalizeData”, and then conducted principal component analysis via “RunPCA”. Finally, we utilized the k-nearest neighbor classification algorithm through “FindNeighbors” and “FindClusters” to classify single-cell samples into distinct clusters. We then scored each cluster based on the normalized expressions of canonical markers, which were listed in Supplementary Table 1 [23].

Statistical analysis

R software (version 4.2.2) was utilized to conduct all the statistical analysis. Figure panels were pieced together by Adobe Illustrator (CC 2019). Statistical significance was ascribed to any difference with a P-value or adjusted P-value falling below 0.05. Spearman’s correlation coefficient was employed to determine the association between continuous variables. Additionally, a t-test or Wilcoxon rank-sum test was performed for statistical comparison between two groups. ^*p < 0.05, ^**p < 0.01, ^***p < 0.001.

Data availability statement

The original contributions presented in the study are included in the article Supplementary Material. Further inquiries can be directed to the corresponding author.

The workflow of this study

The objective of this study was to identify ubiquitination-related biomarkers for TB. Data analysis was conducted using ssGSEA, and machine learning to determine the hub genes that significantly influence the TB response. Subsequently, validation and diagnostic value analyses were performed to assess the potential of hub genes. These findings offer novel perspectives for diagnosing and treating of TB, adding to the existing scientific knowledge in this field. In this work, all mRNA expression-related analyses were achieved by using the gene expression data from the GEO database. The workflow of this study was provided in Figure 1.

Figure 1. The detailed work process of this study.

Ub-related DEGs identified in TB

In total, based on the relevant functional items acquired from the IUUCD database, ssGSEA was conducted to investigate 6 Ub families (E1, E2, E3, DUB, UBD, ULD) in our study and the gene sets in each family were presented in Supplementary Table 2. Among these Ub families, the ssGSEA score of E1, E3, and UBD showed a significant difference in TB and HC in both GSE157657 (Figure 2A) and GSE62525 dataset (Figure 2B). Compared to HC, the abundance of Ub-related genes (E1, E3, and UBD) significantly decreased in TB.

Figure 2. Ub-related DEGs identified in TB. ssGSEA score of E1, E2, E3, DUB, UBD, ULD in GSE157657 (A) and GSE62525 (B) dataset. The volcano plotting of DEGs in GSE157657 (C) and GSE62525 (D) dataset. (E) The overlapping of genes between DEGs and Ub-related genes. DEGs, differentially expressed genes. ssGSEA, single sample gene set enrichment analysis. ^***p < 0.001.

Next, we analyzed GSE157657 and GSE62525 datasets and identified 3448 and 9530 DEGs (|log2 FC| > 0.264 and P < 0.05), respectively (Figure 2C, 2D). Of these, 1874 and 5049 genes were upregulated, while 4481 and 1574 genes were downregulated, respectively. To gain a better understanding of TB pathogenesis, we performed a cross-comparison of gene expression profiles. Ultimately, 96 Ub-related DEGs (UbDEGs) were identified as candidate genes for further analysis (Figure 2E).

Functional and pathway enrichment analysis of UbDEGs

GO and KEGG enrichment pathway analysis were performed on the 96 UbDEGs to investigate their potential roles. As expected, the GO analysis indicated significant enrichment of ubiquitination-related biological processes (such as protein ubiquitination, proteasome-mediated ubiquitin-dependent protein catabolic process, ubiquitin-dependent protein catabolic process). The UbDEGs were also enriched in immune-related biology processes (for example, regulation of immune system process, innate immune response) (Figure 3A). The KEGG analysis identified significant enrichment of the Ubiquitin-mediated proteolysis and Shigellosis signaling pathway among the DEGs. All these results further indicated that these UbDEGs were also closely related to immune function. Additionally, we observed enrichment in cellular compartments such as cytosol, cytoplasm, nucleoplasm, and ubiquitin ligase complex (Figure 3B). Furthermore, in terms of molecular function, identical protein binding, protein binding, carbohydrate binding and transmembrane signaling receptor activity were found to be enriched (Figure 3C).

Figure 3. Enrichment in GO analyses based on UbDEGs. Enrichment in GO BP (A), CC (B) and MF (C) analysis. Abbreviations: GO: Gene Ontology; BP: Biological process; CC: Cellular component; MF: Molecular function.

Identification of the Ub-related hub genes via machine learning

96 UbDEGs as candidate genes were applied for machine learning. We applied the LASSO and selected λ.min to identify 26 Ub-related genes as significant genes (Figure 4A). Then, through SVM with 10-fold cross-validation, 91 Ub-related genes were obtained (Figure 4B). Using Xgboost, and RF algorithms, we ranked the genes based on the importance and extracted the top 51 and top 30 (Figure 4C, 4D), respectively. Ultimately, we identified 11 Ub-related hub genes, which included WDFY1, FBXL15, ZBTB1, VHL, TRIM7, UTP15, EML5, TRIM68, CORO6, ZNF131, and UBA7 (Figure 4E). Moreover, most of the Ub-related hub genes were highly correlated with one another (Figure 4F).

Figure 4. Employment of machine learning techniques in identifying the Ub-related hub gene. (A–D) Construction of Ub-related genes using LASSO, SVM, Xgboost, and RF. (E) The overlapping of hub genes between the four machine learning algorithms illustrated by Venn diagram. (F) The relationship between the 11 Ub-related hub genes. Abbreviations: LASSO: least absolute shrinkage and selection operator; SVM: support vector machine; XGBoost: eXtreme Gradient Boosting; RF: random forest.

To predict the risk of TB, a nomogram prediction model was constructed (Figure 5A). We evaluated the nomogram’s performance using a calibration curve for TB risk prediction, which showed satisfactory agreement in this cohort (Figure 5B). The C-index of the nomogram was 0.962 (95% CI: 0.94338–0.98062), indicating good discrimination by the model. The nomogram demonstrated sound predictive capability based on its apparent performance. The decision curve analysis for the nonadherence nomogram was depicted in Figure 5C, which shows that the threshold probability range of 0.14-1.00 is clinically significant. Additionally, we analyzed the ROC curves of seven gene signatures and noted that the 11 Ub-related hub genes had the highest AUC value of 0.974 (Figure 5D), indicating their excellent diagnostic value. Similar results have also been observed in GSE62525 and GSE83456 (Supplementary Figure 1A, 1B). These findings revealed that all 11 hub genes exhibit outstanding diagnostic accuracy.

Figure 5. (A) Developed nonadherence nomogram. (B) The calibration curves for predicting nonadherence using the nomogram in the cohort are illustrated. (C) Decision curve analysis for the nomogram. (D) ROC curve of Ub-related hub genes in TB diagnosis.

Evaluation of immune cell infiltration

Given that UbDEGs were associated with immune-related biological processes (Figure 3A), we investigated immune cell infiltration to gain further insight into the immune regulation of TB. The boxplot (Figure 6A) showed that TB patients exhibited a higher proportion of monocytes and macrophages M0, but a lower proportion of naïve B cells. We then analyzed the correlation between the 11 Ub-related hub genes and immune cells (Figure 6B). It is worth noting that CD4+memory resting T cells and monocytes cells were significantly correlated with all core genes (P < 0.001). Our analysis revealed a significant correlation between these genes and immune cells, suggesting their pivotal role in the immune process of TB.

Figure 6. Immune cell enrichment in patients with TB. (A) The proportion of all 22 types of immune cells. (B) The relationship between 11Ub-related hub genes and 22 immune cells. (C) The location of the 11 Ub-related genes. (D) Dot plot showing the 11 Ub-related genes in each cell cluster. ^*p < 0.05, ^***p < 0.001.

To further verify the correlation between these hub genes and immune cells, we utilized scRNA-seq on PBMCs derived from public datasets to determine the locations of the 11 hub genes. Our findings indicated that WDFY1, ZBTB1, and ZNF131 were dominantly expressed in monocytes, T, NK, B, and dendritic cells, while the expression levels of FBXL15, TRIM7, UTP15, EML5, TRIM68, and CORO6 were relatively low in all subclusters (Figure 6C, 6D).

Consensus clustering analysis of Ub-related clusters and GSVA of biological pathways

Consensus clustering analysis was conducted using the “Consensus Cluster Plus” package in R software. Based on the expression of 11 Ub-related hub genes to investigate more thoroughly the interactions and consistency among DEGs, we established that k = 2 yielded the most consistent grouping (Figure 7A). After that, the 565 TB samples were divided into two distinct subclusters, cluster A (n = 163) and cluster B (n = 402). The expression levels of 11 hub genes in these two clusters were visualized via a boxplot (Figure 7B). Most Ub-related hub genes showed significantly different expression levels between clusters B and A, except for WDFY1.

Figure 7. Identification of Ub-related subclusters in TB samples and GSVA. (A) Subclusters were performed with 11 Ub-related hub genes. (B) Boxplot showing the expression level of 11 Ub-related hub genes in cluster A and B. Enrichment of pathways based on the KEGG (C), HALLMARK (D), and Reatcome pathway (E). GSVA, gene set variation analysis. ^**p < 0.01; ^***p < 0.001.

By utilizing GSVA analysis, we identified several pathways with differential expression. These pathways were shown by a heatmap (Figures 7C–7E). Compared with cluster B, the expression of KEGG pathways linked with Basal transcription factors, small cell lung cancer, colorectal cancer, T cell receptor signaling pathway, TGF-beta signaling pathway were dramatically reduced in cluster A, while drug metabolism other enzymes, pentose phosphate pathway, and fructose and mannose metabolism were higher in cluster A (Figure 7C). In cluster A, Hallmark activities of DNA repair, oxidative phosphorylation, fatty acid metabolism, myogenesis were higher, while G2M checkpoint, E2F targets, hedgehog signaling were lower (Figure 7D). According to the Reactome-based pathway, GSVA analysis displayed significant enrichment in insertion of tail anchored proteins into the endoplasmic reticulum membrane, post chaperonin tubulin folding pathway, WNT5A dependent internalization of FZD2 FZD5 and ROR2 in cluster A. In contrast, transcriptional activation of mitochondrial biogenesis, synthesis of PIPs at the early endosome membrane, cohesion loading onto chromatin, and meiotic synapsis were enriched in cluster B (Figure 7E).

Functional distinctions between the two Ub-related subclusters

In order to obtain a more profound comprehension of the functional differences between the two subclusters, a differential expression analysis was conducted. A total of 678 DEGs were identified with |log2 fold change| > 0.585 and P < 0.05, out of which 210 were upregulated and 468 were downregulated. The volcano plotting was employed to depict the distribution of DEGs between the two subclusters (Figure 8A). Additionally, GO and KEGG enrichment analysis was conducted on the 678 DEGs to better understand their possible molecular processes and functions. The genes showed enrichment in innate immune response, defense response to virus, antibacterial humoral response, negative regulation of viral genome replication, and innate immune response in mucosa, as was indicated by BP- Gene Ontology analysis (Figure 8B). Moreover, the KEGG enrichment analysis illustrated that these genes were predominantly enriched in transcriptional misregulation in cancer, acute myeloid leukemia, staphylococcus aureus infection, pertussis, and complement and coagulation cascades (Figure 8C). Immune cell infiltration investigation revealed that the abundance of neutrophils, resting memory CD4 T, activated memory CD4 T, activated NK, and activated Dendritic cells were higher in cluster B. In contrast, CD8 T, regulatory T cells, monocytes, and macrophages M0 were lower in cluster B than cluster A (Figure 8D). All these findings were consistent with GSVA results.

Figure 8. Functional and immune cell enrichment analysis between Ub-related subcluster. (A) The volcano plotting of DEGs. (B) Enrichment items in GO BP and (C) KEGG pathway analysis. (D) The proportion of all 22 types of immune cells. ^*p < 0.05; ^***p < 0.001.

Ub-related hub genes in TB-related disease

Numerous publications have already reported on the correlation between TB and various diseases, such as rheumatoid arthritis (RA) [24], chronic obstructive pulmonary disease (COPD) [25], corona virus disease 2019 (COVID-19) [26], and lung adenocarcinomas (LA) [27]. To delve further into the potential involvement of Ub-related genes in these diseases, we conducted an analysis using ssGSEA. Compared to HC, the score of E1 family was significantly changed in all TB-related diseases (Figure 9). The E2, DUB and UBD scores showed differently in RA, COVID-19 and LA. The results showed in Figure 9 suggested that ubiquitination may play a part in the development of these conditions. We observed varying levels of gene expression in the samples of these diseases, with TRIM68 standing out as having a significantly different expression across the range of diseases (Figure 9). That may suggest an essential role for TRIM68 in the association between TB and these diseases.

Figure 9. ssGSEA score of E1, E2, E3, DUB, UBD, ULD and the expression of 11 Ub-related hub genes in RA (A), COPD (B), COVID-19 (C), and LA (D) dataset. Abbreviations: RA: rheumatoid arthritis; COPD: Chronic obstructive pulmonary disease; COVID-19: corona virus disease 2019; LA: lung adenocarcinomas.

Additionally, TRIM68 displayed stronger correlations with a range of immune cells in various diseases (Figure 10), indicating a possible involvement in the immune processes related to these conditions.

Figure 10. The relationship between TRIM68 and 22 immune cells in four TB-related diseases. ^*p < 0.05; ^***p < 0.001.