Repression of enhancer RNA PHLDA1 promotes tumorigenesis and progression of Ewing sarcoma via decreasing infiltrating T‐lymphocytes: A bioinformatic analysis

Background: The molecular mechanisms of EWS-FLI-mediating target genes and downstream pathways may provide a new way in the targeted therapy of Ewing sarcoma. Meanwhile, enhancers transcript non-coding RNAs, known as enhancer RNAs (eRNAs), which may serve as potential diagnosis markers and therapeutic targets in Ewing sarcoma.

Materials and methods: Differentially expressed genes (DEGs) were identified between 85 Ewing sarcoma samples downloaded from the Treehouse database and 3 normal bone samples downloaded from the Sequence Read Archive database. Included in DEGs, differentially expressed eRNAs (DEeRNAs) and target genes corresponding to DEeRNAs (DETGs), as well as the differentially expressed TFs, were annotated. Then, cell type identification by estimating relative subsets of known RNA transcripts (CIBERSORT) was used to infer portions of infiltrating immune cells in Ewing sarcoma and normal bone samples. To evaluate the prognostic value of DEeRNAs and immune function, cross validation, independent prognosis analysis, and Kaplan–Meier survival analysis were implemented using sarcoma samples from the Cancer Genome Atlas database. Next, hallmarks of cancer by gene set variation analysis (GSVA) and immune gene sets by single-sample gene set enrichment analysis (ssGSEA) were identified to be significantly associated with Ewing sarcoma. After screening by co-expression analysis, most significant DEeRNAs, DETGs and DETFs, immune cells, immune gene sets, and hallmarks of cancer were merged to construct a co-expression regulatory network to eventually identify the key DEeRNAs in tumorigenesis of Ewing sarcoma. Moreover, Connectivity Map Analysis was utilized to identify small molecules targeting Ewing sarcoma. External validation based on multidimensional online databases and scRNA-seq analysis were used to verify our key findings.

Results: A six-different-dimension regulatory network was constructed based on 17 DEeRNAs, 29 DETFs, 9 DETGs, 5 immune cells, 24 immune gene sets, and 8 hallmarks of cancer. Four key DEeRNAs (CCR1, CD3D, PHLDA1, and RASD1) showed significant co-expression relationships in the network. Connectivity Map Analysis screened two candidate compounds, MS-275 and pyrvinium, that might target Ewing sarcoma. PHLDA1 (key DEeRNA) was extensively expressed in cancer stem cells of Ewing sarcoma, which might play a critical role in the tumorigenesis of Ewing sarcoma.

Conclusion: PHLDA1 is a key regulator in the tumorigenesis and progression of Ewing sarcoma. PHLDA1 is directly repressed by EWS/FLI1 protein and low expression of FOSL2, resulting in the deregulation of FOX proteins and CC chemokine receptors. The decrease of infiltrating T‐lymphocytes and TNFA signaling may promote tumorigenesis and progression of Ewing sarcoma.

Introduction

Ewing sarcoma is an aggressive tumor, which typically affects bones and soft tissue in children, adolescents, and young adults (Grunewald et al., 2018). With significant racial disparity, the overall incidence for Ewing sarcoma is ∼1.5 cases per million in Europe, and the peak age is 15 years old (Jawad et al., 2009). Ewing sarcoma is also the second common bone cancer (Ferguson and Turner, 2018) and it usually develops in the diaphysis of bones and metastasizes to lungs and bones. Besides, the primary tumor site varies with age, older patients (20–24 years of age) with a higher proportion of pelvic and axial primary tumors, metastatic diseases, and worse outcomes (Worch et al., 2018). Treatment of patients with Ewing sarcoma includes surgery, chemotherapy, and/or radiation therapy and so on (Grunewald et al., 2018). Currently, the 5-year overall survival is 65–75 percent for patients with localized disease. However, patients with metastatic disease have a strikingly lower 5-year overall survival of less than 30 percent, and those with isolated pulmonary metastasis have approximately 50 percent 5-year overall survival (Gaspar et al., 2015).

Ewing sarcoma is driven by a recurrent t (11; 22) (q24; q12) chromosomal translocation (Aurias et al., 1984) that results in the FET–ETS fusions. The most common fusion is EWS–FLI1 (Delattre et al., 1992), which encodes an oncogenic transcription factor (May et al., 1993), regulating different target genes (Cidre-Aranaz and Alonso, 2015) governing the initiation and progression of Ewing sarcoma (Riggi et al., 2014). Therefore, the molecular mechanisms of EWS-FLI–mediating target genes and downstream pathways may provide a new way in the targeted therapy of Ewing sarcoma.

Enhancers are discrete DNA regulatory elements with specific sequence motifs; they interact with target gene promoters and then enhance the transcription of target genes (Blackwood and Kadonaga, 1998). Meanwhile, enhancers also transcript non-coding RNAs, known as enhancer RNAs (eRNAs) (Kim et al., 2010). Recent progress have found that the transcription of active enhancer mostly initiates cell transcription and 40,000–65,000 eRNAs express in human cells (Andersson et al., 2014; Arner et al., 2015; Lee et al., 2020). Besides the direct mechanism, eRNA can also be elicited by tissue-specific transcription factors (TFs). Importantly, activation of tumorigenesis often converges to the destabilization of eRNAs (Zhang et al., 2019; Lee et al., 2020). However, the functional mechanisms of eRNAs in Ewing sarcoma are still unknown. We proposed that eRNAs may serve as potential diagnosis markers and therapeutic targets in Ewing sarcoma.

In this study, based on an integrated bioinformatics analysis, differential expressed eRNAs, transcription factors, hallmark signaling pathways, and immune cells/functions were identified between Ewing sarcoma samples and normal bone samples. Moreover, we also constructed a complete regulatory network to reveal the potential upstream and downstream mechanisms of further exploring the prognostic biomarkers and treatment targets, which provided a basis and reference for the prognostic risk of Ewing sarcoma tumorigenesis.

Materials and methods

Data collection

