Abstract

Background. Tuberculosis (TB) is usually caused by Mycobacterium tuberculosis, which has the highest mortality rate among infectious diseases. This study is designed to identify the key genes affecting the diagnosis and treatment of TB. Methods. GSE54992, which included 39 peripheral blood mononuclear cell (PBMC) samples, was extracted from the Gene Expression Omnibus database. After the samples were classified into type and time groups by limma package, the differentially expressed genes (DEGs) were analyzed using the Analysis of Variance. Using pheatmap package, hierarchical cluster analysis was performed for the DEGs. Then, the key modules correlated with TB were selected using the WGCNA package. Finally, functional and pathway enrichment analyses were carried out using clusterProfiler package. Results. The DEGs in subclusters 3, 6, 7, and 8 were chosen for further analyses. Based on WGCNA analysis, blue and green modules in type group and pink module in time group were selected as key modules. From the key modules, 9 (including BAX and ARPC1B) hub genes in type group and 6 (including DHX36) hub genes in time group were screened. Through pathway enrichment analysis, the TNF signaling pathway was enriched for the green module. Conclusion. DHX36, BAX, and ARPC1B might be key genes acting in the mechanisms of TB. Besides, the TNF signaling pathway might also be critical for the diagnosis and therapy of the disease.

1. Introduction

As an infectious disease, tuberculosis (TB) is mainly induced by Mycobacterium tuberculosis (MTB) and usually affects the lungs [1, 2]. Latent TB has no obvious signs, and approximately 10% of them can develop into active disease [3]. The typical symptoms of active TB include night sweats, coughing blood, weight loss, and fever [4]. Latent TB is not spread, but active TB can be spread via the air when lung TB patients spit, cough, sneeze, or speak [5]. Smokers and people with human immunodeficiency virus (HIV)/Acquired Immune Deficiency Syndrome (AIDS) are at high risk of active TB, and thus early screening and treatment of high-risk people and vaccination are the main methods for preventing TB [68]. A third of the world’s population suffers from TB, and the mortality of TB ranks first among infectious diseases [9]. Active TB affects over 10 million people and leads to 1.3 million death cases in 2016 [10, 11]. Therefore, TB should be deeply investigated to reveal its mechanisms.

The product of the Intracellular Pathogen Resistance 1 (IPR1) gene may function in integrating signals induced by intracellular pathogens through mediating cell death, innate immunity, and pathogenesis, and thus IPR1 is a candidate gene controlling the host resistance to TB [12]. Murine β-defensin–3 (mBD3) and mBD4 expression are induced by mycobacterial infection, which may play roles in controlling mycobacterial growth during TB infection [13]. As an inner membrane transporter of MTB, transmembrane transport protein MmpL3 (MMPL3) participants in the transport of trehalose monomycolate and is a novel target for the treatment of TB patients [14]. Programmed death 1 (PD-1) is implicated in the functions of T cell effector against MTB; therefore, PD-1 can mediate the immune response in hosts during human TB [15, 16]. Through the interleukin (IL) 4Rα signaling pathway, the T helper (Th) 2 response regulates the alternative activation of macrophages and thus contributes to the intracellular persistence of MTB [17]. However, the above researches only report a part of the genes involved in TB, and more studies should be conducted to fully reveal the pathogenesis of the disease.

In 2014, Cai et al. determine the expression pattern of C1q in peripheral blood mononuclear cells (PBMCs) to explore the function of C1q in TB, finding that C1q is closely related to the active disease and disease severity in TB and serves as a diagnostic biomarker for the disease [18]. Nevertheless, more genes correlated with the diagnosis and progression of TB needed to be explored to prevent the deterioration of TB. Through performing comprehensive bioinformatics analyses for the microarray dataset uploaded by Cai et al. [18], the genes playing key roles in the progression of TB were investigated. This study might broaden our understanding of the mechanisms of TB and promote the diagnosis and treatment of the disease.

2. Materials and Methods

2.1. Data Source

The expression profiling data of TB (accession number: GSE54992), which was based on the platform of GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array, was extracted from Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/) database. There were a total of 39 PBMC samples in GSE54992, including 6 samples from healthy donors (HD), 6 samples from patients with latent TB infection (LTBI), 9 samples from TB patients (TB), 9 samples from TB patients after antituberculosis treatment for 3 months (TB3m), and 9 samples from TB patients after antituberculosis treatment for 6 months (TB6m). The samples were from participants recruited at Shenzhen Third People’s Hospital from May 2011 to December 2012. PBMCs were separated from heparinized whole blood as described previously [19]. Cai et al. [18] uploaded GSE54992, and their research was approved by the Institutional Review Board of Shenzhen Third People's Hospital and obtained the informed consent of all participants.

2.2. Data Preprocessing and Differential Expression Analysis

