This article is mentioned in:

Introduction

Ewing sarcoma (ES) is a highly malignant tumor with a low survival rate and high rate of metastasis. In a British cohort of patients with ES, ES accounts for 14% of pediatric bone tumors (1). Comprehensive strategies, including localized surgery, radiotherapy and chemotherapy, have been developed for the treatment of patients with ES (2). However, 30–40% of patients develop recurrence or metastasis after comprehensive therapy (3). The 5-year survival rate of patients with ES was only 55% in a British cohort of patients (1). Therefore, further research on the pathogenesis of ES is required to improve the prognosis of affected patients.

According to a previous study, the Ewing sarcoma breakpoint region 1 (EWS)-E26 transformation-specific (ETS) fusion gene is a major factor in ES (4). EWS-ETS fusion genes may confer upon tumors the capacity for metastasis and invasion by altering RNA transcriptional regulation and epigenetic modification (5,6). Previous studies on the EWS-ETS fusion gene have revealed the unique role of this fusion gene in the development of ES; however, treatment targeting the EWS-ETS fusion gene has been difficult to achieve, and there is increasing evidence demonstrating that EWS-ETS may not be the sole driver in metastatic ES (7–9). In addition, the prognostic value of the EWS-ETS fusion gene in patients with ES remains unclear. A retrospective study identified a significant association between EWS-friend leukemia virus integration 1 (FLI1) transcript subtypes and patient outcomes (10). However, two prospective studies were unable to validate this observation (11,12). Consequently, recognition of new therapeutic or predictive biomarkers to contribute to the treatment of ES is urgently needed to improve the prognosis of these patients.

Weighted gene co-expression network analysis (WGCNA) is an advanced approach to studying the associations between genes and clinical traits. Compared with traditional microarray analysis methods, WGCNA analysis uses a soft threshold instead of the hard threshold of traditional differential gene screening, which facilitates the screening of valuable genes with small fold changes that may be important in the gene regulation cascade (13).

In addition, WGCNA is an efficient method of studying the associations between genes and clinical phenotypes (14). In the present study, WGCNA was used to classify genes with similar expression patterns into different modules. Associations between the co-expression modules and clinical features of patients were analyzed. The most relevant modules to overall vital survival (OVS) and metastasis status of patients with ES were selected for further analysis.

Materials and methods

Data collection and preprocessing

