Purpose:

Sarcoma is the second most common solid tumor type in children and adolescents. The high level of tumor heterogeneity as well as aggressive behavior of sarcomas brings serious difficulties to developing effective therapeutic strategies for clinical application. Therefore, it is of great importance to identify accurate biomarkers for early detection and prognostic prediction of sarcomas.

Experimental Design:

In this study, we characterized three subtypes of sarcomas based on tumor immune infiltration levels (TIIL), and constructed a prognosis-related competing endogenous RNA (ceRNA) network to investigate molecular regulations in the sarcoma tumor microenvironment (TME). We further built a subnetwork consisting of mRNAs and lncRNAs that are targets of key miRNAs and strongly correlated with each other in the ceRNA network. After validation using public data and experiments in vivo and in vitro, we deeply dug the biological role of the miRNAs and lncRNAs in a subnetwork and their impact on TME.

Results:

Altogether, 5 miRNAs (hsa-mir-125b-2, hsa-mir-135a-1, hsa-mir92a-2, hsa-mir-181a-2, and hsa-mir-214), 3 lncRNAs (LINC00641, LINC01146, and LINC00892), and 10 mRNAs (AGO2, CXCL10, CD86, CASP1, IKZF1, CD27, CD247, CD69, CCR2, and CSF2RB) in the subnetwork were identified as vital regulators to shape the TME. On the basis of the systematic network, we identified that trichostatin A, a pan-HDAC inhibitor, could potentially regulate the TME of sarcoma, thereby inhibiting the tumor growth.

Conclusions:

Our study identifies a ceRNA network as a promising biomarker for sarcoma. This system provides a more comprehensive understanding and a novel perspective of how ceRNAs are involved in shaping sarcoma TME.

Translational Relevance

We characterized and verified three subtypes of sarcoma according to the levels of tumor infiltrating immune cells, which are of critical importance in shaping the tumor microenvironment (TME) to promote the development of tumor and affect the cancer therapy. Furthermore, we constructed a prognosis-related ceRNA network containing the key miRNAs and their strongly correlated target mRNAs and lncRNAs based on the subtypes identified. Finally, we identified a sarcoma-targeting drug, a pan-HDCA inhibitor trichostatin A, which can regulate the TME of sarcoma, thereby inhibiting tumor growth. Our results suggest that the prognosis-related ceRNA could serve as a potential biomarker and therapeutic target of sarcoma in clinical situations.

Sarcoma is a heterogeneous group of mesenchymal neoplasms, which can arise from virtually any anatomic site. This gives sarcoma intricate histological classifications. Typically, soft tissue and primary bone sarcomas are the two main clinical types of sarcoma (1). Even though sarcomas are relatively rare in adults, it accounts for 12%–15% of overall pediatric tumor cases in Europe (2). In instances with high histology grade and tumor lager than 5 cm2, the recurrence risk is higher than 50% after surgery. For these reasons, it is typically recommended that patients with high-risk localized tumors receive a combination of surgery, radiotherapy, and chemotherapy (3). But this strategy modestly affects the average survival of metastatic sarcoma that is 12–18 months (4). Chemotherapy was usually applied before surgery to reduce the tumor size in sarcomas that is called neoadjuvant therapy. The typical chemotherapy drugs for soft tissue sarcomas are doxorubicin and ifosfamide. However, the beneficial effect of combination of adjuvant chemotherapy and surgery to high-risks soft tissue sarcoma are still controversial. For example, one study from UK did not find significant improvement in overall survival (OS) for resected STSs with neoadjuvant chemotherapy (5). Targeted drugs are also commonly used to control the development of STS. Presently, drugs targeting VEGF, platelet-derived growth factor, mTOR, IGF-1R, CDK4, c-Kit, and MET are either approved by FDA or undergoing clinical trials for STS treatment. However, the drug resistance and target missing in some types of STS limit the therapeutic effect of targeted drugs (6). In addition, although the application of immunotherapy has received massive of benign outcome in cancer treatment, the therapeutic effect of immune checkpoint inhibitors is not ideal for treating sarcoma. Ipilimumab, an anti-CTLA4 mAb, can neither significantly improve clinical effect nor evoke serological antigen neutralization responses in advanced or metastatic synovial sarcoma even with high-expression level of CTLA-4 (7). The assessment for safety and activity of anti–PD-1 antibodies on patients with advanced STS and bone sarcoma showed that the primary endpoint of overall response was not ideal for either cohort (8). Hence, exploring new therapeutic targets is necessary and urgent for enriching therapeutic approach for sarcoma. It is also important to further understand the underlying regulatory mechanisms in sarcoma immune response.

The difficulty of sarcoma treatment indicates that the interaction of the tumor components in sarcoma tissue is very intricate. Especially the failure of immunotherapy suggests that there might be some other elements regulating the interaction between tumor cells and immune cells, but not simply immune check point. Therefore, the determining factors of immunotherapy may hide in the complicate tumor microenvironment (TME) of sarcoma. Multiple factors, for example, the inherent antigenicity of the tumor, “neoantigens” and other protein mutations in the tumor cells, mutation burden of tumors, and infiltration of immune effector cells in the tumor site, contribute to the response of patients with sarcoma to immunotherapy. Among these factors, the specific pattern of tumor-infiltrating lymphocytes in the TME is closely related to better outcomes of patients with cancer, regardless of the therapy regimen. TME itself is an independent factor to participate in tumor immune response. Explicit evidence has confirmed that recruitment, activation, and reprogramming of immune and stromal cells are regulated by the TME on tumor sites (7). The immunosuppressive TME further promote the immune escape of tumor cells (3). Meanwhile, the regulatory effect of host immune system to tumor is also achieved through TME. The host immune system can alter the components of the TME to affect the cancer development and progression (7). On the basis of these critical roles of TME in tumorigenesis and immune response, researchers highlighted the prognostic value of assessing the immune features of TME. Especially, the immune histopathological and molecular biomarkers of TME can be used to evaluate the patients’ treatment response. Zhu and Hou (8) described that the macrophages are the largest population that infiltrated in the TME of sarcoma among the 22 immune cell types. In addition, a negative relationship was found among CD8+ T cells and the M0/M2 macrophages, but there was a positive relationship between the CD8+ T cells and the M1 macrophages in sarcomas. In addition, the specific types of macrophages or mast cells in TME promote the growth of the tumor by producing growth factors (9) and maintain chronic inflammation (10). However, research on components of TME in sarcomas is still relatively scant. Therefore, the study in respect to the TME components and immune system biomarkers of sarcoma may be conducive to predict the prognosis and therapeutic responses.

The identification of competing endogenous RNA (ceRNA) highlighted that various types of RNA participate in regulating gene expression at the posttranscriptional level harmoniously. mRNAs, the transcripts of pseudogenes, long non-coding RNAs (lncRNAs), and circRNAs, can mediate the translation and stability of a target gene via competitively binding to the miRNA. The malfunction of ceRNA networks contribute to many cancer processes, such as epithelial-to-mesenchymal transition, metastasis, immune infiltration (11), and so on. As expected, dysregulation of ceRNA networks was commonly observed in sarcoma. The characteristic of ceRNA networks in cancer area is distinct from the normal tissue, thus the feature of ceRNA networks can be used to distinguish tumors from normal tissues, and mesenchymal stem cells (12). ceRNA networks can even be used to identify recurrent tumors from primary tumor tissues (13), and to characterize different tumor subtypes no matter base on the molecular component (14), or TME assessment (8).

In this study, we downloaded the public transcriptomic data and corresponding clinical information of sarcomas from The Cancer Genome Atlas (TCGA). Forty-five immune-related signatures were collected to estimate the immune infiltration and classify subtypes of sarcomas. We grouped sarcoma cases into 3 subtypes according to their immune cell infiltration levels (TIILs, tumor immune infiltration levels). Subsequently, we constructed a dysregulated ceRNA network by comparing the sarcoma subtypes we defined and validated the expression of genes in vivo and in vitro. Eventually, we elucidated the roles of the ceRNA network in TME and successfully predicted that a drug could inhibit the tumor growth of sarcoma. This study will construct a systematic view of TME-related signatures in sarcoma and facilitate the development of therapeutic drug targeting TME for sarcoma.

Data downloading and pre-processing

The mRNA, lncRNA, and miRNA expression profiles and related clinical information of human sarcoma cohort data were downloaded from TCGA database (RRID: SCR_003193), fetched on May 28, 2021. Totally, RNA-sequencing (RNA-seq) data of 259 patients with sarcoma (normalized FPKM data) were collected, and 257 of them had miRNA sequencing (normalized RPM data). The RNA-seq and miRNA sequence data were derived from Illumina HiSeq and miRNASeq platform, respectively. The comprehensive gene annotation file was downloaded from GENECODE (www.gencodegenes.org) to identify the protein-coding RNA and lncRNA from the total RNA expression profile. The R package “miRBaseConverter” (RRID: SCR_023873) was used to convert mature miRNAs to hairpin miRNAs.

Identification of immune cell infiltration subtypes

To estimate the infiltration levels of different types of immune cells, immunogenomic pathways and the activity of immune-related biological processes in patients with sarcoma, a total of 45 immune signatures were collected, including 28 immune cell types (15) and 17 immune-related pathways and biological functions (16). The infiltration levels of 45 immune signatures were quantified using single sample gene set enrichment analysis (ssGSEA, RRID: SCR_003199) implemented in the R package “GSVA” (version 1.24.0, RRID: SCR_021058). According to the results of ssGSEA analysis, 259 patients with sarcoma in the TCGA sarcoma cohort were grouped into three clusters and we defined as high, low, and moderate groups of immune infiltration using hierarchical clustering analysis. The proportion of 64 immune and stromal cell types among the subtypes was estimated using the function “deconvolute_xcell” in the R package “immunedeconv” (RRID: SCR_023869). Stromal score, immune score, ESTIMATE score, and tumor purity were calculated using the algorithm “CIBERSORT” (RRID: SCR_016955), which is based on 22 stromal and immune cell types, to evaluate immune cell infiltration among sarcoma subgroups.

Differential expression and functional enrichment analyses

