Identification potential biomarkers in pulmonary tuberculosis and latent infection based on bioinformatics analysis

Background The study aimed to identify the potential biomarkers in pulmonary tuberculosis (TB) and TB latent infection based on bioinformatics analysis. Methods The microarray data of GSE57736 were downloaded from Gene Expression Omnibus database. A total of 7 pulmonary TB and 8 latent infection samples were used to identify the differentially expressed genes (DEGs). The protein-protein interaction (PPI) network was constructed by Cytoscape software. Then network-based neighborhood scoring analysis was performed to identify the important genes. Furthermore, the functional enrichment analysis, correlation analysis and logistic regression analysis for the identified important genes were performed. Results A total of 1084 DEGs were identified, including 565 down- and 519 up-regulated genes. The PPI network was constructed with 446 nodes and 768 edges. Down-regulated genes RIC8 guanine nucleotide exchange factor A (RIC8A), basic leucine zipper transcription factor, ATF-like (BATF) and microtubule associated monooxygenase, calponin LIM domain containing 1 (MICAL1) and up-regulated genes ATPase, Na+/K+ transporting, alpha 4 polypeptide (ATP1A4), histone cluster 1, H3c (HIST1H3C), histone cluster 2, H3d (HIST2H3D), histone cluster 1, H3e (HIST1H3E) and tyrosine kinase 2 (TYK2) were selected as important genes in network-based neighborhood scoring analysis. The functional enrichment analysis results showed that these important DEGs were mainly enriched in regulation of osteoblast differentiation and nucleoside triphosphate biosynthetic process. The gene pairs RIC8A-ATP1A4, HIST1H3C-HIST2H3D, HIST1H3E-BATF and MICAL1-TYK2 were identified with high positive correlations. Besides, these genes were selected as significant feature genes in logistic regression analysis. Conclusions The genes such as RIC8A, ATP1A4, HIST1H3C, HIST2H3D, HIST1H3E, BATF, MICAL1 and TYK2 may be potential biomarkers in pulmonary TB or TB latent infection.


Background
Pulmonary tuberculosis (TB) is a widespread and fatal infectious disease. It is caused by various strains of mycobacteria, usually Mycobacterium tuberculosis [1]. It is estimated that one third of the world's population are infected with M. tuberculosis [2]. More than 90 % of infected individuals remain asymptomatic with a latent infection [3]. With aging or immune system deteriorating, M. tuberculosis can reactivate and cause severe pulmonary TB [4]. Roughly 10 % of the latent infections can progress to active TB. The general signs and symptoms of this disease include fever, chills, night sweats, loss of appetite, weight loss, and fatigue [5]. Approximately, there are 9 million newly diagnosed cases of pulmonary TB and 1.5 million deaths annually, mostly in developing countries [6]. Therefore, uncovering therapeutic biomarkers in pulmonary TB would supply new insights for the diagnosis and treatment of this disease.
Numerous studies have been done to investigate the potential biomarkers for the treatment of pulmonary TB. For example, the serum CA-125 level is found significantly higher in active pulmonary TB than in inactive TB or normal sample, suggesting that CA-125 may be a beneficial parameter in determination of pulmonary TB activity [7]. Pollock et al. [8] suggested that M. tuberculosis Rv1681 protein was a diagnostic marker of active pulmonary TB. Additionally, Chowdhury et al. [9] reported that the serum interleukin (IL)-6 level of the active pulmonary TB patients following anti-tuberculosis drug therapy played an important role in immuneprotection of the host against M. tuberculosis infection. Although many factors have been found, the diagnostic efficiency of pulmonary TB is still unsatisfactory [10]. Therefore, it is necessary to identify novel potential therapeutic biomarkers in pulmonary TB.
In the present study, the microarrays data GSE57736 were downloaded to identify the differentially expressed genes (DEGs) between pulmonary TB and latent tuberculosis infection samples. This dataset is deposited by Guerra-Laso et al. [11], the study of whom demonstrates that IL-26 is a candidate gene for TB susceptibility. In this study, we aimed to use different bioinformatics method to identify the DEGs between the two kinds of samples. Based on the obtained DEGs, we performed proteinprotein interaction (PPI) network construction and network-based neighborhood scoring analysis. Besides, the hierarchical clustering analysis, functional enrichment analysis, correlation analysis and logistic regression analysis of DEGs were performed as well. Findings of this study may help to explore potential targets for the diagnosis and treatment in pulmonary TB.

Affymetrix microarray data
The array data of GSE57736 based on the platform of GPL13497 (Agilent-026652 Whole Human Genome Microarray 4x44K v2) was downloaded from Gene Expression Omnibus database, which was deposited by Guerra-Laso et al. [11]. The dataset available in this analysis contained 15 peripheral blood samples from seven pulmonary TB patients and eight latent tuberculosis infections. Among the seven pulmonary TB patients, there were three men and four women (average 82.7 years) with different clinical conditions: psoriasis (one patient), previous heart failure (one patient), arterial hypertension (two patients), bronchial asthma (one patient), chronic obstructive pulmonary disease (two patients), and prostate cancer (one patient). The eight latent tuberculosis infection samples included six men and two women (average 81.1 years), which had scored a positive result in the QuantiFERON-TB Gold in-tube test (Cellestis, Carnegie, Vic., Australia).

