Bioinformatics analysis of Ewing's sarcoma: Seeking key candidate genes and pathways.

Ewing's sarcoma (ES) is the second most common bone tumor among children and adolescents worldwide. However, the genes and signaling pathways involved in ES tumorigenesis and progression remain unclear. The present study used two gene-expression profile datasets (GSE17674 and GSE31215) to elucidate key potential candidate genes and pathways in ES. Differentially expressed genes (DEGs) were identified and a functional enrichment analysis was performed. A protein-protein interaction (PPI) network was constructed, and the most significant module in the PPI network was selected from the Search Tool for the Retrieval of Interacting Genes/Proteins database. A total of 278 genes were identified by comparing the tumor samples with non-cancerous samples; these included 272 upregulated and 6 downregulated genes. The pathway analysis demonstrated significant enrichment in the positive regulation of transcription in the DEGs coding for RNA polymerase II promoter, plasma membrane and chromatin binding pathways in cancer in general. There were 269 nodes and 292 edges in the PPI network. Finally, MYC, IGF1, OAS1, EZH2 and ISG15 were identified as the hub genes according to the degree levels. The survival analysis revealed that EZH2 is associated with a poor prognosis in patients with ES. In conclusion, the DEGs, associated pathways and hub genes identified in the present study help elucidate the underlying molecular mechanisms of ES carcinogenesis and progression, and provide potential molecular targets and biomarkers for ES. Copyright: © Zhang et al.

Introduction

Ewing's sarcoma (ES) is the second most common primary bone tumor among children and adolescents worldwide (1). In recent years, with the increase in ES research, significant progress has been made towards the diagnosis, treatment, as well as prognosis of this disease. The current standard treatment for ES involves a five-drug chemotherapy regimen (vincristine, doxorubicin, cyclophosphamide, ifosfamide and etoposide), and local resection or radiation therapy, or both (2). Although these treatments improve the survival rate of patients, the high rates of recurrence and metastasis in patients with ES necessitate the urgent development of new diagnostic strategies and therapeutic agents to improve patient prognosis. As the molecular mechanisms of ES tumorigenesis and progression are not yet entirely understood, there remains to be a number of unresolved issues in the diagnosis and treatment of ES. New ES-associated technologies, as well as drugs, have emerged, but they did not meet the clinical standards. Therefore, it is essential to identify new genes and pathways associated with ES tumorigenesis and patient prognosis, in order to help unravel the relevant underlying molecular mechanisms, and to help discover novel diagnostic markers and therapeutic targets.

Bioinformatics is a relatively new interdisciplinary subject that emerged in the late 1980s alongside the launch of the Human Genome Project. It reflects the permeation and integration of biology, computer science, mathematics, physics (3) Through the acquisition, processing, storage, retrieval and analysis of experimental biological data, one achieves the purpose of interpreting the biological meaning and mechanisms underlying the activities in question from within the data. With the rapid development of microarray technology and bioinformatic analyses, microarrays are being used extensively for detecting gene expression levels, particularly in seeking differentially expressed genes (DEGs) (4). The application of microarrays has produced a large amount of data, which have been uploaded and stored in public databases. The use of microarray databases allows the identification of thousands of genes in ES, which can be utilized to screen for more molecular biomarkers. In the present study, two microarray datasets from the Gene Expression Omnibus (GEO) database (5) were downloaded and analyzed in order to investigate the characteristics of the ES genomic expression profiles, as well as to screen for DEGs. DEGs were identified by comparing tumor samples with non-cancerous samples. The potential functions of these DEGs were investigated via Gene Ontology (GO) terms (6) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (7) pathway enrichment analyses. In addition, the present study also constructed a protein-protein interaction (PPI) network to identify the hub genes in ES. A total of 278 genes were identified in the present study, and these findings may contribute to the diagnosis and therapy of ES in the future.

Materials and methods

Microarray data

