Identification of Hub Genes and Typing of Tuberculosis Infections Based on Autophagy-Related Genes

From GEO (https://www.ncbi.nlm.nih.gov/geo) database, TB microarray data were downloaded (GSE54992, normal: 6, disease: 27; Cai et al. 2014). From the Human Autophagy Database (HADb; http://www.autophagy.lu/index.html), autophagy-related genes were obtained (Table SI).

Differential analysis of gene expression matrices (|logFC| > 0.585, adjust_p < 0.05) between healthy samples, as well as TB samples, was performed using the R package limma (Ritchie et al. 2015) to identify TB-associated DEGs (Table SII). Subsequently, TB-related DEGs were intersected with human autophagy-related genes to obtain autophagy-related DEGs in TB. R package heatmap (Zhao et al. 2014) was used to present the heat map of autophagy-related DEGs in TB, and R package ggplot2 (Ito and Murphy 2013) was used to draw volcano plots and box plots of autophagy-related DEGs.

GO and KEGG enrichment analyses (p < 0.05) of autophagy-related DEGs in TB were performed using the R package clusterProfiler (Yu et al. 2012). GO enrichment analysis consists of cellular components (CC), biological processes (BP), and molecular functions (MF).

Using the STRING database (https://string-db.org (Szklarczyk et al. 2023)), PPI analysis was performed on DEGs associated with autophagy in TB. Cytoscape software was used to visualize the PPI network (Shannon et al. 2003), and the topology characteristics of the PPI network were counted using the NetworkAnalyzer module. The correlation analysis of autophagy-related DEGs was performed using the R package corrplot ((https://cran.r-project.org/web/packages/corrplot/index.html (Wei and Simko 2021)).

To investigate whether miRNAs regulate the expression of these genes in TB, we collected miRNAs associated with M. tuberculosis infection from the Human microRNA Disease Database (http://www.cuilab.cn/hmdd (Lu et al. 2008; Li et al. 2014; Huang et al. 2019; Shen et al. 2022)); the experimentally validated miRNA databases associated with human disease), and we further narrowed down the scope based on whether there was a causal relationship (Table SIII). To explore the function of these miRNAs, miRPath online software (http://www.microrna.gr/miRPathv3) was used, and the GO enrichment analysis in terms of BP was performed on mRNAs corresponding to these miRNAs (p < 0.01).

The miRTarBase database (https://mirtarbase.cuhk.edu.cn/~miRTarBase/miRTarBase_2022/php/index.php) (Huang et al. 2022); the experimentally validated miRNA-target gene interaction database) was used to search the above-resulting miRNAs to find the target genes of these miRNAs, and autophagy-related DEGs were intersected with target genes for constructing miRNA-mRNA regulatory networks, and visualized using Cytoscape software.

CytoHubba plug-in of Cytoscape software (Chin et al. 2014) was used to calculate hub genes using five algorithms. The intersection of hub genes obtained from these five algorithms was taken, namely, autophagy-related hub genes in TB. The disease samples were grouped into high- and low-expression groups according to the median value of hub gene expression. GSEA software (Mootha et al. 2003; Subramanian et al. 2005) was used to conduct pathway enrichment analysis for the two groups to explore the signaling pathways affected by hub genes in TB.

R package ConsensusClusterPlus (Wilkerson and Hayes 2010) was used for k-means clustering analysis of the autophagy-related DEGs, and subgroup classification of the samples was conducted. ssGSEA analysis was performed using the R package GSVA (Hänzelmann et al. 2013) to assess the immune microenvironment of each sample in the subtype grouping. Afterward, the R package estimate (Yoshihara et al. 2013) was used to calculate the immune, stromal, and ESTIMATE scores. To further verify the correlation between subgroups and immunity, we measured the immune checkpoint expression levels of the samples. R package CIBERSORT (Newman et al. 2019) was utilized to assess the level of immune cell infiltration in each sample, and the levels of immune infiltration in the subgroups were compared.

In this study, 3,269 DEGs associated with TB were identified in the GEO dataset, of which 1,727 genes were significantly up-regulated, and 1,542 genes were significantly down-regulated. Two hundred twenty-two autophagy-related genes were obtained from the HADb. The DEGs associated with TB were intersected with autophagy-related genes to obtain 47 autophagy-related genes in TB (Fig. 1A and Table SIV). Subsequent differential expression analysis of these 47 genes revealed that 25 genes were significantly up-regulated and 22 were significantly down-regulated (Fig. 1B and 1C). We further studied the expression patterns of these 47 autophagy-related DEGs between normal and TB-infected samples (Fig. 1D and 1E). It was found that ARSA, ULK2, PTK6, ATG16L2, and TNFSF10 were the top five significantly up-regulated autophagy-related genes in TB. IL24, NFKB1, SPHK1, VEGFA, and CCL2 were the top five significantly down-regulated autophagy-related genes in TB. These genes were thought to most likely influence autophagy in TB.

Fig. 1.

To identify interactions between autophagy-related DEGs, PPI analysis was performed. The results showed the interactions between these autophagy-related genes (Fig. 2A). This PPI network's topological parameters (degree) were statistically analyzed, and the degree of the top 20 genes was shown (Fig. 2B). SQSTM1, MAPK8, UVRAG, HSPA5 and GABARAPL1 were the top five genes in PPI degree value and core genes in the PPI network. Correlation analysis was performed to explore the expression correlation of these autophagy-related genes. Fig. 2C showed the relationship between 47 autophagy-related DEGs in the GSE54992 dataset. The significant positive correlations between SQSTM1 and NF-κB; MAPK8 and NF-κB, GABARAPL1, HSPA5; UVRAG and BNIP1; HSPA5 and NF-κB, DNAJB9, GABARAPL1, SQSTM1, CCL2; GABARAPL1 and CCL2, SPHK1, SESN2 could be observed from the figure.

Fig. 2.

GO and KEGG enrichment analyses were performed to analyze the potential biological functions of these autophagy-related DEGs. GO enrichment analysis showed that genes were mainly enriched in cellular response to external stimulus, macroautophagy, and response to nutrient levels. The genes were mainly enriched in CC GO terms such as autophagosome, phagophore assembly site, and autophagosome membrane; genes were mainly enriched in MF GO terms in ubiquitin protein ligase binding (Fig. 3A). In KEGG enrichment analysis, autophagy-related DEGs were mainly involved in apoptosis, autophagy-animal, protein processing in the endoplasmic reticulum, and sphingolipid signaling pathway (Fig. 3B).

Fig. 3.

A search of the HMDD database identified 17 miRNAs associated with M. tuberculosis infection (Table SV), of which six miRNAs had causal relationship (hsa-let-7, hsa-mir-155, hsa-mir-206, hsa-mir-26a, hsa-mir-30a, and hsa-mir-32). GO enrichment analysis uncovered that the mRNAs corresponding to these miRNAs were mainly involved in BP, such as the mitotic cell cycle and cellular nitrogen compound metabolic process (Fig. 4A). Through miRTarbase, 3,095 target genes of these six miRNAs were predicted, and 2,623 target genes of miRNAs were obtained after removing repeat values. These target genes were subsequently intersected with 47 autophagy-related DEGs to obtain 20 autophagy-related genes targeted by miRNAs (Fig. 4B). To investigate the targeting relationship between miRNAs and autophagy-related genes, we constructed a regulatory network map of miRNA-mRNA using the Cytoscape software (Fig. 4C). There was a targeting relationship between these six miRNAs and multiple autophagy-related DEGs. Among them, GABARAPL1, CCL2, and NF-κB were all regulated by hsa-mir-155. MAPK8 and HSPA5 were regulated by hsa-mir-30a. UVRAG was regulated by hsa-mir-32. hsa-let-7 mainly acted on CDKN1A, and hsa-mir-26a mainly regulated the expression of BID and ULK2.

Fig. 4.

Five algorithms, including maximal clique centrality (MCC), the density of maximum neighborhood component (DMNC), maximum neighborhood component (MNC), degree and edge percolated component (EPC) in Cytohubba plug-in were used to calculate the critical coefficients of DEGs, respectively. Then top 10 hub genes (Table I) were selected, from which overlapping genes of these five algorithms were obtained, namely GABARAPL1 and ULK2 (Fig. 5A). Subsequently, to explore the influence of GABARAPL1 and ULK2 on the relevant functions of TB, a GSEA enrichment analysis was conducted. According to the threshold of FDR < 0.25, only GABARAPL1 was found to be significantly enriched. Its main enriched pathways were fatty acid metabolism, peroxisome, pyruvate metabolism, nucleotide excision repair, TCA cycle, and oxidative phosphorylation (Fig. 5B).

Fig. 5.

Autophagy-related DEGs in TB were subjected to consensus clustering genes by using the ConsensusClusterPlus package, and the TB samples were divided into subgroups. According to the area under the cumulative distribution function (CDF) curve, we could see that the internal stability of the cluster was the best when the number of subgroups k = 2 (Fig. 6).

Fig. 6.

Then, to study the immune differences among the autophagy gene-related subgroups, ssGSEA analysis was performed on these two subtypes. The estimate package was used to calculate each sample's immune score, stromal score, and ESTIMATE score. This analysis revealed heterogeneity in the immune microenvironment of different subgroups (Fig. 7A). The stromal, immune, and ESTIMATE scores of subgroup 1 were significantly higher than those of subgroup 2 (Fig. 7B–7D). To verify the accuracy of these results, we analyzed the expression levels of immune checkpoints and the infiltration levels of immune cell populations in different subgroups. The results exhibited that immune checkpoints were significantly upregulated in subgroup 1 compared with subgroup 2, and the level of immune cell infiltration was also increased (Fig. 7E and 7F). The monocytes, dendritic cells resting, B cells memory, and plasma cells were remarkably highly expressed in subgroup 1.

Fig. 7.