ES data were downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) as dataset GSE17679 (15). Probes were mapped to gene symbols using the hgu133plus2.db (version 3.2.3) package (http://bioconductor.org/packages/hgu133plus2.db/) in R. If multiple probes were mapped to gene symbols, the mean value was regarded as the expression value of the gene and the repeated probes were dismissed.

Construction of a weighted co-expression network

The WGCNA package in R is a widely used method of identifying co-expression networks (16). Pearson's correlation coefficient was calculated to construct the correlation matrix. A soft-thresholding function was used to transform the correlation matrix into a weighted adjacency matrix. To acquire a co-expression network with a balance between scale-independence and mean connectivity, the scale-independence and mean connectivity were calculated in different powers by the soft-connectivity algorithm. The adjacency matrix was transformed into a topological overlap matrix (TOM). 1-TOM was used as the distance measurement to cluster genes into co-expression modules with a deep split value of 2, minimum size cutoff of 20 and maximum module size of 5,000. A merge height of 0.3 was set as a criterion to cluster similar modules. To identify the association between a co-expression module and a clinical feature, P-values and correlation coefficients were calculated to produce trait-module heatmaps.

PPI network construction and functional annotation

The Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database (date of access, August 2018; http://string-db.org) was used to construct the protein-protein interaction (PPI) network, which contained experimental and in silico predicted PPI information. The data source was set to only the experiment source. Low-connection proteins were hidden. Hub genes are genes that serve a key role in the regulation of the network. In the present study, the betweenness centrality (BC) was used to measure the importance of a node in the PPI network. Cytohubb in Cytoscape software (version 3.7.1, http://www.cytoscape.org) was used to screen for the hub genes, and the top 50 genes with high BC values were selected for subsequent analysis. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment analysis was performed using the Metascape website (date of access, August 2018; http://metascape.org).

Survival analysis

Since ES is a rare cancer, after carefully searching the GEO database, ArrayExpress (https://www.ebi.ac.uk/arrayexpress/) and TCGA database (https://cancergenome.nih.gov/) with loose screening criteria (including clinical information), only one GEO dataset, GSE63157 (17), met the requirements (detailed clinical information) of the present study. Following annotation, the expression profiles of 85 samples and 17,186 probes were obtained. The Survival and Survminer packages in R were used to analyze and visualize the survival data and P<0.05 was considered to indicate a statistically significant difference. All the R packages used in the present study can be found at http://www.bioconductor.org or https://cran.r-project.org/web/packages.

Results

Construction of weighted gene co-expression network

The expression matrix and clinical data of GSE17679 (15) were downloaded from the GEO database. The scale-independence and mean connectivity were calculated at different thresholds (Fig. 1). Appropriate soft-thresholding has a relatively good balance between scale-independence and mean connectivity of the weighted co-expression network; scale-independence >0.85 and average connectivity <100 were used as the criteria for a suitable soft threshold. As presented in Fig. 1, a soft threshold power of 5 was chosen for the identification of co-expression modules. A cluster dendrogram based on the dissimilarity of the topological overlap matrix was generated (Fig. 2). Acceptable discriminability was present between each module in the similarity heatmap plot (Fig. 3).

Figure 1. Scale-independence and mean connectivity of the network in different soft-threshold powers. The left panel displays the correlation of soft threshold with scale-free fit index. The right panel displays the influence of soft-threshold power on mean connectivity.

Figure 2. Cluster dendrogram based on the dissimilarity of the topological overlap matrix. The different colors correspond to the co-expression modules in Ewing sarcoma.

Figure 3. Associations among modules. (A) Hierarchical clustering of module eigengenes that summarize the modules identified in the weighted gene co-expression network analysis. The branches of the dendrogram grouped together indicate a positive correlation. (B) The eigengene adjacency heatmap. Each row and column in the heatmap corresponds to a module (labeled with colors). Red, high adjacency; blue, low adjacency. Red squares along the diagonal axis represent meta-modules.

Identification of metastasis-related co-expression modules and OVS-related co-expression modules

The first principal component of each module was defined as the module eigengene (ME), which represented the expression levels of all genes in the module. Gene significance (GS) was used to evaluate the correlation between genes and clinical features by linear regression; the module significance (MS) was calculated by averaging the absolute GS values of all genes in the co-expression module. The most statistically significant module (P<0.0001) in each clinical feature (OVS and metastasis) were further analyzed.

As presented in Fig. 4, OVS was significantly correlated with the green module (r2=0.56; P=4×10−4). The royal blue module was correlated with event-free survival (EFS) (r2=0.5; P=0.002), whereas the green-yellow module was most significantly correlated with metastatic state (r2=0.67, P=9×10−6). Subsequent analyses focused on the metastasis-associated co-expression module (yellow-green) and the OVS-associated co-expression module (green).

Figure 4. Heatmap of the correlation between identified modules and clinical features. The x-axis corresponds to clinical features, and the y-axis represents the identified modules. The color scale (blue to red) indicates correlation. Data are presented as the correlation coefficient in the top row and the P-value in the bottom row. OVS, overall vital survival; EFS, event-free survival; STATE, metastatic state.

PPI network in co-expression modules

Following analysis of the STRING database, two PPI networks of the top 50 BC genes in the green and green-yellow modules were constructed. The PPI network of the green module is presented in Fig. 5; the top 10 genes of the green module by BC value were proteasome subunit α4 (PSMA4; BC=12,145), G-protein subunit γ5 (BC=11,152), clathrin light chain A (BC=8,382), ring-box 1 (BC=7,332), ribophorin II (BC=7,288), chaperonin-containing TCP1 subunit 6B (BC=6,931), serine/threonine kinase receptor-associated protein (STRAP; BC=6,736), 3-phophoinositide-dependent protein kinase 1 (BC=6,652), L1 cell adhesion molecule (L1CAM; BC=6,468) and Janus kinase 1 (BC=6,346). The PPI network of the green-yellow module is presented in Fig. 6. The top 10 genes in green-yellow module by BC value were vascular endothelial growth factor D (BC=300), serpin family A member 1 (BC=252), G-protein subunit α transducing 3 (BC=229), glutamate ionotropic receptor AMPA subunit 1 (BC=200), Met proto-oncogene (BC=198), leucine-rich repeat kinase 2 (BC=160), fibroblast growth factor receptor 2 (BC=154), calcium voltage-gated channel auxiliary subunit γ2 (CACNG2; BC=136), phospholipase C β1 (BC=112) and γ-aminobutyric acid type B receptor subunit 2 (GABBR2; BC=94).

Figure 5. Protein-protein interaction network of the top 50 high-BC genes in the overall vital survival-associated module. Nodes of the network represent proteins; the size of each node represents the BC of the node. Lines represent the interactions between the proteins according to the STRING database. BC, betweenness centrality.

Figure 6. Protein-protein interaction network of the top 50 high-BC genes in the metastasis-associated module. Nodes of the network represent proteins; the size of each node represents the BC of the node. Lines represent the interactions between the proteins according to the Search Tool for the Retrieval of Interacting Genes/Proteins database. BC, betweenness centrality.

Functional enrichment of metastasis- and OVS-associated modules

Functional enrichment analysis was used to determine the biological significance of the top 50 high-BC genes in the OVS- and metastasis-associated modules (Table I; Fig. 7). The results demonstrated that the OVS-associated module was mainly enriched in the cell cycle and immune response, including ‘adaptive immune system’ [log10(P)=−10.45], ‘pathways in cancer’ [log10(P)=−7.99], ‘cell cycle, mitotic’ [log10(P)=−7.81], ‘PID TGFBR PATHWAY’ [log10(P)=−6.92], ‘L1CAM interactions’ [log10(P)=−6.67] and ‘cell cycle’ [log10(P)=−6.57].

Figure 7. Heatmap of enriched terms across the top 50 high-betweenness centrality genes across the gene lists of the overall vital survival associated module.

Table I. Top 6 clusters with representative enriched terms (one per cluster).

Table I.

Top 6 clusters with representative enriched terms (one per cluster).

Category	Description	%a	Log10(P)b	Genes
R-HSA-1280218	Adaptive immune system	30	−10.45	CLTA, DYNC1I2, FKBP1A, PDPK1, PRKCB, PSMA4, FBXW4, SKP1, TRAF6, UBE2H, RBX1, HERC5, UBA6, TUBA1C, CDC26
hsa05200	Pathways in cancer	20	−7.99	CCND1, FZD2, GNG5, JAK1, PRKCB, RXRA, TGFBR1, TRAF3, TRAF6, RBX1
R-HSA-69278	Cell cycle, mitotic	22	−7.81	CCND1, CCND3, CCNH, DYNC1I2, PRKCB, PSMA4, SKP1, LPIN2, RBX1, TUBA1C, CDC26
M-286	PID TGFBR PATHWAYc	10	−6.92	FKBP1A, PDPK1, PPP1CA, TGFBR1, STRAP
R-HSA-373760	L1CAM interactions	12	−6.67	ANK1, CLTA, EPHB2, L1CAM, NUMB, TUBA1C
hsa04110	Cell cycle	12	−6.57	CCND1, CCND3, CCNH, SKP1, RBX1, CDC26

a Percentage of all genes of the top 50 betweenness centrality genes of the green module that are identified in each Gene Ontology term; only input genes with at least one ontology term annotation are included in the calculation.

b P-value in log base 10. L1CAM, L1 cell adhesion molecule.

c PID TGFBR PATHWAY TGF-β associated cellular signaling pathway.

The enrichment analysis of the top 50 high-BC genes of the metastasis-associated module was performed using Metascape (Table II; Fig. 8); the genes in this module were mainly associated with cellular metabolism and internal signal transduction, including ‘defective CSF2RB causes pulmonary surfactant metabolism dysfunction 5 (SMDP5)’ [log10(P)=−14.69], ‘PI3K-Akt signaling pathway’ [log10(P)=−8.59], ‘metabolism of xenobiotics by cytochrome P450’ [log10(P)=−7.91], ‘Ras signaling pathway’ [log10(P)=−6.26], ‘retinol metabolism’ [log10(P)=−4.93] and ‘neurotransmitter receptors and postsynaptic signal transmission’ [log10(P)=−4.68].

Figure 8. Heatmap of enriched terms across the gene lists of the metastasis-associated module, colored according to P-value.

Table II. Top 6 clusters with representative enriched terms (one per cluster).

Table II.

Top 6 clusters with representative enriched terms (one per cluster).

GO	Description	%a	Log10(P)b	Genes
R-HSA-5688849	Defective CSF2RB causes pulmonary surfactant metabolism dysfunction 5 (SMDP5)	12	−14.69	SFTPB, SFTPC, SFTPD, SFTA3, SFTPA1, SFTPA2
hsa04151	PI3K-Akt signaling pathway	20	−8.59	COL4A3, COL4A4, EFNA1, FGF5, FGFR2, VEGFD, ITGB6, MET, COL6A6, COL6A5
hsa00980	Metabolism of xenobiotics by cytochrome P450	12	−7.91	ADH1B, CYP1B1, CYP2B6, GSTA1, GSTA2, MGST1
hsa04014	Ras signaling pathway	14	−6.26	EFNA1, FGF5, FGFR2, VEGFD, MET, PLA2G1B, PLA2G10
hsa00830	Retinol metabolism	8	−4.93	ADH1B, CYP2B6, HSD17B6, SDR16C5
R-HSA-112314	Neurotransmitter receptors and postsynaptic signal transmission	10	−4.68	GRIA1, GABBR2, CACNG2, PLCB1, GNAT3

a Percentage of all genes of the top 50 betweenness centrality genes of the green-yellow module that are identified in each Gene Ontology term; only input genes with at least one ontology term annotation are included in the calculation.

b P-value in log base 10. CSF2RB, colony-stimulating factor 2 receptor β common subunit.

Prognostic significance of the identified genes

To determine the association between the identified hub genes and OVS of patients, the GEO GSE63157 dataset was used. All the genes in each module were sorted according to BC value, and R was used to calculate the association between the total survival time of patients with ES, and the top 10 hub genes ranked by BC were assessed. The top hub genes ranked by P-value of survival analysis are shown in Fig. 9. Significant associations were observed between the OVS of patients with ES and the four genes (PSMA4, L1CAM, STRAP and CTLA) (all P<0.05; Fig. 9). The association between the top 10 hub genes ranked by BC (2 genes were excluded as the probeID in the GeneChip® platform was unavailable) of the metastasis-related module and the OVS of the patients were also examined. It was determined that two genes were statistically significant prognostic factors for ES. High expression of CACNG2 and low expression of GABBR2 were favorable prognostic factors for patients with ES (Fig. 10).

Figure 9. Survival curve of the top 4 significant genes in the overall vital survival-associated module. (A) L1CAM, L1 cell adhesion molecule; (B) STRAP, serine/threonine kinase receptor-associated protein; (C) PSMA4, proteasome subunit α4; (D) CTLA, cytotoxic T-lymphocyte-associated protein 4.

Figure 10. Survival curves of (A) CACNG2 and (B) GABBR2. Genes in the metastasis-associated module were also associated with the overall vital survival of patients with Ewing sarcoma. CACNG2, calcium voltage-gated channel auxiliary subunit γ2; GABR2, γ-aminobutyric acid type B receptor subunit 2.

Discussion

ES is the second most common pediatric bone tumor (18). Although a comprehensive treatment plan for ES already exists, it remains an invasive tumor with a high recurrence rate and a low 5-year survival rate (2). Previous studies have revealed that the EWS-ETS fusion gene may be the driving gene for ES; however, targeted therapy for the EWS-ETS fusion gene is still difficult to achieve (2). Therefore, it is necessary to identify new target genes, which may lead to the improvement of prognosis in patients with ES.

In the present study, two modules associated with OVS and metastatic state were identified by weighted co-expression analysis. The GO and KEGG pathway enrichment analysis of the OVS-associated module revealed that the genes in this module were enriched in the immune response and cell cycle, including terms such as ‘adaptive immune system’ and ‘pathways in cancer’. The results of the WGCNA analysis demonstrated that PSMA4 served a hub role in the PPI network of the OVS-associated module. Previous studies have reported that mutations in PSMA4 may be associated with familial lung cancer (19,20), which suggests that PSMA4 may be a potential proto-oncogene. In addition, PSMA4 serves a role in chemotherapy resistance and immune response (21). PSMA4 forms a complex with proteasome activator subunits 3 and 4, which is required in the processing of major histocompatibility complex (MHC) class I-presented antigenic peptides (22–24). MHC class I serve a crucial role in the activation of cytotoxic T-cells, which is responsible for the tumor-associated immune response (25,26). Therefore, PSMA4 may be associated with tumor immune disorders, and further research is needed to elucidate the role of PSMA4 in the carcinogenesis of ES.

The protein encoded by L1CAM is an axonal glycoprotein that belongs to the immunoglobulin supergene family. L1CAM is a cell adhesion molecule associated with the prognosis of ovarian cancer, melanoma and endometrial cancer (27–29), and multivariate Cox survival analysis has suggested that L1CAM may be an independent prognostic variable in a number of different types of cancer, although there were no studies examining it potential as a prognostic variable in ES (30). These results suggest that L1CAM may serve a crucial role in tumor progression. In addition to that in cancerous tissues, abnormal L1CAM expression has also been identified in precancerous lesions. Geismann et al (31) revealed the role of L1CAM in the malignant transformation of pancreatic cancer, and an association with precancerous lesions was also identified in inflammatory bowel disease (32) and endometriosis (33). Taken together, previous studies have demonstrated the crucial role of L1CAM in the prognosis of carcinomas originating from epithelial tissue. The results of the present study revealed that L1CAM may also be a potential prognostic biomarker of ES.

In a previous study, STRAP was described to negatively regulate the transforming growth factor b (TGF-β) signaling pathway by binding to Smad7 (34). STRAP regulates the mitogen-activated protein kinase/extracellular signal-regulated kinase signaling pathway in a unmber ofdifferent types of tumor (35). Overexpression of STRAP has been identified in osteosarcoma (36). In addition, downstream of STRAP, the TGF-β signaling pathway is a crucial osteoblastic signaling pathway in the early stage of bone formation (37), and has been identified to be overexpressed in high-grade osteosarcoma tissue (38). Since STRAP serves osteogenic and oncogenic roles in bone tumors, further evidence is needed to determine the role of STRAP in the tumorigenesis of ES.

In the present study, KEGG pathway and GO functional enrichment showed that the top 50 BC genes in the metastasis-associated module were mostly enriched in terms associated with metabolism and signaling pathways. According to the literature, CACNG2 encodes a protein involved in the trafficking of glutamatergic AMPA receptors (39). Previous studies have revealed the key role of CACNG2 in mental and neuropathic disorders, including bipolar disorder, schizophrenia (40) and chronic pain (39). GABBR2 is a protein involved in neurological disorders, such as epilepsy (41) and Huntington's disease (42). However, the roles of CACNG2 and GABBR2 in ES are still unclear. The results of the present study suggested that these two genes may serve as potential biomarkers for metastasis or prognosis in patients with ES. Further research is needed to elucidate their role in this disease.

In conclusion, the present study identified several potential molecules involved in the metastasis and prognosis of ES. The WGCNA analysis identified certain co-expression modules in ES, which were associated with clinical features of the samples, to screen for clinical feature-associated genes. The hub genes identified in the present study require further research due to the lack of associated studies. The PPI network and functional enrichment of the hub genes may provide a view into the gene co-expression interaction framework of ES.

Acknowledgements

Not applicable.

Funding

No funding was received.

Availability of data and materials

The datasets generated and/or analyzed during the current study are available in the GEO repository,(accession nos. GSE63157 and GSE17679).

Authors' contributions

BW, JL and XL analyzed the data and wrote the manuscript. YO designed the experiment and revised the manuscript. The final version of the manuscript has been approved by all authors.

Ethics approval and consent to participate

Not applicable.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References