Bioinformatics Analysis of the Association between Ewing’s Sarcoma and Tuberculosis Comorbidity ()
1. Background
Tuberculosis (TB) is a multisystem disease caused by infection with Mycobacterium tuberculosis. TB has variable clinical manifestations, and the organs commonly involved include respiratory, gastrointestinal, lymphatic, skin, central nervous system and other systems [1] [2], among which nearly 75% of patients are mainly affected in the lungs [3], making lung disease the most common manifestation. Tuberculosis is one of the top 10 causes of death worldwide, and the World Health Organization (WHO) estimates that 1.8 billion people (about a quarter of the global population) are infected with Mycobacterium tuberculosis [4]. In 2020, 101 million people developed TB, compared to 106 million in 2021, and 1.5 million people died from TB in 2020, compared to 1.6 million in 2021. In addition, there was a 3.6% increase in TB incidence in 2021 compared to 2020, which is a reversal of the nearly 2% annual decline over the past 20 years [5]. It is well known that adolescence is a period of increased risk and burden of tuberculosis, and both the prevalence of mycobacterium tuberculosis infection and the incidence of tuberculosis have increased significantly [6], the reasons for which are not fully understood and may be related to the increased susceptibility of adolescents to mycobacterium tuberculosis.
Ewing sarcoma (ES) is the second most common primary malignant bone tumor in children, adolescents and young adults after osteosarcoma, accounting for 10% - 15% of all osteosarcomas [7]. Ewing’s sarcoma is highly aggressive and characterized by rapid tumor growth and active metastasis. Most Ewing sarcomas are bony and most often involve the bones, especially the long bones, pelvis, and ribs, with the most common occurrence occurring in the shaft of the long tubular bone [8] [9]. However, 20% to 30% of ES may be extrasolar and can occur in soft tissues, such as the thoracic wall or pleural cavity [10]. Translocations are commonly present in tumor cells, and somatic chromosomal translocations that fuse the Ewing sarcoma Breakpoint Region 1 (EWSR1) gene with members of the ETS transcription factor family have always been observed in Ewing sarcoma. EWSR1 fuses the Friend leukemia Integration 1 (FLI1) transcription factor to produce EWSR1-FLI1 fusion protein (also known as EWS-FLI1). Many studies have shown that EWS-FLI1 t (11; 22) (q24; q12) is the most common translocation, seen in 85% of patients with Ewing’s sarcoma (11; 22) (q24; q12) Translocation fuses the 5’ end of EWSR1 to the 3’ end of FLI1, producing the fusion protein EWSR1-fli1 [8] [11] [12]. EWSR1-FLI1 binds to chromatin and can change the epigenetic state, thereby altering gene expression and leading to the occurrence of tumors [13]. Bone tumor microenvironment (TME) plays an important role in tumor genesis, development and metastasis, and its cell composition and physical properties can affect the biological behavior of cancer cells. In the process of tumor growth, hypoxia is an important regulator of angiogenesis [14], and the transcription factor hypoxia-inducible factor 1 (HIF-1) is a major regulator of hypoxia. For example, in osteosarcoma and Ewing’s sarcoma [15], under hypoxic TME, HIF-1 directly activates the transcription of pro-angiogenesis genes, including VEGF, and promotes angiogenesis [16] [17]. In addition to regulating angiogenesis, hypoxia is directly related to the progression, metastasis, and patient prognosis of osteosarcoma. Thus, hypoxia enhances the expression and transcriptional function of the EWS/FLI1 fusion protein in a HIF-1-dependent manner in Ewing’s sarcoma cells, thereby increasing its aggressiveness [18]. In addition, HIF-1 can also activate signaling pathways such as Akt and STAT3 signaling pathways to promote the proliferation and migration of osteosarcoma cells [19] [20]. Studies have shown that the expression of HIF-1α remains elevated in granulomas of patients with active tuberculosis [21]. In conclusion, TB can be used as a risk factor for ES metastasis, but whether the two are related at the genetic level remains unclear. Therefore, it is urgent to explore new therapeutic targets to provide theoretical basis for clinical treatment of ES complicated with TB.
In summary, through bioinformatics analysis of ES and TB collected in Gene Expression Omnibus (GEO) and their control population gene chip data, this study studied the common key genes and their molecular networks and explored the possible molecular biological functions involved. It provides some theoretical reference for revealing the potential correlation and molecular mechanism of the two.
2. Materials
2.1. Obtaining Data
We searched for “Ewing’s Sarcoma” and “Tuberculosis” in gene expression omnibus (GEO) of the National Center for Biotechnology Information, respectively. The original data of ES patient data chip GSE17674 and TB patient data chip GSE57736 were obtained. Differential expression analysis of ES and TB disease mRNA microarray data was conducted based on R-language limma package screening for differentially expressed genes [22]. The logarithmic absolute value of the test statistic P < 0.05 and the fold change (FC) |log FC| > 2 were set as screening conditions. They were screened for differentially expressed mRNAs (DEmRNAs). The intersection of the differential genes of ES and NFH diseases was used to obtain the common target genes of both
2.2. GO Function and KEGG Signal Pathway Enrichment Analysis
The R clusterProfiler package was used for GO and KEGG enrichment analysis, and the GO function of the two common DEmRNAs was enriched to analyze the biological processes involved in them, and the KEGG signaling pathway was enriched at the same time.
2.3. PPI and Core Module Analysis
Degs were analyzed using the search tool for the retrieval of interacting genes/proteins (STRING) 11.0 online tool. The calculation results of STRING were imported into Cytoscape 3.9.1 software, and the plug-ins CytoNCA, Cytohubba and MCODE were used as the protein interaction network diagram to co-express core genes and modules.
2.4. Hub Gene Localization
Tissue localization of HUB gene was performed using online database BioGPS (http://biogps.org/#goto=welcome). The highest intertissue expression level was greater than twice the second expression level, indicating that HUB gene can be used as a common expression marker gene for ES and OP.
3. Results
3.1. Identification of Comorbidity mRNA
The R language limma package was used to analyze the differential expression of gene chips for the two diseases respectively, and it was found that the expression levels of most genes were basically consistent, indicating that the data was suitable for further analysis, as shown in Figure 1.
Figure 1. Box map of gene expression. Note: A is GSE57736 data chip box diagram, bGSE17674 data chip box diagram; Blue is the normal group, red is the experimental group, and the black lines in the figure are basically the same straight line, indicating that the total gene expression is consistent, which can be analyzed in the next step.
By analyzing the above data set, a total of 300 differential genes were identified between TB patients and normal individuals, and a total of 3029 differential genes were identified between ES patients and normal individuals (Figure 2). We identified the intersection of these two datasets and obtained a total of 55 co-expressed genes, shown in Venn diagram (Figure 2).
Figure 2. Wayne diagram of ES and NFH intersection genes.
3.2. Intersection Gene GO and KEGG Enrichment Analysis
GO and KEGG pathway enrichment analysis of ES and NFH common target genes was performed by using R cluster Profiler package. The signaling pathways in Biological Processes (BP), Cellular Components (CC), Molecular Functions (MF) and KEGG were analyzed respectively. As can be seen:
1) BP function is mainly enriched in the positive regulation of interleukin−8 production, viral life cycle, viral process, the biological process involved in interaction with host and other aspects; 2) CC function is mainly concentrated in the intrinsic component of endoplasmic reticulum membrane, an integral component of endoplasmic reticulum membrane, etc. 3) MF is mainly concentrated in nuclear retinoid X receptor binding, nuclear receptor binding, nuclear receptor coactivator activity, and nuclear retinoic acid receptor binding, etc. 4) There was no significant enrichment of the KEGG pathway, as shown in Figure 3.
3.3. Construction of Common Expression PPI and Prediction of Core Genes
Using STRING online tool, 55 common DEmRNAs were analyzed, and the minimum interaction score was set as 0.4, and the interaction diagram consisting of 15 nodes and 62 connections was obtained. The PPI interaction diagram was obtained after the calculation result of STRING was imported into Cytoscape 3.9.1 software, as shown in Figure 4. Generic expression gene modules were obtained by MOCDE algorithm in Cytoscape plugin (Figure 5). Three Hub genes were obtained through 10 algorithms including MCC, DMNC and MNC in CytoHubba plug-in (Figure 5), namely CD3D, LCK and KLRB1.
Figure 3. GO functional enrichment analysis and KEGG pathway analysis. Note: (a): BP, (b): CC, (c): MF.
Figure 4. PPI network of differential genes.
Figure 5. Core module of Hub gene and core genes.
4. Discussion
Immune cell infiltration is an indicator of host immune response to cancer antigens [23]. In addition, functional enrichment analysis showed that CD3D plays an important role in inflammation and immune responses. High expression of CD3D is significantly enriched in antigen processing and presentation, cell adhesion molecules, cytokine-cytokine receptor interactions, and B-receptor signaling pathways. Antigen processing and presentation are essential for a successful humoral response [24]. Increased antigen presentation can lead to adaptive immune responses, including antibody production [25]. In addition, adhesion molecules are critical for immune cell homing to inflamed tissues and lymphatic organs and play an important role in immune homeostasis in both health and disease [26]. In addition, cytokine-cytokine receptor interactions and B-cell receptor signaling pathways are closely related to the tumor immune microenvironment and are crucial for tumor immune response [27]. In addition, the JAK/STAT pathway was enriched in the group with high CD3D expression. The JAK/STAT cytokine signaling pathway is associated with proliferation, immune and inflammatory responses [28]. Most cytokines involved in the immune response use the JAK/STAT signaling pathway [29]. These findings confirm that elevated levels of CD3D expression may be caused by an anti-tumor immune response that recruits and activates various types of immune cells.
Lck expression has been detected in a variety of solid cancers, including breast cancer [30]-[33], colon cancer [34]-[36] and lung cancer [37] [38]. These observations lead to the hypothesis that Lck may have a cancer-promoting function and may therefore represent a potential diagnostic biomarker and therapeutic target for solid cancers. In fact, Lck inhibitors are not only used to treat leukemia, but also various solid cancers [39].
The expression of Lck in patients with cholangiocarcinoma is associated with tumor recurrence [40]. Treatment of bile duct cancer cells with the SFK inhibitor dasatinib showed tumor inhibition in both in vitro and in vivo models. In this form of cancer, the molecular pathologic mechanisms of Lck in tumor development have been studied more closely. Transcriptional coactivator Yes associated protein (YAP), which can be induced by PDGF receptor signaling, plays an important role in the pathogenesis of cholangiocarcinoma [41] and has been found to be regulated by Lck. Figure 2) [40]. Lck phosphorylation of YAP (at Y357), dasatinib and SirNA-targeted Lck knockout blocked YAP tyrosine phosphorylation. In addition, treatment of bile duct cancer cells with dasatinib induced the redistribution of YAP from the nucleus to the cytoplasm and down-regulated the expression of YAP target genes involved in carcinogenesis [40].
Lck also appears to play a role in cancer stem cells (CSC) in endometrioid cancer models [40]. In this form of cancer, the cell surface complement inhibitor CD55 regulates self-renewal and cisplatin resistance in a complement-independent manner. Interestingly, CD55 is signaled by the transmembrane adaptor protein LIME (Lck interacting Transmembrane adaptor) (Figure 2) [40], which we have determined interacts with Lck [41]. Inhibition of LIME expression in CSC leads to reduced levels of active Lck and impaired CD55-mediated signaling. In addition, Lck plays a role in the upregulation of genes associated with DNA repair, including MLH1 and BRCA1, and in the resistance of CSC to cisplatin. In fact, treating these cells with the Fyn/Lck inhibitor saracatinib or inhibiting Lck with shRNA can make CSC sensitive to cisplatin [40]. Taken together, these data suggest that Lck represents a new and important target for tumor therapy by modulating migration, tumor growth, and cancer dryness.
KLRB1 is down-regulated at both gene and protein levels in most tumors, which is consistent with previous findings [42]. Interestingly, KLRB1 expression in kidney cancer tissues (KIRC and KIRP) showed the opposite result. A previous study reported that KLRB1 is generally associated with a good prognosis [43]. Our results also confirmed that patients with higher KLRB1 expression achieved longer survival in different cancers such as breast cancer, melanoma and thyroid cancer. Unlike these cancers, patients with low-grade gliomas with high KLRB1 expression have a poorer prognosis. For patients with hepatocellular carcinoma, co-expression of CD161 and IL-7R can enhance the expression of IL-2, TNF-αα and perforin, which can improve prognosis [40]. In addition, we found that patients with liver cancer expressing KLRB1 had longer PFI. Recent studies on oropharyngeal squamous cell carcinoma have shown that down-regulation of CD161 can lead to immune escape of cancer cells [44], while IL-17 and IFN-γ produced by CD161 + T cells can reduce tumor load and improve OS [45]. In addition, another study reported that KLRB1 was also associated with favorable outcomes in non-small cell lung cancer [46], which is consistent with our findings. We also investigated the relationship between patient clinical phenotype and KLRB1 expression. Notably, KLRB1 expression was significantly associated with age in 11 types of cancer, which may help guide clinical treatment protocols and drug selection. Our study also found that the higher the tumor stage, the lower the KLRB1 expression level in LUAD, SKCM, THCA and TGCT. These results all suggest that KLRB1 can be used as prognostic-related predictors of different cancers and may play a functional role in a variety of tumors.
In this study, bioinformatics methods were used to select suitable gene data chips through GEO database, and R language was used to screen differentially expressed genes related to the incidence of ES and ONFH. The differentially expressed genes in the two datasets were mapped with volcano maps, and 668 common differentially expressed genes were obtained by intersection of the co-expressed genes with Venn map. GO, KEGG and KEGG enrichment analysis of differential genes showed that differential genes were mainly enriched in muscle structure development and actin cytoskeleton tissue, and were related to striated muscle contraction pathway and RUNX3 regulation of CDKN1A transcription. As the therapeutic target of ES, EWS/FLI is currently only found to be expressed in ES, and RUNX protein is closely related to the Runt domain of EWS/FLI. There are three subtypes of RUNX, RUNX1 and RUNX3 can promote the occurrence of various tumors, and RUNX3 can be detected in all ES tumors. In this study, it was found that the signaling pathways involved in DEGs mainly focused on RUNX3’s regulation of CDKN1A transcriptional signal transduction pathways and other aspects, which is consistent with the above-mentioned skeletal immune mechanisms and provides a direction for further research on which factors in the immune system are involved.
CD5 + B cells (also known as B1 cells) are a group of congenital B cells responsible for the production of natural antibodies, rapid humoral immunity and immune regulation, which can produce IL-10 and inhibit T cell response, and CD5 cells are significantly increased in TB patients [47]-[49]. However, CXCR6 was not significantly associated with TB and ES, Shahd Ezzeldin et al. also reported that CD5 was significantly expressed in the metastasis of Ewing’s sarcoma in early childhood, and the Hub gene we screened also contained CD5 genes, indicating that CD5 can be used as an early landmark molecule of TB and ES, but how the two are interconnected through CD5 remains to be further studied.
5. Summary and Outlook
All of these genes are involved in the occurrence and development of ES and ONFH to varying degrees, especially many of them are related to inflammation and immune function abnormalities, which provides new ideas for further research on the role and correlation of the immune system in the pathogenesis of the two.
In summary, the target genes related to both diseases were explored through data mining, which on the one hand provided a new early warning indicator for ES combined with ONFH, and on the other hand, provided a reference for the next step to reveal the mechanism of the association between the two.