Screening of Potential Key Biomarkers for Ewing Sarcoma: Evidence from Gene Array Analysis

Introduction

ES is a poorly differentiated and aggressive tumor, which is called the three most common bone tumors together with osteosarcoma and chondrosarcoma.^1–4 The origin of the tumor is not clear. Pathologically, a tumor composed of small and round cells is considered to be original in the neuroectoderm through electron microscopy and immunohistochemical analysis. Mesenchymal-stem cells, which are deemed to originate from primitive bone marrow, usually have the following features: Modified by fusion transcripts involving EWS-ERG or EWS-FLI1 genes.^5–8 ES’s clinical manifestations are not specific, mainly manifested as local masses, pain, etc.; X-rays manifested as osmotic destruction of the backbone, with onion skin-like and needle-like periosteal reactions. There are ethnic differences in its incidence, that European and American populations is higher than that in Asian populations. There are also differences in treatment effects. Although researchers have worked hard to improve the treatment and diagnosis plan, people with metastases or relapses at the time of diagnosis still show poor prognosis.⁹ The wide application of bioinformatics analysis and microarray technology at the genomic level, especially the screening of gene heredity and variation, is conducive to our screening of functional pathways and differentially expressed genes involved in the ES mechanism. However, the independence and false positives of microarray analysis make it difficult for us to obtain reliable results. Therefore, to screen potential biomarkers of ES, we downloaded and analyzed 3 mRNA data sets from the Gene Expression Omnibus (GEO) database to obtain the DEGs between non-tumor tissue and tumor tissue. Subsequently, a series of analyses were carried out, such as protein interaction (PPI) network, gene ontology (GO), pathway enrichment analysis, immune analysis, survival analysis, etc., to help us understand the molecular mechanism of carcinogenesis. Finally, 629 DEGs and 10 hub genes were screened out.

Materials and Methods

Gene Chip Data

GEO¹⁰ is a public genomics database with a high storage capacity in the entire microarrays dataset and gene expression data. Three gene datasets were downloaded from GEO, which are GSE68776,¹¹ GSE45544^12,13 and GSE17674.^14–16 Transform the platform probe information into corresponding gene symbols. The GSE68776 dataset contains 32 tumor biopsy specimens and 33 normal adult tissue samples. GSE45544 contains 8 Ewing sarcoma patient samples (MuET-x), 6 Ewing sarcoma patient samples (TUMK00XX) and 22 normal tissues (PBMC, spleen, thymus, stomach, uterus, fetal brain, fetal liver) sample. GSE17674 contains the Ewing sarcoma family from 44 tumor patients and 18 normal muscle samples.

DEG Recognition

GEO2R is an interactive web tool. Use GEO2R to compare multiple sets of samples to identify DEG across experimental conditions. The false discovery rate of Benjamini and Hochberg and the adjusted P-value (adj.P) can be used to limit false positives and find statistically significant genes. Adjust the probe set, delete the probe set without the corresponding gene symbol, and delete the gene with multiple probe sets. LogFC <-1 or logFC (fold change)> 1, adj.P <0.01, the above three conditions are considered to be statistically significant.

DEG’s GO and KEGG Enrichment Analysis