The R package “limma” (3.46.0, RRID: SCR_010943) was applied for the analysis of identifying differential expression genes with criteria of adjusted P value of <0.05 and |fold change| >2 among different sarcoma subtypes. The R package “clusterProfiler” (3.18.1, RRID: SCR_016884) was used to perform gene ontology (GO) enrichment and Kyoto encyclopedia of genes and genomes (KEGG) analysis. The FDR with BH adjusted <0.05 was considered as significant enrichment.

Survival analysis and ceRNA network construction

First, the interactions between miRNAs and mRNAs were acquired from three commonly used databases: miRDB (RRID: SCR_010848), TargetScan (RRID: SCR_010845), and miRWalk (RRID: SCR_016509). Moreover, the mRNA–miRNA interaction pairs shared among three databases were selected for further construction. The interactions between miRNAs and lncRNAs were collected from miRCode (RRID: SCR_023870). Subsequently, the identified differentially expressed (DE) miRNAs, DE mRNAs, and DE lncRNAs, between high and low immune infiltration subgroups, were filtered out to construct the global ceRNA network. The Cytoscape (3.8.2, RRID: SCR_003032) was used for visualizing the ceRNA network. Finally, to build a prognosis-related ceRNA network, we used the R package “survival” (RRID: SCR_021137) to perform Kaplan-Meier (KM) curve analysis of DE miRNAs in the global ceRNA network. The results with a P value of <0.05 were retained as key miRNAs for the construction of prognosis-related ceRNA network. Besides, univariate and multivariate Cox regression analyses were conducted to identify the independent prognostic indicators for OS.

Construction of PPI network

PPI network was built using the Search Tool for the Retrieval of Interacting Genes (STRING; RRID: SCR_005223) to investigate the interaction of DE mRNAs in the prognosis-related ceRNA network. The protein interactions were filtered by the median confidence score (0.4) and visualized by Cytoscape (3.8.2, RRID: SCR_003032). The hub genes in the PPI network were identified by the plugin cytoHubba (RRID: SCR_017677) according to the degree of the connectivity and top 30 hub genes were identified through the degree ranking method. The PPI network were analyzed using MCODE (RRID: SCR_015828) plugin and the top module was selected for further analysis.

Correlation analysis of lncRNA and mRNA

In the prognosis-related ceRNA network, the targets lncRNAs and mRNA of individual key miRNA were selected to perform Pearson correlation analysis. Pearson Correlation Coefficient (PCC) >0.4 and P < 0.05 were used to yield the interaction pairs with high confidence, and the candidate pairs were used to construct a ceRNA subnetwork to investigate the potential mechanisms.

Validation of gene expression

The online database, Expression Atlas (RRID: SCR_007989), which supports the gene and protein expression across species and biological conditions, was used to screen the candidate datasets to validate the expression of genes in the ceRNA subnetwork. Finally, GSE16779, GSE31019, GSE35493, GSE66354, GSE21050, GSE41293, and GSE43045 from Gene Expression Omnibus (GEO, RRID:SCR_005012) and PRJEB36314 from European Bioinformatics Institute (RRID:SCR_004856) databases were selected to further verification. In addition, a rat model of rhabdomyosarcoma induced by military heavy metals was also used to further validation. Briefly, one of eight pure metal pellets was surgically implanted individually into the gastrocnemius muscle of 3-month-old male Sprague Dawley rats (n = 8/group). Tantalum (Ta), which is inert both in vivo and in vitro, with biocompatibility in devices, lants, and stents, was used as a control. RNA-seq was performed according to the standard procedures. The data were deposited in GEO under the accession number of GSE157929. According to the results, the rats implanted cobalt (Co) and nickel (Ni) were identified rhabdomyosarcoma, therefore, we just collected the groups of Ta, Co, and Ni implant for further validation.

Gene set enrichment analysis

To investigate the biological pathways of miRNAs and lncRNAs in the subnetwork, GSEA was conducted on the basis of TCGA the sarcoma data. According to the cutoff value, the median value of the expression, the samples were divided into low- and high-expression groups. The permutations of the gene set were performed 1,000 times per analysis and gene sets with |normalized enrichment score | > 1 and P < 0.05 were considered as significant enrichment.

Evaluation of immune cell infiltration

The algorithm “CIBERSORT” (RRID: SCR_016955) was used to estimate the relative abundance of 22 immune cell types. The PCCs among the relative abundance of 22 immune cell types was calculated to identify the patterns of the immune cell infiltration. PCC was also calculated to evaluate the correlation between the expression of lncRNAs and miRNAs, and immune cell infiltration. Besides, sarcoma samples from TCGA databases were divided into low- and high-expression groups, then the Wilcoxon rank-sum test was performed to compare the differences in immune cell infiltration between the low- and high-expression groups.

Identification of potential drug

Connectivity Map (CMap) analysis, which uses a reference database containing drug-specific gene expression profiles and compares it with a disease-specific gene signatures, was adopted to identify the potential drugs. In this study, the R package “Dr. Insight” (RRID: SCR_023871) was applied to performing CMap analysis. The results of the t test statistic scores between low and high immune infiltration subtypes from the “limma” analysis were used as input data. The drugs with an FDR <0.1 were considered as key targets for the therapy of patients with sarcoma. Moreover, “oncoPredict” (RRID: SCR_023872) R package was used to predict the sensitivity of the drug.

Methylation analysis

The Human Disease Methylation Database Version 3.0 RRID: SCR_017488, which is an online website and supports the abnormal methylome of human diseases, was used to analyze the methylation of genes in the ceRNA subnetwork. The t test with P value of <0.05 was considered as significant difference methylation between healthy people and patients with sarcoma.

Mouse model construction

ICR mice were purchased from Silaike Experiment Animal Co., Ltd. All animal experiment protocols were ethically approved by University of Macau Animal Experiment Ethics Committee (UMARE–041–2020). A total of 1×106 of S180 cells were planted in 4-week-old male ICR mice by subcutaneous injection to construct sarcoma mouse models. Supplements of TSA were start when the size of tumor reached 0.104 cm3. 18 ICR mice were divided into 3 groups: The control, TSA high-dose, and TSA low-dose groups. Mice were intraperitoneally receiving 0.5 mg/kg of TSA and 2 mg/kg of TSA every day for 2 weeks in low-dose and high-dose group, respectively. Mice in the control group only received equal volumes of DMSO. The mice were sacrificed, and the tumor tissues were harvested for following experiment.

3D microsphere cell culture

The human liposarcoma cell line SW872, a gift from Dr. Henry Hang Fai Kwok (University of Macau), was cultured in DMEM (Gibco Life Technologies, 11965–092) containing 10% FBS (Cell Bio, FSP500) with 100 IU/mL penicillin, and 100 μg/mL streptomycin (PS; Gibco, 15140–122). For three-dimensional (3D) culture, the SW872 cells were dissociated using Trypsin-EDTA (Gibco, 25200–056) and seeded in the low-attachment plate (Costar, 7007) at a total of 1×105 per well culturing of 48 hours to allow aggregation to dense spheres, then formed spheres were resuspended in Matrigel solution (50% DMEM medium, and 50% Matrigel matrix, Corning 354234), and seeded in pre-warmed 24-well culture plates at 30 μL per drop. Once cell-Matrigel drops were solidified at 37°C, 500 μL/well SW872 culture medium was added in to initiate 3D culture (17).

TSA treatment and PBMC co-culture

The 3D cultures were divided into seven groups, high-level peripheral blood mononuclear cell (PBMC high), low-level PBMC (PBMC low), 1 or 0.2 μmol/L TSA-treated cells co-cultured with high-level or low-level PBMC, and a control group without any treatment. Each group contained 5-well plates. TSA treatment groups were pretreated with RPMI-1640 medium (Gibco Life Technologies, A4192301) containing Trichostatin A (TSA) 1 or 0.2 μmol/L, 10% FBS, and 1%PS for 2 days incubation at 37°C before adding PBMCs. The stock TSA solution (1 mmol/L) was prepared with pure DMSO following the manufacturer's instructions. Meanwhile, the control and PBMC-only groups were cultured with the same medium but change the TSA solution to 1:1,000 pure DMSO. The PBMCs were isolated from the whole blood provided by donors as previously described. The PBMC high groups received a total of 5×105 PBMC cells per well, and low-level groups received a total of 1×105 PBMC cells per well. PBMCs and the SW872 sphere were cocultured in RPMI-1640 medium containing 100 IU/ml of IL2, 10% FBS, and 1%PS, at 37°C. After 2 days of incubation, the supernatant was removed, and the cells were washed with PBS. 5-well plates in the same group were dissolved in TRizol solution (Invitrogen Life Technologies) and mixed in one tube. The mixtures were stored at −80°C for future experiments.

ELISA

ELISA assays were performed with commercial kits according to the individual manufacturer's instructions. The following ELISA kit was used: The Mouse HDAC ELISA kit from ML BIO (ml037192). Tumor tissues were homogenized with PBS then centrifuged for 3,000 rpm, 20 minutes. 10 μL of supernatant were mixed with working buffer in 98-well plate. Incubate the plate in 37°C, 15 minutes, following by adding 50 μL stop buffer. The plate was read at 450 nm by spectrophotometry.

RNA purification and RT-qPCR

Total RNA was extracted by TRizol-chloroform method as previously described. RNA was reversely transcribed using the iScript cDNA Synthesis Kit (1708890, Bio-Rad). The primers used for 3D culture samples are shown in Supplementary Table S1. The primers used for mice tissue were listed in Supplementary Table S2. GAPDH gene was set as internal reference. All PCR reactions were performed with iTaq Universal SYBR green supermix (1725124, Bio-Rad) in a total volume of 10 μL, with the following condition: Activated polymerase at 95°C for 30 seconds, then 42 cycles of denaturation at 95°C for 5 seconds, annealing at 60°C for 45 seconds, and extension at 72°C for 30 seconds. Data were analyzed using GraphPad Prism 6.

Flow cytometry