Data preprocessing and differential expression analysis
The probe IDs were converted into corresponding gene symbols based on the annotation information on the platform. When multiple probes corresponded to a same gene, the average expression value was calculated to represent the gene expression level. The limma package [12] in R was used to identify DEGs between pulmonary TB and TB latent infection samples. The Benjamin and Hochberg (BH) [13] method was used to adjust the raw p-values [false discovery rate (FDR)]. Then, log 2 -fold change (log 2 FC) was calculated. Only genes with |log 2 FC| > 1.0 and FDR < 0.05 were selected as DEGs.

PPI network construction
Human Protein Reference Database (HPRD, http:// www.hprd.org/) [14] is a database of curated proteomic information pertaining to human proteins. Search Tool for the Retrieval of Interacting Genes (STRING, http:// string.embl.de/) [15] is an online database which collects comprehensive information of proteins. In our study, the DEGs were mapped into STRING and HPRD databases to identify significant protein pairs with confidence score > 0.4. Then the PPI network was constructed based on these protein pairs using Cytoscape software [16].

Network-based neighborhood scoring
Neighborhood scoring [17] is a local method for prioritizing candidates based on the distribution of DEGs in the network. Gene in PPI network was assigned a score, which was based on its FC and the FC of its neighbors. The score of each node in PPI network was calculated with the neighborhood scoring method [18]. When the hub node and its neighborhood nodes were significantly highly expressed, the score > 0; When the hub node and its neighborhood nodes were significantly lowly expressed, the score < 0. Therefore, the top 50 nodes with higher scores and the last 50 nodes with lower scores were identified as important genes.
In order to confirm the efficiency of these important genes differentiating pulmonary TB and TB latent infection samples, hierarchical clustering analysis of the important genes was performed using cluster software [19]. The results were presented by TreeView software [20]. The expression profile data were filtered and normalized using cluster software. In detail, genes that were expressed in at least 80 % of the samples were selected. Besides, the genes and samples were normalized with median center method [21].

Functional enrichment analysis
Gene Ontology (GO, http://www.geneontology.org) database [22] is a collection of a large number of gene annotation terms. Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.genome.ad.jp/kegg/) knowledge database [23] is applied to identify the functional and metabolic pathway. Database for Annotation, Visualization and Integrated Discovery (DAVID, https://david.ncifcrf.gov/) [24] is a tool that provides a comprehensive set of functional annotation for large list of genes. In this study, the important genes were performed GO and KEGG pathway enrichment analyses with DAVID. With the enrichment threshold of p-value < 0.05, the DEGs enrichment results in GO terms and KEGG pathways were obtained.

Correlation analysis
The immune system dysfunction has been suggested to play an important role in the occurrence of pulmonary tuberculosis [25]. Therefore, it is necessary to analyze the correlations of genes associated with immune system. In the present study, the Pearson correlation coefficient (PCC) was calculated across 15 samples to investigate potential regulatory relationships between important genes. The gene pairs with |PCC| > 0.5 were selected for further analysis.

Logistic regression analysis
In order to identify the risk biomarkers of pulmonary TB, we performed multivariate logistic regression analysis for the gene pairs with significant correlations (|PCC| > 0.5) using SPSS 19.0 software (SPSS Inc., Chicago, Illinois, USA) [26]. The genes with p-value < 0.05 were selected as feature genes.

Identification of DEGs
Based on the thresholds of |log 2 FC| > 1.0 and FDR < 0.05, a total of 1084 DEGs were identified between pulmonary TB and TB latent infection samples, including 565 down-regulated genes and 519 up-regulated genes. The result was shown in volcano plot (Fig. 1).

Network-based neighborhood scoring
The top 50 nodes with higher scores and the last 50 nodes with lower scores were selected and the top 5 and last 5 genes were shown in Table 1. For instance, sirtuin 5 (SIRT5), tyrosyl-tRNA synthetase (YARS), sphingomyelin  Furthermore, hierarchical clustering analysis for these important genes showed that these 100 important genes could differentiate the pulmonary TB samples and TB latent infection samples (Fig. 3).

Functional enrichment analysis
GO enrichment analysis was carried out for the important genes. The significant (p < 0.05) GO biological process (BP) terms of up-and down-regulated genes were shown in Table 2 (p-values in ascending order). The down-regulated genes were significantly enriched in regulation of osteoblast differentiation (p = 0.00816) and positive regulation of hydrolase activity (p = 0.018904). Besides, the up-regulated genes were mainly related to nucleoside triphosphate biosynthetic process (p = 0.003512), mRNA export from nucleus (p = 0.004389) and purine nucleoside triphosphate metabolic process (p = 0.005796).
In addition, 2 pathways were enriched by the upregulated important genes (Table 2), including adipocytokine signaling pathway and cardiac muscle contraction pathway. However, the down-regulated genes were not enriched in any pathways.

Correlation analysis
A total of 950 gene pairs were identified in Pearson correlation analysis. The top 10 highly correlated gene pairs were shown in Table 3. Specially, the expression levels of the top 4 correlated gene pairs (PCC > 0.9) were shown in Fig. 4, that was, RIC8A-ATP1A4, HIST1H3C-HIST2H3D, HIST1H3E-BATF and MICAL1-TYK2, besides, all of them showed positive correlations.

Discussion
In this study, a total of 1084 DEGs including 565 down-and 519 up-regulated genes were selected. The up-regulated genes were mainly related to nucleoside triphosphate biosynthetic process. The down-regulated genes were significantly enriched in regulation of osteoblast differentiation. The gene pairs RIC8A-ATP1A4, HIST1H3C-HIST2H3D, HIST1H3E-BATF and MICAL1-TYK2 were identified with highly positive correlations. Besides, they were selected as feature genes in logistic regression analysis.
RIC8A encoding protein interacts with guanine nucleotide binding protein (G protein) [27]. It has been reported that RIC8A controls Drosophila neural progenitor asymmetric division by regulating heterotrimeric G proteins [28]. G protein is an important signal transducing molecule in cells [29], which activates MAP kinase signaling [30]. Elkington et al. [31] have reported that active pulmonary TB can be mediated by MAP kinase signaling pathway. In this study, RIC8A was down-regulated in pulmonary TB, suggesting that RIC8A may be associated with the pulmonary TB development through regulating MAP kinase signaling pathway with G proteins.
ATP1A4 is a member of P-type cation transport ATPase family and belongs to Na, K-ATPase subfamily. The P-type ATPases remove Ca 2+ against very large concentration gradients in eukaryotic cells and play an important role in intracellular calcium homeostasis [32]. Importantly, calcium homeostasis involves in apoptosis and regulates important cellular events triggered upon infection of macrophages with pathogenic mycobacteria [33]. It has been reported that the M. tuberculosis blocks the delivery of the Na, K-ATPases [34]. Additionally, Rao et al. [35] have shown that de novo ATP synthesis is essential for the viability of hypoxic nonreplicating mycobacteria, requiring the cytoplasmic membrane to be fully energized. Interestingly, ATP1A4 was found enriched in nucleoside triphosphate biosynthetic process in this study. Therefore, we speculated that ATP1A4 may play a vital role in the occurrence of pulmonary TB by controlling ATP synthesis.  HIST1H3C, HIST2H3D and HIST1H3E belong to histone H3 family, which are responsible for controlling the dynamics of the chromosomal fiber in eukaryotes by regulating histone acetylation. Importantly, this process is essential in modulating gene transcription through chromatin organization, and perturbation of this process can result in aberrant gene transcription and cause some diseases, including lung diseases [36,37]. Additionally, histones play a central role in DNA repair and DNA replication [38]. Boshoff et al. [39] reported that DNAdamaging agents were rich in vivo produced by host cells due to an effort to eradicate the M. tuberculosis. BATF belongs to the adaptor-related protein 1 (AP-1)/ activating transcription factor (ATF) superfamily of transcription factors. AP-1 family transcription factors control the differentiation of lymphocyte cells in immune system [40]. Lymphocytes are crucial in the immune defence against M. tuberculosis, which can secrete interferons (ITFs) in response to M. tuberculosis infection [41]. It has been reported that ITF-γ, a product of T lymphocytes, contributes to protective immunity against M. tuberculosis by activating macrophages in pulmonary TB [42]. Taken together, although the role of BATF in pulmonary TB has not been studied, we speculate that BATF may be involved in the occurrence of TB via immune system.
TYK2 encodes a tyrosine kinase belonging to Janus kinases (JAKs) family. It has been reported that JAKs are activated following interactions between cytokines and their cognate receptors on cell surface [43]. TYK2 negatively regulates adaptive Th1 immunity by mediating IL-10 signaling and promoting ITF-γ-dependent IL-10 reactivation [43]. Redford et al. [44] have reported that IL-10 can suppress the functions of macrophage and dendritic cell, which were required for the capture, control and initiation of immune responses to M. tuberculosis. Therefore, TYK2 may involve in the pathogenesis of pulmonary TB via regulating IL-10. For MICAL1, it can act as a cytoskeletal regulator [45]. Specially, cell migration and phagocytosis are critically dependent on cytoskeletal rearrangements [46]. It has been reported that cell migration and phagocytosis are important for resistance against pulmonary TB [47]. Therefore, the down-regulation of MICAL1 may be related to pulmonary TB via controlling cell migration and phagocytosis.

Conclusions
In conclusion, the present study identified several key gene pairs (RIC8A-ATP1A4, HIST1H3C-HIST2H3D, HIST1H3E-BATF, MICAL1-TYK2) associated with pulmonary TB or TB latent infection by comprehensive bioinformatics methods, which may provide new insights for the diagnosis and treatment of this disease. However, this study had some limitations. On the one hand, the sample size was small which might cause a high rate of false positive results. Secondly, there was no experimental verification. Therefore, further genetic and experimental studies with larger sample sizes are still needed to confirm the findings in this study.