RNA-sequencing (RNA-seq) data of 85 Ewing sarcoma samples were downloaded from Treehouse database (https://treehousegenomics.soe.ucsc.edu/public-data/#datasets), an RNA database of children’s tumors. RNA-seq data of 3 normal bone samples were downloaded from SRA database (https://www.ncbi.nlm.nih.gov/sra/). For validation, we also obtained gene expression profiles of 256 sarcoma samples from TCGA database (https://tcga-data.nci.nih.gov). Also single-cell RNA sequencing (scRNA-seq) data of GSE146221 were downloaded from Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE146221) to verify our results. Batch effects of these RNA-seq data were reduced using normalization and batch-effect correction methods.

Next, the eRNA expression profiles of Ewing sarcoma and the target gene list corresponding to eRNAs were downloaded from eRNA in cancer (eRic) database (https://hanlab.uth.edu/eRic/) (Zhang et al., 2019), which benefits researchers to obtain eRNA expression profile, as well as the target genes and drug response of eRNA across TCGA samples. Besides, based on the gene location in hg38 genome, ChIP seeker package was utilized to annotate the official gene symbol of each eRNA (Yu et al., 2015).

Moreover, expression profiles of 318 transcription factors (TFs) were downloaded from Cistrome database (http://cistrome.org/) (Zheng et al., 2019). 50 hallmarks of cancer and 29 immune gene sets were obtained from the Molecular Signatures Database (MSigDB) (http://software.broadinstitute.org/gsea/msigdb) (Liberzon et al., 2015). This study was approved by the Ethics Committee of Tongji University School of Medicine.

Differential expression analysis

First off, differential expression analysis was conducted to identify differentially expressed genes (DEGs) between Ewing sarcoma samples and normal bone samples by utilizing the Linear Models for Microarray Data (limma) package (Smyth, 2004). Specifically, DEGs were distinguished according to | Log2 fold-change (FC) | > 1 and false discovery rate (FDR) < 0.05. Also, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were conducted to reveal the biological function of DEGs. Likewise, differentially expressed eRNAs (DEeRNAs) and target genes (DETGs), as well as the differentially expressed TFs between Ewing sarcoma samples and normal bone samples were also identified based on the same criteria.

Construction of prognostic prediction model in sarcoma

To evaluate the prognostic value of DEeRNAs, we used sarcoma samples from TCGA database to conduct cross-validation and independent prognosis analysis. The sarcoma samples were randomly assigned into training set (156 samples) and testing set (100 samples). Training set was utilized to construct the prognostic prediction model, while testing set was utilized to evaluate the prediction model.

Before constructing the prognostic prediction model, lasso regression was applied to avoid overfitting. Then, univariate Cox regression analysis was performed to select DEeRNAs in relation to prognosis. The DEeRNAs independently associated with prognosis, screened by multivariate Cox regression once again, were eventually integrated into the prognostic prediction model. Thus, the risk score of each sarcoma sample was calculated according to the following formula:

$\mathbf{R}\mathbf{i}\mathbf{s}\mathbf{k}\mathbf{S}\mathbf{c}\mathbf{o}\mathbf{r}{\mathbf{e}}_{\mathbf{i}}={\mathbf{\beta }}_{\mathbf{1}}\times \mathbf{g}\mathbf{e}\mathbf{n}{\mathbf{e}}_{\mathbf{1}}+{\mathbf{\beta }}_{\mathbf{2}}\times \mathbf{g}\mathbf{e}\mathbf{n}{\mathbf{e}}_{\mathbf{2}}+{\mathbf{\beta }}_{\mathbf{3}}\times \mathbf{g}\mathbf{e}\mathbf{n}{\mathbf{e}}_{\mathbf{3}}+\mathbf{...}+{\mathbf{\beta }}_{\mathbf{j}}\times \mathbf{g}\mathbf{e}\mathbf{n}{\mathbf{e}}_{\mathbf{j}}.$

Among the formula, “i” was the order number of sarcoma samples, while “j” was the quantity of DEeRNAs in this model. “β” was the regression coefficient of corresponding DEeRNAs. Each sarcoma sample was given a risk score, and based on the mean of risk scores, 256 sarcoma samples was classified as high-risk group and low-risk group. The same was true in training set and testing set. Through ranking risk score of each sarcoma sample, scatter dot plot and heatmap were delineated to display the survival time and the expression of independent prognostic factors in high-risk group and low-risk group. Additionally, receiver operator characteristic (ROC) curve was conducted to evaluate the efficiency of the prediction model. In high-risk group and low-risk group, function enrichment analysis was also conducted using GO and KEGG analysis, as well as the hallmark of cancer gene sets.

Validation of immune clustering among DEeRNAs

To infer portions of infiltrating immune cells in Ewing sarcoma and normal bone samples, expression of DEeRNAs between Ewing sarcoma samples and normal bone samples as well as the correlation in 8 immune cell types were identified using cell-type identification by estimating the relative subsets of RNA transcripts (CIBERSORT). CIBERSORT was performed with 1,000 permutations, where a threshold <0.05 was recommended. Also, correlation analysis was applied to infer the associations between different types of immune cells.

To validate the prognostic value of immune proportions, CIBERSORT was also implemented in 221 sarcoma samples in TCGA database after removing the missing data. Defined by the primary gene signature file LM22 of CIBERSORT, 22 types of immune cells were identified. Through single-sample gene set enrichment analysis (ssGSEA), the immune infiltration degrees of 29 types of immune cells were detected using 29 immune gene sets from MSigDB. Eventually, Kaplan–Meier survival analysis was utilized to display the correlation between survival and immune proportions in sarcoma samples.

Identification of differentially expressed hallmarks of cancer and immune gene sets

Gene Set Variation Analysis (GSVA) (Hanzelmann et al., 2013) was conducted to detect the expression of hallmarks of cancer in Ewing sarcoma and normal bone samples. Then, differential expression patterns of 50 hallmarks of cancer between Ewing sarcoma and normal bone samples were determined by differential expression analysis using limma R package (Smyth, 2004). The immune infiltration degrees of 29 types of immune cells in Ewing sarcoma and normal bone samples were detected using ssGSEA based on their specific surface markers (Barbie et al., 2009).

Construction of DEeRNA regulatory network for Ewing sarcoma oncogenesis

First of all, DEeRNAs and DETGs annotated by eRic database, as well as DETFs were retrieved from the above screening. Then, differentially expressed hallmarks of cancer were quantified as continuous variables by GSVA, and immune cells and gene sets were separately obtained from CIBERSORT and ssGSEA. Subsequently, co-expression analysis was conducted among the aforementioned factors, which were illustrated in different colors. Purple indicated the immune cell types by CIBERSORT, blue indicated the hallmarks of cancer by GSVA, indigo blue indicated the immune gene sets by ssGSEA, yellow indicated potential upstream DETFs of DEeRNAs, and pink indicated potential DETGs of DEeRNAs. The interaction pairs between DEeRNAs and DETFs, DETGs, immune cell types by CIBERSORT, hallmarks of cancer by GSVA, and immune gene sets by ssGSEA were utilized to construct the regulatory network for Ewing sarcoma oncogenesis. In the network, we set thresholds as cor. (correlation coefficient) > 0.85 and p < 0.05 between DEeRNAs and DETGs; cor. > 0.70 and p < 0.05 between DEeRNAs and DETFs; cor. > 0.50 and p < 0.05 between DEeRNAs and infiltrating immune cells; cor. > 0.50 and p < 0.05 between DEeRNAs and immune gene sets; cor. > 0.60 and p < 0.05 between DEeRNAs and hallmarks of cancer. Besides, the Pearson co-expression analysis was also utilized to estimate the correlation between the six components in the regulatory network.

Identification of candidate small-molecule drugs

In the Connectivity Map (CMap) database (https://portals.broadinstitute.org/cmap/) (Lamb et al., 2006), DEG maps were utilized to predict the associations between small molecule drugs and various diseases. The positive score was the same as the reference gene expression profile, whereas the negative score may be the opposite. Here, CMap was used to determine small molecule drugs that may target Ewing sarcoma based on the expression profiles. Specifically, the database was utilized to screen enrichment fractions < -0.85 and p < 0.05, and small molecule drugs with negative scores were considered as candidate therapeutic molecules.

ATAC-seq validation of key DEeRNAs

Assay for Transposase Accessible Chromatin with high-throughput sequencing (ATAC-seq) data of key DEeRNAs were obtained from chromatin accessibility landscape of primary human cancers (https://gdc.cancer.gov/about-data/publications/ATACseq-AWG), which were then used to identify the chromatin accessibility in the location of these DEeRNAs (Corces et al., 2018).

External validation

To further demonstrate the reliability of our findings, multidimensional external validation was conducted based on multiple online databases. First off, the Human Protein Atlas (Uhlen et al., 2015), cBioportal (Cerami et al., 2012), and Oncomine (Rhodes et al., 2004) databases were used to compare the expression of DEeRNAs between normal and pathological tissues. Besides, Encyclopedia of Cancer Cell Lines (CCLE) (Ghandi et al., 2019) was used to show the expression of DEeRNAs across various different cancer cell lines. Also, Gene Expression Profiling Interactive Analysis (GEPIA) was a web-based tool to conduct survival analysis of single gene (Tang et al., 2017; Li et al., 2021).

Moreover, CR Cistrome database (http://cistrome.org/db/#/) (Wang Q et al., 2014)was applied to elucidate the interaction between DETFs and DEeRNAs in the chromatin level, based on chromatin-immunoprecipitation followed by sequencing (ChIP-seq) for Histone 3 Lysine 27 acetylation (H3K27ac). Furthermore, eRic database (Zhang et al., 2019) (https://hanlab.uth.edu/eRic/) was utilized to validate the expression, clinical relevance, target genes, and drug response of DEeRNAs.

Single-cell RNA sequencing transcriptome analysis

The single-cell RNA sequencing (scRNA-seq) data of GSE146221 were downloaded from Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE146221), which included Ewing sarcoma cell lines CHLA9, CHLA10, and TC71 (Miller et al., 2020). All data were integrated by “IntegrateData” function and analyzed by the R toolkit Seurat (http://satijalab.org/seurat/). Those cells were extracted for the following analysis which had more than 100,000 transcripts expressing. After the top 2,000 variable genes were filtered via “vst” method, “FindConservedMarkers,” and “FindMarkers” function, the marker genes of each cell type were identified. The MKI67, CD44, CD24, and PROM1, markers of tumor stem cells, were also utilized to determine the tumor stem cells. Data dimensionality were reduced by principal component analysis (PCA) and the top 20 principal components (PCs) were extracted for the next clustering analysis and Uniform Manifold Approximation and Projection for Dimension Reduction (UMAP) analysis. “CellCycleScoring” function and markers of phases were utilized to visualize the cell cycle stage. At last, “iTALK” package (Wang et al., 2019) was used to identify the ligand and receptor pairs in different cell types, and the “edgebundleR” package (https://github.com/garthtarr/edgebundleR) was used to visualize the intercellular communication.

Statistical analysis

All statistical analyses of this study were conducted by R version 3.6.1 and two-tailed p < 0.05 was required for statistical significance.

Results

Identification of DEGs and functional enrichment analysis

The analysis process of this study was presented in Supplementary Figure S1. A total of 4,941 DEGs were identified between 85 Ewing sarcoma patients and 3 normal bone samples, the expression of which was illustrated in the heatmap (Figure 1A). The volcano plot of the DEGs was illustrated in Figure 1B. GO and KEGG enrichment analyses were conducted using R’s cluster Profiler software package. The most significant GO items of biological processes (BPs), cellular components (CCs), and molecular functions (MFs) were skeletal system development, extracellular matrix, and positive regulation of cell migration, respectively (Figure 1C). Cytokine–cytokine receptor interaction, proteoglycans in cancer, and transcriptional mis-regulation in cancer were the most critical KEGG pathways, in which most DEGs were enriched (Figure 1D). Furthermore, a total of 669 eRNAs were defined as DEeRNAs between Ewing sarcoma patients and normal bone samples from 5,100 eRNAs, which were illustrated by the heatmap (Figure 1E) and volcano plot (Figure 1F). The heatmap and volcano plot of 664 DETFs identified between Ewing sarcoma patients and normal bone samples were shown in Figures 1G,H.

FIGURE 1

FIGURE 1. Identification of DEGs, DEeRNAs, and DETFs. (A) The DEG analysis between 85 Ewing sarcoma patients and 3 normal samples. (B) The volcano plot of a total of 4,941 DEGs identified between Ewing sarcoma samples and normal samples. (C) The GO enrichment analysis of DEGs. (D) The KEGG enrichment analysis of DEGs. (E) The heatmap of 669 DEeRNAs identified between Ewing sarcoma samples and normal bone samples. (F) The volcano plot of DEeRNAs. (G) The heatmap of 664 DETFs identified between Ewing sarcoma samples and normal bone samples. (H)The volcano plot of DETFs.

Construction of a prognostic prediction model in sarcoma

To evaluate the prognostic value, we used 256 sarcoma samples from TCGA database to conduct cross-validation and independent prognosis analysis. The sarcoma samples were randomly assigned into training set (156 samples) and testing set (100 samples). First off, lasso regression was applied to screen DEeRNAs to avoid overfitting (Supplementary Figures S2A,B). The univariate Cox regression analysis and multivariate Cox regression analysis of training set were displayed in Figures 2B,C. Eventually, 13 DEeRNAs significantly related to prognosis were merged into the prediction model. According to the computing formula of risk score, the sarcoma samples were sorted into high-risk group and low-risk group (Supplementary Figures S2D,E). The differential expression and heatmap of 13 prognostic-related DEeRNAs were illustrated in Figure 2A and Supplementary Figure S2C. Ranking by the risk score of each sample, the scatter dot plot, and the distribution curve were shown in Supplementary Figures S2D,E. The area under the curve (AUC) of ROC curve was 0.725 in all sets, 0.728 both in training set and testing set. Also, in all sets, the AUC of ROC curve was 0.777 at 1-year, 0.733 at 2-years, and 0.768 at 3-years. These values showed a good predictability of the prognostic prediction model in sarcoma.

FIGURE 2

FIGURE 2. Construction of prognostic prediction model in sarcoma. (A) The differential expression of 13 DEeRNAs in the prognostic prediction model. (B) The univariate Cox regression analysis of training set. (C) The multivariate Cox regression analysis of training set. (D) The summary of clinical information of 256 sarcoma samples in TCGA database. (E) The classification of sarcoma samples based on the risk score and censor group. (F) The ROC curves of the prognostic prediction model at 1-year, 2-years and 3-years (G) The GO analysis of low-risk group and high-risk group. (H) The KEGG analysis of low-risk group and high-risk group. (I) The hallmarks of cancer identified by GSVA in low-risk group and high-risk group.

The function enrichment analysis was also conducted in both high-risk group and low-risk group. In low-risk group, B cell activation and some other adaptive immune response gene sets were enriched in Go analysis; cytokine–cytokine receptor interaction and toll-like receptor signaling pathway were enriched in KEGG analysis, representatively, inflammatory response and myogenesis were enriched in hallmarks of cancer. In high-risk group, embryonic organ development and embryonic morphogenesis gene sets were detected by Go analysis, hedgehog signaling pathway and TGF beta signaling pathway were detected in KEGG analysis as well as hallmarks of cancer (Figures 2G–I).

CIBERSORT analysis and co-expression analysis of Ewing sarcoma

We explored the relationship between DEeRNA expression and cancer-infiltrating immune cells, and depicted a summary of the cell compositions in Ewing sarcoma samples and normal bone samples by CIBERSORT algorithm. The proportions of 8 immune cells in 85 Ewing sarcoma patients and 3 normal samples were presented by the bar plot, encompassing B cells, cancer-associated fibroblasts, CD4⁺ T cells, CD8⁺ T cells, endothelial cells, macrophages, NK cells, and uncharacterized cells (Figure 3A). Compared to normal bone tissues, infiltration of endothelial cells (p < 0.01) and B cells (p < 0.01) was increased, whereas infiltration of cancer-associated fibroblasts (p < 0.05), NK cells (p < 0.05), and CD4⁺ T cells (p < 0.05) was decreased in Ewing sarcoma samples, which suggested that these immune cells had a significantly prognostic value for Ewing sarcoma (Figure 3B). While, the heatmap showed the co-expression patterns between CD4⁺ T cells and CD8⁺ T cells (R = 0.53); CD4⁺ T cells and macrophages (R = 0.51); endothelial cells and macrophages (R = 0.47); CD8⁺ T cells and macrophages (R = 0.72), indicating a strong correlation between these immune cells in Ewing sarcoma (Figure 3C).

FIGURE 3

FIGURE 3. CIBERSORT analysis and co-expression analysis. (A) The proportions of 8 immune cells in 85 Ewing sarcoma patients and 3 normal samples explored by CIBERSORT analysis. (B) The immune infiltration of Ewing sarcoma and normal bone tissues. (C) The co-expression patterns of immune cells in Ewing sarcoma.

Validation of immune clustering in sarcoma

To validate the prognostic value of immune proportions, CIBERSORT was also implemented in 221 sarcoma samples in TCGA database after removing the missing data. The 22 immune fractions of sarcoma samples were displayed in Figure 4A, classified by risk score. The immune subtypes of high-risk group and low risk group sarcoma samples were displayed in Figure 4B. Specifically, T cells CD8, T cells CD4 memory resting, and macrophages M2 comprised a large proportion of immune cells (Figure 4C). On the other hand, 29 immune gene sets were quantified by ssGSEA, which was displayed in Figure 4D. Compared to low-risk group, high-risk group had lower immune function, significant in CCR, check-point, cytolytic activity, DCs, HLA, inflammation promoting, mast cells, neutrophils, NK cells, parainflammation, pDCs, T cells, TIL, Treg, and IFN response. The Kaplan–Meier survival curves below also typically displayed good correlation between immune function and survival in sarcoma samples (Figure 4E).

FIGURE 4

FIGURE 4. Validation of immune clustering in sarcoma. (A) The 22 immune fractions of sarcoma samples identified by CIBERSORT. (B) The classification of 221 sarcoma samples based on the risk score and immune subtypes. (C) The 22 immune fractions of sarcoma samples classified by risk score. (D) The function score of 29 immune gene sets identified by ssGSEA. (E) The Kaplan–Meier survival curves of immune gene sets, representatively.

Identification of DETFs, differential hallmarks of cancer, and immune gene sets

Representatively, the heatmap and volcano plot of 68 DETFs in Ewing sarcoma samples and normal bone samples were shown in Figures 5A,B. A total of 21 differential hallmarks of cancer were identified from 50 hallmark pathways between Ewing sarcoma samples and normal bone samples, which were shown in the heatmap and volcano plot (Figures 5C,D). Besides, the correlation of GSVA score of hallmark pathways and Ewing sarcoma was investigated (Figure 5E). Immune cell infiltration status was evaluated using ssGSEA to validate the associations between the Ewing sarcoma samples and normal bone samples with tumor immune characteristics. Specifically, 29 immune-related terms, or immune functions, were quantified in the heatmap to unravel the abundance of diverse immune cell types in Ewing sarcoma samples and normal bone samples (Figure 5F).

FIGURE 5

FIGURE 5. Identification of DETFs, differentially expressed hallmarks of cancer, and immune gene sets co-expressed with DEeRNAs. (A) The heatmap of 68 DETFs in Ewing sarcoma samples and normal bone samples. (B) The volcano plot of DETFs. (C) The heatmap of 21 differentially expressed hallmarks of cancer identified between Ewing sarcoma samples and normal bone samples. (D) The volcano plot of differentially expressed hallmarks of cancer. (E) The correlation of GSVA score of hallmark pathways and Ewing sarcoma. (F) The heatmap of 29 immune-related terms evaluated by ssGSEA between the Ewing sarcoma samples and normal bone samples.

The network construction and Connectivity Map Analysis

After the co-expressed analysis, the heatmap showed the expression of most significant DEeRNAs, DETFs, and DETGs in Figure 6A. A total of six different dimension regulatory network was constructed with 17 DEeRNAs, 29 DETFs, 9 DETGs, 5 immune cells by CIBERSORT, 24 immune gene sets by ssGSEA, and 8 hallmarks of cancer by GSVA, which showed the potential regulatory relationships across these factors (Figure 6B). Four key DEeRNAs (CCR1, CD3D, PHLDA1, and RASD1) showed significant co-expression relationships in the six different dimension regulatory network. We supposed that these DEeRNAs may play crucial roles in the tumorigenesis of Ewing sarcoma. Furthermore, the interaction coefficients among these components were illustrated by the heatmap by Pearson correlation analysis (Figure 6C).

FIGURE 6

FIGURE 6. The network construction and Connectivity Map Analysis. (A) The heatmap of DEeRNAs, DETFs, and DETGs. (B) The six different dimension regulatory network, encompassing 17 DEeRNAs, 29 DETFs, 9 DETGs, 5 immune cells by CIBERSORT, 24 immune gene sets by ssGSEA, and 8 hallmarks of cancer by GSVA. (C) The cor-expression heatmap of these components above. (D) The CMap analysis of Ewing sarcoma as well as other 33 cancer types.

The heatmap depicted the enrichment score of each compound analyzed by CMap in Ewing sarcoma, as well in other 33 cancer types (Malta et al., 2018). Importantly, MS-275 and pyrvinium with the highest specificity and the lowest p value were considered as the best compounds that might target Ewing sarcoma (Figure 6D).

ATAC-seq and external validation

Figure 7 depicted the accessible chromatin sites at the key DEeRNAs, including CCR1, CD3D, PHLDA1, and RASD1 (Figures 7A–D). Furthermore, we analyzed data from public databases to assess the prognostic effects of key DEeRNAs and potential regulatory mechanisms in Ewing sarcoma. Based on the human protein atlas database, we examined the expression level of CD3D, MAZ, and PHLDA1 by immunohistochemical (IHC) staining assay and observed that there was medium expression of CD3D, MAZ, and PHLDA1 in normal tissue. Representative IHC images were presented in Supplementary Figure S3. Based on Oncomine database, we identified that expression of CCR1, CD3D, MAZ, PHLDA1, and RASD1 was higher in tumor than that in normal tissue at the pan-cancer level (Supplementary Figures S4A–E). Additionally, expression of CCR1, CD3D, MAZ, PHLDA1, and RASD1 in various different tissues was determined based on CCLE database (Supplementary Figures S5A–E). Taken together, the expression of CCR1, CD3D, PHLDA1, and RASD1 in multiple databases were summarized in Supplementary Table S1. In cBioPortal database, the correlation of mutation count and overall survival of key DEeRNAs were shown in Supplementary Figures S6A–E. Also, the survival analysis of PHLDA1 in GEPIA database was displayed in Supplementary Figure S6G.

FIGURE 7

FIGURE 7. ATAC-seq analysis of key DEeRNAs. (A) The accessible chromatin sites of CCR1 analyzed by ATAC-seq. (B) The accessible chromatin sites of CD3D analyzed by ATAC-seq. (C) The accessible chromatin sites of PHLDA1 analyzed by ATAC-seq. (D)The accessible chromatin sites of RASD1 analyzed by ATAC-seq.

To explore the role of enhancer-specific histone in modifications of eRNA transcription, ChIP-seq data of H3K27ac were downloaded and analyzed. The UCSC Genome Browser tracks showed enrichment of H3K27ac on multiple loci in the DEeRNAs (CCR1, CD3D, PHLDA1, and RASD1) (Supplementary Figures S7–S10). Results of external validation in eRic database showed the specific chromatin localization and target genes of the four key DEeRNAs (CCR1, CD3D, PHLDA1, and RASD1), as well as potential drugs that may target these DEeRNAs in different cancers (Supplementary Table S2). We further investigated the expression of key DEeRNAs between tumor and normal samples among different cancer types and identified that the expression level of CCR1, CD3D, PHLDA1, and RASD1 was significantly up-regulated in tumor tissue, as compared with normal tissue. Additionally, prognostic effect of key DEeRNAs was displayed between high expression group and low expression group (Supplementary Figures S11–S14).

Single-cell RNA-seq transcriptome analysis

Unsupervised clustering clearly identified 12 cell clusters (Figure 8A, left). By utilizing the expression of differentially expressed marker genes, we attributed these clusters to 3 Ewing sarcoma cell lines (CHLA9, CHLA10, and TC71) based on hierarchical similarities (Figure 8A). The heatmap displayed the up- or down-regulated genes in the 12 clusters (Figure 8B). The dot plots showed the proportion of cells expressing tumor stemness-related gene markers (CD44 and MKI67) and key DEeRNAs (PHLDA1 and RASD1) and their scaled relative expression level in 12 cell clusters (Figure 8C). Specifically, MKI67 (a known nuclear marker of proliferation) was highly expressed in all cell clusters, indicating high cellular proliferative activities in these cancer cells. The cell number and proportions of 3 main cell subtypes were quite diverse among the 12 cell clusters (Figure 8D). As demonstrated in the top DEGs, specifically cancer stem lineage clusters expressed high levels of stemness feature genes (MKI67, CD44, CD24, and PROM1) and key DEeRNAs (PHLDA1 and RASD1) (Figure 9A). The cell cycle distribution of 12 cell clusters was shown in the UMAP plot (Figure 9B). Cells within cluster 5 were mainly in G2 phase while cells in cluster 3 were mainly in S phase. The ligand-receptor plot displayed ligand-receptor pairs among those clusters (Figure 9C). All these results showed that PHLDA1 and RASD1 (key DEeRNAs) were extensively expressed in cancer stem cells of Ewing sarcoma, which were potential targets for tumor treatment.

FIGURE 8

FIGURE 8. Single-cell transcriptomic analysis of Ewing sarcoma cell lines. (A) The distribution of 12 clusters in 3 Ewing sarcoma cell lines (CHLA9, CHLA10, and TC71). (B) Gene co-expression of top 5 genes in clusters. (C) Significantly up- or down-regulated genes in clusters. (D) Cell number and cell proportion of 12 clusters in 3 cell lines.

FIGURE 9

FIGURE 9. The key biomarkers extensively expressed in Ewing sarcoma stem cell. (A) The distribution of marker genes, and MKI67, PHLDA1 and RASD1 were extensively high expressed in all clusters. (B) The cell cycle distribution and the cell cycle score in 12 clusters. (C) The ligand-receptor pairs among 12 clusters.

Discussion

Ewing sarcoma is the second common bone cancer, with strikingly low 5-year overall survival after metastasis (Gaspar et al., 2015). Ewing sarcoma is characteristic with a recurrent chromosomal translocation and the EWS-FLI fusion may provide a new way in the targeted therapy of Ewing sarcoma (Cidre-Aranaz and Alonso, 2015). eRNAs are generated during the transcription of active enhancer (Zhang et al., 2019). In human cancers, eRNAs are specific to tumor types (Lee et al., 2020). Various eRNAs have been demonstrated to be differentially expressed in prostate cancer (Zhao et al., 2016). In breast cancer cells, estrogen-induced transcription of eRNAs was identified to be significantly upregulated (Crudele et al., 2020). Conversely, a recent study showed that expression of eRNAs was significantly decreased in throat cancer (Zhang et al., 2019). Collectively, activation of oncogenes or oncogenic pathways was associated with aberrant generation of eRNAs in human cancers, and eRNAs may play a broad role in the pathophysiology of Ewing sarcoma.

To the best of our knowledge, this is the first study to show DEeRNAs which are potentially engaged in the cellular transition from the normal cells into malignant cells and their potential regulatory relationships in Ewing sarcoma. Herein, an integrated bioinformatics analysis was performed to determine differential eRNA and target gene expression between Ewing sarcoma and normal samples. Differentially infiltrating immune cells were detected by CIBERSORT between Ewing sarcoma samples and normal bone samples. To verify the prognostic power of DEeRNAs and immune proportions, we used sarcoma samples from TCGA database into cross-validation and independent prognosis analysis, as well as Kaplan–Meier survival analysis. In addition, we constructed a DEeRNA co-expressed regulatory network of Ewing sarcoma, encompassing 17 DEeRNAs, 29 DETFs, 9 DETGs, 5 immune cells by CIBERSORT, 24 immune gene sets by ssGSEA, and 8 hallmarks of cancer by GSVA. Importantly, four DEeRNAs (CCR1, CD3D, PHLDA1, and RASD1) were considered to have significant co-expression relationships in the six different dimension regulatory networks. Moreover, Connectivity Map Analysis was applied to pursue small molecules targeting Ewing sarcoma. ATAC-seq data were utilized to provide information on chromatin accessibility of key DEeRNAs. In the end, external validation based on multidimensional online databases and scRNA-seq analysis were used to verify our key findings, which showed that the screened DEeRNAs play a critical role in the tumorigenesis of Ewing sarcoma and could be utilized as important reference markers for future research.

The signal axes of four key eRNAs (CCR1, CD3D, PHLDA1, and RASD1) were as follows: BATF-CCR1-complemen; BATF-CD3D-allograft rejection; FOSL2-RASD1-tnfa signaling via NFKB; and FOSL2-PHLDA1-FOXC1-tnfa signaling via NFKB. Importantly, signal axis FOSL2-PHLDA1-FOXC1-TNFA signaling via NFkB was extracted for the subsequent analyses by theoretical basis and literature review, which will be explained in detail in the following sections as potential mechanism related to the tumorigenesis of Ewing sarcoma. The correlation coefficient between FOSL2 and PHLDA1 was 0.89 (p < 0.001); between PHLDA1 and FOXC1 was 0.87 (p < 0.001); between PHLDA1 and TNFA signaling via NFkB was 0.70605164 (p < 0.001). In the interaction and correlation network, PHLDA1 was also related to cancer-associated fibroblasts (R = 0.67; p < 0.001) and CCR (R = 0.53; p < 0.001).

In Ewing sarcoma, EWS-FLI fusions encode oncogenic proteins functioning as a transcription factor regulating abnormal transcription (Sanchez et al., 2008). Well, a number of studies have described target genes mediated by EWS/ETS proteins. In particular, PHLDA1 has been reported to be few target genes that are directly repressed by the binding of EWS/FLI1 through meta-analysis and experiments in vitro (Boro et al., 2012). PHLDA1 (pleckstrin homology-like domain family, member 1) gene is one of the members of the PHLDA gene family (Frank et al., 1999), which has been reported to suppress tumorigenesis (Chen et al., 2018). To be specific, PHLDA1 may repress tumorigenesis by inducing apoptosis and inhibiting cell growth (Neef et al., 2002; Chen et al., 2018). In melanoma (Neef et al., 2002), breast cancer (Nagai et al., 2007), oral cancer (Coutinho-Camillo et al., 2013), and stomach cancers (Zhao et al., 2015), the reduced expression of PHLDA1 has already been described. Moreover, PHLDA1 is not only a tumor suppressor, but also a new targeted therapy to re-sensitize drug-resistant cancer cells (Fearon et al., 2018).

FOS-like antigen 2 (FOSL2) is a member of activator protein-1 (AP-1) transcription factor family (Tulchinsky, 2000), which is involved in cell proliferation, transformation, and death (Shaulian and Karin, 2002). FOSL2 plays a key role in bone development (Bozec et al., 2013). FOSL2 is expressed in stromal cells of human chondroblastic and osteoblastic osteosarcomas, and the deficiency of FOSL2 induces a differentiation defect in osteoblasts both in vivo and in vitro experiments (Bozec et al., 2010). In addition, FOSL2 has been reported to exert a specific function of mediating TGF-β pathway in extracellular matrix (ECM) remodeling (Busnadiego et al., 2013) and in non-small cell lung cancer (Wang J et al., 2014). In adult T-cell leukemia, aberrantly expressed FOSL2 has been demonstrated to induce CCR4 expression MDM2 (Nakayama et al., 2008).

The Forkhead box C1 (FOXC1) is a member of the Forkhead box (FOX) family, a group of transcription factors and the Fox family are involved in cellular proliferation, differentiation, and death (Lehmann et al., 2003). As a consequence, the deregulation of FOX proteins is able to promote tumorigenesis and cancer progression (M yatt and Lam, 2007). Recently, FOXC1 is demonstrated to be a critical transcriptional regulator for the development and maintenance of hematopoietic stem and progenitor cells (HSPCs) in bone marrow (Omatsu and Nagasawa, 2015). FOXC1 is preferentially expressed to maintain haematopoietic stem and progenitor cells in the adipoosteogenic progenitor CAR cells of developing adult bone marrow (Omatsu et al., 2014). FOXC1 is also able to inhibit CAR cell differentiation into adipocytes, by upregulating CXCL12 and stem cell factor (SCF) (Omatsu et al., 2014). On the other hand, FOXC1 is responsible for governing quiescence by the nuclear factor of activated T-cells 1 (NFATC1) and BMP signaling in stem cells (Wang et al., 2016). Similarly, in basal-like breast cancer (BLBC), FOXC1 may increase cancer stem cell (CSC) properties by cellular mechanisms (Han et al., 2015).

Tumor necrosis factor alpha (TNF‐α) is a cytokine produced by activated macrophages, T lymphocytes, and natural killer (NK) cells, and exerts a wide function in cellular apoptosis and survival, as well as inflammation and immunity. Also, TNF-α is now used in isolated limb perfusion for treatment of soft tissue sarcoma (STS) and other large tumors (Eggermont et al., 2003). Through the activation of nuclear transcription factors, such as NFkB (nuclear factor kappa B) and AP-1, TNF-α is able to modulate the expression of a majority of different genes (Schutze et al., 1992). However, NFkB plays a critical role in preventing cell death induced by TNF-α (Beg and Baltimore, 1996). Aberrant NF-kB expression has been described in many human cancers and tips apoptosis–proliferation balance toward malignant growth (Lin and Karin, 2003).

CC chemokine receptors include CCR 1-10 and CC chemokines are ligands to CCR1-10. The movement of immune cells is driven by CC chemokine receptors and CC chemokines (Hughes and Nibbs, 2018). In cancer, the expression of CC chemokine receptors promotes metastasis and may provide new targets for cancer immunotherapy (Mollica Poeta et al., 2019). CCR7 mediates lymphocyte migration, and CCR9 is involved in rare metastases to the small intestine in melanoma (Zlotnik et al., 2011). In Ewing sarcoma, the expression of CCR5-ligand, CCL5, is positively related to the number of infiltrating CD8⁺ T-lymphocyte and patients with high numbers of infiltrating T-lymphocytes have an overall survival advantage (Berghuis et al., 2011).

Signal axis, FOSL2-PHLDA1-FOXC1-TNFA signaling via NFkB, is first reported to be associated with tumorigenesis and progression of Ewing sarcoma. All that being said, the limitations of bioinformatics study are easy to see and unavoidable. First, the sample size of Ewing sarcoma and normal bone samples in our study was limited. Although our results were validated by sarcoma samples, scRNA-seq analysis and multidimensional online databases, larger sample size and more comprehensive data are needed to get more reliable and more accurate results. Second, the direct regulating mechanism of FOSL2-PHLDA1-FOXC1-TNFA signaling via NFkB in Ewing sarcoma is unclear. Laboratory-based experiment and clinical study are tremendously needed to explore the direct-action mechanism of the signal axis in Ewing sarcoma. Our hypothesis may just provide a new way for the treatment of Ewing sarcoma.

Conclusion

In summary, we presume PHLDA1 is a key regulator in the tumorigenesis and progression of Ewing sarcoma. PHLDA1 is directly repressed by the binding of EWS/FLI1 protein and low expression of FOSL2, resulting in the deregulation of FOX proteins and CC chemokine receptors. T-lymphocytes expressing less CC chemokine receptors may not migrate to the tumor site. Inhibition of infiltrating T-lymphocytes and TNFA signaling may promote tumorigenesis and progression of Ewing sarcoma.

Data availability statement

Publicly available datasets were analyzed in this study. These data can be found at: Treehouse database (https://treehousegenomics.soe.ucsc.edu/public-data/#datasets), SRA database (https://www.ncbi.nlm.nih.gov/sra/), TCGA database (https://tcga-data.nci.nih.gov) and GEO database (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE146221).

Ethics statement

The study was approved by the Ethics Committee of Shanghai General Hospital, School of Medicine, Shanghai Jiaotong University.

Author contributions

Conception/design: RH, DH, SW, SX, YL, MJ, XZ, SC, XY, WZ, JL, HL, ZH, HZ, and HY. Collection and/or assembly of data: RH, DH, SW, SX, YL, MJ, XZ, SC, XY, WZ, JL, HL, ZH, HZ, and HY. Data analysis and interpretation: RH, DH, SW, SX, YL, MJ, XZ, SC, XY, WZ, JL, HL, ZH, HZ, and HY. Manuscript writing: RH, DH, SW, SX, YL, MJ, XZ, SC, XY, WZ, JL, HL, ZH, HZ, and HY. Final approval of the manuscript: RH, DH, SW, SX, YL, MJ, XZ, SC, XY, WZ, JL, HL, ZH, HZ, and HY.

Funding

This study was supported in part by the National Natural Science Foundation of China (Grant No. 81772856); Youth Fund of Shanghai Municipal Health Planning Commission (No.2017YQ054; 2017Y0117); Interdisciplinary Program of Shanghai Jiao Tong University (No.YG2017MS26); Shanghai Talent Development Fund (No.2018094); Henan medical Science and Technology Research Project (No. 201602031); and Key Project of Provincial and Ministerial Co-construction of Henan Medical Science and Technology (No. SBGJ202002031). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Acknowledgments

The authors appreciate the Treehouse database team, SRA database team, and TCGA database team for providing an open-access resource for the study.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The handling editor WX declared a shared parent affiliation with the author HZ at the time of review.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2022.952162/full#supplementary-material

Supplementary Figure S1 | The schematic diagram of analytic processing.

Supplementary Figure S2 | Cross-validation and model diagnosis of sarcoma. (A) Coefficient curves of variables in the prognostic prediction model. (B) The correlation between λ and partial likehood deviance. At the variable number of 13, the partial likehood deviance was lowest. (C) The heatmap of 13 DEeRNAs in all sets, training set, and testing set. (D) The scatter dot plot of all sets, training set, and testing set. (E) The distribution curve of risk score in all sets, training set, and testing set. (F) The area under the curve (AUC) of ROC curve was 0.725 in all sets, 0.728 both in training set and testing set.

Supplementary Figure S3 | Human protein atlas database. The medium expression of CD3D (A), MAZ (B), and PHLDA1 (C) in normal tissue by immunohistochemical staining assay.

Supplementary Figure S4 | Oncomine database. The expression of CCR1 (A), CD3D (B), MAZ (C), PHLDA1 (D), and RASD1 (E) in tumor and normal tissues at the pan-cancer level.

Supplementary Figure S5 | CCLE database. The expression of CCR1 (A), CD3D (B), MAZ (C), PHLDA1 (D), and RASD1 (E) in various different tissues.

Supplementary Figure S6 | cBioPortal database and GEPIA database. The correlation of mutation count of CCR1 (A), CD3D (B), RASD1 (C), PHLDA1 (D), and MAZ (E) and overall survival. (F) The summary of patients’ clinical information and genetic alteration of DEeRNAs. (G) The survival analysis of PHLDA1.

Supplementary Figure S7 | The ChIP-seq data of CCR1 analyzed by UCSC Genome Browser.

Supplementary Figure S8 | The ChIP-seq data of CD3D analyzed by UCSC Genome Browser.

Supplementary Figure S9 | The ChIP-seq data of PHLDA1 analyzed by UCSC Genome Browser.

Supplementary Figure S10 | The ChIP-seq data of RASD1 analyzed by UCSC Genome Browser.

Supplementary Figure S11 | The expression and prognostic effect of CCR1 between tumor and normal samples among different cancer types on eRic database.

Supplementary Figure S12 | The expression and prognostic effect of CD3D between tumor and normal samples among different cancer types on eRic database.

Supplementary Figure S13 | The expression and prognostic effect of PHLDA1 between tumor and normal samples among different cancer types on eRic database.

Supplementary Figure S14 | The expression and prognostic effect of RASD1 between tumor and normal samples among different cancer types on eRic database.

References

Andersson, R., Gebhard, C., Miguel-Escalada, I., Hoof, I., Bornholdt, J., Boyd, M., et al. (2014). An atlas of active enhancers across human cell types and tissues. Nature 507 (7493), 455–461. doi:10.1038/nature12787

PubMed Abstract | CrossRef Full Text | Google Scholar

Arner, E., Daub, C. O., Vitting-Seerup, K., Andersson, R., Lilje, B., Drablos, F., et al. (2015). Transcribed enhancers lead waves of coordinated transcription in transitioning mammalian cells. Science 347 (6225), 1010–1014. doi:10.1126/science.1259418

PubMed Abstract | CrossRef Full Text | Google Scholar

Aurias, A., Buffe, D., Zucker, J. M., and Mazabraud, A. (1984). Translocation involving chromosome 22 in Ewing's sarcoma. a cytogenetic study of four fresh tumors. Cancer Genet. Cytogenet. 12 (1), 21–25. doi:10.1016/0165-4608(84)90003-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Barbie, D. A., Tamayo, P., Boehm, J. S., Kim, S. Y., Moody, S. E., Dunn, I. F., et al. (2009). Systematic RNA interference reveals that oncogenic KRAS-driven cancers require TBK1. Nature 462 (7269), 108–112. doi:10.1038/nature08460

PubMed Abstract | CrossRef Full Text | Google Scholar

Beg, A. A., and Baltimore, D. (1996). An essential role for NF-kappaB in preventing TNF-alpha-induced cell death. Science 274 (5288), 782–784. doi:10.1126/science.274.5288.782

PubMed Abstract | CrossRef Full Text | Google Scholar

Berghuis, D., Santos, S. J., Baelde, H. J., Taminiau, A. H., Egeler, R. M., Schilham, M. W., et al. (2011). Pro-inflammatory chemokine-chemokine receptor interactions within the Ewing sarcoma microenvironment determine CD8(+) T-lymphocyte infiltration and affect tumour progression. J. Pathol. 223 (3), 347–357. doi:10.1002/path.2819

PubMed Abstract | CrossRef Full Text | Google Scholar

Blackwood, E. M., and Kadonaga, J. T. (1998). Going the distance: a current view of enhancer action. Science 281 (5373), 60–63. doi:10.1126/science.281.5373.60

PubMed Abstract | CrossRef Full Text | Google Scholar

Boro, A., Pretre, K., Rechfeld, F., Thalhammer, V., Oesch, S., Wachtel, M., et al. (2012). Small-molecule screen identifies modulators of EWS/FLI1 target gene expression and cell survival in Ewing's sarcoma. Int. J. Cancer 131 (9), 2153–2164. doi:10.1002/ijc.27472

PubMed Abstract | CrossRef Full Text | Google Scholar

Bozec, A., Bakiri, L., Jimenez, M., Schinke, T., Amling, M., and Wagner, E. F. (2010). Fra-2/AP-1 controls bone formation by regulating osteoblast differentiation and collagen production. J. Cell Biol. 190 (6), 1093–1106. doi:10.1083/jcb.201002111

PubMed Abstract | CrossRef Full Text | Google Scholar

Bozec, A., Bakiri, L., Jimenez, M., Rosen, E. D., Catala-Lehnen, P., Schinke, T., et al. (2013). Osteoblast-specific expression of Fra-2/AP-1 controls adiponectin and osteocalcin expression and affects metabolism. J. Cell Sci. 126 (Pt 23), 5432–5440. doi:10.1242/jcs.134510

PubMed Abstract | CrossRef Full Text | Google Scholar

Busnadiego, O., Gonzalez-Santamaria, J., Lagares, D., Guinea-Viniegra, J., Pichol-Thievend, C., Muller, L., et al. (2013). LOXL4 is induced by transforming growth factor β1 through Smad and JunB/Fra2 and contributes to vascular matrix remodeling. Mol. Cell. Biol. 33 (12), 2388–2401. doi:10.1128/MCB.00036-13

PubMed Abstract | CrossRef Full Text | Google Scholar

Cerami, E., Gao, J., Dogrusoz, U., Gross, B. E., Sumer, S. O., Aksoy, B. A., et al. (2012). The cBio cancer genomics portal: an open platform for exploring multidimensional cancer genomics data. Cancer Discov. 2 (5), 401–404. doi:10.1158/2159-8290.CD-12-0095

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, Y., Takikawa, M., Tsutsumi, S., Yamaguchi, Y., Okabe, A., Shimada, M., et al. (2018). PHLDA1, another PHLDA family protein that inhibits Akt. Cancer Sci. 109 (11), 3532–3542. doi:10.1111/cas.13796

PubMed Abstract | CrossRef Full Text | Google Scholar

Cidre-Aranaz, F., and Alonso, J. (2015). EWS/FLI1 target genes and therapeutic opportunities in ewing sarcoma. Front. Oncol. 5, 162. doi:10.3389/fonc.2015.00162

PubMed Abstract | CrossRef Full Text | Google Scholar

Corces, M. R., Granja, J. M., Shams, S., Louie, B. H., Seoane, J. A., Zhou, W., et al. (2018). The chromatin accessibility landscape of primary human cancers. Science 362 (6413), eaav1898. doi:10.1126/science.aav1898

PubMed Abstract | CrossRef Full Text | Google Scholar

Coutinho-Camillo, C. M., Lourenco, S. V., Nonogaki, S., Vartanian, J. G., Nagai, M. A., Kowalski, L. P., et al. (2013). Expression of PAR-4 and PHLDA1 is prognostic for overall and disease-free survival in oral squamous cell carcinomas. Virchows Arch. 463 (1), 31–39. doi:10.1007/s00428-013-1438-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Crudele, F., Bianchi, N., Reali, E., Galasso, M., Agnoletto, C., and Volinia, S. (2020). The network of non-coding RNAs and their molecular targets in breast cancer. Mol. Cancer 19 (1), 61. doi:10.1186/s12943-020-01181-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Delattre, O., Zucman, J., Plougastel, B., Desmaze, C., Melot, T., PeterM., , et al. (1992). Gene fusion with an ETS DNA-binding domain caused by chromosome translocation in human tumours. Nature 359 (6391), 162–165. doi:10.1038/359162a0

PubMed Abstract | CrossRef Full Text | Google Scholar

Eggermont, A. M., de Wilt, J. H., and Ten Hagen, T. L. (2003). Current uses of isolated limb perfusion in the clinic and a model system for new strategies. Lancet. Oncol. 4 (7), 429–437. doi:10.1016/s1470-2045(03)01141-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Fearon, A. E., Carter, E. P., Clayton, N. S., Wilkes, E. H., Baker, A. M., Kapitonova, E., et al. (2018). PHLDA1 mediates drug resistance in receptor tyrosine kinase-driven cancer. Cell Rep. 22 (9), 2469–2481. doi:10.1016/j.celrep.2018.02.028

PubMed Abstract | CrossRef Full Text | Google Scholar

Ferguson, J. L., and Turner, S. P. (2018). Bone cancer: diagnosis and treatment principles. Am. Fam. Physician 98 (4), 205–213.

PubMed Abstract | Google Scholar

Frank, D., Mendelsohn, C. L., CicconE, E., Svensson, K., Ohlsson, R., and Tycko, B. (1999). A novel pleckstrin homology-related gene family defined by ipl/tssc3, TDAG51, and Tih1: tissue-specific expression, chromosomal location, and parental imprinting. Mamm. Genome 10 (12), 1150–1159. doi:10.1007/s003359901182

PubMed Abstract | CrossRef Full Text | Google Scholar

Gaspar, N., Hawkins, D. S., Dirksen, U., Lewis, I. J., Ferrari, S., Le Deley, M. C., et al. (2015). Ewing sarcoma: current management and future approaches through collaboration. J. Clin. Oncol. 33 (27), 3036–3046. doi:10.1200/JCO.2014.59.5256

PubMed Abstract | CrossRef Full Text | Google Scholar

Ghandi, M., Huang, F. W., Jane-Valbuena, J., Kryukov, G. V., Lo, C. C., McDonald, E. R., et al. (2019). Next-generation characterization of the cancer cell line Encyclopedia. Nature 569 (7757), 503–508. doi:10.1038/s41586-019-1186-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Grunewald, T. G. P., Cidre-Aranaz, F., Surdez, D., Tomazou, E. M., de Alava, E., Kovar, H., et al. (2018). Ewing sarcoma. Nat. Rev. Dis. Prim. 4 (1), 5. doi:10.1038/s41572-018-0003-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Han, B., Qu, Y., Jin, Y., Yu, Y., Deng, N., Wawrowsky, K., et al. (2015). FOXC1 activates smoothened-independent hedgehog signaling in basal-like breast cancer. Cell Rep. 13 (5), 1046–1058. doi:10.1016/j.celrep.2015.09.063

PubMed Abstract | CrossRef Full Text | Google Scholar

Hanzelmann, S., Castelo, R., and Guinney, J. (2013). GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinform. 14, 7. doi:10.1186/1471-2105-14-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Hughes, C. E., and Nibbs, R. J. B. (2018). A guide to chemokines and their receptors. FEBS J. 285 (16), 2944–2971. doi:10.1111/febs.14466

PubMed Abstract | CrossRef Full Text | Google Scholar

Jawad, M. U., Cheung, M. C., Min, E. S., Schneiderbauer, M. M., Koniaris, L. G., and Scully, S. P. (2009). Ewing sarcoma demonstrates racial disparities in incidence-related and sex-related differences in outcome: an analysis of 1631 cases from the SEER database, 1973-2005. Cancer 115 (15), 3526–3536. doi:10.1002/cncr.24388

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, T. K., Hemberg, M., Gray, J. M., Costa, A. M., Bear, D. M., Wu, J., et al. (2010). Widespread transcription at neuronal activity-regulated enhancers. Nature 465 (7295), 182–187. doi:10.1038/nature09033

PubMed Abstract | CrossRef Full Text | Google Scholar

Lamb, J., Crawford, E. D., Peck, D., Modell, J. W., Blat, I. C., Wrobel, M. J., et al. (2006). The connectivity map: using gene-expression signatures to connect small molecules, genes, and disease. Science 313 (5795), 1929–1935. doi:10.1126/science.1132939

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, J. H., Xiong, F., and Li, W. (2020). Enhancer RNAs in cancer: regulation, mechanisms and therapeutic potential. RNA Biol. 17 (11), 1550–1559. doi:10.1080/15476286.2020.1712895

PubMed Abstract | CrossRef Full Text | Google Scholar

Lehmann, O. J., Sowden, J. C., Carlsson, P., Jordan, T., and Bhattacharya, S. S. (2003). Fox's in development and disease. Trends Genet. 19 (6), 339–344. doi:10.1016/S0168-9525(03)00111-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, C., Tang, Z., Zhang, W., Ye, Z., and Liu, F. (2021). GEPIA2021: integrating multiple deconvolution-based analysis into GEPIA. Nucleic Acids Res. 49 (W1), W242–W246. doi:10.1093/nar/gkab418

PubMed Abstract | CrossRef Full Text | Google Scholar

Liberzon, A., Birger, C., Thorvaldsdottir, H., Ghandi, M., Mesirov, J. P., and Tamayo, P. (2015). The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell Syst. 1 (6), 417–425. doi:10.1016/j.cels.2015.12.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, A., and Karin, M. (2003). NF-kappaB in cancer: a marked target. Semin. Cancer Biol. 13 (2), 107–114. doi:10.1016/s1044-579x(02)00128-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Myatt, S. S., and Lam, E. W. (2007). The emerging roles of forkhead box (Fox) proteins in cancer. Nat. Rev. Cancer 7 (11), 847–859. doi:10.1038/nrc2223

PubMed Abstract | CrossRef Full Text | Google Scholar

Malta, T. M., Sokolov, A., Gentles, A. J., Burzykowski, T., Poisson, L., Weinstein, J. N., et al. (2018). Machine learning identifies stemness features associated with oncogenic dedifferentiation. Cell 173 (2), 338–354. doi:10.1016/j.cell.2018.03.034

PubMed Abstract | CrossRef Full Text | Google Scholar

May, W. A., Gishizky, M. L., Lessnick, S. L., Lunsford, L. B., Lewis, B. C., Delattre, O., et al. (1993). Ewing sarcoma 11;22 translocation produces a chimeric transcription factor that requires the DNA-binding domain encoded by FLI1 for transformation. Proc. Natl. Acad. Sci. U. S. A. 90 (12), 5752–5756. doi:10.1073/pnas.90.12.5752

PubMed Abstract | CrossRef Full Text | Google Scholar

Miller, H. E., Gorthi, A., Bassani, N., Lawrence, L. A., Iskra, B. S., and Bishop, A. J. R. (2020). Reconstruction of ewing sarcoma developmental context from mass-scale transcriptomics reveals characteristics of EWSR1-FLI1 permissibility. Cancers (Basel) 12 (4), E948. doi:10.3390/cancers12040948

PubMed Abstract | CrossRef Full Text | Google Scholar

Mollica Poeta, V., Massara, M., Capucetti, A., and Bonecchi, R. (2019). Chemokines and chemokine receptors: new targets for cancer immunotherapy. Front. Immunol. 10, 379. doi:10.3389/fimmu.2019.00379

PubMed Abstract | CrossRef Full Text | Google Scholar

Nagai, M. A., Fregnani, J. H. T. G., Netto, M. M., Brentani, M. M., and Soares, F. A. (2007). Down-regulation of PHLDA1 gene expression is associated with breast cancer progression. Breast Cancer Res. Treat. 106 (1), 49–56. doi:10.1007/s10549-006-9475-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Nakayama, T., Hieshima, K., Arao, T., Jin, Z., Nagakubo, D., Shirakawa, A. K., et al. (2008). Aberrant expression of Fra-2 promotes CCR4 expression and cell proliferation in adult T-cell leukemia. Oncogene 27 (23), 3221–3232. doi:10.1038/sj.onc.1210984

PubMed Abstract | CrossRef Full Text | Google Scholar

Neef, R., Kuske, M. A., Prols, E., and Johnson, J. P. (2002). Identification of the human PHLDA1/TDAG51 gene: down-regulation in metastatic melanoma contributes to apoptosis resistance and growth deregulation. Cancer Res. 62 (20), 5920–5929.

PubMed Abstract | Google Scholar

Omatsu, Y., and Nagasawa, T. (2015). The critical and specific transcriptional regulator of the microenvironmental niche for hematopoietic stem and progenitor cells. Curr. Opin. Hematol. 22 (4), 330–336. doi:10.1097/MOH.0000000000000153

PubMed Abstract | CrossRef Full Text | Google Scholar

Omatsu, Y., Seike, M., Sugiyama, T., Kume, T., and Nagasawa, T. (2014). Foxc1 is a critical regulator of haematopoietic stem/progenitor cell niche formation. Nature 508 (7497), 536–540. doi:10.1038/nature13071

PubMed Abstract | CrossRef Full Text | Google Scholar

Rhodes, D. R., Yu, J., ShanKer, K., Deshpande, N., Varambally, R., Ghosh, D., et al. (2004). ONCOMINE: a cancer microarray database and integrated data-mining platform. Neoplasia 6 (1), 1–6. doi:10.1016/s1476-5586(04)80047-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Riggi, N., Knoechel, B., Gillespie, S. M., Rheinbay, E., Boulay, G., Suva, M. L., et al. (2014). EWS-FLI1 utilizes divergent chromatin remodeling mechanisms to directly activate or repress enhancer elements in Ewing sarcoma. Cancer Cell 26 (5), 668–681. doi:10.1016/j.ccell.2014.10.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Sanchez, G., Bittencourt, D., Laud, K., Barbier, J., Delattre, O., Auboeuf, D., et al. (2008). Alteration of cyclin D1 transcript elongation by a mutated transcription factor up-regulates the oncogenic D1b splice isoform in cancer. Proc. Natl. Acad. Sci. U. S. A. 105 (16), 6004–6009. doi:10.1073/pnas.0710748105

PubMed Abstract | CrossRef Full Text | Google Scholar

Schutze, S., Machleidt, T., and Kronke, M. (1992). Mechanisms of tumor necrosis factor action. Semin. Oncol. 19 (2 Suppl. 4), 16–24.

PubMed Abstract | Google Scholar

Shaulian, E., and Karin, M. (2002). AP-1 as a regulator of cell life and death. Nat. Cell Biol. 4 (5), E131–E136. doi:10.1038/ncb0502-e131

PubMed Abstract | CrossRef Full Text | Google Scholar

Smyth, G. K. (2004). Linear models and empirical bayes methods for assessing differential expression in microarray experiments. Stat. Appl. Genet. Mol. Biol. 3, Article3. doi:10.2202/1544-6115.1027

PubMed Abstract | CrossRef Full Text | Google Scholar

Tang, Z., Li, C., Kang, B., Gao, G., Li, C., and Zhang, Z. (2017). GEPIA: a web server for cancer and normal gene expression profiling and interactive analyses. Nucleic Acids Res. 45 (W1), W98–W102. doi:10.1093/nar/gkx247

PubMed Abstract | CrossRef Full Text | Google Scholar

Tulchinsky, E. (2000). Fos family members: regulation, structure and role in oncogenic transformation. Histol. Histopathol. 15 (3), 921–928. doi:10.14670/HH-15.921

PubMed Abstract | CrossRef Full Text | Google Scholar

Uhlen, M., Fagerberg, L., Hallstrom, B. M., Lindskog, C., Oksvold, P., Mardinoglu, A., et al. (2015). Proteomics. tissue-based map of the human proteome. Science 347 (6220), 1260419. doi:10.1126/science.1260419

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, J., Sun, D., Wang, Y., Ren, F., Pang, S., Wang, D., et al. (2014). FOSL2 positively regulates TGF-β1 signalling in non-small cell lung cancer. PLoS One 9 (11), e112150. doi:10.1371/journal.pone.0112150

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Q., Huang, J., Sun, H., Liu, J., Wang, J., Qin, Q., et al. (2014). CR cistrome: a ChIP-seq database for chromatin regulators and histone modification linkages in human and mouse. Nucleic Acids Res. 42 (Database issue), D450–D458. doi:10.1093/nar/gkt1151

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, L., Siegenthaler, J. A., Dowell, R. D., and Yi, R. (2016). Foxc1 reinforces quiescence in self-renewing hair follicle stem cells. Science 351 (6273), 613–617. doi:10.1126/science.aad5440

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Y., Wang, R., Zhang, S., Song, S., Jiang, C., Han, G., et al. (2019). iTALK: an R Package to characterize and illustrate intercellular communication. bioRxiv. doi:10.1101/507871

CrossRef Full Text | Google Scholar

Worch, J., Ranft, A., DuBois, S. G., Paulussen, M., Juergens, H., and Dirksen, U. (2018). Age dependency of primary tumor sites and metastases in patients with Ewing sarcoma. Pediatr. Blood Cancer 65 (9), e27251. doi:10.1002/pbc.27251

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, G., Wang, L. G., and He, Q. Y. (2015). ChIPseeker: an R/bioconductor package for ChIP peak annotation, comparison and visualization. Bioinformatics 31 (14), 2382–2383. doi:10.1093/bioinformatics/btv145

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Z., Lee, J. H., Ruan, H., Ye, Y., Krakowiak, J., Hu, Q., et al. (2019). Transcriptional landscape and clinical utility of enhancer RNAs for eRNA-targeted therapy in cancer. Nat. Commun. 10 (1), 4562. doi:10.1038/s41467-019-12543-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, P., Lu, Y., and Liu, L. (2015). Correlation of decreased expression of PHLDA1 protein with malignant phenotype of gastric adenocarcinoma. Int. J. Clin. Exp. Pathol. 8 (5), 5230–5235.

PubMed Abstract | Google Scholar

Zhao, J., Zhao, Y., Wang, L., Zhang, J., Karnes, R. J., Kohli, M., et al. (2016). Alterations of androgen receptor-regulated enhancer RNAs (eRNAs) contribute to enzalutamide resistance in castration-resistant prostate cancer. Oncotarget 7 (25), 38551–38565. doi:10.18632/oncotarget.9535

PubMed Abstract | CrossRef Full Text | Google Scholar

Zheng, R., Wan, C., Mei, S., Qin, Q., Wu, Q., Sun, H., et al. (2019). Cistrome data browser: expanded datasets and new tools for gene regulatory analysis. Nucleic Acids Res. 47 (D1), D729–D735. doi:10.1093/nar/gky1094

PubMed Abstract | CrossRef Full Text | Google Scholar

Zlotnik, A., Burkhardt, A. M., and Homey, B. (2011). Homeostatic chemokine receptors and organ-specific metastasis. Nat. Rev. Immunol. 11 (9), 597–606. doi:10.1038/nri3049

PubMed Abstract | CrossRef Full Text | Google Scholar