1. Introduction
Tuberculosis (TB) is one of the most common infectious diseases, caused by the pathogen Mycobacterium tuberculosis (MTB). It was linked to a high rate of infection and a long-term disease course [Reference Dheda1]. TB is one of the top 10 causes of death according to the World Health Organization (WHO)'s global TB report for 2020. Each year, about ten million people were infected with tuberculosis. TB was identified as a repetitive immune reaction, and a subset of patients with tuberculosis will grow into active tuberculosis. It has a close connection with poverty. Around 90% of tuberculosis patients stayed in underdeveloped countries [Reference Strausz2]. There are many diagnostic tools available today, including the tuberculin skin examination (TST), the interferon gamma release assay (IGRA), and imaging procedures. They all lack precision and are more expensive and highly technical [Reference Cavusoglu3, Reference Zhou, Li and Zhou4].
microRNAs (miRNAs) are noncoding RNA molecules that have 22 to 23 nucleotides [Reference Huang5]. They regulate gene expression at the posttranscriptional stage by facilitating mRNA degradation and preventing mRNA translation [Reference Correia6]. When a patient is infected with tuberculosis, major miRNAs are released into the bloodstream. MiRNAs have shown to play an important role in a variety of pathological and physiological mechanisms in tuberculosis [Reference Backes, Meese and Keller7]. Microarray analysis, miRNA, and gene have been extensively utilized in the quest for biomarkers and therapeutic targets in recent years. The discovery of novel DE miRNAs and DEGs provides useful and accurate perspectives for future study [Reference Zheng, Qing and Song8].
To validate the significant genes and miRNAs, two mRNA expression profile datasets, GSE62147 [Reference Tientcheu, Maertzdorf and Weiner9] and GSE34608-GPL6480 [Reference Maertzdorf, Weiner and Mollenkopf10], and three miRNA expression profile datasets, GSE34608-GPL7731 [Reference Maertzdorf, Weiner and Mollenkopf10], GSE29190 [Reference Wang, Yang and Sun11], and GSE49951 [Reference Siddle, Deschamps and Tailleux12], were used. The 250 DEGs in speech profile datasets were, respectively, identified. The enrichment analysis of DEGs showed 5 KEGG pathways were related to the development and progress of TB, including ‘hepatocellular carcinoma,’ ‘Kaposi sarcoma-associated herpesvirus infection,’ ‘phosphatidylinositol signaling system,’ ‘circadian entrainment,’ and ‘apelin signaling pathway.’ A protein-protein interaction (PPI) network of DEGs and DEGs-DE miRNAs was created.10 hub genes were then chosen.
The present study established two DE miRNAs, hsa-mir-7 and hsa-mir-451. The gene-encoding hsa-mir-7 and hsa-mir-451 were separately located in human chromosomal region 9q21.32 and 17qll.2. The recent documents presented that Hsa-mir-7 and hsa-mir-451were involved in cancer-related biological processes and innate immune response [Reference Bai and Wu13, Reference Alamdari-Palangi, Karami, Karami and Baazm14]. But, their future roles in TB remain unknown.
2. Materials and Methods
2.1. Microarray Data Screening
The bioinformatics analysis was carried out in accordance with the protocol depicted in Figure 1. The term 'Tuberculosis” or “microRNA and Tuberculosis” was searched in the Gene Expression Omnibus database (GEO). There are some database criteria that are helpful with more analysis: I) clinical study pairs must have both healthy and tuberculosis patients, and II) mRNA and miRNA data from plasma (Affymetrix miRNA 4.0) were retrieved. As a result, two gene expression profile datasets, GSE62147 and GSE34608-GPL6480, were used in the primal mRNA datasets. GSE62147 had 14 TB-positive donors and 14 healthy donors. GSE34608-GPL6480 had eight TB donors and eighteen healthy donors. Six donors with tuberculosis and three healthy donors contained GSE29190. GSE49951 had 71 TB-positive donors and 71 healthy donors.
Figure 1 Flow diagram of bioinformatics analysis.
2.2. Selecting DEGs and DE miRNAs and Constructing a Volcanic Map
To diagnose DEGs and DE miRNAs in TB patients, the GEO2R online research method (https://www.ncbi.nlm.nih.gov/geo/geo2r/) was used. Exact cutoff values (P < 0.05 and |log fold change (FC)|≥1) were established. The two gene and three miRNA datasets were submitted to VENN's online tool (http://bioinformatics.psb.ugent.be/webtools/Venn/). The DEGs and DE miRNAs that overlap were filtered out. The volcanic chart of overlapping DEGs was constructed using the online method of bioinformatics (http://www.bioinformatics.com.cn/).
2.3. GO and KEGG Pathway analysis
Gene Ontology (Go) was applied to perform the category biological process (BP), molecular function (MF), and cellular component (CC) enrichment analysis [Reference Ashburner, Ball and Blake15]. The KEGG was a set of databases that provides information regarding biological mechanisms, cellular processes, chemical substances, and diseases [Reference Kanehisa and Goto16]. The R software package clusterProfiler was used in functional gene annotations and KEGG enrichment analysis [Reference Zhang, Zhu and Pu17]. A cutoff criterion (P < 0.05 and FDR<0.05) was set. Also, the top 10 pathways with the maximum number of genes for the corresponding term were chosen. Significant items of GO and KEGG were submitted to R software package ggplot2 to visualize and merge enriching analysis. The enrichment dot bubble method was used to build the bubble plot.
2.4. Construction of Network and Identification of Top Modules and Hub Genes
Permanent DEGs were uploaded to the STRING online tool (https://string-db.org/) [Reference Szklarczyk, Gable and Lyon18]. A cutoff point (interaction score >0.4) was established. The PPI network was constructed using the CytoScape program [Reference Shannon19]. The MCODE plug-in was used to determine the PPI network's top module. Strict cutoff conditions (degree cutoff = 2, node score cutoff = 0.2, k-core = 2, and maxdepth = 100) were established. Meanwhile, the top 10 hub genes were identified using the CytoHubba plug-in.
2.5. Identification of the Target Gene and Construction of the DE miRNA-DEGs Regulatory Network
To predict miRNA target genes, the online tools TargetScan (http://www.targetscan.org/vert71) [Reference Peterson, Thompson, Ufkin, Sathyanarayana, Liaw and Congdon20], miRDB (http://www.mirdb.org/) [Reference Weng, Zhang and Huang21], miRWalk (http://mirwalk.uni-hd.de/) [Reference Sticht, De La Torre, Parveen and Gretz22], and miRTarBase (http://mirtarbase.cuhk.edu.cn/) [Reference Weng, Zhang and Huang21] were used. The criterion for goal gene selection is cumulative weighted context++ score online > -0.5. The target gene that met the criteria was chosen from three databases. The overlapping target gene of miRNA was discovered using the VENN online method. CytoScape was used to upload miRNA and DEGs. To visualize and merge analysis, a DE miRNA-DEGs regulatory network was developed.
3. Results
3.1. Screening of DEGs and DE miRNAs
Under the precise cutoff criteria employed (P < 0.05 and |log fold change (FC)|≥1), there were 614 and 4211 DEGs extracted from GSE6214 and GSE34608-GPL6480, respectively. Following that, DE miRNAs were extracted from GSE34608-GPL7731, GSE29190, and GSE49951, respectively, yielding 174, 36, and 23 DE miRNAs. Two gene expression profile datasets yielded a total of 181 overlapping DGEs (135 downregulated genes and 46 upregulated genes) (Figure 2(b) and Table 1). Two overlapping DE miRNAs (two downregulated and zero upregulated) were derived from three miRNA expression profile datasets (Figure 2(c) and Table 2). The volcano plot of three gene expression profile datasets revealed a substantial difference between regular and tuberculosis patients (Figure 2(a)).
Figure 2 (a) Volcano plot showing downregulated and upregulated differentially expressed genes. (b) Venn diagram displaying the number of common differentially expressed genes (DEGs) between normal people and patients with TB. (c) Venn diagram for common differentially expressed miRNA (DE miRNA).
Table 1 Differentially expressed genes between the normal people and patients with TB.
DEGs, differentially expressed genes; TB, tuberculosis
Table 2 Differentially expressed miRNA between the normal people and patients with TB.
3.2. GO Functional and KEGG Pathway Enrichment analysis
DEGs that were consistent were submitted to R software for GO functional and KEGG pathway enrichment analysis. GO-CC showed that upregulated DEGs were mainly associated with ‘sarcoplasm’ and downregulated DEGs were mainly associated with ‘specific granule,’ ‘tertiary granule,’ ‘secretory granule lumen,’ ‘cytoplasmic vesicle lumen,’ ‘vesicle lumen,’ ‘secretory granule membrane,’ ‘specific granule membrane,’ ‘specific granule lumen,’ ‘tertiary granule membrane,’ ‘primary lysosome,’ and ‘azurophil granule.’ GO-BP showed that downregulated DEGs were primarily involved in ‘neutrophil-mediated immunity,’ ‘neutrophil activation,’ ‘neutrophil degranulation,’ ‘neutrophil activation involved in immune response,’ ‘response to lipopolysaccharide,’ ‘response to a molecule of bacterial origin,’ ‘defense response to a bacterium,’ ‘regulation of response to biotic stimulus,’ ‘humoral immune response,’ and ‘positive regulation of anion transport.’ GO-MF showed that downregulated DEGs were mainly related with ‘cysteine-type endopeptidase inhibitor activity,’ ‘virus receptor activity,’ and ‘exogenous protein binding’ (Figure 3(a)–3(d) and Table 3).
Figure 3 The gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of differentially expressed genes. The biological process (BP), bellular component (CC), and molecular function (MF) consist of GO enrichment analysis commonly. (a) The enrichment dot bubble of GO-BP analysis, (b) the enrichment dot bubble of GO-CC analysis, and (c) the enrichment dot bubble of GO-MF analysis derived from downregulated DEGs. (d) The enrichment dot bubble of GO analysis derived from upregulated DEGs. (e) The enrichment dot bubble of KEGG analysis derived from DEGs.
Table 3 Gene ontology analysis of common upregulated and downregulated DEGs.
GO, Gene Ontology; DEGs, differentially expressed genes; BP, biological process; CC, cellular component; MF, molecular function; Count, number of DEGs.
The result represented that upregulated DEGs were significantly associated with 5 KEGG pathways, including ‘hepatocellular carcinoma,’ ‘Kaposi sarcoma-associated herpesvirus infection,’ ‘phosphatidylinositol signaling system,’ ‘circadian entrainment,’ and ‘apelin signaling pathway’ (Figure 3(e), Table 4).
Table 4 KEGG pathway analysis of common upregulated DEGs.
KEGG, Kyoto Encyclopedia of Genes and Genomes; DEGs, differentially expressed genes; Count, number of DEGs.
3.3. Hub Gene Identification Based on the DEG PPI Network and Module Analysis
The STRING online tool was updated with 181 common DEGs (135 downregulated and 46 upregulated). Significant DEGs were identified in 135 of 181 DEGs. Significant DEGs were visualized in detail using the CytoScape software. The PPI network was constructed with 122 nodes and 241 edges (Figure 4(a)). Additionally, 59 of 181 DEGs were not part of the PPI network. The hub genes contained Cathelicidin Antimicrobial Peptide (CAMP), CEA Cell Adhesion Molecule 8 (CEACAM8), CD59 molecule (CD59), C-Type Lectin Domain Containing 5A (CLEC5A), Stomatin (STOM), Membrane-Spanning 4-Domains A3 (MS4A3), G Protein-Coupled Receptor 84 (GPR84), Late Endosomal/Lysosomal Adaptor, MAPK and MTOR Activator 3 (LAMTOR3), Defensin Alpha 4 (DEFA4), and CEA Cell Adhesion Molecule 1 (CEACAM1) (Figure 4(b)). The top significant module was defined from the PPI network using the MCODE plug-in based on its degree value. The top module featured fifteen nodes and forty-five corners (Figure 4(c)).
Figure 4 (a) Circular nodes show the DEGs. The blue nodes stand for downregulated DEGs, and red nodes stand for upregulated DEGs. (b) The 10 hub genes identified from the PPI network. (c) The top module extracted from the PPI network. The blue nodes stand for downregulated DEGs, and red nodes stand for upregulated DEGs.
3.4. Construction of the DEGs-DE miRNA Network
Because of the shortcomings of each dataset, hsa-mir-7 and hsa-mir-451 were each submitted to four accurate datasets for target gene prediction. Using the VENN online research website, miRNA target genes were integrated with DEGs. The identified data revealed that hsa-target mir-7 genes shared 19 common genes with DEGs (Figure 5(a)) and the hsa-goal mir-451 gene shared four genes with DEGs (Figure 5(b)). Meanwhile, CytoScape developed the DEGs-DE miRNA PPI network (Figure 6).
Figure 5 (a) VENN diagram identified the consistent genes between the DEGs and the target genes of miRNA-7. (b) VENN diagram identified the consistent genes between the DEGs and the target genes of miRNA-451.
Figure 6 The PPI network of DEGs and DE miRNA. The blue nodes indicate the downregulated genes. The red nodes indicated the upregulated genes. The blue V-shape nodes stand for downregulated miRNA.
4. Discussion
In recent years, a number of studies have been carried out to reveal the potential mechanisms of tuberculosis. The prevalence of tuberculosis has been steadily rising. Traditional studies have two flaws: a single genetic case and a limited cohort [Reference Huang, Sherman and Lempicki23]. Bioinformatics analysis of gene expression profile datasets of TB patients is now being used to screen more reliable data. A total of 181 common DEGs (135 downregulated genes and 46 upregulated genes) were extracted from two gene expression profile datasets in this study. The hsa-mir-7 and hsa-mir-451, two DE miRNAs, were extracted from four gene expression profile datasets. Both hsa-mir-7 and hsa-mir-451 were reduced in expression. The 181 DEGs underwent GO and KEGG enrichment analysis, allowing them to be classified into BP, CC, MF, and KEGG groups. Finally, a DEGs-DE miRNA PPI network was built. The PPI network was used to screen the top module and ten hub genes.
To determinate the underlying molecular mechanisms in the TB process, the most enriched BP, CC, and MF pathways were combined with downregulated and upregulated DEGs for comprehensive analysis, separately. GO analysis presented that upregulated DEGs were mainly related with sarcoplasm and downregulated DEGs were mainly related with neutrophil, antibacterial activities, primary lysosome, and exogenous protein binding. The comprehensive analysis of GO and hub genes demonstrated that hub genes mainly associated with neutrophil and primary lysosome signal pathways. KEGG analysis presented that upregulated DEGs were mostly involved in hepatocellular carcinoma, Kaposi sarcoma-associated herpesvirus infection, the phosphatidylinositol signaling system, circadian entrainment, and the apelin signaling pathway. Most of the abovementioned KEGG pathways play an important role in the immune response and apoptosis of Mycobacterium tuberculosis. Kaposi sarcoma-associated herpesvirus is well known to be involved in antiapoptosis, enhancement of cytokine production, and cell proliferation [Reference Gonçalves, Uldrick and Yarchoan24]. The previous study reported that MTB enhances bacterial virulence by avoiding host cell death [Reference Briken and Miller25]. Phosphatidylinositol is a lipid anchor, which as virulence factor and modulate host immune response in MTB [Reference Belcher Dufrisne, Jorge and Timóteo26]. Also, the apelin decreases mitochondrial apoptosis, mitochondrial ROS-triggered oxidative damage, and NF-κB activation to inhibit acute lung injury (ALI) and acute respiratory distress syndrome (ARDS) [Reference Yan, Wang, Cao and Chen27]. Circadian entrainment was also involved in innate immune response and photoperiod significant impairment and enhancing immune function [Reference Onishi, Maneval and Cable28].
According to a recent study, hsa-mir-7 impaired NF-κB and AKT transcriptional activity. The NF-κB and AKT pathways are important inflammatory-associated pathways and inhibited host cell autophagy in tuberculosis [Reference Zhou, Li and Zhou4, Reference He, Sun and Zhang29, Reference Ye, Yang and Huang30]. The hsa-mir-7 was thought to be involved in innate immune responses in a previous study [Reference Zhao, Chen and Guo31]. Meanwhile, hsa-mir-451 protects against cell death caused by ischemia/reperfusion injury, cancer, and myocardial I/R injury [Reference Gilfillan, Das, Shah, Alam and Bhandari32–Reference Kong, Feng and Yang34]. CAMP, CEACAM8, CD59, CLEC5A, STOM, MS4A3, GPR84, LAMTOR3, CEACAM1, and DEFA4 are among the ten hub genes discovered. They act as TB immune mediators and, additionally, as an active participant in the control of the inflammatory response [Reference Zhang, Liu, Liao and Zhao35–Reference Mukherjee, Sarkar-Roy, Wagener and Majumder40]. Meanwhile, macrophages are innate immune cells that serve as a guardian in the immune system's protection and response to microbial infection [Reference Hmama, Peña-Díaz, Joseph and Av-Gay41]. GPR84 is a member of the family of metabolic G protein-coupled receptors [Reference Wittenberger, Schaller and Hellebrand42]. GPR84 activation enhances bacterial adhesion and phagocytosis in macrophages. According to a recent study, elevated glucose levels promote GPR84 expression [Reference Recio, Lucy and Purvis43]. STOM is a member of the family of integral membrane proteins [Reference Hiebl-Dirschmied, Entler, Glotzmann, Maurer-Fogy, Stratowa and Prohaska44]. As previously mentioned, STOM's primary role is to regulate glucose transporter type 1 operation. Furthermore, elevated glucose levels augment macrophage anti-inflammatory function [Reference An, au and Ma45]. GPR84 and STOM are closely associated with the adaptive immune response, as determined by experimentation. CAMP is an antimicrobial peptide that is synthesized at the C-terminus of proteins [Reference Gudmundsson, Magnusson, Chowdhary, Johansson, Andersson and Boman46]. It eradicated MTB and slowed the progression of tuberculosis. CAMP is a vitamin D-related gene that would be triggered in response to vitamin D. It is essential for macrophages infected with MTB to benefit from the antimicrobial reaction induced by vitamin D. CAMP activation enhances the effectiveness of immune responses in tuberculosis [Reference Li, Østerhus and Johnsen47]. CEACAM1 and CEACAM8 are members of the class of heavily glycosylated carcinoembryonic antigens [Reference Berling, Kolbinger and Grunert48, Reference Kuroki, Arakawa, Matsuo, Oikawa, Nakazato and Matsuoka49]. The primary role of CEACM1 is antiapoptotic. It is involved in granulocyte survival [Reference Singer, Klaile and Scheffrahn50]. It is classified as an immunoregulatory checkpoint regulator. Previous research demonstrated that inadequate CEACAM1 resulted in inflammatory exacerbation [Reference Horst, Najjar, Wagener and Tiegs51]. CEACM8 is a responsive granulocyte biomarker. It performs two functions: it recognizes neutrophils and degrades extracellular matrix, thus stimulating the immune response [Reference Huang, Pan and Ma52]. MTB has been shown to control the immune system by controlling the macrophage cell cycle. It prevents macrophages from entering the interphase and gap phase 1 phases [Reference Tipgomut, Wongprommoon, Takeo, Ittiudomrak, Puthong and Chanchao53]. MS4A3 is a member of the family of membrane-spanning 4A genes [Reference Adra, Lelias and Kobayashi54]. MS4A3 is intimately linked to the cell cycle. As a consequence, the activation of the kinase-associated phosphatase (KAP) leads to cell cycle arrest. CLEC5A is a member of the superfamily C-type lectin/C-type lectin-like domain (CTL/CTLD) [Reference Adra, Lelias and Kobayashi54]. CLEC5A controls cell development by inducing apoptosis and arresting the cell cycle. Additionally, it plays a critical function in the production of inflammation and serves as a critical gene for the clinical management of pulmonary inflammation [Reference Teng, Chen and Hsu55]. DEFA4 is a part of the lipocalin family, which is involved in the transportation of vitamins, lipids, and steroid hormones [Reference Palfree, Sadro and Solomon56]. It inhibits bacterial growth by attaching to pathogenic bacterial siderophores. It plays a significant role in immune response and protection. In recent years, LAMTOR3 and CD59 have been extensively used to inhibit the development of a variety of cancers [Reference Weinstock, Anliker and von Zabern57, Reference Lunin, Munger, Wagner, Ye, Cygler and Sacher58]. CD59's primary role is to regulate immune cell activation throughout the tumor microenvironment [Reference Chen, Ding and Li59]. Additionally, the top module was defined through the PPI network of DEGs. The top module included a large number of genes associated with macrophages and the innate immune response. Previously published studies established that macrophages act as a host cell for MTB and are active in the innate immune response. MTB destroyed the macrophage's main immune systems, antigen introduction and intracellular killing [Reference Hmama, Peña-Díaz, Joseph and Av-Gay41].
The important DEGs were established using integrated bioinformation analysis. GPR84, STOM, CAMP, CEACM8A, MS4A3, LAMTOR3, DEFA4, CLEC5A, and CD59 were the ten hub genes that were filtered out. Two DE miRNAs, hsa-mir-7 and hsa-mir-451, were discovered. Some important pathways, including neutrophil, antibacterial activities, primary lysosome, exogenous protein binding, the apelin signaling pathway, Kaposi sarcoma-associated herpesvirus infection, and the phosphatidylinositol signaling system, are involved in TB. DEGs and DE miRNAs may be used as sensitive biomarkers in the clinical diagnosis of tuberculosis. The network of DE miRNA-DEGs revealed the molecular mechanism causing tuberculosis. More accurate clinical samples and experiments are needed to validate the cause of tuberculosis. In the future, bioinformatics research could be able to identify novel genes and pathways for use in genomic therapy for tuberculosis.
Data Availability
There are no files that need to be deposited in public repositories among the data in the manuscript. The manuscript contains all necessary information about the data collected.
Conflicts of Interest
The authors declare no conflicts of interest.
Authors’ Contributions
Shejie Shen contributed equally to this work.
Open access