Mice tumor tissues were homogenized and filtrated by 200-mesh strainers. Then, isolated cells were collected by centrifugation at 800 rpm for 10 minutes and resuspended in PBS. A total of 1×106 cells in 400-μL PBS were incubated with antibodies at 4°C for 30 minutes avoiding light. The following antibodies were used: Anti-mouse CD4 antibody (BioLegend 100509 0.25 μg/test), anti-mouse CD8 (BioLegend 100707 0.25 μg/test), anti-mouse CD3 (BioLegend 100235 0.5 μg/test), anti-mouse CD206 (BioLegend 141705 0.5 μg/test), anti-mouse CD86 (BioLegend 105011 0.25 μg/test), anti-mouse F4/80 (BioLegend 123107 0.25 μg/test) to examinate the percentage of T cells. Flow cytometry was performed using Cytoflex 5 (Beckman) and analyzed using the Flowjo software (BD Biosciences)

IHC

IHC staining was performed as per standard protocol. Briefly, the formalin-fixed, paraffin-embedded mice tumor tissues were cut into 5-μm-thick sections. The sections were deparaffinized in xylene and rehydrated in a graded series of alcohols. The sections were washed with PBS and then boiled in R-Buffer-A (Electron Microscopy Sciences). Then, sections were washed with PBS followed by incubation with quenching solution (10 mL of 30% H2O2 to 90 mL of absolute methanol) for 20 minutes. Sections were incubated with blocking solution (50% 3% BSA and 50% Animal-Free Blocker, Vector Laboratories) for at 1 hour at room temperature (RT) and incubated with primary antibody against CD3, CD4, or CD8 at 4°C overnight. Sections were then incubated with Biotinylated goat anti-mouse and rabbit (Abcam) for 10 minutes at RT after washing with PBS for 3 times. Then Streptavidin peroxidase was added to incubate for 10 minutes at RT. DAB solution (DAB substrate kit peroxidase, Vector Laboratories) was prepared before use (1 drop reagent 1, 2 drops reagent 2, and 1 drop reagent 3 in 2.5 mL ddH2O) and added to the sections. Once the color was changed, sections were washed with tap water immediately for 5 minutes and counterstained in hematoxylin for 1 minute and washed with tap water for 5 minutes. Slides then were dehydrated in a graded series of alcohols and put into xylene for 5 minutes. Slides were placed to dry and DPX mounting was covered onto the tissues and a coverslip was placed on each slide. Images were taken using Carl Zeiss LSM 880 super-resolution microscopes. A 4-point scale was used to define the different positive rate of CD3+, CD4+, and CD8+ in tumor tissue of mice bearing with S180 according to the IHC results. +, ++, +++, and ++++ indicated positive rate located in <30%, 30%∼40%, 40%∼50%, >50%, respectively.

Statistical analysis

All statistical analyses were conducted using R software (version 4.2.0, RRID: SCR_001905). The R packages “ggplot2” (version 2.2.1, RRID: SCR_014601) and “pheatmap” (1.0.12, RRID: SCR_016418) were applied to plotting heat maps and other visualized graphs. The univariate and multivariate Cox regression analyses were used to identify the prognostic variables that were significantly correlated with OS and evaluate the hazard ratio (HR), 95% confident interval (CI) of HR and a P value. The OS curves were established by a KM survival curve function “ggsurvplot” (RRID: SCR_021094), as implemented in R package “survminer” (version 0.4.2, RRID: SCR_021094), and the difference in survival distributions between sarcoma subgroups was estimated by the two-sided log-rank test. For comparisons between two groups, Wilcoxon analysis was performed, whereas for comparisons among three or more groups, Kruskal–Wallis analysis was applied.

Data availability

The authors declare that all data supporting the findings of this study are available from the corresponding authors on reasonable request. The source code for main analysis is deposited at “https://github.com/lva85/Prognostic_ceRNA_network_in_sarcoma.git.”

Identification, validation, and characterization of sarcoma subtypes

To assess the TME in patients with sarcoma, we collected the public data of sarcomas from the TCGA database, including clinical information, RNA-seq, and miRNA-seq data. As the pipeline (Fig. 1) shows, first, we performed the TIILs analysis. The immune metagenes for 28 immune cell types are determined using 37 microarray data from a number of different studies (15). The 17 immune-related pathways and biological functions are generated by constructing network of a given system and mining networks with graph theory (18). On the basis of the infiltration levels of 45 immune-related gene signatures for each sample, we conducted hierarchical clustering and identified three apparent clusters (i.e., C1, C2, and C3 in Supplementary Fig. S1A). The rationality and reliability of the clustering was further confirmed by investigating the expression of the immune-related cluster of differentiation (CD) and IL molecules. We found that lots of the CD and IL family genes showed significantly different expression among the three clusters (Supplementary Fig. S1B and S1C). In addition, we estimated the proportion of 64 immune and stromal cell types. We found that the proportion of some myeloid cell types is distinct among the three clusters, including monocytes, neutrophils, macrophages, macrophage M1, macrophage M2, dendritic cells (DC), activated DCs, conventional DCs, plasmacytoid DCs, immature DCs, and basophils. Meanwhile, lymphoid cells also showed significant difference among three subtypes, such as CD4+ memory T cells, CD4+ naïve T cells, CD4+ T cells, CD4+ Tcm (central memory) cells, CD4+ Tem (effector memory) cells, CD8+ T cells, CD8+ naïve T cells, CD8+ Tcm cells, CD8+ Tem cells, gamma delta T cells (Tgd cells), natural killer cells (NK cells), NK T cells, B cells, class-switched memory B cells, and plasma cells (Supplementary Fig. S1D). We further found that there were 43 signatures that showed significant differences between C1 and C2, 29 between C1 and C3, and 23 between C2 and C3 (Supplementary Fig. S1E). The comparison between C1 and C2 covered all different signatures that also appeared in the comparison between C1 versus c3 and C2 versus C3. Finally, after removing redundancy, we defined three subtypes of sarcoma: Low TIILs (Low, represented by C1), high TIILs (High, represented by C2), and moderate TIILs (Moderate, represented by C3), according to the infiltration levels of the 43 signatures (Fig. 2A).

Figure 1.

Diagram depicting the pipeline of the study design. In total, there are seven steps for the analysis: 1. Identification of subtypes; 2. Validation of the subtypes; 3. Characterizations of the subtypes; 4. Construction of the ceRNA network; 5. Validation of the expression of the genes in the subnetwork; 6. Identification of the ability of the subnetwork to reshape the TME; and 7. Prediction of drugs.

Figure 1.

Diagram depicting the pipeline of the study design. In total, there are seven steps for the analysis: 1. Identification of subtypes; 2. Validation of the subtypes; 3. Characterizations of the subtypes; 4. Construction of the ceRNA network; 5. Validation of the expression of the genes in the subnetwork; 6. Identification of the ability of the subnetwork to reshape the TME; and 7. Prediction of drugs.

Close modal
Figure 2.

Identification, validation, and characterization of subtypes of sarcoma. A, Definition of sarcoma subtypes: C1 represents low TIILs, C2 represents high TIILs, and C3 represents moderate TIILs. B, Clustering confirmation using Consensus Clustering Analysis. C1, C2, and C3 represent sarcoma subtypes with low, high, and moderate TIILs, respectively. C, Comparison for the clustering between methods of Hierarchical Clustering (HC) and Consensus Clustering (CC). D, Schemes of experiments with tumor-bearing mice. E, Flow cytometry analysis of CD4+ T cells for tumor tissues of S180-bearing mice. F, Flow cytometry analysis of CD8+ T cells for tumor tissues of S180-bearing mice. G, Comparison of tumor purity among subtypes. H, Comparison of tumor microenvironment score among subtypes. I, Overall survival analysis for the subtypes. J, Distribution of the niduses of patients with sarcoma. The P values were calculated using the Wilcoxon rank-sum test; *, P < 0.05; ****, P < 0.0001.

Figure 2.

Identification, validation, and characterization of subtypes of sarcoma. A, Definition of sarcoma subtypes: C1 represents low TIILs, C2 represents high TIILs, and C3 represents moderate TIILs. B, Clustering confirmation using Consensus Clustering Analysis. C1, C2, and C3 represent sarcoma subtypes with low, high, and moderate TIILs, respectively. C, Comparison for the clustering between methods of Hierarchical Clustering (HC) and Consensus Clustering (CC). D, Schemes of experiments with tumor-bearing mice. E, Flow cytometry analysis of CD4+ T cells for tumor tissues of S180-bearing mice. F, Flow cytometry analysis of CD8+ T cells for tumor tissues of S180-bearing mice. G, Comparison of tumor purity among subtypes. H, Comparison of tumor microenvironment score among subtypes. I, Overall survival analysis for the subtypes. J, Distribution of the niduses of patients with sarcoma. The P values were calculated using the Wilcoxon rank-sum test; *, P < 0.05; ****, P < 0.0001.

Close modal

To validate the sarcoma subtypes, we reanalyzed the data using Consensus Clustering approach (RRID: SCR_016954). From the consensus cumulative distribution function plot (Supplementary Fig. S2A), delta area plot (Supplementary Fig. S2B), and cluster-consensus value plot (Supplementary Fig. S2C), we found when the number of clusters was 3 (Fig. 2B), the clustering showed high stability. The heat map visualization for the infiltration levels of immune-related signatures is also obviously distinct among these three clusters and presents low (C1), high (C2), and moderate (C3) TIILs pattern (Supplementary Fig. S2D). The comparison for the clustering results between these two approaches also showed high consistence, the overlap of the samples highly reaches to over 80% for every corresponding subtype (Fig. 2C). Furthermore, we developed an S180-bearing mouse sarcoma model to validate the sarcoma subtypes (Fig. 2D). On the basis of the flow cytometry analysis of CD4+ and CD8+ T cells from a cohort of 10 mice, we found that three mice could be classified into low T cells infiltration, three mice infiltrated with high levels of T cells, and remaining four mice could be grouped into moderate T cells infiltration (Fig. 2E and F).