DAVID (The Database for Annotation, Visualization and Integrated Discovery) (v6.8)¹⁷ helps researchers understand the biological significance behind many genes through a comprehensive set of functional annotation tools. KEGG (https://www.kegg.jp) is a resource database, molecular data sets generated by genome sequencing and other high-throughput experimental technologies, using molecular-level large-scale information to understand biological systems.¹⁸ The gene ontology (GO) knowledge base is the world’s largest knowledge base of gene functions and information sources. The information in this knowledge base is the basis for the computational analysis of genetics and molecular biology experiments in biomedicine, that is readable by both humans and machines.¹⁹ DAVID online database was used for biological analysis of DEG, then, P<0.05 was considered statistically significant.

Network (PPI) and Module Analysis

STRING (V11.0) (http://string-db.org) is a database of predicted protein–protein interactions. The interactions include physical and functional associations, that originate from the transfer of information between predictive and computational organisms, as well as the interaction of other (main) database aggregations.²⁰ Using STRING to analyze the function of proteins, that may provide information about the mechanism of related diseases. In our study, the STRING database was used to construct the DEG protein interaction network, and the interactions with a comprehensive score (supported by data)>0.4 were considered to be statistically significant. Cytoscape (v3.8.0) is an open-source application software used to visualize complex networks in bioinformatics and integrate network attribute data.²¹ Cytoscape’s plug-in MCODE (version v2.0.0) is used to cluster topology on a given network to find highly interconnected areas.²² CytoHubba (v0.1) is another plug-in in Cytoscape. It uses a variety of topology algorithms to explore and predict important nodes and subnets in the network.²³ Use Cytoscape to make a protein interaction network (PPI), and then use MCODE to analyze the important modules. Select parameters: MCODE score>5, k score=2, degree cutoff=2, maximum depth=100, node score cutoff=0.2. Then filter out the key sub-networks, namely hub genes (a total of 10) through the cytoHubba in the module. Finally, use DAVID to analyze the hub genes by KEGG and GO.

Selection and Analysis of Key Genes

Oncomine is a tumor gene chip database that can be used to analyze gene expression differences and predict co-expressed genes. It can also be classified according to clinical information such as issue type, tumor stage and grade.²⁴ Using the Oncomine database, under the Ewing’s sarcoma classification, the top ten co-expressed genes of each hub gene were selected. After being classified, a total of 85 genes including genes were obtained. Draw the co-expression network through String and Cytoscape. BiNGO (Biological Network Gene Oncology) (v3.0.3) is a web tool used to find that certain gene ontology (GO) are statistically overrepresented in a certain set of genes or biological networks.²⁵ We used BiNGO, a plug-in of Cytoscape, to analyze the biological process of the hub gene and made a network interaction map. UCSC Xena is a multi-omics online exploration tool for clinical or phenotypic data.²⁶ The hierarchical clusters and sample types of hub genes were constructed by UCSC Xena. TIMER is a broad resource tool for systematic analysis of cancer immune infiltration.^27–29 TIMER (v2.0) was used to analyze the immune infiltration of hub genes, then, showed the correlation between gene expression and immune analysis, the correlation between gene expression and clinical. The hub gene sample code needs to be downloaded from the TCGA (The Cancer Genome Atlas) (https://portal.gdc.cancer.gov/) database for Immune infiltration analysis. OncoLnc (http://www.oncolnc.org) is a tool for interactively exploring the relevance of survival.³⁰ Oncolnc was used to analyze the overall survival of hub genes. Seven genes were obtained through the results of survival analysis and immune correlation analysis. Then use Oncomine to analyze and compare the expression of 8 genes in different sarcoma types, and 3 key genes were screened out. Finally, analyze the correlation between the key gene’s expression and Ewing’s sarcoma tumor grade, stage, EWSR1-FLI1 fusion and TP53 mutation on Oncomine.

Results

Identify DEGs in Microarrays

After standardizing the microarray results, DEG was determined (8207 in GSE68776, 2115 in GSE45544, and 28893 in GSE17674) (Figure 1A). As is shown, the picture of Venn (Figure 1B), the overlap between these 3 datasets contains 629 genes, consisting of 423 down-regulated and 206 up-regulated genes between ES tissue and non-cancer tissue.

Figure 1 (A) The volcano plot shows the relationship between the fold change and the P value in each group of microarrays. Red is an up-regulated gene, and blue is a down-regulated gene. Adjust P value <0.05. (B) Venn diagram, select DEG fold change> 2, P value <0.01, 3 microarray data sets show that 629 genes overlap. (C) 629DEG’s PPI network.

GO and KEGG Enrichment Analysis of DEG

The functions and pathways of DEGs were analyzed by DAVID. GO analysis results showed that the biological process (BP) changes of DEGs were significantly obvious in protein phosphorylation, regulation of related pathways by p53 mediator, intercellular adhesion, and mitotic cell cycle. The changes in molecular function (MF) mainly includes protein kinase activity, carbohydrate binding, cadherin binding involved in cell-cell adhesion, oxidoreductase activity, and protein serine/threonine kinase activity. The cell composition (CC) changes of DEGs mainly includes sarcoma, spindle pole, and cell adhesion junctions. KEGG pathway analysis shows that DEG is mainly rich in cell cycle, glioma, cancer pathway and P53 pathway (Table 1).

Table 1 Enrichment Analysis (629DEG’s)

PPI Network and Module Construction

The PPI network (Figure 1C) and the most important modules (Figure 2A) were built through Cytoscape. CytoHubba was then used to screen out the genes (Figure 2B), of which CCNB2 is an up-regulated gene, and the rest are the down-regulated genes. The functional analysis of hub genes by DAVID showed that the up-regulated genes in this module are mainly enriched in cell division and cycle, mitosis, protein binding, etc., and the down-regulated genes are mainly enriched in protein kinase activity, protein localization to kinetochore, protein phosphorylation, chromosome centromere regions, etc. (Table 2).

Table 2 Enrichment Analysis (10DEG’s)

Figure 2 (A) The important module is obtained from a PPI network with 48 nodes and 782 edges. Red is an up-regulated gene; light blue is a down-regulated gene. (B) Using cytohubba’s MCC algorithm to screen out the top ten genes in the module, we can see the connection between them. The weight of these genes in the module network becomes higher as the color darkens.

Selection and Analysis of Key Genes

We screened out 10 hub genes. The names and functions of these genes are shown in Table 3. Oncomine, String and Cytoscape were used to plot the co-expression network of these genes (Figure 3A) and biological process analysis diagram (Figure 3B). Through hierarchical clustering, we can see that pivot genes can basically distinguish ES samples from non-tumor samples (Figure 4A). The correlation between gene expression and clinical results (Figure 4B) showed that among 260 patients with sarcoma, the expression levels of AURKB, BUB1, CDCA8, CDK1, MAD2L1, and MCOPH were higher. The correlation between immune infiltration and gene expression is shown in Figure 5. Survival analysis showed that the overall survival rate of patients with CDCA8, NCAPH, CDK1, MAD2L1, CCNB2, and FANCI was relatively high (>30%) (Figure 6). Combining the results obtained in Figures 4B, 6 and 8 genes (AURKB, BUB1, CDCA8, CDK1, MAD2L1, FANCI, NCAPH and CCNB2) were screened out. Through the screening of 8 genes, the key genes are CDCA8, MAD2L1 and FANCI (Figure 7). Finally, the analysis lists several expression patterns of the three genes that are relatively important in ES (Figure 8). Based on the above analysis, it is suggested that CDCA8, MAD2L1 and FANCI may have a great importance in the occurrence or progression of ES.

Table 3 Annotation (10DEG’s)

Figure 3 PPI-network and Bioanalysis chart. (A) Co-expression of hub genes. The blue nodes are co-expressed genes, and the red nodes are hub genes. (B) Some biological process diagrams of hub genes. The size of the node refers to the ontology of the gene. The color depth of the node refers to the adjusted P value.

Figure 4 (A) Hierarchical clustering: down-regulated genes are marked in blue, and up-regulated genes are marked in red. Sample type: Different colors in the picture correspond to different sample types. The sarcoma sample data (n=271) comes from the TCGA database. P<0.05 is considered statistically significant. (B) The correlation between gene expression and clinical outcome. (SARC, n=260).

Figure 5 The relationship between gene expression and immune infiltration.

Figure 6 Survival analysis. CDCA8, NACPH, CDK1, MAD2L1, CCNB2 and FANCI have a relatively high overall survival rate (>30%). P<0.05 was considered statistically significant.

Figure 7 Comparison of genes in multiple analyses. FANCI: Cancer Type: Sarcoma. 1. Neale Multi-cancer. 2. Neale Multi-cancer 2; MAD2L1: Cancer Type: Sarcoma. 1. Barretina CellLine 2. 2. Neale Multi-cancer. 3. Rothenberg CellLine; CDCA8: Cancer Type: Sarcoma. 1. Barretina CellLine 2. 2. Neale Multi-cancer. 3. Neale Multi-cancer 2. 4. Rothenberg CellLine; NCAPH: Sarcoma Type: Ewing’s Sarcoma. 1. Baird Sarcoma. Cancer Type: Sarcoma. 2. Neale Multi-cancer. 3. Neale Multi-cancer 2; AURKB: Cancer Type: Sarcoma. 1. Barretina CellLine 2. 2. Bittner Multi-cancer. 3. Neale Multi-cancer. 4. Wooster CellLine; CCNB2: Sarcoma Type: Ewing’s Sarcoma. 1. Bittner Sarcoma. Cancer Type: Sarcoma. 2. Neale Multi-cancer. 3. Rothenberg CellLine; BUB1: Cancer Type: Sarcoma. 1. Bittner Multi-cancer. Sarcoma Type: Ewing’s Sarcoma. 2. Henderson Sarcoma; CDK1: Sarcoma Type: Ewing’s Sarcoma. 1. Baird Sarcoma. Cancer Type: Sarcoma. 2. Barretina CellLine 2. The rank for a gene is the median rank for that gene across each of the analyses. The p-value for a gene is its p-value for the median-ranked analysis. In the heat map of 8 genes, CDCA8, MAD2L1 and FANCI have more or higher expressions.

Figure 8 The relationship between the expression of genes (CDCA8, MAD2L1 and FANCI) and clinical Ewing’s sarcoma grade, staging, EWSR1-FLI1 pathway and tp53 mutation pathway. MAD2L1 has higher expression in the four cases, while CDCA8 and FANCI have higher expression in the grades of ES.

Discussion

In the past few years, constantly updated technologies have made contributions to obtain more accurate diagnosis and finer classification of diseases, that have highlighted the importance of molecular heterogeneity between and within tumors and secondary genetic changes. Some new large-scale sequencing technologies have helped exploration the genomic pattern of Ewing’s sarcoma, proving that ES is a cancer related to gene fusion, and belongs to the transcription factor family of EWSR1 and ETS (especially FLI1) repeated translocations.^31–33 Similar to other cancers caused by gene fusions, ESFT has a low somatic mutation load, indicating that the fusion of EWSR1-ETS has an advantage as a driver. We also need to pay attention to mutations that lead to cancer recurrence, including several tumor suppressors such as STAG2, TP53 and CDKN2A. More notably, ESFT shows obvious mutations in genes involved in the kinase signaling pathway.³⁴ Targeted fusion is achieved through changes in the related transcriptome or epigenome, It has always been an active and challenging field.

In this study, we analyzed 3 mRNA microarray data sets and screened out 629 differentially expressed genes (423 down-regulated and 206 up-regulated genes). Then the hub gene was identified for functional enrichment analysis to explore the interaction between DEGs. The down-regulated genes are mainly enriched in kinase activity, protein localization to kinetochore, protein phosphorylation, chromosome centromeric regions, etc. The up-regulated genes are mainly enriched in protein binding, cell cycle and cell division, etc. We all know that mitosis and cell cycle processes have a great influence on tumor dysregulation.^35–37 Many tumor processes also involve protein kinases. For instance, in various human malignant tumors, it is very common for Protein kinases to facilitate signal transduction of MAPK/ERK, making it the key molecule in the dysregulation of this pathway. Consequently, the MAPK becomes a key signaling pathway involved in tumorigenesis and regulating the uptake of glucose by malignant cells.³⁸ Phosphorylation is a relatively common post-translational modification (PTM), which regulates many cell functions under healthy conditions. However, changes in the phosphorylation pathway can lead to serious diseases, especially cancer.³⁹ The instability of chromosomes is considered a portent of bad tumors; and the strict regulation of the region where the kinetochore and chromosome centromere are combined also reveals the possible role of RNA on the centromere.⁴⁴ In short, these theories and our results finally coincide.

Among the 10 genes we screened, CDCA8, MAD2L1 and FANCI are of great significance. Some studies have shown that overexpression of CDCA8 contributes to the proliferation of tumor cells, such as colorectal cancer and lung cancer cells.^40,41 In addition, high CDCA8 expression has also been found to indicate a poor prognosis for gastric cancer.⁴² However, there is currently no evidence suggesting that CDCA8 is a specific marker for the above tumors. Based on RNA sequencing, bioinformatics analysis, and protein levels, CDCA8 silencing can down-regulate the levels of cyclin B1 and p-cdc2 and explain how it induces G2/M arrest.⁴³ Selective suppression of the CDCA8-AURKB pathway may also become an effective way to treat cancer.⁴⁰ According to our research, CDCA8 may be a potential treatment target for ES, but it has yet to be confirmed. MAD2L is associated with mitochondrial checkpoints. Defects in the control of mitotic checkpoints are thought to contribute to chromosomal instability and aneuploidy.⁴⁴ Some studies suggest that MAD2L1 is associated with endometrial cancer, gastric cancer, and liver cancer. In the case of endometrial cancer with lymphatic metastasis, MAD2L1 is significantly overexpressed.⁴⁵ In gastric cancer, the up-regulation of miR-30a-3p can reduce the expression of MAD2L1, which in turn inhibits its cell proliferation.⁴⁶ In addition, certain changes in HCC cells (such as proliferation and migration)⁴⁷ can be controlled by restraining MAD2L1 and miR-200c-5p. Therefore, defects in mitotic checkpoints may help to increase the sensitivity of certain tumor cells to mitotic spindle inhibitors. We highly suspect that MAD2L may affect the mitosis process of ES cells, according to the analysis of MAD2L1 in Figure 8. In the study of FANCI gene expression, most of them are related to malignant hematological diseases of children, but reports related to lung cancer and breast cancer are also not rare. Children’s FANCI gene was detected in 94 genes related to hematological malignancies. The sequence analysis and testing of this gene was done by the US INVITAE laboratory.⁴⁸ In lung adenocarcinoma tissues, compared with neighboring normal tissues, the mRNA and protein of FANCI are overexpressed.⁴⁹ According to a recent study, FANCI may promote cell metabolism when DNA repair is not required.⁵⁰ Importantly, a study verified that there is a certain relationship between FANCI and sarcoma. They observed that 9 of 66 patients (13.6%) had 10 types (including FANCI and TP53) cancer-related genes contain at least one germline mutation gene that can predict disease. They also successfully verified that the two variants (FANCE and FANCI) did not show a loss of heterozygosity, but structural data indicated that the positions of these two genes were related to important protein interactions.^51–53 This finding further supports our belief that FANCI may be one of the potential key biomarkers of Ewing’s sarcoma.

Document retrieval results indicate that the interaction between Ewing’s sarcoma and hub genes (BUB1, CDK1, CCNB2, AURKB, PLK4, CENPF, and NCAPH) has not been widely reported. Overexpression of BUB1 plays an important role in breast cancer.⁵⁴ Dysregulated mitotic kinases are often associated with cancer. In addition, miR-10b, which is highly related to the BUB1 gene, was down-regulated in osteosarcoma samples.⁵⁵ Some scholars proposed that CDK1 may be related to breast cancer and explained that circMETTL3 affects breast cancer through the circMETTL3/miR-31-5p/CDK1 axis.⁵⁶ In pancreatic cancer cells driven by K-Ras, knocking out cell cycle regulators CDK1 (or CDK2) or transcription regulators CDK7 (or CDK9) is as effective as knocking out K-Ras.⁵⁷ The high expression of CCNB2 has some effects on bladder cancer, lung cancer and colorectal cancer.^58,59 The expression of AURKB is often associated with lung cancer and liver cancer. They believe that the interaction between RB1 and AURKB may be related to the small cell lung cancer driver.⁶⁰ PLK4 is associated with lung cancer and tumors of the reproductive system, nervous system, digestive system, etc. According to data, Plk4 is widely over-expressed in tumor samples from cancer patients.⁶¹ The abnormal expression of CENPF is related to many malignant tumors, including liver cancer and breast cancer.^62–65 Studies have shown that the interaction between LANA and CENPF has an impact on the genome of Kaposi’s sarcoma-associated herpes virus (KSHV).⁶⁶ The expression of NCAPH is related to breast cancer, lung cancer and colon cancer. In patients with hormone receptor-related breast cancer, the upregulation of NCAPH indicates a poor prognosis.⁶⁷ The expression of NCAPH and Mcl-1 suggests that patients with non-small cell lung cancer may have a poor prognosis.⁶⁸ Many NCAPH mutations have been found in patients with colon cancer.⁶⁹

In summary, this study aims to screen and identify DEGs that may be involved in the carcinogenesis or progression of ES. 629 DEGs and 10 hub genes were screened using the public gene expression comprehensive dataset. Furthermore, survival analysis suggested that CDCA8, NACPH, CDK1, MAD2L1, CCNB2 and FANCI were significantly associated with patient prognosis. Especially CDCA8, MAD2L1 and FANCI which may provide new clues for studying the mechanism of ES from the perspective of bioinformatics. However, further experiments and exploration are still needed to verify these results.