Identification of aberrantly methylated-differentially expressed genes and potential agents for Ewing sarcoma.

BACKGROUND: Human DNA methylation is a common epigenetic regulatory mechanism, and it plays a critical role in various diseases. However, the potential role of DNA methylation in Ewing sarcoma (ES) is not clear. This study aimed to explore the regulatory roles of DNA methylation in ES.
METHODS: The microarray data of gene expression and methylation were downloaded from the Gene Expression Omnibus (GEO) database, and analyzed via GEO2R. Venn analysis was then applied to identify aberrantly methylated-differentially expressed genes (DEGs). Subsequently, function and pathway enrichment analysis was conducted, a protein-protein interaction (PPI) network was constructed, and hub genes were determined. Besides, a connectivity map (CMap) analysis was performed to screen bioactive compounds for ES treatment.
RESULTS: A total of 135 hypomethylated high expression genes and 523 hypermethylated low expression genes were identified. The hypomethylated high expression genes were enriched in signal transduction and the apoptosis process. Meanwhile, hypermethylated low expression genes were related to DNA replication and transcription regulation. The PPI network analysis indicated C3, TF, and TCEB1 might serve as diagnostic and therapeutic targets of ES. Furthermore, CMap analysis revealed 6 chemicals as potential options for ES treatment.
CONCLUSIONS: The introduction of DNA methylation characteristics over DEGs is helpful to understand the pathogenesis of ES. The identified hub aberrantly methylated DEGs and chemicals might provide some novel insights on ES treatment. 2021 Annals of Translational Medicine. All rights reserved.

Introduction

Ewing sarcoma (ES) is a malignant tumor that occurs in bone or soft tissue and has a high incidence rate in children and adolescents (1). It is highly aggressive and easy to relapse after resection (2). Despite advanced multidisciplinary therapies, the prognosis of ES patients remains poor (3). A chromosomal translocation between the ESR1 gene on chromosome 22 and the FLI1 gene on chromosome 11 is present in approximately 85% of ES cases (4). Epigenetic alterations of tumor suppressor genes and protein-coding proto-oncogenes have greatly influenced the development of ES (5). For all these reasons, research on the epigenetics of ES is urgently needed for the improvement of clinical outcomes.

Epigenetics refers to the study of changes in gene function but does not involve changes in gene sequences (6). Epigenetic regulation contributes to fundamentally transcriptional regulation through several factors, including DNA methylation, non-coding RNA, and histones modification. EWS/FLI1 fusion and few recurrent mutations are not enough to explain the disease progression of ES, thus aberrant DNA methylation has been the subject of intensive research (7). Aberrant DNA methylation is common in most tumors and there emerging widespread hypomethylation leading to chromosomal instability while hypermethylation in particular CpG islands of some gene promoter related to gene silencing in ES cells (8), including DNA hypomethylation at enhancers regulated by the EWS-FLI1 fusion protein (9), the comprehensive profile of the aberrantly methylated-differentially expressed genes (DEGs) and their roles in the development of ES remain unclear.

Until now, little research has simultaneously analyzed gene expression and gene methylation data and their link to ES development, as the sequencing platform and related databases are updated, novel therapeutic targets and pathological mechanism might be found. In the present study, the gene expression profile (GSE45544) and gene methylation profile (GSE118872) were downloaded from the Gene Expression Omnibus (GEO) database, and a series of bioinformatics analyses were subsequently performed. We identified aberrantly methylated-DEGs and their correlative signaling pathways of ES. The hub genes that might be therapeutic targets were also determined. In the end, we performed a connectivity map (CMap) analysis and found several bioactive compounds that might be potential therapeutic agents for ES.