We further investigated the characteristics of each subtype. We found that stromal score, immune score, ESTIMATE score (Supplementary Fig. S3A), tumor purity (Fig. 2G), and TME score (Fig. 2H) were also stratified among three sarcoma subtypes and the scores were consistent with the TIILs of subtypes. We then compared the 10-year OS among the subtypes and the results came out with statistical differences (log-rank P = 0.035). Importantly, the OS rates increased as the TIILs increased (Fig. 2I). The results of differential transcriptomic expression analysis demonstrated that the most DE transcriptome was identified in the comparison between C2 (High) and C1 (Low). Notably, 1,442 genes were downregulated from C1 (Low) to C2 (High), and only 15 genes were downregulated from C1 (Low) to C3 (Moderate). Moreover, almost all of the upregulated genes from C3 (Moderate) to C2 (High; 378 out of 380) were also upregulated from C1 (Low) to C2 (High; Supplementary Fig. S3B and S3C). The detailed DE analysis results are organized in Supplementary Tables S3–S5. A lot of immune-related GO terms and pathways were enriched using the DE genes between C1 (Low) and C2 (High) subtypes (Supplementary Tables S6 and S7). Besides, we found that the majority of niduses of most patients were derived from the lower extremity and retroperitoneum/upper abdominal. However, the subtypes we defined almost appear in all different nidus sites, implying the common characteristic of soft tissue sarcoma (Fig. 2J).

Construction of prognosis-related ceRNA network

According to the characteristics of sarcoma subtypes, especially the moderate TIILs of the C3 (Moderate) subtype, which appears to be close to C1 (Low), therefore, we focused on the C1 (Low) and C2 (High). We annotated the biotypes of DE transcriptome between C1(Low) and C2(High) using the reference annotation file from GENCODE database and identified 127 lncRNAs and 1,645 mRNAs (Supplementary Fig. S4A; Supplementary Table S8). Next, we identified 171 DE pre-miRNAs between C1(Low) and C2(High) using the limma package with the criteria of BH adjusted P value < 0.05 (Supplementary Fig. S4B; Supplementary Table S9). 24 of DE pre-miRNAs were significantly related to patients’ prognosis according to the KM survival curves (Supplementary Table S10). Subsequently, 94 out of 127 DE lncRNAs were filtered to simulate the interaction with 55 out of 171 DE pre-miRNAs based on the miRcode database. Then we further identified 54 out of 55 pre-miRNAs targeted 993 out of 1,645 DE mRNAs from miRDB, TargetScan, miRWalk databases. Consequently, we constructed a global ceRNA interaction network, including 5,407 edges, 993 mRNAs, 94 lncRNAs, and 54 pre-miRNAs (Supplementary Fig. S4C).

The KM survival curves were applied to assessing the association of 54 pre-miRNAs in the initial ceRNA network with patient prognosis. The results indicated that 9 pre-miRNAs were significantly negatively correlated with OS, that is, hsa-mir-125b-2, hsa-mir-130a, hsa-mir-135a-1, has-mir-181a-2, has-mir-214, has-mir-301a, has-mir-9–2, has-mir-92a-1, and has-mir-92a-2 (Fig. 3A). These 9 pre-miRNAs were considered as key miRNAs. Ultimately, we built a prognosis-related ceRNA network with nine key miRNAs targeting 74 DE lncRNAs, 493 DE mRNAs and containing 576 nodes and 957 edges, totally (Fig. 3B). The whole-miRNA recognition elements between miRNAs and lncRNAs, miRNAs and mRNAs were deposited in Supplementary Tables S11 and S12, respectively. All the nodes in the prognosis-related ceRNA network were included in the univariate Cox regression analysis. 6 miRNAs, 11 lncRNAs, and 83 mRNAs with a P value of <0.05 were subjected to multivariate Cox regression analysis. Finally, only 1 miRNA and 15 mRNA were identified as independent prognostic factors for OS in patients with sarcoma (Supplementary Table S13).

Figure 3.

Construction of prognosis-related ceRNA network. A, The KM curves of the dysregulated miRNAs with significant association with OS in the global ceRNA network. B, Prognosis-related ceRNA network construction. C, Bubble plot for the GO and KEGG enrichment analysis using the genes in prognosis-related ceRNA network. D, The subnetwork construction using the pairs with strong co-expression between mRNAs and lncRNAs and targeting miRNAs. E, Heat map for the correlation between nodes in the subnetwork and clinic features.

Figure 3.

Construction of prognosis-related ceRNA network. A, The KM curves of the dysregulated miRNAs with significant association with OS in the global ceRNA network. B, Prognosis-related ceRNA network construction. C, Bubble plot for the GO and KEGG enrichment analysis using the genes in prognosis-related ceRNA network. D, The subnetwork construction using the pairs with strong co-expression between mRNAs and lncRNAs and targeting miRNAs. E, Heat map for the correlation between nodes in the subnetwork and clinic features.

Close modal

To investigate the biological processes and pathways that the prognosis-related ceRNA network involved, we performed GO and KEGG enrichment analysis of 493 mRNAs with R package “clusterProfiler.” GO functional enrichment analysis notably indicated that the biological function of ceRNA network was enriched under immune response groups, including “immune response-regulating signaling pathway,” “lymphocyte differentiation,” “mononuclear cell differentiation,” “positive regulation of IL2 production,” and “positive regulation of leukocyte cell–cell adhesion.” Surprisingly, besides immune response, the ceRNA network was also dominant in epigenetically related functions, such as “chromatin organization,” “histone modification,” and “protein deacylation or deacetylation.” KEGG pathway analysis strongly demonstrated that ceRNA network was strongly related to cancer-related pathways, such as the “TNF signaling pathway” and “transcriptional misregulation in cancer.” The significant enrichment of “Herpes simplex virus 1 infection” and “Pertussis” in patients with sarcoma are similar to molecular response induced by virus infection (Fig. 3C).

To further depict the biological functions of key miRNAs in the prognosis-related ceRNA network, we used GSEA analysis to predict the enriched KEGG pathways between key miRNAs’ high- and low-expression groups. The high expression of hsa-mir-135a-1 and hsa-mir-181a-2 showed a similar positive correlation to axon guidance and regulation of the actin cytoskeleton. Besides, it was also found that hsa-mir-135a-1 was involved in the WNT signaling pathway and the leukocyte transendothelial migration pathway. hsa-mir-181a-2 participated in the pathways of spliceosome, NK cell-mediated cytotoxicity and Fc gamma R-mediated phagocytosis. The high expression of hsa-mir-214 was positively correlated with the T-cell receptor signaling pathway but not very statistically convinced. The high expression of hsa-mir-92a-2 was positively related to NK cell-mediated cytotoxicity and primary immunodeficiency but negatively related to autoimmune thyroid disease, complement and coagulation cascades, immune network for IgA production, allograft rejection, leishmania infection, type Ⅰ diabetes mellitus, viral myocarditis, and systemic lupus erythematosus. In addition, the high expression of hsa-mir-125b-2 was positively related to apoptosis, neurotrophin signaling pathway, p53 signaling pathway, and toll-like receptor signaling (Supplementary Fig. S5A).

To access the cellular functions of the prognosis-related ceRNA network, we further constructed a PPI network using the mRNAs from the network based on the STRING database. The PPI network was visualized by Cytoscape and contained 394 nodes and 1,306 edges (Supplementary Fig. S5B). The top 30 hub genes with the highest degree were identified by cytoHubba plugin (Supplementary Fig. S5C; Supplementary Table S14). Moreover, the top module containing 18 nodes and 121 edges (Supplementary Fig. S5D) was identified by the MCODE plugin. GO enrichment analysis showed that mainly four categories of biological processes were affected by the top module, that is, positive regulation of leukocyte proliferation, positive regulation of response to cytokine stimulus, positive regulation of NIK/NF-kappaB signaling and morphogenesis of an endothelium (Supplementary Fig. S5E). KEGG pathway enrichment analysis for the top module demonstrated some concrete immune-related pathways, which were more apparently associated with cancer, and some immunodeficiency virus infection pathways (Supplementary Fig. S5F). Subsequently, we evaluated the association of top 30 hub genes, and genes in the top module with OS using KM analysis. We found that LY75, CD86, TNFSF10, CSF2RB, CD69, CCR2, CASP1, CD27, IKZF1, IRF8, CD247, and CXCL10 were positively correlated with OS, whereas DMT3A and AGO2 were negatively correlated with OS (Supplementary Fig. S6A). Similarly, we performed the KM analysis for the lncRNAs in the prognosis-related ceRNA network, and found that high expression of LINC00641, PRKCZ-AS1, LINC01355, and KIZ-AS1 brings benefit to OS. Although high expression of PIK3CD-AS1, LINC01678, LINC00892, LINC01146, PSMB8-AS1, LINC01907, and CKMT2-AS1 reduced the OS (Supplementary Fig. S6B).

To validate the hypothesis that lncRNA positively regulates mRNA expression through interacting with miRNA in the ceRNA network, we calculated the PCC between the target lncRNAs and target mRNAs of nine key miRNAs, individually (Supplementary Table S15). The interaction pairs with |PCC| > 0.4 and a P value of <0.05 were considered as co-expression pairs. Finally, we identified LINC00641 positively regulated AGO2 through hsa-mir-125b-2; LINC01146 positively regulated CD86, CXCL10, and CASP1 through hsa-mir-135a-1; LINC00892 positively regulated CD69 through hsa-mir-181a-2 and CD27 through hsa-mir-214; LINC00892 and LINC01146 simultaneously positively regulated CD247, CCR2, and CSF2RB through hsa-mir-214; LINC01146 positively regulated IKZF1 through hsa-mir-92a-2 (Supplementary Fig. S6C). The subnetwork was visualized in Fig. 3D. In addition, the relationships among clinical variables and the expression of nodes in the subnetwork were assessed to distinguish the sarcoma by subtype, sex, histological diagnosis, residual tumor, and tumor total necrosis (Fig. 3E).

Validation of gene expression in vivo and in vitro

We screened the Expression Atlas database to verify differences in gene expression between sarcoma and normal samples both in vivo and in vitro (Table 1). Specifically, in the soft tissue sarcoma cell line HT1080, we found that both CXCL10 and CASP1 were upregulated in the cell treated with IFNα and/or U0126, a MEK pathway inhibitor, for 6 and 12 hours, respectively (19). Furthermore, CXCL10 was also found to be upregulated in atypical teratoid/rhabdoid tumor (20, 21). CSF2RB was found to be downregulated in myxosarcoma (22). However, when the 93T449 cell line, a well-differentiated liposarcoma from the retroperitoneum, was exposed by hypoxia, CXCL10 and CASP1 were found to be downregulated comparing with normoxia condition (23). In addition, in a mouse model with undifferentiated pleomorphic sarcoma induced by Kras and p53 mutations, Cd86, Cxcl10, Casp1, Ccr2, and Csf2rb were found to be upregulated except Cd27 was downregulated (24). Furthermore, in mouse model with clear cell sarcoma induced by TAT-Cre, Cxcl10, Ccr2, and Csf2rb were identified to be upregulated, but Cd69, Cd27, and Ikzf1 were downregulated (25). Similar results were also found in clear cell sarcoma mouse model induced by Rosa26CreER, except Ccr2 was reversely dysregulated (25). In addition, Cxcl10 and Csf2rb were also upregulated in synovial sarcoma, whereas Cd69, Cd27, and Ikzf1 were downregulated (25).