The present study downloaded two gene expression datasets (GSE17674 (8) and GSE31215 (9) from the GEO database based on Affymetrix GPL570 platform data (Affymetrix Human Genome U133 Plus 2.0 Array). The GSE17674 dataset contained 62 tissue samples, including 44 ES samples and 18 normal skeletal muscle samples; the GSE31215 dataset contained 4 ES samples and 4 non-cancerous human pediatric mesenchymal stem cell samples.

Expression analysis of DEGs

The bioinformatics analysis involved downloading the original Celestia files and the platform-probe annotation information file. The raw data were first processed using the robust multi-array average (RMA) function (10) and later transformed into expression values utilizing the affy package (11) (version 3.5.0) in R language (12). The probe IDs were changed to the corresponding gene symbol according to the annotation information in the platform. The missing parts of the data were filled with the algorithm K-nearest neighbors (KNNimpute) (13). The limma package (14) of R was used for detecting the DEGs between the tumor samples and non-cancerous samples. A logFC (|log2-fold change|) value >1 and P<0.05 were considered to indicate a statistically significant difference.

GO and KEGG enrichment analysis

The Database for Annotation, Visualization and Integrated Discovery (DAVID; version 6.8) (15) is an online biological information database that provides a series of functional annotation tools to investigate the biological function of genes. In order to analyze the function of DEGs, the present study inputted the DEGs into the DAVID for the GO (6) and KEGG (16) enrichment analyses. GO annotates and classifies genes via biological pathways, molecular functions and cellular location, whereas KEGG identifies the relevant pathways in which these differential genes may be involved. Gene counts >5 and P<0.05 were selected as the cut-off criteria for the enrichment analysis.

PPI network and hub genes

The online Search Tool for the Retrieval of Interacting Genes (STRING) (17) is an online database that provides critical assessment and integration of PPIs. Interactions with a combined score >0.9 was considered statistically significant in the present study. The PPI network was visualized using Cytoscape (version 3.6.1) (18), and the Cytoscape plugin Molecular Complex Detection (MCODE) (19) was used for identifying the most significant module in the PPI network. Subsequently, MCODE scores >5, degree cutoff=2, node score cutoff=0.2, Max depth=100 and k-core=2 were used as filtering criteria. KEGG and GO analyses were performed for genes in this module. The top 5 genes sorted by the degree levels in the Cytoscape plugin, cytoHubba (20) were defined as hub genes, and these 5 genes were utilized in the co-expression gene network and mutation survival analysis, using Kaplan-Meier curves generated with the cBioPortal (http://www.cbioportal.org) online platform (21).

Results

Identification of DEGs

The overlap among the two datasets contained 361 genes (10,212 in the GSE17674 dataset, and 611 in the GSE31215 dataset), as presented in the Venn diagram (Fig. 1). The results demonstrated the upregulation of 272 genes in the two datasets, and the downregulation of 6 genes. Genes that were not simultaneously up- or downregulated in the two datasets were excluded, and a total of 278 genes were identified as DEGs (Table I).

Figure 1.

Venn diagram demonstrating the genes selected with a fold change >2 and a P-value <0.05 from the RNA expression profiling sets GSE17674 and GSE31215.

Table I.

Compared with non-cancerous samples, 278 DEGs were identified from the datasets analyzed, 272 of which were upregulated and 6 of which were downregulated in tumor samples.

Status	DEGs
Downregulated	SYNM, GPAT3, LINC01119, PPL, LOC101929897, HAS3
Upregulated	RIDA, EWSAT1, KHDRBS3, CDH8, SAMD9L, HOXD9, LPAR4, OGDHL, SV2B, SLC5A6, WBSCR17, MYO10, CAPN15, GFRA2, GADD45B, NR2F1, CCDC155, LRRC4C, COL21A1, ZNF57, CA11, DLX3, CD99P1, CDS1, HSD17B2, PARM1, RASSF2, TMEM176A, SOX8, RPP25, POLA1, CHGBPIP5K1C, PCDHB5, ZNF711, CHST1, PDGFD, TGIF2, HERC5, TMEM91, ADGRL1, PLSCR1, CLCA1, LINC01503, C5, ALDH16A1, ATP6V0B, ITGB2-AS1, CYTL1, TBCCD1, CLSTN2, PALM2, COL9A3, ASAP1, PRKAR2B, SULT1A1, DDIT4, PNP, FGFR3, IFITM3, ISL1, HLA-E, LAPTM5, ARTN, ADARB1, MYL9, TRIM69, MAPT, ADRA1D, SPATA13, MND1, TLE2, GMFG, FBXO15, C1orf226, TDRD3, FAM213A, LAMB3, FMOD, ZNF195, INTS6L, HIST1, H2BD, KCNE4, IFIT5, SP110, ASGR1, SLCO5A1, STK32B, GALNT14, IFITM1, IGSF21, ZNF703, ITGB2, OAS3, ST6GAL1, NYNRIN, ERFE, CLDN1, YPEL5, SYT4, OAS1, FREM1, CITED2, IL1RAP, MBD4, STMN4, NR0B1, LMO2, ERF, PTPN22, NCKAP1L, CSRNP2, BHLHE22, LINS1, PEG3-AS1, CAV2, RAP1GAP, ALG6, ADGRE5, INTU, L3MBTL3, HIST1H2BE, TOX3, EFNA1, EPHA3, SORD, TRIM35, MEIS1, NME7, ODF2L, KIAA1462, IGF1, ZNF468, ZNF432, TRHDE, NCKAP5, UNC5A, PROS1, TMEM71, PLXDC2, EMILIN1, LBH, HIST1H2BH, CADPS2, ZNF426 TCF4, LOC102724275, SSBP4, FAM107B, TEAD2, MAFB, SMARCC1, SH2B3, CPNE4, LINGO1, GOT1L1, CRIP1, RGS2, POU3F1, NPY, CRIP2, APELA, CARD16, ENPEP, ISG15, LECT1, GDF10, ADRB3, FAT3, MSX1, HIST1H2AC, TRIB2, SYT1, OLFML3, AKAP7, CD14, HS3ST4, PLSCR4, CLIP2, SEC11C, RAB11FIP1, S1PR3, LOXHD1, OLFM3, HOXD13, RUNX3, PGLYRP2, UTS2, KIT, DNAJC12, LGALS8, H2BFS, CPVL, CALCB, OLFM1, RCOR1, MFAP4, CSRP1, ZDHHC21, BGN, STEAP2, DLG2, IRS2, DUSP6, TFAP2B, TRPM4, GPR137B, IVNS1ABP, GIMAP2, AMER2, MXRA5, PCDH17, RARRES2, HOOK1, JAK1, TPBG, ITM2A, MYC, SMA4, APOE, SLAIN1, CCK, ABHD6, JAKMIP2, MMP9, GRP, RGL1, SCG2, EGR2, GALNT7, NPTXR, FCGRT, CDH11, NUDT11, NELL2, NAIP, MARCKSL1, ALDH7A1, LIPI, TOX2, TNC, SLC26A2, GLCE, APCDD1, FAT4, EDNRA, KMO, TMEFF2, TSPAN13, CCDC171, LY96, CYP26B1, ID4, DAPK1, PAX7, ATP1A1, RBM11, DKK2, EZH2, KDSR, CKS2, PRSS35, ADGRG2, PCDH8, ID2, PTPN13, STEAP1, HMCN1, FAM84B, BCL11B, PRKCB, NKX2-2

GO enrichment and KEGG pathway analysis

The present study imported the DEGs into DAVID in order to classify the functions with the GO and KEGG pathway analyses. The GO analysis revealed an enrichment of DEGs in 37 functional clusters (Table II). With respect to the biological process terms, the DEGs were primarily enriched in the positive regulation of, ‘positive regulation of transcription from RNA polymerase II promoter’, ‘positive regulation of cell proliferation’, ‘transcription from RNA polymerase II promoter’, ‘positive regulation of transcription, DNA-template’, ‘inflammatory response’, ‘cell adhesion’ and ‘cell proliferation’ (Fig. 2A). As for the cell component terms, the DEGs were primarily enriched in, ‘plasma membrane’, ‘extracellular exosome’, ‘integral component of plasma membrane’, ‘extracellular region’ and ‘extracellular space’ (Fig. 2B). In regard to the molecular function terms, the DEGs were primarily enriched in ‘chromatin binding’, ‘sequence-specific DNA binding’, ‘transcriptional activator activity, RNA polymerase II core promoter proximal region sequence-specific binding’ and ‘RNA polymerase II core promoter proximal region sequence-specific DNA binding’ (Fig. 2C). The KEGG pathway analysis demonstrated significantly enriched DEGs in certain cancer-associated pathways, including ‘phosphoinositide 3-kinase-protein kinase B (PI3K-Akt) signaling pathway’, ‘Rap1 signaling pathway’, ‘signaling pathways regulating the pluripotency of stem cells’, ‘cyclic guanosine monophosphate-protein kinase G signaling pathway’ and ‘transcriptional misregulation in cancer’ (Fig. 2D).

Table II.

GO and KEGG pathway enrichment analysis of DEGs in Ewing sarcoma samples.

Category	ID	Term	Count
BP	GO:0045071	Negative regulation of viral genome replication	6
BP	GO:0060337	Type I interferon signaling pathway	7
BP	GO:0051607	Defense response to virus	10
BP	GO:0045665	Negative regulation of neuron differentiation	6
BP	GO:0008284	Positive regulation of cell proliferation	17
BP	GO:0045944	Positive regulation of transcription from RNA polymerase II promoter	27
BP	GO:0001501	Skeletal system development	8
BP	GO:0001503	Ossification	6
BP	GO:0050830	Defense response to Gram-positive bacterium	6
BP	GO:0007156	Homophilic cell adhesion via plasma membrane adhesion molecules	8
BP	GO:0045893	Positive regulation of transcription, DNA-templated	16
BP	GO:0006954	Inflammatory response	13
BP	GO:0010628	Positive regulation of gene expression	10
BP	GO:0000165	MAPK cascade	10
BP	GO:0008283	Cell proliferation	12
BP	GO:0007267	Cell-cell signaling	9
BP	GO:0007155	Cell adhesion	13
BP	GO:0006366	Transcription from RNA polymerase II promoter	14
CC	GO:0005615	Extracellular space	36
CC	GO:0005576	Extracellular region	40
CC	GO:0005623	Cell	7
CC	GO:0031012	Extracellular matrix	11
CC	GO:0070062	Extracellular exosome	55
CC	GO:0005887	Integral component of plasma membrane	31
CC	GO:0043005	Neuron projection	9
CC	GO:0005667	Transcription factor complex	8
CC	GO:0005794	Golgi apparatus	20
CC	GO:0000139	Golgi membrane	15
CC	GO:0005886	Plasma membrane	71
MF	GO:0003682	Chromatin binding	16
MF	GO:0001077	Transcriptional activator activity, RNA polymerase II core promoter proximal region sequence-specific binding	12
MF	GO:0005179	Hormone activity	6
MF	GO:0030246	Carbohydrate binding	8
MF	GO:0005088	Ras guanyl-nucleotide exchange factor activity	6
MF	GO:0003714	Transcription co-repressor activity	8
MF	GO:0000978	RNA polymerase II core promoter proximal region sequence-specific DNA binding	11
MF	GO:0043565	Sequence-specific DNA binding	14
KEGG	hsa04550	Signaling pathways regulating pluripotency of stem cells	8
KEGG	hsa04015	Rap1 signaling pathway	9
KEGG	hsa05219	Bladder cancer	4
KEGG	hsa04666	Fc gamma R-mediated phagocytosis	5
KEGG	hsa05200	Pathways in cancer	12
KEGG	hsa04151	PI3K-Akt signaling pathway	11
KEGG	hsa04022	Cgmp-PKG signaling pathway	7
KEGG	hsa05202	Transcriptional misregulation in cancer	7

GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; BP, biological process; CC, cellular component; MF, molecular function.

Figure 2.

Functional enrichment analysis of DEGs in Ewing's sarcoma. GO analysis demonstrated significantly enriched DEGs in (A) biological process terms and (B) cell component terms and (C) molecular function terms. (D) Significantly enriched terms obtained from Kyoto Encyclopedia of Genes and Genomes pathway analysis. DEGs, differentially expressed genes; GO, Gene Ontology.

The present study subsequently constructed a PPI network utilizing STRING and Cytoscape in order to depict the interactions between the identified DEGs (Fig. 3A). There were 269 nodes and 292 edges in the network. A significant module was built with 8 nodes and 28 edges, which had the highest MCODE score (Fig. 3B). The functional analysis of genes involved in this module revealed that the genes were enriched primarily in multi-organism processes, immune system processes and 2′-5′-oligoadenylate synthetase activity.

Figure 3.

(A) PPI network of differentially expressed genes was constructed using Cytoscape. (B) The most significant module with the highest Molecular Complex Detection score was obtained from the PPI network. PPI, protein-protein interaction.

Hub gene analysis

The top 5 genes identified as hub genes were MYC, IGF1, OAS1, EZH2 and ISG15. The network of the hub genes and their co-expression genes was analyzed using the cBioPortal online platform (Fig. 4A). The mutation survival analysis was only available for EZH2 from the database, and the results revealed that the EZH2 mutation was significantly associated with worse disease-free survival time (Fig. 4B and C).

Figure 4.

(A) Hub genes and their co-expression genes were analyzed in cBioPortal. Nodes with a bold black outline represent the hub genes. Nodes with a thin black outline represent the co-expression genes. (B) Disease-free survival time and (C) overall survival time expression analyses of EZH2 were performed using the cBioPortal online platform. P<0.05 was considered to indicate a statistically significant difference.

Discussion

ES is primarily observed in bones and muscles, as well as in soft tissues. It is the second most common primary malignant bone tumor in children and adolescents, second only to osteosarcoma (1). The ES family of tumors is an aggressive childhood cancer that includes classic ES, Askin's tumor and peripheral primitive neuroectodermal tumor (1). Previous studies on ES have reported the characteristic t(11;22) chromosomal translocation, and the resulting EWS-FLI1 gene fusion was revealed to have an impact on diagnosis (22); despite this, the specific pathogenesis of ES is not yet understood. Microarray technology enables the investigation of genetic alterations in cancer and is a useful approach for identifying new biomarkers in a number of different diseases (23). Advances in microarray technology enabled the elucidation of the molecular mechanism underlying ES development, and also in the identification of novel diagnostic biomarkers.

In the present study, two gene profile datasets were selected from the GEO for the bioinformatics analysis. A total of 278 genes were identified through comparing tumor samples with non-cancerous samples. The KEGG enrichment analysis in the present study revealed that the identified DEGs were enriched primarily in cancer-associated pathways and the PI3K-Akt signaling pathway. Earlier studies have reported that the PI3K-Akt signaling pathway participates in the proliferation, invasion and migration of tumor cells, and these results support those of the present study (24,25). A PPI network was constructed for the identified DEGs, and the most significant module was subsequently extracted from the PPI network. A total of five key genes were defined based on the degree rank, and included two of the genes in this module.

As a proto-oncogene and a nuclear phosphoprotein coding gene, MYC plays a role in cell cycle progression, apoptosis and cellular transformation (26). Previous studies have reported the involvement of MYC amplification in a multitude of different types of human cancer (27,28). Diseases associated with MYC include Burkitt's lymphoma and multiple myeloma (29). MYC participation is evident in multiple pathways from the KEGG pathway analysis of DEGs. In an earlier report, Kawano et al (25) reported increased MYC expression by microRNA-20b, which promotes cell proliferation in ES cells. A previous study has demonstrated that MYC-driven cancer cells promote tumorigenesis through immune escape, suggesting that MYC-induced tumors may be particularly sensitive to immuno-oncological interventions (30).

Insulin growth factor 1 (IGF1) encodes a protein involved in mediating growth and development. As a ligand for the IGF1-receptor (IGF1R), IGF1 binds to the α subunit of IGF1R, resulting in the activation of intrinsic tyrosine kinase activity, thereby activating the PI3K/AKT and the RAS/mitogen-activated protein kinase pathways (31). Earlier evidence suggested that IGF1 pathway plays a significant role in the pathogenesis of ES, and that IGF1R serves as a potential therapeutic target both in vitro and in vivo (32). The IGF1 pathway [IGF1R, 2 ligands (IGF1 and IGF2), and 6 IGF-binding proteins (IGF-BP1-6)] serves as a biomarker for disease progression in patients with ES; in these patients, high circulating levels of IGF-1 and IGF-BP3 at baseline are associated with a favorable outcome following treatment with chemotherapy (33). In addition, the differential expression of IGF-I transcripts has been demonstrated in bladder cancer and colorectal carcinoma (34).

The 2′-5′-oligoadenylate synthetase 1 (OAS1) belongs to the OAS family, which inhibits cellular protein synthesis and resistance to viral infection (35). Lee et al (36) demonstrated that OAS1 negatively regulates the expression of chemokines and interferon-responsive genes in human macrophages. Literature search results did not reveal any association between OAS and the carcinogenesis or progression of tumors; thus suggesting that ES occurrence may be associated with OAS1, but further investigations are required.

EZH2 is a protein-coding gene that encodes a histone methyltransferase that serves as the catalytic subunit of PRC2. It has been reported that mutations in EZH2 are associated with a number of different types of cancer, such as gastric cancer and hematological malignancies, and high EZH2 expression is associated with metastatic disease in pediatric soft tissue sarcomas (37). Richter et al (38) demonstrated that EZH2 promotes tumorigenesis by shaping the oncogenicity and stem cell phenotype in ES pathology (38). Furthermore, EZH2 inhibitors have been used as targeted therapies in the treatment of a variety of malignancies (39). EZH2 is expressed in numerous types of immune cell and has distinct functions. In recent years, significant advances have been made by using EZH2 as a target of immunotherapy in various different types of cancer (39). An increasing number of EZH2-targeted drugs, such as tazemetostat, CPI-1205, GSK2816126 and SHR2554, are being tested using in vitro biochemical and cell analyses, in vivo preclinical models and patient clinical trials, and have demonstrated antitumor effects in various types of malignant tumor. Reports have demonstrated that treatment with EZH2 inhibitor sensitizes ES cells to effective cytolysis by G-specific chimeric antigen receptor gene-modified T cells (40).

Interferon-stimulated gene 15 (ISG15) encodes an ubiquitin-like protein that is conjugated to intracellular target proteins upon activation by IFN-α and IFN-β (41). Previous studies have reported that ISG15 plays a vital role in antiviral and antitumor mechanisms (42–44). The expression of ISG15 has implications in a wide range of human tumors, such as digestive system tumor, breast cancer and lung cancer (45–47). Upregulation of ISG15 has been reported in breast carcinoma cells and was proposed as a novel breast cancer marker with prognostic significance. Reports have demonstrated that strategies for manipulating ISG15 levels through the IFN-α/ISG15/p53 axis may prove to be effective in the treatment of cervical cancer (40). Hence, further investigation of the role of ISG15 would help in not only the identification of novel drug targets, but would also provide new opportunities for intervention in cancer progression.

In conclusion, the present study analyzed two ES-associated cohort profile datasets and identified candidate genes and pathways that may participate in ES tumorigenesis and progression. These results may not only contribute to an improved understanding of the underlying molecular mechanisms, but also provide a series of potential biomarkers and therapeutic targets. However, further experiments are required to verify the findings of the present study.