We present the following article in accordance with the STREGA reporting checklist (available at https://dx.doi.org/10.21037/atm-21-4972).

Methods

Microarray data

The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013). The gene expression profiling datasets (GSE45544) and gene methylation profiling datasets (GSE118872) were downloaded from the GEO (https://www.ncbi.nlm.nih.gov/geo/) of the National Center for Biotechnology Information (NCBI). A total of 21 ES and 22 normal samples were included in GSE45544 (platform: GPL6244 Affymetrix Human Gene 1.0 ST Array). A total of 21 ES and 3 normal samples were included in GSE118872 (platform: GPL8490 Illumina Human Methylation 27 Bead Chip).

Data processing

We used GEO2R, an R-based tool provided in GEO, to analyze original submitted data from microarrays. The DEGs and differentially methylated genes (DMGs) were identified with the criteria of P<0.05 and |t|>2. Subsequently, we obtained hypomethylated-upregulated genes by overlapping hypomethylated and upregulated genes, and by taking the intersection of hypermethylated and downregulated genes, we obtained hypermethylated-downregulated genes.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis

Functional and pathway enrichment analysis was conducted with hypergeometric distribution through the Database for Annotation, Visualization, and Integrated Discovery (DAVID, https://david.ncifcrf.gov/) (10). GO elucidated the roles of the selected genes in 3 aspects: molecular function (MF), biological process (BP), and cellular component (CC) (11). KEGG pathway enrichment analysis of selected genes was performed to identify important pathways. P value <0.05 was considered significant.

Protein-protein interaction (PPI) network construction and module analysis

To demonstrate the interaction between identified genes, a PPI network was built in the Search Tool for the Retrieval of Interacting Genes (STRING) database with the criteria of an interaction score ≥0.4. For further analysis, the network was visualized in Cytoscape, an open-source software platform for visualizing complex networks (https://cytoscape.org/). Molecular Complex Detection (MCODE), an app for cluster analysis in Cytoscape, was used to find important functional modules of selected genes. Degree cut off =2, node score cut off =0.2, k-score =2, and max depth =100 were regarded as significant. CytoHubba, an app in Cytoscape, was used to identify hub genes with the Matthews Correlation Coefficient (MCC) algorithm, and node =5 was set as the criterion.

Identification of small-molecule compounds by CMap analysis

The CMap, a gene-expression-based drug development database, makes it efficient for researchers to find potential agents that might reverse transcriptional profiles of relevant diseases (12). For CMap analysis, the hub genes were converted into probes and subsequently submitted to the “query” section. A P value <0.01 and enrichment score <0 were set as screening criteria.

Statistical analysis

The DEGs and DMGs were analyzed and obtained using the limma package of the R software containing empirical Bayes test statistics, with the threshold of |t|>2 and P≤0.05 as statistically significant.

Results

Identification of DMGs and DEGs in ES

To find DEGs or DMGs, we used GEO2R to analyze the corresponding profiling datasets. We then identified 1,876 upregulated genes and 6,805 downregulated genes from the gene expression microarray (GSE45544). We also obtained 843 hypomethylated genes as well as 1,503 hypermethylated genes from gene methylation profiling microarrays (GSE118872). Subsequently, by overlapping 1,876 upregulated DEGs and the 843 hypomethylated DMGs, a total of 135 hypomethylated-upregulated genes were obtained. Similarly, by taking the intersection of 6,805 downregulated DEGs and 1,503 hypermethylated DMGs, we obtained 523 hypermethylated-downregulated genes ().

Figure 1

Venn diagrams showing the aberrantly methylated-differently expressed genes from gene expression microarrays (GSE45544) and gene methylation microarrays (GSE118872) (A) Hypomethylated and upregulated genes. (B) Hypermethylated and downregulated genes.

GO and KEGG enrichment analysis

To visualize the attributes of the genes found above, we performed GO and KEGG enrichment analysis through DAVID. For hypomethylated-upregulated genes, we found that the MFs of those genes were enriched in calcium ion binding, oxidoreductase activity, and serine-type endopeptidase inhibitor activity. As for BPs, these genes indicated enrichment in intracellular signal transduction and positive regulation of the apoptotic process. For CCs, the analysis suggested enrichment mainly focused on extracellular regions. It indicated hypomethylated-upregulated genes might function in the tumor microenvironment of ES. Hypomethylated-upregulated genes might significantly act in some pathways, including metabolic pathways, phagosomes, and peroxisome proliferators-activated receptors (PPAR) signaling pathways ().

Figure 2

GO and KEGG enrichment analysis of hypomethylated and upregulated genes in ES. GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; ES, Ewing sarcoma.

For hypermethylated-downregulated genes, MF enrichment indicated protein binding, DNA binding, and structural constituents of the ribosome. Enriched BPs included cell division, apoptotic process, and transcription from the RNA polymerase II promoter. Additionally, enriched cell components included the nucleus and nucleoplasm, and cytosol and cytoplasm. For hypermethylated-downregulated genes, the KEGG analysis suggested enrichment of those genes in pyrimidine metabolism, cell cycle, purine metabolism, and ubiquitin-mediated proteolysis ().

Figure 3

GO and KEGG enrichment analysis of hypermethylated and downregulated genes in ES. GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; ES, Ewing sarcoma.

PPI network construction, module analysis, and hub gene selection

To identify the relationships between hypomethylated-upregulated genes, we conducted PPI network analysis according to the STRING database, and then performed an MCODE analysis in Cytoscape. This yielded 3 important modules which were displayed in . We analyzed hypermethylated-downregulated genes in the same way, and the 3 important modules were illustrated in . The selected hub genes were shown in .

Figure 4

Three modules of hypomethylated-upregulated genes. (A) Top 1 module. (B) Top 2 module. (C) Top 3 module.

Figure 5

Three modules of hypermethylated-downregulated genes. (A) Top 1 module. (B) Top 2 module. (C) Top 3 module.

Figure 6

Hub genes of aberrantly methylated-differently expressed genes in ES. (A) Five hub genes of hypomethylated-upregulated genes. (B) Five hub genes of hypermethylated-downregulated genes. ES, Ewing sarcoma.

CMap analysis

Using the CMap database, we found a total of 6 bioactive compounds that might exert therapeutic effects on ES, including iloprost, prednisolone, diphenhydramine, butacaine, benzbromarone, and fendiline, as shown in . Among them, iloprost, benzbromarone, and fendiline were reported to have an obvious anti-tumor effect. The others, including prednisolone and diphenhydramine, were usually combined with other drugs for tumor clinical treatment.

Table 1

Potential therapeutic agents for ES

Rank	CMap name	Mean	N	Enrichment	P value
1	Iloprost	−0.823	3	−0.964	0.00012
2	Prednisolone	−0.371	5	−0.773	0.00106
3	Diphenhydramine	−0.444	5	−0.747	0.00192
4	Butacaine	−0.35	4	−0.776	0.00515
5	Benzbromarone	−0.642	3	−0.834	0.00899
6	Fendiline	−0.605	3	−0.831	0.00973

ES, Ewing sarcoma; CMap, connectivity map.

Discussion

Despite progress in the treatment of ES, the survival rate remains poor (13). Most patients diagnosed with ES have a chromosomal translocation with rare recurrent mutations (14), which were not enough to explain cancer evolution. Some researchers have reported that chromosomal translocation can facilitate widespread DNA methylation alterations and further influence disease progression, such as PTRF silencing caused by hypermethylation of PTRF contributes to the tumorigenic behavior of ES (15), Therefore, it is essential to elucidate the underlying DNA methylation changes in ES.

For hypomethylated-upregulated genes, GO enrichment analysis demonstrated that the BPs were associated with intracellular signal transduction, proteolysis, and positive regulation of apoptotic processes. This is reasonable because previous studies have demonstrated that both altered signal transduction cascades and evasion of apoptosis are hallmarks of cancers, including ES (16,17). Poly(ADP-ribose) polymerase (PARP) proteolysis has been detected in the process of apoptosis in ES (18). The KEGG pathway enrichment analysis showed several significantly enriched pathways including arginine and proline metabolism, the peroxisome proliferator-activated receptor (PPAR) signaling pathway, and pertussis. The PPAR signaling pathway has been reported to be involved in a variety of cancers, but its specific function is unclear (19). Some researchers believe that PPAR pathway activation could promote the initiation or development of tumors, while others hold the opposite view (20). Arginine and proline are 2 important nonessential amino acids, and their synthesis could promote tumor growth (21). Besides, pertussis toxin can influence the stimulation and inhibition of adenylate cyclase, which plays an important role in ES activities (22).

The PPI network illustrated the functional connections between gene expression and gene methylation in ES. Through Cytoscape, the top 5 hub genes were presented: C3, TF, PNPLA2, VWA1, and FUCA2. The C3 gene encodes complement C3, which is a central protein associated with most complement activation pathways, and the C3 gene is a direct target of PPAR receptor-α, which could function in certain functional pathways of ES metabolism (23,24). The TF gene encodes transferrin, an iron-binding transport protein that is usually overexpressed on the surface of tumor cells, including ES cells (25,26), some researchers found that transferrin and its receptor were also closely associated with tumor cell growth and proliferation (27). The PNPLA2 gene encodes a protein named adipose triglyceride lipase, which is expressed in human tumors and related to tumor aggressiveness (28). The VWA1 gene encodes von Willebrand factor A-domain-related protein (WARP), the function of which is still unclear (29). However, there is translational upregulation of VWA1 in hepatocellular carcinoma (HCC), and cell proliferation can be inhibited after knockdown of VWA1 (30). The expression of the FUCA2 can be upregulated by interferon-gamma (IFN-γ), which is closely associated with the initiation and development of tumors (31); however, little is known about its function in ES. The driving factor EWS/FLI1 fusion reprograms DNA methylation. Regulating the methylation level of hub genes in disease progression may play a role in delaying disease progression or curing, which could be achieved by regulating methylases and demethylases, and there are already epigenetic modifiers, including DNA methylation inhibitors and histone deacetylase inhibitors, approved for cancer treatment.

For hypermethylated-downregulated genes, GO enrichment analysis demonstrated the BPs to be cell division, apoptotic process, transcription from the RNA polymerase II promoter, and DNA replication. MF enrichment included protein binding, DNA binding, poly(A) RNA binding, and structural constituent of ribosome. The KEGG analysis showed that the enriched pathways were pyrimidine metabolism, cell cycle, purine metabolism, and ubiquitin-mediated proteolysis. These results implied the important role of DNA replication and transcription dysregulation in the development of ES, which might be regulated by DNA methylation (32).

After the construction of the PPI network, the top 5 hub genes were identified through Cytoscape: CDC20, TCEB1, CUL2, ASB6, and KLHL13. The CDC20 gene is important in the regulation of the cell cycle, and its expression is upregulated in a variety of malignant tumors (33,34). Hypoxia-inducible factor-1α (HIF-1α) was found to directly bind to the ES-FLI-1 promoter and could be closely related to the invasiveness of ES, while TCEB1 could increase HIF stabilization (35,36). Thus, a TCEB1-targeted drug might be helpful for ES treatment. The CUL2 gene is believed to be a candidate tumor suppressor gene, and the protein it encodes is a member of the cullin family, which participates in the ubiquitination and degradation of certain cellular proteins (37,38). The ASB6 gene encodes a kind of ankyrin repeat protein that contains a C-terminal suppressor of cytokines signaling (SOCS) box motif, and the SOCS functions in ubiquitination and proteasomal degradation of certain proteins (39). The Kelch domain-containing protein encoded by KLHL13 was found to be associated with cytoskeletal functions (40). The role of ASB6 and KLHL13 in the development of ES requires further exploration. Regulating the methylation level of hub genes may be an effective treatment.

The CMap database is commonly used for drug development. In this study, based on CMap analysis, we determined 6 small-molecule compounds that might be potential treatment options for ES. Iloprost was shown to be of the highest significance among our results. As a prostaglandin analog, iloprost has antithrombotic, vasodilatation, and cytoprotection effects (41). In some malignant tumors, increased platelets are frequently observed, and researchers have speculated that those platelets played a critical role in tumor invasion and metastasis (42). Some preclinical studies have demonstrated that iloprost could reduce the adhesion of tumor cells to platelets, inhibit tumor angiogenesis, and hematogenous spread (43). Thus, iloprost might be a potential tumor inhibitor for clinical use. Benzbromarone was also reported to have an anti-tumor effect. It might influence the motility of vascular endothelial cells, and subsequently inhibit tumor angiogenesis (44). Recent studies have shown that the FDA approved drug, fendiline, could influence multiple tumor phenotypes, including proliferation, migration, invasion, and independent growth (45).

In this study, several limitations should be acknowledged. First, we only analyzed 2 datasets, and a larger sample size is required to verify our results. Also, in vitro and in vivo experiments are required to further examine the effects of the filtered hub genes and associated pathways on ES. Analysis of clinical parameters and prognosis is also lacking due to the limitations of data. Therefore, further experiments are indispensable to verify the specific function of methylation and hub genes in ES.

Conclusions

In summary, this study identified aberrantly methylated-DEGs and their roles in ES using bioinformatics tools. We identified 10 hub genes, among them, C3, TF, and TCEB1 might be potential therapeutic targets for ES. The 6 identified bioactive compounds might provide some novel insights for ES treatment.

The article’s supplementary files as