Table 1.

Validation of the expression of genes in the subnetwork through Expression Atlas database.

GenePadjLog2(FC)Control/CaseSpeciesGEOPMID
Cd86 9.547e−10 2.3 4/17 Mus musculus GSE16779 19956606 
CXCL10 2.213e−6 4.4 3/3 Homo sapiens GSE31019 22970192 
 1.132e−7 5.5 3/3    
 4.611e−4 1.7 3/3    
 6.742e−3 1.2 3/3    
 4.489e−6 2.5 3/3    
 6.136e−3 1.5 7/20  GSE35493 23382118 
 5.142e−4 1.1 13/17  GSE66354 25968456 
 9.474e−3 −1.7 3/3  GSE21050 29423096 
Cxcl10 6.920e−3 2.2 6/4 Mus musculus GSE41293 & GSE43045 23410975 
 8.102e−5 2.1 6/6    
 1.722e−5 4.1 4/17  GSE16779 19956606 
CASP1 4.295e−3 1.1 3/3 Homo sapiens GSE31019 22970192 
 5.871e−4 1.3 3/3    
 1.318e−2 1.2 3/3    
 7.113e−4 1.3 3/3    
 4.054e−2 −1.3 3/3  GSE21050 29423096 
Casp1 6.832e−7 1.9 4/17 Mus musculus GSE16779 19956606 
Cd69 3.674e−3 −2.6 6/4 Mus musculus GSE41293 & GSE43045 23410975 
 2.914e−3 −2.5 6/6    
 1.628e−6 −3.9 6/5    
CD27 3.549e−2 1.6 3/3 Homo sapiens GSE21050 29423096 
Cd27 4.781e−5 −3.2 6/6 Mus musculus GSE41293 & GSE43045 23410975 
 1.863e−7 −4.6 6/5    
 4.494e−5 −3.7 6/4    
 3.137e−6 −1.4 4/17  GSE16779 19956606 
Cd247 2.764e−9 −5.3 6/5 Mus musculus GSE41293 & GSE43045 23410975 
 5.487e−3 −3.1 6/4    
 3.316e−4 −2.9 6/6    
Ccr2 7.983e−10 4.6 4/17 Mus musculus GSE16779 19956606 
 8.511e−3 1.1 6/4  GSE41293 & GSE43045 23410975 
 1.386e−2 −1.1 6/5    
CSF2RB 6.559e−9 −3.1 9/6 Homo sapiens PRJEB36314 32330934 
Csf2rb 1.434e−10 2.4 6/4 Mus musculus GSE41293 & GSE43045 23410975 
 1.025e−5 1.2 6/5    
 2.238e−2 6/6    
 3.808e−6 1.9 4/17  GSE16779 19956606 
Ikzf1 1.490e−3 −2 6/4 Mus musculus GSE41293 & GSE43045 23410975 
 1.610e−3 −2 6/6    
 2.347e−3 −1.8 6/5    
GenePadjLog2(FC)Control/CaseSpeciesGEOPMID
Cd86 9.547e−10 2.3 4/17 Mus musculus GSE16779 19956606 
CXCL10 2.213e−6 4.4 3/3 Homo sapiens GSE31019 22970192 
 1.132e−7 5.5 3/3    
 4.611e−4 1.7 3/3    
 6.742e−3 1.2 3/3    
 4.489e−6 2.5 3/3    
 6.136e−3 1.5 7/20  GSE35493 23382118 
 5.142e−4 1.1 13/17  GSE66354 25968456 
 9.474e−3 −1.7 3/3  GSE21050 29423096 
Cxcl10 6.920e−3 2.2 6/4 Mus musculus GSE41293 & GSE43045 23410975 
 8.102e−5 2.1 6/6    
 1.722e−5 4.1 4/17  GSE16779 19956606 
CASP1 4.295e−3 1.1 3/3 Homo sapiens GSE31019 22970192 
 5.871e−4 1.3 3/3    
 1.318e−2 1.2 3/3    
 7.113e−4 1.3 3/3    
 4.054e−2 −1.3 3/3  GSE21050 29423096 
Casp1 6.832e−7 1.9 4/17 Mus musculus GSE16779 19956606 
Cd69 3.674e−3 −2.6 6/4 Mus musculus GSE41293 & GSE43045 23410975 
 2.914e−3 −2.5 6/6    
 1.628e−6 −3.9 6/5    
CD27 3.549e−2 1.6 3/3 Homo sapiens GSE21050 29423096 
Cd27 4.781e−5 −3.2 6/6 Mus musculus GSE41293 & GSE43045 23410975 
 1.863e−7 −4.6 6/5    
 4.494e−5 −3.7 6/4    
 3.137e−6 −1.4 4/17  GSE16779 19956606 
Cd247 2.764e−9 −5.3 6/5 Mus musculus GSE41293 & GSE43045 23410975 
 5.487e−3 −3.1 6/4    
 3.316e−4 −2.9 6/6    
Ccr2 7.983e−10 4.6 4/17 Mus musculus GSE16779 19956606 
 8.511e−3 1.1 6/4  GSE41293 & GSE43045 23410975 
 1.386e−2 −1.1 6/5    
CSF2RB 6.559e−9 −3.1 9/6 Homo sapiens PRJEB36314 32330934 
Csf2rb 1.434e−10 2.4 6/4 Mus musculus GSE41293 & GSE43045 23410975 
 1.025e−5 1.2 6/5    
 2.238e−2 6/6    
 3.808e−6 1.9 4/17  GSE16779 19956606 
Ikzf1 1.490e−3 −2 6/4 Mus musculus GSE41293 & GSE43045 23410975 
 1.610e−3 −2 6/6    
 2.347e−3 −1.8 6/5    

In addition, we applied RNA-seq data from rhabdomyosarcoma rat model for further validation. In this study, all rats in the 12-month Ni (nickel)-implanted group and two of the 12-month Co (cobalt)-embedded rats developed tumors, commonly rhabdomyosarcoma, whereas the other 12-month Co-embedded rats only developed spindle cell mesenchymal tumor (26). Cd247 and Cd27 were identified to be upregulated in the first and third month, whereas Casp1, Cd86, Cxcl10, and Ikzf1 were only found to be upregulated in third month after Co-embedded compared with Ta (Tantalum)-embedded rats. In addition, in the first month, Csf2rb was upregulated after Co-embedded whereas Cxcl10 was downregulated after Ni-embedded. However, mRNA level of Ago2 was constant at any time points and in all metal-embedded groups included in split of the rest genes were also not significantly dysregulated in the sixth and twelfth month in the rats with metal implants (Fig. 4A).

Figure 4.

Validation for the expression of genes in the subnetwork. A, Expression of the genes in the subnetwork using RNA-seq data from rhabdomyosarcoma (n = 8 per group). The P values were calculated by the Wilcoxon rank-sum test. Ta, Co, and Ni represent the rats embedded with Ta (Tantalum, control group), Co (cobalt), and Ni (nickel), respectively. B, Expression of the genes in the subnetwork using S180-bearing mice model (n = 3 per group). C, Schematic representation of the experiments in vitro examining the expression of the genes between low and high TIIL subtypes using the SW872 cell line. D, Expression of marker genes of T cells (n = 3 per group). E, Expression of genes in low and high TIIL subtypes (n = 3 per group). The P values were calculated by the t test; *, P < 0.05; **, P < 0.01; ***, P < 0.0001.

Figure 4.

Validation for the expression of genes in the subnetwork. A, Expression of the genes in the subnetwork using RNA-seq data from rhabdomyosarcoma (n = 8 per group). The P values were calculated by the Wilcoxon rank-sum test. Ta, Co, and Ni represent the rats embedded with Ta (Tantalum, control group), Co (cobalt), and Ni (nickel), respectively. B, Expression of the genes in the subnetwork using S180-bearing mice model (n = 3 per group). C, Schematic representation of the experiments in vitro examining the expression of the genes between low and high TIIL subtypes using the SW872 cell line. D, Expression of marker genes of T cells (n = 3 per group). E, Expression of genes in low and high TIIL subtypes (n = 3 per group). The P values were calculated by the t test; *, P < 0.05; **, P < 0.01; ***, P < 0.0001.

Close modal

To further verify the expression of the genes in the subnetwork, we sacrificed 10 mice bearing S180 tumor cells. And on the basis of the subtype's confirmation, we found that CXCL10, IKZF1, CD27, CD69, CCR2, and CSF2RB showed significant upregulation from low to high TIILs subtypes. Though the expression of the remaining four genes has no statistical significance, the direction of their dysregulation is consistent with our analysis (Fig. 4B). In addition, we also designed experiments in vitro to further confirmation the expression of the genes (Fig. 4C). As the experiment designed, we co-cultured SW872 cell lines with low and high PBMC infiltration, respectively (Fig. 4D). The RT-qPCR results also consistently demonstrated the dysregulation of the genes in the subnetwork except for AGO2, but the expression of AGO2 also showed decreased from low to high TIILs subtypes (Fig. 4E).

Exploration of potential drug targets