The original data in GSE54992 were preprocessed using the R packages Affy [20, 21] (version 1.52.0, http://bioconductor.org/packages/release/bioc/html/affy.html) and limma [22] (version 3.32.2, https://bioconductor.org/packages/release/bioc/html/limma.html). The preprocessing processes were background correction, normalization, log2 conversion, and probe annotation. The probes without matching gene symbols were filtered out. Afterward, the expression value of the gene corresponding to several probes was acquired by calculating the mean value of the probes.

Using the R package limma [22], the data were conducted with standardization analysis, and the samples were classified into type groups (including HD, LTBI, and TB samples) and time groups (including TB, TB3m, and TB6m samples). Based on the Analysis of Variance (ANOVA) [23], the differentially expressed genes (DEGs) between type and time groups were analyzed. The DEGs were defined as genes with value ≤ 0.05 and |fold change (FC)| ≥ 2.

2.3. Hierarchical Cluster Analysis

To identify the genes with similar expression patterns, hierarchical cluster analysis was performed for the DEGs using the R package pheatmap [24] (version 1.0.2, https://cran.r-project.org/web/packages/pheatmap/index.html). The distance calculation algorithm, genetic clustering method, and the clustering method for gene clusters separately were Euclidean, kmeans, and hcluster. To screen the targets that could be used for the diagnosis and treatment of TB, the genes significantly dysregulated between TB and HD/LTBI groups and that were near the expression in HD group along with TB-TB3m-TB6m treatments were defined as TB-specific genes and utilized for the subsequent analyses.

2.4. Weighted Gene Coexpression Network Analysis (WGCNA)

WGCNA is an algorithm developed for investigating module information from high-throughput data [25]. The R package WGCNA (version 1.61, https://cran.r-project.org/web/packages/WGCNA/index.html) [25] was applied for analyzing the DEGs and the expression data of the DEGs were utilized as the input data for building the coexpression network. The main processes of WGCNA were coexpression network construction and module identification, the identification of disease-associated modules, enrichment analysis for key modules and protein network construction, and the identification and enrichment analysis of hub nodes in the key modules.

2.5. Functional and Pathway Enrichment Analysis

Using the R package clusterProfiler (version 3.4.4, https://bioconductor.org/packages/release/bioc/html/clusterProfiler.html) [26], Gene Ontology (GO) [27] functional and Kyoto Encyclopedia of Genes and Genomes (KEGG) [28] pathway enrichment analyses for the DEGs were conducted. The Benjamini and Hochberg (BH) method [29] was used for adjusting the values, and the adjusted value <0.01 was set as the threshold for significant results.

3. Results

3.1. Differential Expression Analysis

After the 54676 probes in GSE54992 were preprocessed, 23520 genes were obtained. Principal component analysis (PCA) for the samples showed that HD and LTBI samples had little differences, and TB and TB3m samples had similar expression patterns (Figure 1).

There were 520, 2931, and 2887 DEGs separately in LTBI vs. HD, TB vs. HD, and TB vs. LTBI type comparison groups. Besides, a total of 462, 1502, and 741 DEGs separately were screened in TB3m vs. TB, TB6m vs. TB, and TB6m vs. TB3m time comparison groups (Table 1). The cluster heatmap for all DEGs showed that the DEGs could separate the samples in different groups very well (Figure 2).

3.2. Hierarchical Cluster Analysis

With kmeans = 9, the DEGs were performed with hierarchical cluster analysis (Figure 3). Subsequently, the genes in each cluster were compared, and the genes (a total of 3108 genes) in subclusters 3, 6, 7, and 8 were chosen for the following analyses.

3.3. WGCNA Analysis

The expression matrixes of the 3108 DEGs in the key subclusters were extracted and taken as the input data for constructing coexpression network. Coexpression network should have the characteristics of a scale-free network. Therefore, the weighting parameter β (soft threshold power) needed to be negatively correlated with the square of the correlation coefficients between log (k) and log (p(k)). The higher the square was, the closer the coexpression network was to scale-free network. When the square firstly approached 0.85, the corresponding β = 10 was suitable for building coexpression network (Figure 4(a)). Besides, the mean connectivity was −1.35 when the β value was 10 (Figure 4(b)).

After the system clustering tree was obtained for the genes, 12 network modules (at least 30 genes were involved in each module) were identified (Figure 4(c)). Then, cluster analysis for the modules was conducted, and a total of 8 modules (black module, involving 82 genes; blue module, involving 953 genes; brown module, involving 473 genes; green module, involving 1267 genes; grey module, involving 103 genes; pink module, involving 75 genes; purple module, involving 47 genes; red module, involving 108 genes) were finally obtained after merging the closely clustered modules (height was set at 0.1) (Figure 4(c)).

Gene significance (GS) is defined as the mediated value of each gene in the linear regression between gene expression and the sample traits. Module significance (MS) was defined as the average GS within modules and was calculated to measure the correlation between modules and sample traits. If GS and MS are highly correlated, it means that genes are the most important elements of modules and are highly significantly associated with the trait. According to the absolute value of the correlation coefficient between each module and disease state, blue, green, and red modules were the top 3 modules in the type group. Meanwhile, pink, blue, and green modules were the top 3 modules in the time group (Figure 5(a)). Based on the absolute value of the GS in each module, the key modules in type (blue and green modules) (Figure 5(b)) and time (brown and pink modules) (Figure 5(c)) groups separately were selected. As a result, blue and green modules in type group and pink module in time group were selected by both of the two methods and thus used for screening the hub genes related to the disease.

There separately were 874 and 27 hub nodes in type and time groups with the cutoff of MS > 0.8 ( value < 0.01) and GS > 0.2 ( value < 0.01). Besides, the top 10 genes were selected as candidate hub genes with networkScreening function in the WGCNA package. In addition, the hub nodes were intersected with these candidate hub genes, and the overlapped genes were redefined as the hub genes. Finally, 9 hub genes (including BCL2-associated X protein, BAX; and Actin-Related Protein 2/3 Complex, Subunit 1B, ARPC1B) in type group and 6 hub genes (including DEAH (Asp-Glu-Ala-His) box polypeptide 36, DHX36) in time group were screened. The expression diagrams of the 15 hub genes are shown in Figure 6. The 15 hub genes were mainly involved in pink and blue modules, among which 6 hub genes in time group were specifically downregulated expressed in TB, and 9 hub genes in type group were specifically upregulated expressed in TB.

3.4. Enrichment Analysis for Key Modules and Hub Genes

Functional (Figure 7(a)) and pathway (Figure 7(b)) enrichment analyses for the key modules showed that no significant functional term and pathway were enriched for the pink module. Besides, the tumor necrosis factor (TNF) signaling pathway was enriched for the green module. Moreover, the hub genes in type and time groups were also conducted with enrichment analysis. The results showed that the hub genes in the time group were implicated in the functional term of 7-methylguanosine mRNA capping (Figure 7(c)). However, the hub genes had no significantly enriched pathways.

4. Discussion

In this study, the DEGs in type and time comparison groups separately were screened. After performing hierarchical cluster analysis, the DEGs in subclusters 3, 6, 7, and 8 were chosen for further analyses. WGCNA analysis indicated that blue and green modules in the type group and pink module in the time group were key modules. Subsequently, 9 (including BAX and ARPC1B) and 6 (including DHX36) hub genes separately were identified in type group and time group. Pathway enrichment analysis showed that the TNF signaling pathway was enriched for the green module.

Increased B-cell CLL/lymphoma 2 (BCL2) and decreased BAX are detected in macrophages, and BCL2 overexpression in macrophages carrying MTB may be related to its intracellular survival [30]. Immunohistochemical staining shows that overexpressed BAX, P53, and Fas cell surface death receptor (FAS) have correlations with reduced BCL2 in TB granulomas [31]. MTB infection causes apoptosis of human neutrophils through inducing reactive oxygen species- (ROS-) dependent expression change of Bax/Bcl-x(L) and caspase-3 activation [32]. The recombinant Bacille Calmette-Guerin (rBCG): BAX strain contributes to the induction of Th1 protective immune responses, which may be a promising vaccine candidate for TB [33]. ARPC1B deficiency can lead to severe combined immunodeficiency with signs of mild bleeding and immune disorder [34]. A previous study considers that DHX36 is correlated with the DNA biosensor of MTB [35]. These indicated that DHX36, BAX, and ARPC1B might be related to the mechanisms of TB.

Proteinase-activated receptor-2 (PAR2), TNF, and galectin 9 (GAL9) pathways are essential for limiting MTB growth, and lipoarabinomannan (LAM) can reduce their activation to promote the intracellular growth of MTB [36]. The interactions between MTB grown in the condition of hypoxia and host macrophages can induce the TNF signaling pathway, DNA-damage stress response, and activation of apoptosis, especially, MTB-H bacilli which are sensitive to TNF-governed killing [37, 38]. TNF plays a role in MTB-induced macrophage apoptosis, which is correlated with the TNF- and c-Cbl-dependent FLIP(S)-degradation pathway [39]. Thus, the genes in the green module might function in TB via the TNF signaling pathway.

In conclusion, DHX36, BAX, and ARPC1B might be involved in the diagnosis and treatment of TB. In addition, the TNF signaling pathway might also be important for the development of TB. However, experimental researches should be carried out in the future to support our results.

Abbreviations

PCA:Principal component analysis
HD:Healthy donors
LTBI:Latent TB infection
TB:Tuberculosis
DEGs:Differentially expressed genes
MTB:Mycobacterium tuberculosis
HIV:Human immunodeficiency virus
AIDS:Acquired Immune Deficiency Syndrome
IPR1:Intracellular Pathogen Resistance 1
PBMCs:Peripheral blood mononuclear cells
TNF:Tumor necrosis factor
LAM:Lipoarabinomannan.

Data Availability

The datasets analyzed for this study can be found in the Gene Expression Omnibus (GEO) database (GSE54992) (http://www.ncbi.nlm.nih.gov/geo/).

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Authors’ Contributions

YZ was responsible for conception and design of the research. YL and HL were responsible for acquisition of data. QL performed analysis and interpretation of data. WW performed statistical analysis. YZ and YL drafted the manuscript. ZJ and WL revised the manuscript for important intellectual content. All authors read and approved the final manuscript.