Next, CMap analysis was applied to identifying potentially effective drugs on patients with sarcoma. The results indicated that TSA, a histone deacetylases (HDAC) inhibitor, was a critical drug that might have therapeutic value (Fig. 5A), and the IC50 value of TSA varied from subtypes (Fig. 5B). In addition, the IC50 value of TSA was significantly affected by high or low expression of nodes in the subnetwork. Higher IC50 value was associated with high-expression levels of CASP1, CCR2, CD247, CD27, CD69, CD86, CSF2RB, CXCL10, IKZF1, LINC00892, and LINC01146. Meanwhile, it is also related to low expression of AGO2, LINC00641, hsa-mir-135a-1, and hsa-mir-181a-2 (Supplementary Fig. S7A and S7B). Because TSA is a pan-HDAC inhibitor, we tried to find some reasonable epigenetic modification evidence of TSA by analyzing the methylation levels of the genes in the subnetwork in healthy people and patients with sarcoma. There was a significant difference in the methylation levels of AGO2, CCR2, CD247, CD86, CD69, IKZF1, and CSF2RB between healthy people and patients with sarcoma. But CASP1, CXCL10, LINC00641, and CD27 showed no significant difference between these two groups (Supplementary Fig. S7C).

Figure 5.

Key drugs that had potential therapeutic effects on patients with sarcoma. A, Identification of key drugs through CMap analysis. B, Comparison of drug sensitivity among subtypes of sarcoma. The P values were calculated by the Wilcoxon rank-sum test. C, Schematic diagram showing the schedule of treatment of S180-bearing mice. D, TSA treatment reduced HDAC activity changes among mice bearing S180. E, Tumor size among different groups. F, Visualization for the correlation between the nodes in the subnetwork and IC50 value of TSA. G, UMAP visualizations of single-cell RNA-seq from metastatic lung samples of mice bearing M3–9-M rhabdomyosarcoma tumor cells, colored by cell types. H, UMAP plots showing expression of the genes in the subnetwork. I, Schematic diagram showing the schedule of treatment of the SW872 cell line. J, The expression of markers of CD4+/CD8+ T cells. K, The expression comparison of the genes in the subnetwork between low TIIL SW872 cells with and without TSA treatment. L, The expression comparison of the genes in the subnetwork between high TIIL SW872 cells with and without TSA treatment. TSA0.2 and TSA1 mean the concentration of TSA is 0.2 and 1.0 μmol/L, respectively. The P values were calculated by the t test; *, P < 0.05; **, P < 0.01; ***, P < 0.0001.

Figure 5.

Key drugs that had potential therapeutic effects on patients with sarcoma. A, Identification of key drugs through CMap analysis. B, Comparison of drug sensitivity among subtypes of sarcoma. The P values were calculated by the Wilcoxon rank-sum test. C, Schematic diagram showing the schedule of treatment of S180-bearing mice. D, TSA treatment reduced HDAC activity changes among mice bearing S180. E, Tumor size among different groups. F, Visualization for the correlation between the nodes in the subnetwork and IC50 value of TSA. G, UMAP visualizations of single-cell RNA-seq from metastatic lung samples of mice bearing M3–9-M rhabdomyosarcoma tumor cells, colored by cell types. H, UMAP plots showing expression of the genes in the subnetwork. I, Schematic diagram showing the schedule of treatment of the SW872 cell line. J, The expression of markers of CD4+/CD8+ T cells. K, The expression comparison of the genes in the subnetwork between low TIIL SW872 cells with and without TSA treatment. L, The expression comparison of the genes in the subnetwork between high TIIL SW872 cells with and without TSA treatment. TSA0.2 and TSA1 mean the concentration of TSA is 0.2 and 1.0 μmol/L, respectively. The P values were calculated by the t test; *, P < 0.05; **, P < 0.01; ***, P < 0.0001.

Close modal

To validate the therapeutic value of TSA, subcutaneous tumor models were first constructed by injecting S180 cells. When the tumor volume reached about 100 mm3, the inoculated mice were then randomly allocated to the following groups: DMSO, TSA 0.5 mg/kg, TSA 2 mg/kg (Fig. 5C). After two weeks administration, we detected the HDAC activity in tumor tissues. We found that TSA significantly reduced the HDAC activity and high-dose administration of TSA contributed to stronger inhibition of HDAC (Fig. 5D). We further found that TSA could obviously inhibit the growth of tumor (Fig. 5E; Supplementary Fig. S7D and S7E) without significant effects on the body weight (Supplementary Fig. S7F).

In addition, the IC50 value of TSA was identified to be positively correlated with the expression of mRNAs and lncRNAs in the subnetwork except for AGO2 and LINC00641. All miRNAs in the subnetwork were weakly and negatively correlated with the IC50 value of TSA (Fig. 5F). We further identified the co-expression patterns among 22 tumor-infiltrating immune cell types through correlation analysis (Supplementary Fig. S8A). In addition, five miRNAs and three lncRNAs in the subnetwork were significantly correlated with infiltration levels of 19 immune cells except activated DCs, gamma delta T cells and regulatory T cells (Supplementary Fig. S8B). We also found that the infiltration of many immune cells was significantly different between low- and high-expression groups of miRNAs and lncRNAs (Supplementary Fig. S8C). These results gave us some clues that TSA might regulate TME through targeting the nodes in the subnetwork. We then downloaded scRNA-seq data from metastatic lung of mice bearing M3–9-M rhabdomyosarcoma tumor cells for further confirmation. After determining the parameters of the features: Cell barcodes with >500 and <4,500 genes detected and <20% mitochondrial gene expression (Supplementary Fig. S8D), we detected 17 cell clusters (Supplementary Fig. S8E) and the DE genes in each cluster (Supplementary Table S16). Cell types were annotated with the reference of MCA3.0 (27), CellMarker2.0 (RRID: SCR_018503), SingleR.MouseRNAseqData (28), and SingleR.ImmGenData (ref. 29; Supplementary Table S16). Finally, we determined the final cell types and visualized in UMAP plot (Fig. 5G). The expression of specific marker genes was visualized in bubble plot (Supplementary Fig. S8F). Subsequently, to figure out whether the genes are expressed specifically in any cell type, we investigated the expression the genes in the subnetwork via this scRNA-seq data. We found that Cd27, Cd247, and Cd69 showed specific expression in B cells and/or Cd8+ T cells, Ccr2 was specifically expressed among monocytes, NK cells, and Cd8+ T cells, Cd86 was overexpressed in B cells and monocytes, Csf2rb showed relatively high expression in neutrophil, eosinophil, and monocytes. In contrary to Ago2, Casp1, Ikzf1, those were almost high expressed among all cell types, Cxcl10 showed relatively conservative expression among all cell types (Fig. 5H).

We further demonstrated that TME was altered within S180-bearing mice injected with different dose of TSA. We found that the infiltration of CD4+/CD8+ T cells was enhanced by TSA (Supplementary Fig. S9; Supplementary Table S17). Meanwhile, we investigated the therapeutic effects of TSA in the SW872 cell line (Fig. 5I), we found that for high TILLs sarcoma subtype, TSA could significantly increase the infiltration of CD4+/CD8+ T cells. However, for low TIILs sarcoma subtype, TSA slightly changed the infiltration of CD4+/CD8+ T cells without significance (Fig. 5J). Furthermore, we found that AGO2, CASP1, CXCL10, CD86, and CSF2RB showed significant dysregulation after adding TSA (1 μmol/L) to low TIILs SW872 cells with upregulation for first two genes and downregulation for the rest three genes (Fig. 5K). For high TIILs SW872 cells, the expressions of CXCL10, CD86, CASP1, and CD69 were significantly inhibited after adding TSA (1 μmol/L), but CD27 and CSF2RB were significantly enhanced (Fig. 5L).

Although sarcoma constitutes only approximately 1% of all human malignancies cases, it presents the second most common solid tumor in children and adolescents (30). Currently, there are more than 100 histological subtypes characterized, and the number of subtypes is continuously increasing due to newly identified molecular profiling. Therefore, the novelty with more precise classification of sarcoma is still imperative. Clues from the previous study showed the significant immunogenic response occurred in rhabdomyosarcoma rats with early stages (26). Furthermore, the infiltration of immune cells in pleomorphic sarcoma is associated with the tumor morphology and anatomical location, and closely associated with OS of patients with sarcoma (31). All the evidence contributes to our focus on the immune infiltration of TME. In our study, we integrated the 45 immune-related signatures with more classical immune cells and curated immunogenetic pathways, and performed ssGSVA analysis to estimate the infiltration for every single patient. The results are obviously, significantly, and robustly showed there are three subtypes among patients with sarcoma, that is, C1 (Low), C2 (High), and C3 (Moderate). Differing from the general classification based on the different estimated scores of 22 popular immune cell types (8, 32), we collected more comprehensive signatures, which can depict the complete status of the immune infiltration. The approach to scoring each signature for each sarcoma patient using more original algorithm, that is, ssGSVA benefits for the more precise classification of sarcoma and promotes to the complete description of the characteristics of patients with sarcoma, such as the higher the levels of immune infiltration, the smaller the proportion of tumor cells and the more complicated the TME in the patient. And the consistent classification was also identified using Consensus Clustering analysis. Besides, from the number of patients with different tumor sites, we realized that the TME stratification was a common characteristic for soft tissue sarcoma. These results were further confirmed through a small cohort of S180-bearing mice model. Moreover, we infer that a large difference exists between C1 (Low) and C2 (High), and slight difference between C1 (Low) and C3 (Moderate) from the OS of subtypes and DE genes among subtypes.

With the classification of sarcoma, our study further revealed that dysregulated ceRNA network was involved in the progression of sarcoma and identified the potential prognosis biomarker. In our study, we first developed an S180-bearing mouse model to validate the sarcoma subtypes, and further validated the expression of the genes in the subnetwork in vivo and in vitro. Moreover, according to the OS of three subtypes, we found that the lower the TIILs, the lower survival rates in patients with sarcoma. Consistently, the five key miRNAs in the subnetwork are all dangerous factors according to the univariate Cox regression analysis (HR > 1) and upregulated in patients with sarcoma with low TIILs (C1 subtype) compared with patients with sarcoma with high TIILs (C2 subtype). Meanwhile, the individual KM analysis of these five miRNAs also demonstrated the poor survival rates in their high expression group of them. Accordingly, in the subnetwork, the target mRNAs, including CASP1, CCR2, CSF2RB, and IKZF1, and target lncRNAs, including LINC01146 and LINC00892, are found to be protective factors, downregulated in C1 subtype, and lower survival rates in their low expression groups. Only target lncRNA LINC00641 in the subnetwork is found to be a dangerous factor, upregulated in C1 subtype and lower survival rates in its high-expression group. Besides, AGO2 is upregulated in C1 subtype and shows lower survival rates in its high-expression group. However, CD86, CXCL10, CD69, CD27, and CD247 are identified to be downregulated in C1 subtype and show lower survival rates in their low expression groups in spite of no statistical difference identified in their univariate Cox regression analysis. Actually, studies have revealed the important role of these molecules in sarcoma. MiR-125b was found to inhibit cell biological progression of Ewing's sarcoma by suppressing the PI3K/Akt signaling pathway (33), and develop chemoresistance in Ewing's sarcoma/primitive neuroectodermal (34). Upregulation of miR-181a/miR-212 were found to improve myogenic commitment in murine fusion-negative rhabdomyosarcoma (35). MiR-214–3p was also found to be commonly downregulated by EWSFI1 and by CD99 and its restoration limits Ewing sarcoma aggressiveness (36). MiR-92a was also identified to modulate proliferation, apoptosis, migration, and invasion of osteosarcoma cell lines by targeting Dickkopf-related protein 3 (37). These studies experimentally demonstrated the mediated regulatory mechanism of these five key miRNAs in the occurrence, development, and metastasis of sarcoma. Furthermore, AGO2 has been proved to be a key effector of miRNA-induced silencing complexes and assembles with miRNA to form the complexes (38). It has been verified that the expression of costimulatory molecules CD80 and/or CD86 by a Kaposi's sarcoma tumor cell line induces differential T-cell activation and proliferation (39). In the tumor-inflammatory microenvironment, CASP1 and its processed cytokines, such as IL1β, play an important role in the occurrence and development of cancer (40). Besides, reduced CCR2 is identified to improve the prognosis of sarcoma by remodeling the TME (41). Genomic alterations of IKZF1 could negatively affect immunogenicity and tumor response to immune checkpoint blockade (42). All of these studies are the evidence that the construction of the ceRNA in our research is of reliability and robustness. Moreover, the molecules in the subnetwork are full of potential to be prognosis biomarkers.

Notably, our study also identified that epigenetic modification is closely associated with sarcoma development, though a lot of immune-related pathways are significantly enriched through the DE genes between C1 and C2 subtypes, genes in the prognosis-related ceRNA network, and parent genes of proteins in the PPI network, respectively. Concretely, histone modification, especially deacetylation, was enriched for the genes in the prognosis-related ceRNA network. HDACs can promote deacetylation of histones and tighten their interaction with DNA, resulting in a closed chromatin structure and the inhibition of gene transcription (43). Studies have shown that HDACs influence diverse cellular processes and contribute to sarcoma growth and progression by multiple mechanisms (44). Consistently, we identified TSA, which is pan-HDAC inhibitor of class I and II HDACs, had a potential therapeutic effect on the sarcoma subtypes in our study. In various cancer cells, the shift to an increased acetylation/deacetylation ratio by HDAC inhibitors (HDI) was found to have a substantial effect on their fate (45). HDIs are found to inhibit the proliferation of a variety of transformed cells in vitro and tumor progression in several solid tumors and hematological malignancies (46), induce cell-cycle arrest, differentiation, cell death, and modulate the immune response and decrease angiogenesis (47). Therefore, HDIs seem to be promising anticancer drugs particularly in the combination with other anticancer drugs and/or radiotherapy. HDIs have been found to upregulate tumor-suppressor genes, downregulate oncogenes, induce apoptosis and cell-cycle arrest, decrease invasion, metastasis and angiogenesis, inhibit tumor growth through regulating autophagy, induce reactive oxygen species generation, and induce tumor cell differentiation in sarcomas (48). In our study, TSA was identified to be therapeutic potential to the subtypes of sarcoma, and potentially reverse the low immune infiltration in patients with sarcoma thereby improving the survival of patients. In our in vivo experiments, we confirmed the therapeutic value of TSA, and further found that higher dose of TSA could significantly inhibit the growth of tumors. Except for the dysregulation of the genes in the subnetwork demonstrated in vivo and in vitro, we further demonstrated that some of them showed inconsistent regulation between low and high TIILs subtypes after adding TSA, such as CASP1, CD27, CD247, CCR2, and CSF2RB. This result may imply the complicated regulation of these effectors or the complicated therapeutic effects of TSA. However, this study offers new insights into the relationship between HDI and sarcoma.

Nevertheless, there are still some limitations in this study. First, the main results are analyzed on the basis of TCGA data, which is relatively small and lack of large number of paired samples. Second, the validation experiments are relatively simple, the more exact molecular mechanisms involved the nodes and their corresponding interactions in the subnetwork are worth of deeper investigation.

No disclosures were reported.

The content is solely the responsibility of the authors and does not necessarily represent the official views of the NIH.

D. Leng: Conceptualization, resources, data curation, software, formal analysis, validation, investigation, visualization, methodology, writing–original draft, writing–review and editing. Z. Yang: Validation, investigation, methodology, writing–review and editing. H. Sun: Validation, investigation, methodology, writing–review and editing. C. Song: Formal analysis, visualization, methodology, writing–review and editing. C. Huang: Investigation, methodology, writing–review and editing. K.U. Ip: Validation, methodology. G. Chen: Conceptualization, investigation, methodology, writing–review and editing. C.-X. Deng: Conceptualization, supervision, funding acquisition, investigation, writing–review and editing. X.D. Zhang: Conceptualization, supervision, funding acquisition, investigation, methodology, writing–review and editing. Q. Zhao: Conceptualization, resources, supervision, funding acquisition, methodology, project administration, writing–review and editing.

This study was supported by the National Key R&D Program of China (2019YFA0904400), the Science and Technology Development Fund, Macau SAR (FDCT/0043/2021/A1, FDCT/0004/2019/AFJ and FDCT/0065/2021/A), the Natural Science Foundation of Guangdong Province (2023A1515010549), Shenzhen Science and Technology Project (SGDX2020110309280301), the University of Macau (MYRG2022–00143-FHS), the Ministry of Education Frontiers Science Centre for Precision Oncology, University of Macau (SP2023–00001-FSCPO), Dr. Stanley Ho Medical Development Foundation (SHMDF-VSEP/2022/002), and by U.S. National Institutes of Health (through Grants UL1TR001998, 1U01DK135111 and OT2HL161847) and by the DRC at Washington University (grant no. P30 DK020579). The authors thank the funding supported by the Science and Technology Development Fund of Macau (FDCT/0043/2021/A1 and FDCT/0002/2021/AKP) and Zhongnanshan Medical Foundation of Guangdong Province (ZNSA-2021016). The authors also thank the Office of Scientific Writing at University of Kentucky's College of Public Health for assistance preparing this article as well as Prof. Hubing Shi at Sichuan University.

The publication costs of this article were defrayed in part by the payment of publication fees. Therefore, and solely to indicate this fact, this article is hereby marked “advertisement” in accordance with 18 USC section 1734.

Note: Supplementary data for this article are available at Clinical Cancer Research Online (http://clincancerres.aacrjournals.org/).

1.
Skubitz
KM
,
D'Adamo
DR
.
Sarcoma
.
Mayo Clin Proc
2007
;
82
:
1409
32
.
2.
Stiller
CA
,
Trama
A
,
Serraino
D
,
Rossi
S
,
Navarro
C
,
Chirlaque
MD
, et al
.
Descriptive epidemiology of sarcomas in Europe: report from the RARECARE project
.
Eur J Cancer
2013
;
49
:
684
95
.
3.
Raj
S
,
Miller
LD
,
Triozzi
PL
.
Addressing the adult soft tissue sarcoma microenvironment with intratumoral immunotherapy
.
Sarcoma
2018
;
2018
:
9305294
.
4.
Ayodele
O
,
Razak
ARA
.
Immunotherapy in soft-tissue sarcoma
.
Current Oncol
2020
;
27
:
17
23
.
5.
Woll
PJ
,
Reichardt
P
,
Le Cesne
A
,
Bonvalot
S
,
Azzarelli
A
,
Hoekstra
HJ
, et al
.
Adjuvant chemotherapy with doxorubicin, ifosfamide, and lenograstim for resected soft-tissue sarcoma (EORTC 62931): a multicentre randomised controlled trial
.
Lancet Oncol
2012
;
13
:
1045
54
.
6.
Yuan
J
,
Li
X
,
Yu
S
.
Molecular targeted therapy for advanced or metastatic soft tissue sarcoma
.
Cancer Control
2021
;
28
:
10732748211038424
.
7.
Watnick
RS
.
The role of the tumor microenvironment in regulating angiogenesis
.
Cold Spring Harb Perspect Med
2012
;
2
:
a006676
.
8.
Zhu
N
,
Hou
J
.
Assessing immune infiltration and the tumor microenvironment for the diagnosis and prognosis of sarcoma
.
Cancer Cell Int
2020
;
20
:
577
.
9.
Dunn
GP
,
Old
LJ
,
Schreiber
RD
.
The three Es of cancer immunoediting
.
Annu Rev Immunol
2004
;
22
:
329
60
.
10.
Coussens
LM
,
Werb
Z
.
Inflammation and cancer
.
Nature
2002
;
420
:
860
7
.
11.
Zhang
K
,
Zhang
L
,
Mi
Y
,
Tang
YC
,
Ren
FF
,
Liu
B
, et al
.
A ceRNA network and a potential regulatory axis in gastric cancer with different degrees of immune cell infiltration
.
Cancer Sci
2020
;
111
:
4041
50
.
12.
Wang
Y
,
Gao
Y
,
Guo
S
,
Chen
Z
.
Integrated analysis of lncRNA-associated ceRNA network identified potential regulatory interactions in osteosarcoma
.
Genet Mol Biol
2020
;
43
:
e20190090
.
13.
Chen
DL
,
Lu
YX
,
Zhang
JX
,
Wei
XL
,
Wang
F
,
Zeng
ZL
, et al
.
Long non-coding RNA UICLM promotes colorectal cancer liver metastasis by acting as a ceRNA for microRNA-215 to regulate ZEB2 expression
.
Theranostics
2017
;
7
:
4836
49
.
14.
Zhou
C
,
Chen
Z
,
Xiao
B
,
Xiang
C
,
Li
A
,
Zhao
Z
, et al
.
Comprehensive analysis of GINS subunits prognostic value and ceRNA network in sarcoma
.
Front Cell Dev Biol
2022
;
10
:
951363
.
15.
Charoentong
P
,
Finotello
F
,
Angelova
M
,
Mayer
C
,
Efremova
M
,
Rieder
D
, et al
.
Pan-cancer immunogenomic analyses reveal genotype-immunophenotype relationships and predictors of response to checkpoint blockade
.
Cell Rep
2017
;
18
:
248
62
.
16.
Bhattacharya
S
,
Dunn
P
,
Thomas
CG
,
Smith
B
,
Schaefer
H
,
Chen
JM
, et al
.
ImmPort, toward repurposing of open access immunological assay data for translational and clinical research
.
Sci Data
2018
;
5
:
180015
.
17.
Ding
RB
,
Chen
P
,
Rajendran
BK
,
Lyu
X
,
Wang
H
,
Bao
J
, et al
.
Molecular landscape and subtype-specific therapeutic response of nasopharyngeal carcinoma revealed by integrative pharmacogenomics
.
Nat Commun
2021
;
12
:
3046
.
18.
Chaussabel
D
,
Baldwin
N
.
Democratizing systems immunology with modular transcriptional repertoire analyses
.
Nat Rev Immunol
2014
;
14
:
271
80
.
19.
Christian
SL
,
Zu
D
,
Licursi
M
,
Komatsu
Y
,
Pongnopparat
T
,
Codner
DA
, et al
.
Suppression of IFN-induced transcription underlies IFN defects generated by activated Ras/MEK in human cancer cells
.
PLoS ONE
2012
;
7
:
e44267
.
20.
Birks
DK
,
Donson
AM
,
Patel
PR
,
Sufit
A
,
Algar
EM
,
Dunham
C
, et al
.
Pediatric rhabdoid tumors of kidney and brain show many differences in gene expression but share dysregulation of cell cycle and epigenetic effector genes
.
Pediatr Blood Cancer
2013
;
60
:
1095
102
.
21.
Griesinger
AM
,
Josephson
RJ
,
Donson
AM
,
Levy
JMM
,
Amani
V
,
Birks
DK
, et al
.
Interleukin-6/STAT3 pathway signaling drives an inflammatory phenotype in group a ependymoma
.
Cancer Immunol Res
2015
;
3
:
1165
74
.
22.
Scalise
M
,
Torella
M
,
Marino
F
,
Ravo
M
,
Giurato
G
,
Vicinanza
C
, et al
.
Atrial myxomas arise from multipotent cardiac stem cells
.
Eur Heart J
2020
;
41
:
4332
45
.
23.
Yang
L
,
Forker
L
,
Irlam
JJ
,
Pillay
N
,
Choudhury
A
,
West
CML
.
Validation of a hypoxia related gene signature in multiple soft tissue sarcoma cohorts
.
Oncotarget
2018
;
9
:
3946
55
.
24.
Mito
JK
,
Riedel
RF
,
Dodd
L
,
Lahat
G
,
Lazar
AJ
,
Dodd
RD
, et al
.
Cross species genomic analysis identifies a mouse model as undifferentiated pleomorphic sarcoma/malignant fibrous histiocytoma
.
PLoS ONE
2009
;
4
:
e8075
.
25.
Straessler
KM
,
Jones
KB
,
Hu
H
,
Jin
HF
,
van de Rijn
M
,
Capecchi
MR
.
Modeling clear cell sarcomagenesis in the mouse: cell of origin differentiation state impacts tumor characteristics
.
Cancer Cell
2013
;
23
:
215
27
.
26.
Wen
Y
,
Vechetti
IJ
,
Alimov
AP
,
Hoffman
JF
,
Vergara
VB
,
Kalinich
JF
, et al
.
Time-course analysis of the effect of embedded metal on skeletal muscle gene expression
.
Physiol Genomics
2020
;
52
:
575
87
.
27.
Han
X
,
Wang
R
,
Zhou
Y
,
Fei
L
,
Sun
H
,
Lai
S
, et al
.
Mapping the mouse cell atlas by microwell-seq
.
Cell
2018
;
172
:
1091
107
.
28.
Benayoun
BA
,
Pollina
EA
,
Singh
PP
,
Mahmoudi
S
,
Harel
I
,
Casey
KM
, et al
.
Remodeling of epigenome and transcriptome landscapes with aging in mice reveals widespread induction of inflammatory responses
.
Genome Res
2019
;
29
:
697
709
.
29.
Heng
TS
,
Painter
MW
,
Immunological Genome Project
C
.
The Immunological Genome Project: networks of gene expression in immune cells
.
Nat Immunol
2008
;
9
:
1091
4
.
30.
Grunewald
TG
,
Alonso
M
,
Avnet
S
,
Banito
A
,
Burdach
S
,
Cidre-Aranaz
F
, et al
.
Sarcoma treatment in the era of molecular medicine
.
EMBO Mol Med
2020
;
12
:
e11131
.
31.
Wustrack
RL
,
Shao
E
,
Sheridan
J
,
Zimel
M
,
Cho
SJ
,
Horvai
AE
, et al
.
Tumor morphology and location associate with immune cell composition in pleomorphic sarcoma
.
Cancer Immunol Immunother
2021
;
70
:
3031
40
.
32.
Deng
J
,
Zeng
W
,
Kong
W
,
Shi
Y
,
Mou
X
.
The study of sarcoma microenvironment heterogeneity associated with prognosis based on an immunogenomic landscape analysis
.
Front Bioeng Biotechnol
2020
;
8
:
1003
.
33.
Li
J
,
You
T
,
Jing
J
.
MiR-125b inhibits cell biological progression of Ewing's sarcoma by suppressing the PI3K/Akt signalling pathway
.
Cell Prolif
2014
;
47
:
152
60
.
34.
Iida
K
,
Fukushi
J
,
Matsumoto
Y
,
Oda
Y
,
Takahashi
Y
,
Fujiwara
T
, et al
.
miR-125b develops chemoresistance in Ewing sarcoma/primitive neuroectodermal tumor
.
Cancer Cell Int
2013
;
13
:
21
.
35.
Pozzo
E
,
Giarratana
N
,
Sassi
G
,
Elmastas
M
,
Killian
T
,
Wang
CC
, et al
.
Upregulation of miR181a/miR212 improves myogenic commitment in murine fusion-negative rhabdomyosarcoma
.
Front Physiol
2021
;
12
:
701354
.
36.
De Feo
A
,
Pazzaglia
L
,
Ciuffarin
L
,
Mangiagli
F
,
Pasello
M
,
Simonetti
E
, et al
.
miR-214–3p is commonly downregulated by EWS-FLI1 and by CD99 and its restoration limits Ewing sarcoma aggressiveness
.
Cancers
2022
;
14
:
1762
.
37.
Yu
HY
,
Song
H
,
Liu
L
,
Hu
S
,
Liao
YX
,
Li
G
, et al
.
MiR-92a modulates proliferation, apoptosis, migration, and invasion of osteosarcoma cell lines by targeting Dickkopf-related protein 3
.
Biosci Rep
2019
;
39
:
BSR20190410
.
38.
Zhang
HL
,
Zhao
X
,
Guo
YM
,
Chen
R
,
He
JF
,
Li
L
, et al
.
Hypoxia regulates overall mRNA homeostasis by inducing Met(1)-linked linear ubiquitination of AGO2 in cancer cells
.
Nat Commun
2021
;
12
:
5416
.
39.
Foreman
KE
,
Wrone-Smith
T
,
Krueger
AE
,
Nickoloff
BJ
.
Expression of costimulatory molecules CD80 and/or CD86 by a Kaposi's sarcoma tumor cell line induces differential T-cell activation and proliferation
.
Clin Immunol
1999
;
91
:
345
53
.
40.
Jin
H
,
Jin
X
,
Cao
BR
,
Wang
WB
.
Berberine affects osteosarcoma via downregulating the caspase-1/IL1 beta signaling axis
.
Oncol Rep
2017
;
37
:
729
36
.
41.
Wei
BX
,
Feng
H
,
Wu
H
.
Reduced CCR2 can improve the prognosis of sarcoma by remodeling the tumor microenvironment
.
Int J Gen Med
2022
;
15
:
3043
53
.
42.
Chen
JC
,
Perez-Lorenzo
R
,
Saenger
YM
,
Drake
CG
,
Christiano
AM
.
IKZF1 enhances immune infiltrate recruitment in solid tumors and susceptibility to immunotherapy
.
Cell Syst
2018
;
7
:
92
103
.
43.
Grunstein
M
.
Histone acetylation in chromatin structure and transcription
.
Nature
1997
;
389
:
349
52
.
44.
Schmidt
O
,
Nehls
N
,
Prexler
C
,
von Heyking
K
,
Groll
T
,
Pardon
K
, et al
.
Class I histone deacetylases (HDAC) critically contribute to Ewing sarcoma pathogenesis
.
J Exp Clin Cancer Res
2021
;
40
:
322
.
45.
Hull
EE
,
Montgomery
MR
,
Leyva
KJ
.
HDAC inhibitors as epigenetic regulators of the immune system: impacts on cancer therapy and inflammatory diseases
.
Biomed Res Int
2016
;
2016
:
8797206
.
46.
Miyanaga
A
,
Gemma
A
,
Noro
R
,
Kataoka
K
,
Matsuda
K
,
Nara
M
, et al
.
Antitumor activity of histone deacetylase inhibitors in non–small cell lung cancer cells: development of a molecular predictive model
.
Mol Cancer Ther
2008
;
7
:
1923
30
.
47.
Shanmugam
G
,
Rakshit
S
,
Sarkar
K
.
HDAC inhibitors: targets for tumor therapy, immune modulation, and lung diseases
.
Transl Oncol
2022
;
16
:
101312
.
48.
Xie
C
,
Wu
B
,
Chen
B
,
Shi
Q
,
Guo
J
,
Fan
Z
, et al
.
Histone deacetylase inhibitor sodium butyrate suppresses proliferation and promotes apoptosis in osteosarcoma cells by regulation of the MDM2-p53 signaling
.
Onco Targets Ther
2016
;
9
:
4005
13
.
This open access article is distributed under the Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International (CC BY-NC-ND 4.0) license.

Supplementary data