Huling Li
Dandan Lin
Xiaoyan Wang1
Kai Wang- 1School of Public Health, Xinjiang Medical University, Urumqi, China
- 2School of Continuing Education, Xinjiang Medical University, Urumqi, China
- 3Department of Medical Engineering and Technology, Xinjiang Medical University, Urumqi, China
Background: N6 methyladenosine (m6A)-related noncoding RNAs (including lncRNAs and miRNAs) are closely related to the development of cancer. However, the gene signature and prognostic value of m6A regulators and m6A-associated RNAs in regulating sarcoma (SARC) development and progression remain largely unexplored. Therefore, further research is required.
Methods: We obtained expression data for RNA sequencing (RNA-seq) and miRNAs of SARC from The Cancer Genome Atlas (TCGA) datasets. Correlation analysis and two target gene prediction databases (miRTarBase and LncBase v.2) were used to deduce m6A-related miRNAs and lncRNAs, and Cytoscape software was used to construct ceRNA-regulating networks. Based on univariate Cox regression and least absolute shrinkage and selection operator (LASSO) Cox regression analyses, an m6A-associated RNA risk signature (m6Ascore) model was established. Prognostic differences between subgroups were explored using Kaplan–Meier (KM) analysis. Risk score-related biological phenotypes were analyzed in terms of functional enrichment, tumor immune signature, and tumor mutation signature. Finally, potential immunotherapy features and drug sensitivity predictions for this model were also discussed.
Results: A total of 16 miRNAs, 104 lncRNAs, and 11 mRNAs were incorporated into the ceRNA network. The risk score was obtained based on RP11-283I3.6, hsa-miR-455-3p, and CBLL1. Patients were divided into two risk groups using the risk score, with patients in the low-risk group having longer overall survival (OS) than those in the high-risk group. The receiver operating characteristic (ROC) curves indicated that risk characteristic performed well in predicting the prognosis of patients with SARC. In addition, lower m6Ascore was also positively correlated with the abundance of immune cells such as monocytes and mast cells activated, and several immune checkpoint genes were highly expressed in the low-m6Ascore group. According to our analysis, lower m6Ascore may lead to better immunotherapy response and OS outcomes. The risk signature was significantly associated with the chemosensitivity of SARC. Finally, a nomogram was constructed to predict the OS in patients with SARC. The concordance index (C-index) for the nomogram was 0.744 (95% CI: 0.707–0.784). The decision curve analysis (DCA), calibration plot, and ROC curve all showed that this nomogram had good predictive performance.
Conclusion: This m6Ascore risk model based on m6A RNA methylation regulator-related RNAs may be promising for clinical prediction of prognosis and might contain potential biomarkers for treatment response prediction for SARC patients.
Introduction
Sarcomas (SARC) are comprised of an extensive, heterogeneous, and biologically diverse group of malignant tumors all of which are derived from mesenchymal cells (Hui, 2016). It has more than 100 different subtypes that can occur in any part of the body. Although there are more than 100 subtypes, sarcomas can be divided into two main types including soft-tissue sarcomas (STSs, accounting for 80%) and three major subtypes of bone sarcomas (BSs, including the osteosarcomas, chondrosarcomas, and Ewing’s (EW) sarcomas) (Doyle, 2014; Hui, 2016; Alavi et al., 2018; Ferguson and Turner, 2018; Hatina et al., 2019). SARC are characterized by a low incidence rate (approximately 1% of all malignancies in adults and 10–15% of all malignancies in pediatric cancers) but a poor prognosis in most cases (Hui, 2016; Siegel et al., 2019). Despite advances in the treatment of STS and osteosarcoma, with comprehensive treatment strategies including surgery, chemotherapy, radiotherapy, and targeted therapy, the survival rate of patients with advanced STSs and BSs needs to be improved (Zhu et al., 2019a; Zhang et al., 2021a; Eaton et al., 2021). Although it is difficult to determine the survival for each subtype of sarcoma due to the heterogeneity of the disease, the 5-year survival rates for patients with STS and for patients with bone sarcomas are generally larger than (>) 80% and about 70%, respectively. However, 5-year survival rates of the patients developing advanced-stage STS or various bone sarcomas are less than (<) 20% and between 22% and 57%, respectively (Bone Cancer, 2021; SarcomasTissue, 2021; SEER, 2021). In addition, about 50% of STS patients and about 90% of BS patients end up with distant metastasis, which remains a major cause of death and a barrier to effective treatment (Bacci et al., 1998; Italiano et al., 2011). Sarcomas have a poor prognosis due to their aggressive growth and high risk of metastasis (Aung et al., 2019). Therefore, further investigations into the pathogenesis of sarcomas and identification of the novel prognostic biomarkers that facilitate the improvement of therapy and prognosis of sarcomas are urgently needed.
N6-methyladenosine (m6A) is the prevalent modification in eukaryotic mRNAs, whose reversible methylation may have a profound impact on gene expression regulation (Niu et al., 2013), and plays a critical role in RNA processing. m6A RNA methylation is dynamically regulated by the corresponding m6A regulators, which are well classified into three subtypes, namely, “methylases—writer,” “demethylases—erasers,” or “m6A binding proteins—readers” depending on different functions (Yang et al., 2018; Chen et al., 2019a; Barbieri and Kouzarides, 2020). Increasing evidence has revealed the correlation between m6A and human cancers, including ovarian cancer (Li et al., 2021a), hepatocellular carcinoma (Li et al., 2021b), and colorectal cancer (Ji et al., 2020). The literature has also indicated that m6A-related genes might serve as novel prognostic biomarkers for different cancers, and that the m6A-related risk score may assist in risk assessment and prognostic stratification (Li et al., 2021c; Xu et al., 2021). The m6A regulators are also linked with the progression and prognosis of SARC. In recent years, there have been some studies on the role of m6A methylation regulators in sarcomas. For instance, WTAP, which is highly expressed in osteosarcoma tissues and associated with the worse prognosis of osteosarcoma patients, was found to potentially promote osteosarcoma progression by inhibiting HMBOX1 in an m6A-dependent manner in vitro and in vivo (Chen et al., 2020). Hou et al. (2020) evaluated the relationship between copy number variations (CNVs) and mutations of m6A regulatory factors and the prognosis of patients with soft-tissue sarcomas by using The Cancer Genome Atlas (TCGA) database. Their analysis indicated that CNVs and mutations of KIAA1429, YTHDF3, and IGF2BP1 were independent risk factors predicting OS and DFS. Zhang et al. (2021a) constructed and validated a signature based on m6A-related lncRNAs which could function as independent prognosis-specific predictors in STS, thereby providing novel insights into the roles of m6A-related lncRNAs in STS. Noncoding RNAs (ncRNAs), such as long noncoding RNAs (lncRNAs) and microRNAs (miRNAs), function as key regulators of gene expression, and the competitive endogenous RNA (ceRNA) network is a transcriptional regulatory network at the molecular level, consisting of lncRNA, miRNA, and mRNA (Salmena et al., 2011; Zhu et al., 2019b; Cava et al., 2019). More and more studies have shown that the competing endogenous RNAs (ceRNAs) have been suggested to be involved in essential biological processes and play crucial roles in the initiation and development of neoplasms, and they potentially serve as diagnostic and prognostic markers or therapeutic targets (Qi et al., 2015). At present, the mRNA–miRNA–lncRNA network is mostly studied in the ceRNA network. Recent studies have shown novel roles of the ceRNA network in thyroid cancer, colon adenocarcinoma, gastric cancer, hepatocellular carcinoma, and lung adenocarcinoma (Chang et al., 2020; Hu et al., 2020; Nie et al., 2020; Yang et al., 2020; Zhang et al., 2021b). In addition, a previous study has constructed a ceRNA network and predicted the prognosis of soft-tissue sarcoma recurrence (Huang et al., 2019). Gao et al. (2021) identified differentially expressed mRNAs (DEGs), lncRNAs (DELs), and miRNAs (DEMs) in sarcomas by comparing the gene expression profiles between sarcoma and normal muscle samples in Gene Expression Omnibus (GEO) datasets. Through target gene prediction, a lncRNA–miRNA–mRNA–ceRNA network that contained 113 mRNAs, 69 lncRNAs, and 29 miRNAs was constructed, which might provide insights into further research on the molecular mechanism and potential prognosis biomarkers. Yang et al. (2022) applied starBase and Cytoscape to construct a competing endogenous RNA (ceRNA) network based on m6A-related prognostic lncRNA signature and revealed the prognostic role of m6A-related lncRNAs in osteosarcoma and identified them as potential biomarkers for predicting the prognosis of patients with osteosarcoma. Although most studies have established prognostic models for patients with sarcoma based on the relationship between m6A and noncoding RNA, the biomarkers included in the models are relatively single, and the starting point and perspective of each study are different. So far as we know, the m6A-related ceRNA network and gene signature including three types of potential biomarkers of m6A regulator- and m6A-related noncoding RNAs (lncRNAs and miRNAs) in the regulation of soft-tissue sarcoma (STS) development and progression and prognostic values are largely unexplored. This study is based on the direction to launch a series of explorations on sarcoma.
In this study, based on transcript, somatic mutation, and clinical data obtained from The Cancer Genome Atlas (TCGA) and cBioPortal databases, we conducted extensive analysis. First, the correlation between the expression of 21 widely reported key m6A RNA methylation regulators and the prognosis of SARC was analyzed. Second, a regulatory network of lncRNA (predicted by the LncBase v.2 database)–miRNA (predicted by miRTarBase)–m6A regulators, including 16 miRNAs, 11 m6A regulators, and 104 lncRNAs, was constructed. In addition, using LASSO Cox regression analysis, a risk score model based on miRNA and lncRNA and m6A regulators in the ceRNA network was established to predict the prognosis, immune landscape, and chemosensitivity of SARC patients. Then, the risk score model was evaluated through the ROC curve. These may be potential biomarkers or therapeutic targets of SARC in the future. To further explore the potential relationship between m6Ascore and clinicopathological data, we developed a clinical m6Ascore nomogram to predict the prognosis of sarcoma patients.
Materials and methods
Acquisition of information of patients with SARC
All datasets used in this study were publicly available, and the detailed workflow for risk model construction and subsequent analyses is shown in Figure 1. The RNA sequencing (RNA-seq) transcriptome information (including mRNA, lncRNA, and miRNA expression data), patient clinical information, and somatic mutation status data were obtained from The Cancer Genome Atlas (TCGA)(query) (https://portal.gdc.cancer.gov/) and cBioPortal websites (http://www.cbioportal.org/). Patients with unclear survival time, survival status, and clinicopathological characteristics were excluded. Only patients with overall survival (OS) times of more than or equal to 30 days were included in the dataset.
FIGURE 1. Flow chart of this study.
Selection of m6A genes and m6A-related lncRNAs and miRNAs
According to other published studies, 21 m6A genes were investigated in this study, including writers (METTL3, METTL14, CBLL1, VIRMA [KIAA1429], RBM15, RBM15B, ZC3H13, and WTAP), readers (YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, IGF2BP1, HNRNPA2B1, HNRNPC, FMR1, LRPPRC, and ELAVL1), and erasers (ALKBH5 and FTO). Using the STRING (version 11.0, http://string-db.org/) database, we retrieved the interactions of the 21 m6A RNA methylation regulators and visualized the interactions. The Pearson correlation analysis was used to elucidate the correlation between different m6A RNA methylation regulators. We downloaded the profiles of lncRNAs, miRNA, and m6A genes from TCGA database. Thereafter, the Pearson correlation analysis was conducted to screen the m6A gene-related lncRNAs and miRNAs in SARC samples. LncRNAs with correlation coefficients >0.4 and p < 0.001 were regarded as m6A-related lncRNAs, while miRNAs with correlation coefficients <–0.15 and p < 0.05 were regarded as m6A-related miRNAs.
Construction of the m6A-associated ceRNA network
All miRNAs that were negatively correlated with the 21 selected m6A regulators (METTL3, METTL14, CBLL1, VIRMA [KIAA1429], RBM15, RBM15B, ZC3H13, WTAP, YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, IGF2BP1, HNRNPA2B1, HNRNPC, FMR1, LRPPRC, ELAVL1, ALKBH5, and FTO) were obtained from TCGA-SARC dataset (Pearson correlation coefficient < –0.15 and p < 0.05). Then, the related miRNAs possessing potential binding sites with the 21 selected m6A regulators were simultaneously predicted using the miRTarBase database (http://mirtarbase.cuhk.edu.cn/), and the intersection with the results from the analysis in the previous step was collected. Meanwhile, all lncRNAs were obtained from TCGA-SARC dataset and further screened according to their Pearson correlation coefficients (Pearson correlation coefficient >0.4, p < 0.05) with 21 selected m6A regulators. Next, the related lncRNAs possessing potential binding sites with the interactive m6A-related miRNAs were predicted using the LncBase v.2’s experimental module (http://carolina.imis.athena-innovation.gr/), and the intersection with the results from the analysis in the previous step was obtained. Finally, the lncRNA–miRNA–m6A regulator ceRNA network based on the foundation of the interactions between m6A-related lncRNA and m6A-related miRNA and between m6A-related miRNA and m6A regulators was created and visualized using the “ggalluvial” R package.
Establishment of the m6A-associated risk model
First, by means of the “survival” package in R, univariate Cox regression analysis was utilized to explore the correlation between the genes in the ceRNA network and overall survival (OS) of SARC patients to determine the prognostic-related biomarkers. Then, we used the “glmnet” R package to perform LASSO regression analysis on the genes screened in the univariate Cox regression analysis. Finally, the genes that can be used as independent prognostic factors of OS were screened using multivariate Cox regression analysis, and their regression coefficients (β) were calculated. The following formula was used to calculate the prognostic risk score of each patient:
Risk score = β (mRNA1) * expr (mRNA1) + β (mRNA2) * expr (mRNA2) + … + β (mRNAn) * expr (mRNAn) + β (miRNA1) * expr (miRNA1) + β (miRNA2) * expr (miRNA2) + … + β (miRNAn) * expr (miRNAn) + β (lncRNA1) * expr (lncRNA1) + β (lncRNA2) * expr (lncRNA2) + … + β (lncRNAn) * expr (lncRNAn).
At the same time, the cut-off point of risk score was picked using the “survminer” R package which divided patients into high- and low-risk groups. Subsequently, the Kaplan–Meier survival analysis and log-rank test were used to evaluate the difference in the overall survival (OS) between the risk score groups. Finally, to reflect the prediction ability of the risk score model, we generated the area under the curve (AUC) of the time-dependent receiver operating characteristic (tROC) curves (“riskRegression” package in R) and calculated the area under the curve (AUC) for 1-year, 3-year, and 5-year overall survival (OS).
Analysis of tumor immune signatures and function enrichment for m6Ascore
Tumor immune signatures were evaluated in two aspects: 1) immune checkpoints and 2) the levels of infiltrating immune and immune and stromal scores were calculated using CIBERSORT, TIMER, ssGSEA, and ESTIMATE algorithms. Parts of the results are available at the Genomic Data Commons (GDC, https://gdc.cancer.gov/) and Tumor IMmune Estimation Resource (TIMER2.0, http://timer.cistrome.org/). The gene set enrichment analysis (GSEA) was utilized to understand the biological processes involved in the high- and low-risk groups. Hallmarks in GSEA were used to identify predefined gene sets. A pathway with a |normalized enrichment score (NES)| >1.5, a p-value < 0.05, and a false discovery rate (FDR) < 0.05 was considered to be significant, as described in the Results section.
Analysis of the tumor mutation status in the low- and high-risk groups
The information of somatic mutations in TCGA samples was downloaded from cBioPortal database and TCGA database. The tumor mutational burden (TMB) is defined as the total number of gene mutations per million bases, including the total number of gene-coding errors, base substitutions, and gene insertions or deletions. The TMB value of each patient was calculated using the Perl programming language. Significantly mutated genes between the low- and high-m6Ascore groups and the interaction effect of gene mutations were analyzed using “maftools” R packages. The statistical test for the proportion of mutation was evaluated by the two-side Chi-squared test, and p < 0.05 was considered to be significant.
Prediction of therapeutic sensitivity in patients with different m6Ascore
We studied the predictive capacity of m6Ascore in responding immunotherapy and 138 drugs of chemotherapies/targeted therapies. Based on the public pharmacogenomics database, Genomics of Drug Sensitivity in Cancer (GDSC, https://www.cancerrxgene.org), the 50% inhibiting concentration (IC50) values of the 138 drugs were calculated using the “pRRophetic” R packages. The potential response of patients to immunotherapy was inferred by the tumor immune dysfunction and exclusion (TIDE) algorithm. Generally, a lower TIDE score predicts a better response to immunotherapy. The results of TIDE module analysis of patients with SARC from TCGA dataset were downloaded from the TIDE website (http://tide.dfci.harvard.edu/).
Construction of a predictive nomogram
First, we developed univariate and multivariate Cox regression analyses of the m6Ascore risk signature and other clinicopathological characteristics to confirm the independence of the m6Ascore risk signature. Then, we used the aforementioned factors to establish a nomogram using the “rms” and “regplot” R packages to predict the prognosis of patients with sarcoma. Finally, ROC, C-index, calibration curve analysis, and DCA curves were used to determine whether our established nomogram was suitable for clinical use.
Statistical analysis
The continuous data are expressed as the mean (standard deviation, SD). The categorical data are expressed as frequency and percentage. The relativity between m6Ascore risk signature and immune checkpoint molecules, and TMB were analyzed using the Spearman or Pearson correlation analysis. The Chi-squared test was used to compare different proportions. The differences in proportions of the immune-infiltrating cells, immune checkpoint gene expression, TMB, IC50, and TIDE-related signatures between high- and low-risk groups were compared using the Wilcoxon rank-sum test (Mann–Whitney U test). The Kaplan–Meier survival analysis and log-rank tests were used to analyze the differences in OS between different risk score groups. The univariate and multivariate Cox regression analyses were performed to screen the independent predictors for OS. In all statistical results, except for the special instructions, a two-tailed p-value less than 0.05 indicated statistical significance. All analyses were performed using R software (version 3.6.2, https://www.r-project.org/).
Results
Prognostic analysis of m6A-associated genes in SARC
A total of 255 patients with primary sarcoma had completed clinical data, including 117 (45.9%) males and 138 (54.1%) females, with the mean age being 60.7 ± 14.7. The detailed demographic and clinicopathological data of these SARC patients are shown in Table 1. With the STRING database used to understand the interaction between the 21 m6A RNA methylation regulators, a PPI network was obtained (Figure 2A), which indicated that all enrolled 21 studied genes exhibited gene–gene interactions. In addition, we examined the relationship between 21 regulators through the Pearson correlation test and observed that all m6A RNA methylation regulators were generally positively correlated (Figure 2B). YTHDC1 showed a positive correlation with ZC3H13 (correlation coefficient: 0.8) and METTL14 (correlation coefficient: 0.8). YTHDF3 showed a positive correlation with KIAA1429 (correlation coefficient: 0.8). HNRNPA2B1 showed a positive correlation with ELAVL1 (correlation coefficient: 0.8). Next, according to the aforementioned criteria, we analyzed the prognostic role of the 21 m6A-associated genes in SARC in a total of 245 patients with OS and RNA-seq data. We performed univariate and multivariate Cox regression analyses to explore whether these genes were associated with the prognosis of SARC patients. The results of univariate Cox regression revealed that METTL3 (p = 0.032) and CBLL1 (p = 0.048) were risk genes for SARC (Figure 3A). The multivariate Cox regression results revealed that both CBLL1 (p = 0.033) and RBM15 (p = 0.023) were risk factors for OS. It also revealed that ALKBH5 (HR (hazard ratio) = 0.669, p = 0.021, and 95% CI (confidence interval) [0.475–0.942]) may be a protective gene for OS (Figure 3B).
TABLE 1. Main clinicopathological characteristics of the 255 SARC patients.
FIGURE 2. Interactions and correlations among m6A RNA methylation regulators in SARC. (A) PPI network was constructed to evaluate the interactions among m6A RNA methylation regulators; and (B) Pearson correlation analysis was used to analyze the correlations among m6A RNA methylation regulators.
FIGURE 3. Correlation between the expression levels of m6A-related genes and overall survival (OS) rates in SARC patients (N = 245). (A) Univariate analysis of m6A-related genes associated with OS; (B) and multivariate analysis of m6A-related genes associated with OS.
Identification of key upstream miRNAs of 21 m6A genes
The potential upstream miRNAs of 21 m6A genes were predicted using the experimentally verified microRNA–target gene interaction database miRTarBase. Ultimately, it found that 720 miRNAs interacted with 20 m6A genes (writers: METTL3, METTL14, CBLL1, KIAA1429 [VIRMA], RBM15B, ZC3H13, and WTAP; readers: YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, IGF2BP1, HNRNPA2B1, HNRNPC, FMR1, LRPPRC, and ELAVL; and erasers: ALKBH5 and FTO), while no upstream miRNA interacted with one m6A gene (RBM15) in the database (Supplementary Table S1). Next, the matrix expressions of 21 m6A genes and 2,064 miRNAs were abstracted from TCGA database. We defined the miRNAs that were significantly related to one of the 21 m6A genes (Pearson R < -0.15 and p < 0.05) as m6A-related miRNAs. Finally, 144 miRNAs were discerned as m6A-related miRNAs (Supplementary Table S2). Also, the predicted miRNAs were then intersected with these m6A-related miRNAs to select the 27 m6A-related miRNAs (key miRNAs) that interacted with 12 m6A genes. Cytoscape was used to construct an miRNA–mRNA regulatory network, including 30 miRNA–mRNA pairs (Figure 4).
FIGURE 4. MiRNA–mRNA interaction network of m6A genes.
Identification of upstream lncRNAs of key miRNAs in patients with SARC
Based on the ceRNA hypothesis that lncRNAs can compete with mRNAs for miRNA binding, thus playing a vital role in the process of tumor development and progression (Nie et al., 2020), this study further used the experimental module of LncBase database (http://carolina.imis.athena-innovation.gr/diana_tools/web/index.php?r = Lncbasev2%2Findex-experimental), an online database, to predict the upstream lncRNAs that may interact with the 27 key miRNAs (hsa-miR-124-3p, hsa-miR-455-3p, hsa-miR-34a-5p, hsa-miR-214-5p, hsa-miR-124-5p, hsa-miR-302a-5p, hsa-miR-3202, hsa-miR-125a-3p, hsa-miR-224-3p, hsa-miR-4999-5p, hsa-miR-6766-5p, hsa-miR-6808-3p, hsa-miR-3133, hsa-miR-6885-3p, hsa-miR-1298-3p, hsa-miR-5692c, hsa-miR-378f, hsa-miR-378b, hsa-miR-5693, hsa-miR-3612, hsa-miR-6515-3p, hsa-miR-361-5p, hsa-miR-5008-3p, hsa-miR-6865-3p, hsa-miR-193b-3p, hsa-miR-23c, and hsa-miR-196b-5p), and a total of 1,552 lncRNAs were obtained (Supplementary Table S3). Then, the matrix expressions of 21 m6A genes and 14,789 lncRNAs were abstracted from TCGA database. We defined the lncRNAs that were significantly related to one of the 21 m6A genes (Pearson R > 0.4 and p < 0.001) as m6A-related lncRNAs. Finally, 1,028 lncRNAs were discerned as the m6A-related lncRNAs (Supplementary Table S4). Also, the predicted lncRNAs were then intersected with these m6A-related lncRNAs to select the m6A-related lncRNAs (n = 104) that interacted with 16 m6A-related miRNAs. Cytoscape was used to construct an lncRNA–miRNA regulatory network, including 189 lncRNA–miRNA pairs (Figure 5).
FIGURE 5. MiRNA–lncRNA regulatory network associated with m6A genes in SARC.
Construction of an m6A-related competing endogenous RNA (ceRNA) network
Based on the results mentioned previously, the miRNA–mRNA and lncRNA–miRNA were selected to establish an m6A-associated ceRNA network that contained 104 lncRNAs, 16 miRNAs, and 11 mRNAs (Figure 6). We searched the lncRNAs and miRNAs in our study in the m6A-Atlas database (Tang et al., 2021), RMBase v2.0 (Xuan et al., 2018) and WHISTLE (Chen et al., 2019b) database, which were used to predict m6A methylation sites in these RNAs. The search results (Supplementary Tables S6, S7) showed that most lncRNAs in this study had potential methylation sites. However, we did not retrieve the potential m6A methylation sites of miRNAs in these libraries. Detailed information of the ceRNA network is listed in Supplementary Table S5.
FIGURE 6. Construction of a network of lncRNA–miRNA–m6A regulators.
Construction of a risk model according to m6A-related RNAs in SARC patients
A total of 244 patient samples were included in the whole dataset. A total of 16 miRNAs, 104 lncRNAs, and 11 m6A regulators (mRNAs) in the ceRNA network were selected as candidate biomarkers for the following step analysis. The result of univariate survival analysis showed that five m6A-related lncRNAs (RP11-46C24.7, RP11-283I3.6, SLC25A21-AS1, RP11-81A1.6, and RP11-346C20.4), two m6A-related miRNA (hsa-miR-455-3p and hsa-miR-124-3p), and one m6A regulator (CBLL1) were associated with SARC prognosis (p < 0.1). Among these biomarkers, hsa-miR-455-3p (miRNA), hsa-miR-124-3p (miRNA), and CBLL1 (mRNA) were considered as risk biomarkers (HR > 1). Subsequently, to obtain the most useful predictive features, LASSO Cox regression analysis was performed on eight genes (Figures 7A,B). Based on the final Cox regression model results (Table 2), we selected 3 genes (hsa-miR-455-3p, CBLL1, and RP11-283I3.6) with p < 0.1 to construct the risk score model. Furthermore, we calculated the risk score (m6Ascore) of each sample based on the risk score model, and the formula was shown as follows:
TABLE 2. Cox proportional hazard regression analysis and LASSO of eight genes.
FIGURE 7. Risk model from m6A-related genes. (A-B) LASSO Cox regression analysis of eight m6A-related genes; (C) overall survival analysis for patients in high/low-risk groups; (D) distributions of risk score, different patterns of survival status, and survival time between the high- and low-risk groups; (E-G) and time-dependent ROC curves at 1, 3, and 5 years based on the m6A-related gene signature.
Risk score = (0.09964) * hsa-miR-455-3p + (0.84723) * CBLL1 + (−0.35103) * RP11-283I3.6.
Patients with SARC were divided into the low-risk or high-risk groups with the optimal cutoff (2.492404) of the risk score. The Kaplan–Meier (KM) curve analysis result showed that the low-risk group had a better prognosis than the high-risk group (p = 0.00038) (Figure 7C). The distribution of risk score between low-risk and high-risk groups is depicted in Figure 7D, and a survival status map was plotted to demonstrate the status for each sample using TCGA sarcoma dataset values. A time-dependent ROC curve was performed to evaluate the sensitivity and specificity of the risk score. The m6A-related signature’s AUC values were 0.638 (0.515, 0.761), 0.663 (0.578, 0.749), and 0.679 (0.585, 0.773), respectively, for an OS of 1, 3, and 5 years (Figures 7E–G). This indicates that our risk score model could be used to predict SARC patient survival.
m6Ascore was associated with the SARC immune landscape
Our study explored the relationship between m6Ascore and the biological process of SARC, for which we conducted GSEA. GSEA using TCGA data of the hallmark gene sets indicated that in the two cohorts, DNA repair, E2F targets, G2M checkpoint, mitotic spindle, mTORC1 signaling, MYC targets v1, MYC targets v2, oxidative phosphorylation, protein secretion, and unfold protein response were significantly enriched in the high-risk group (the top ten pathways are shown in Figure 8A), while the angiogenesis was significantly enriched in the low-risk group (the one pathway is shown in Figure 8A; red line). This finding provided insights into the potential biological processes and signaling pathways modulated by m6A-related RNAs in SARC.
FIGURE 8. Function enrichment analysis for m6Ascore and immune landscape of different m6Ascore subgroups. (A) GSEA of samples with high or low m6Ascore. Top ten significant pathways associated with high m6Ascore (p < 0.05 and FDR-adjusted q < 0.05, NES > 1.5). Significant pathways in the red module are associated with low m6Ascore (nominal p < 0.05 and FDR adjusted q < 0.05, NES < 1.5); (B) comparison of immune checkpoint gene expression levels between the low-m6Ascore group and the high-m6Ascore group. *p < 0.05; **p < 0.01; ***p < 0.001; (C) correlation chord chart shows the mutual correlation between m6Ascore and several prominent immune-checkpoint-relevant genes (CD44, VTCN1, TMIGD2, TNFRSF18, TNFRSF25, TNFRSF4, TNFRSF8, and VTCN1); (D) heatmap for immune responses based on CIBERSORT, ESTIMATE, ssGSEA, and TIMER algorithms among the high- and low-risk group; and (E) beeswarm plot for infiltration levels of immune cells in high- and low-risk samples through the CIBERSORT. The statistical difference between the two groups was compared by the Wilcoxon rank-sum test. *p < 0.05; **p < 0.01; ***p < 0.001.
The expression of immune checkpoints was used to predict immunotherapeutic benefits in multiple malignancies. Next, the study was aimed to explore whether m6Ascore could predict immunotherapeutic benefits in SARC patients. We analyzed the correlation between the high- and low-m6Ascore groups and 50 immune checkpoints. The expression of CD44 and VTCN1 was increased in the high-m6Ascore group, while the expression of TMIGD2, TNFRSF18, TNFRSF25, TNFRSF4, TNFRSF8, and NRP1 was increased in the low-m6Ascore group (Figure 8B). A chord chart was used to display the correlation between m6Ascore and immune checkpoint molecules. The results showed that m6Ascore was negatively correlated with the expression of TMIGD2, TNFRSF18, TNFRSF25, TNFRSF4, TNFRSF8, and NRP1 (Figure 8C), indicating that the poor prognosis of high-m6Ascore patients might be due to the tumor immunosuppressive microenvironment. We also investigated the differences in the distribution of infiltrating immune cells between the high-risk and low-risk groups using CIBERSORT, ESTIMATE, ssGSEA, and TIMER algorithms. ESTIMATE score and immune cell types, which were differentially infiltrated between the low- and high-risk groups, are presented in Figure 8D. To analyze the composition of immune cells in different m6Ascore subgroups, we used the Wilcoxon rank-sum test to compare the distribution of immune cells in different m6Ascore subgroups. The results showed that only in CIBERSORT algorithm, the abundance of immune cells including macrophages M0 (p < 0.05) and T-cell CD4 memory activated (p < 0.05) was highly infiltrated in the high-m6Ascore subgroup, while monocytes (p < 0.05) and mast cells activated (p < 0.05) were more abundant in the low-m6Ascore subgroup (Figure 8E). Totally, these findings suggest that there are significant associations between m6Ascore and the tumor immune landscape in SARC.
Mutation status in SARC patients in the high- and low-m6Ascore groups
To investigate m6Ascore-related mechanisms in SARC, somatic mutations data were also analyzed. The frequency of mutations in top 20 genes between the two groups is shown in Figure 9A. A significant mutually exclusive phenomenon was observed among mutations of these genes (Figure 9B). Subsequently, differentially mutated genes between the two groups were detected, and the “maftools” package analysis result showed that CSMD1 was the only one significant differentially mutated gene detected between the high- and low-m6Ascore cohorts (Figure 9C). The CSMD1 mutation burden in the low-m6Ascore subtype was significantly higher than that in the high-m6Ascore subtype.
FIGURE 9. Mutation patterns of SARC patients. (A) Mutational patterns of 169 SARC patients in the low- and high-m6Ascore groups displayed by the oncoplot; (B) interaction effect of genes mutating differentially in patients in the low- and high-m6Ascore groups; and (C) forest plot of genes mutating differentially in patients of the low- and high-m6Ascore groups.
We further found that there was no significant correlation between m6Ascore and TMB (Figure 10A). Also, there was no significant difference in TMB between the patients with high m6Ascore and those with low m6Ascore (Figure 10B). However, we found that low TMB was associated with good OS (Figure 10C, log–rank test, p = 0.0045). Subsequently, we explored whether the combination of m6Ascore and TMB could be a more powerful predictive biomarker for prognosis. We integrated m6Ascore and TMB to stratify all the samples into the high-TMB/low-m6Ascore, low-TMB/low-m6Ascore, high-TMB/high-m6Ascore, and low-TMB/high-m6Ascore groups. As shown in Figure 10D, significant differences were found among all groups (log-rank test, p < 0.0001), and patients in the high-TMB/high-m6Ascore group exhibited poor OS. These results together strongly demonstrated that the risk score was positively correlated with tumor malignancy. Next, the “maftools” R package was used to analyze and summarize the mutation data. The top 20 driver genes with the highest alteration frequency between the aforementioned subgroups are shown in Figure 10E.
FIGURE 10. Relationship of the m6Ascore signature with TMB. (A) correlation between the m6Ascore signature and TMB depicted by scatter plots; (B) comparison of TMB between the high- and low-m6Ascore groups; (C) KM survival curves for the high- and low-TMB groups stratified at the optimal cutoff in TCGA-SARC cohorts (log-rank test, p = 0.0045); (D) Kaplan–Meier survival analysis for four groups stratified by combining the TMB and the m6Ascore signature in TCGA-SARC cohort; and (E) waterfall plot of tumor somatic mutation established by stratifying with m6Ascore and TMB. Each column represents an individual patient. Group1, high TMB /high m6Ascore, green; Group2, high TMB/low m6Ascore, red; Group3, low TMB/high m6Ascore, purple; and Group4, low TMB/low m6Ascore, yellow. The mutational types include frame shift del, frame shift ins, in-frame del, in-frame ins, missense mutation, multi-hit, nonsense mutation, and splice site.
m6Ascore prediction of response to chemotherapy and immunotherapy
To find the potency of m6Ascore as a biomarker for predicting the response of SARC patients to drugs (including chemotherapy, targeted therapy, and immunotherapy), we used the “pRRophetic” algorithm to estimate the therapeutic response based on the half-maximal inhibitory concentration (IC50) available in the Genomics of Drug Sensitivity in Cancer (GDSC) database for each sample. We inferred the IC50 values of the 138 drugs in TCGA-SARC patients. Finally, we found that 69 compounds were screened out for significant differences in the estimated IC50 values between the two groups, and the high-risk group was more sensitive to these compounds including ATRA, cyclopamine, JNK inhibitor, PD173074, and QS11,while the low-risk group was more sensitive to the remaining compounds (Supplementary Figures S1–S3). Figures 11A-N display the top 14 compounds that might be used for further analysis in patients with SARC. In terms of response to immunotherapy, we used the TIDE algorithm to assess the potential clinical efficacy of immunotherapy in different m6Ascore subgroups. The TIDE algorithm assessed expression signatures of T-cell dysfunction and T-cell exclusion to assess tumor immune evasion and integrated them into the TIDE total score. In addition, the TIDE module was used to analyze multiple signatures to estimate tumor immune evasion, such as MDSC, M2 TAM, or CAF signatures. In our results (Figure 11O), the high-m6Ascore subgroup had a relatively higher T-cell exclusion score (p < 0.05), but there was no difference in the T-cell dysfunction scores between the two subgroups (p > 0.05). For other features produced by TIDE, although there was no significant difference in M2 TAM and CAF features between the two groups, the low-risk group had a relatively lower trend, and a lower proportion of MDSC was associated with the low-risk group (p < 0.05). These results suggest that SARC patients in the high-risk group may have a higher potential for the immunosuppressive TME status and be less responsive to ICI therapy. These findings identified the promising role of this risk signature as a predictor for chemotherapy and immunotherapy efficacy in the treatment of SARC patients.
FIGURE 11. Correlation between m6Ascore and therapy. (A–N) IC50 of some chemotherapeutic drugs are in the high- and low-risk patients. *p < 0.05; **p < 0.01; ***p < 0.001; (O) TIDE prediction difference in the high- and low-risk patients. *p < 0.05; **p < 0.01; ***p < 0.001.
Construction and evaluation of the nomogram predicting OS based on the m6A-related signature
To confirm whether the m6A-related signature for OS was an independent prognostic factor, univariate and multivariate Cox regression analyses were performed (Table 3). As the results showed, in the univariate Cox regression analysis, m6Ascore, age, and cancer status were significantly associated with the OS of sarcoma patients. Then, m6Ascore, age, and cancer status were identified as independent prognostic factors of sarcomas via multivariate Cox regression analysis. All these independent factors were combined to establish a nomogram for predicting the 1-, 3-, and 5-year OS (Figure 12A). As shown in Figure 12A, m6Ascore contributes more to the total score than other variables. The 1-, 3-, and 5-year OS rates of patients declined as the total score increased. The C-index of the nomogram model reached 0.744 (95% CI: 0.707–0.784). The calibration plots showed that the nomogram model predicted the overall survival of patients with SARC well (Figure 12B). We compared the clinical net benefit of the nomogram through DCA curves. The nomogram demonstrated a larger net benefit within most of the threshold probabilities (Figures 12C–E), indicating that the nomogram had high potential clinical utility for predicting prognosis in patients with SARC. Finally, Figures 12F–H show the predictive potential of the nomogram using time-dependent ROC curves. The area under the ROC curve (AUC) of the nomogram model for OS was 0.693 at 1 year, 0.772 at 3 years, and 0.834 at 5 years.
TABLE 3. Univariate and multivariate Cox regression analysis of the overall survival (OS) of patients with SARC.
FIGURE 12. Construction and validation of the nomogram. (A) Nomogram based on the multivariate Cox regression model of SARC patients; (B) calibration plots for the internal validation of the current nomogram; the x-axis represents the nomogram-predicted overall survival, and the y-axis represents the actual overall survival of patients with SARC; (C–E) DCA of the nomogram based on m6Ascore for 1-, 3- and 5-year OS prediction; and (F–H) time-dependent ROC curves of nomogram, m6Ascore, cancer status, and age at 1, 3, and 5 years.
Discussion
In the present study, the prognostic significance of 21 m6A RNA methylation regulators in SARC was first explored. The study showed that only CBLL1 was found to be an independent risk factor for OS by univariate and multivariate Cox regression analyses. Numerous studies have shown that m6A regulator-related signatures and m6A regulator-related miRNA or lncRNA signatures may serve as prognostic biomarkers in patients with various types of cancer (Li et al., 2021c; Jin et al., 2021; Xu et al., 2021). However, the expression and functional roles of m6A RNA methylation regulators and their associated regulatory networks in the occurrence and progression of SARC have not been widely discussed. Therefore, in this study, bioinformatics was used to analyze the existing sequencing datasets of integrated sarcoma, and miRNAs and lncRNAs related to m6A regulators were obtained through correlation analysis, and the upstream miRNAs and lncRNAs of m6A regulators were determined using miRTarBase and LncBase v.2 databases, respectively. The two parts of the results of miRNAs and lncRNAs were intersected to key miRNAs and lncRNAs, respectively. Finally, a lncRNA–miRNA–m6A regulator ceRNA network was constructed, which contained 104 lncRNAs, 16 miRNAs, and 11 m6A regulators. We should pay more attention to screening out key lncRNAs, miRNAs, and m6A regulators that can predict OS. Therefore, considering all these genes in the ceRNA network, we performed univariate Cox regression analysis to identify genes associated with clinical prognosis in sarcoma patients. Eight genes were found to be significantly associated with clinical outcomes in sarcoma. Finally, LASSO regression analysis and multivariate Cox regression analysis were performed, and three genes (RP11-283I3.6, hsa-miR-455-3p, and CBLL1) were identified as prognostic biomarkers for SARC, and the three genes were included in the risk scoring model for predicting OS for SARC.
By searching for these genes in the PubMed database, it is found that the mechanisms of hsa-miR-455-3p and CBLL1 in tumor progression or studies related to tumor progression have been reported. The aberrant expression of hsa-miR-455-3p and its prognostic value have been widely reported in various types of human cancers. In our study, hsa-miR-455-3p served as a risk biomarker for SARC, consistent with the results of bioinformatics analysis by Wang et al. (2019), which showed that the differential expression level of miR-455-3p was the most significant in gliomas. Subsequently, its expression in glioma patients was examined. Consistent with the results of bioinformatics analysis, the expression level of miR-455-3p was significantly upregulated in glioma tissues compared with normal tissues. Cox regression analysis further identified miR-455-3p as an independent prognostic indicator of overall survival in glioma patients. hsa-miR-455-3p is overexpressed in skin basal cell carcinoma (BCC) (Sand et al., 2012). However, the expression of miR-455-3p is obviously decreased in the tissues and cells of hepatocellular carcinoma (HCC). This miRNA can impair HCC cell malignancy via suppression of insulin growth factor receptor expression, thereby disrupting glycolysis (Hu et al., 2019; Jiang et al., 2019). Gao et al. (2018), Yang et al. (2017), and Shang et al. (2021) reported that low expression of miR-455-3p in non-small-cell lung cancer, esophageal squamous cell carcinoma (ESCC), and pancreatic cancer (PAAD) tissues was strongly associated with poor prognosis. miR-455-3p acts as a tumor suppressor in esophageal squamous cell carcinoma (ESCC) and inhibits cell proliferation and invasion by targeting FAM83F. miR-455-3p is involved in increasing the expression of HOXC4, promoting transcriptional disorders in cancer. The miR-455-3p–HOXC4 axis is expected to be closely related to the metastasis and prognosis of human pancreatic cancer. Yi et al. (2020) and Sun et al. (2020) together showed that the expression of miR-455-3p was decreased in osteosarcoma tissues and cell lines. Patients with high miR-455-3p expression had satisfactory survival rates. miR-455-3p is a potential clinical therapeutic target and prognostic biomarker inhibiting proliferation, migration, and invasion and enhancing apoptosis. However, Hisaoka et al. (2011) showed that, compared to normal skeletal muscle, hsa-miR-455-3p was significantly upregulated in synovial sarcomas, suggesting that the molecule has a potential oncogenic role, which calls for further investigation to develop a better understanding of the oncogenic mechanisms. Our finding is consistent with this study. CBLL1 plays an important role in tumorigenesis. Hui et al. (2019) found that CBLL1 was upregulated in non-small-cell lung cancer (NSCLC) tissues compared to the adjacent nontumor tissues, and the high expression of CBLL1 was associated with the tumor size in NSCLC tissues. Their results confirmed that CBLL1 promoted the proliferation by promoting G1/S cell cycle transition in NSCLC cells. Moreover, CBLL1 knockdown inhibited cell invasion via increased E-cadherin protein expression and decreased expression of MMP2 and MMP9 in NSCLC cell lines. Previous studies found that CBLL1, an E3 ubiquitin ligase, inhibits ER pathway activity by binding to an ER co-activator and then further inhibits the proliferation and differentiation of BC cells (Makdissi et al., 2009). Zheng et al. (2021) confirmed that higher CBLL1 expression was associated with a better prognosis in BC than lower CBLL1 expression. Functional analysis showed that CBLL1 was related to the ESR1-related pathway, apoptosis-related pathway, cell cycle pathway and immune-related pathway in BC. Although there is no report on the correlation between the expression of RP11-283I3.6 and CBLL1 and SARC progression, based on our results, we speculate that they may be potential biomarkers for sarcoma prognosis. The mechanism of action of these genes in sarcoma is unclear and requires further follow-up studies. Subsequently, based on this risk score feature, we divided all SARC patients into the high-risk and low-risk groups, and KM survival analysis showed that patients in the low-risk group had significantly higher OS than those in the high-risk group. The time-dependent AUC indicated that the risk scoring model had good predictive performance for OS.
The molecular heterogeneity features between high- and low-risk patients were further analyzed. GSEA showed that DNA repair, E2F targets, G2M checkpoint, mitotic spindle, mTORC1 signaling, MYC targets v1, MYC targets v2, oxidative phosphorylation, protein secretion, and unfold protein response were significantly enriched in the high-risk specimens. Meanwhile, activation of angiogenesis was detected in low-risk specimens. The tumor immune microenvironment comprising stromal cells and immune cells correlates with immunotherapy response (Zhou et al., 2020). Components of the immune microenvironment are key determinants of prognosis and response to immunotherapy (Das et al., 2020). Immunotherapy is an emerging new approach to treating a variety of cancers, including sarcomas. Exploring which patient can benefit from immunotherapy remains a great challenge. Here, the association between immune cell infiltration and this risk score was comprehensively analyzed using the CIBERSORT, ESTIMATE, ssGSEA, and TIMER algorithms in the present study. Compared with patients in the high-risk group, the low-risk patients had higher levels of monocytes and mast cells activated and decreased levels of macrophages M0 and T-cell CD4 memory activated infiltration. The infiltration of M0 macrophages is positively related to poor clinical outcomes in human malignancies, including STS, which is in accordance with the findings of our study (Zhu and Hou, 2020). Thus, the worse clinical outcomes of the high-risk group may be associated with infiltrating immune cellular populations. In addition, expression levels of several immune checkpoints, including TMIGD2, TNFRSF18, TNFRSF25, TNFRSF4, TNFRSF8, and NRP1, were also higher in the low-risk group. These data were indicative of this risk signature being closely related to immunotherapy. m6Ascore may be used as an indicator independent of TMB expression to predict the efficacy of ICB. Tumor mutational burden (TMB) has been identified as a biomarker of immunotherapy response (Hellmann et al., 2018a; Hellmann et al., 2018b), where higher TMB predicts higher benefits from immunotherapy (Hellmann et al., 2018a). However, the prognostic value of TMB varies across cancer types according to a pan-cancer study (Ding et al., 2018). In bladder urothelial carcinoma (BLCA), stomach adenocarcinoma (STAD), and uterine corpus endometrial carcinoma (UCEC), high TMB is associated with longer overall survival (OS). In HNSCC, kidney renal clear cell carcinoma (KIRC), and low-grade glioma (LGG), high TMB is associated with shorter OS. The Cox regression result of TMB in SARC patients showed that HR (95% CI) was 1.25 (0.63–2.49). In the present study, patients with lower TMB had better prognosis than SARC patients with higher TMB, which is consistent with the trend of pan-cancer research. The significance of the results may be that the selected threshold is different. However, we did not observe a significant correlation between m6Ascore and TMB, nor did we observe a significant difference in the distribution of TMB between patients in the two different risk groups. The main reason is that with a threshold of p < 0.05, using Fisher’s exact test, CSMD1 was the only one significant differentially mutated gene detected between the high- and low-m6Ascore cohorts, and CSMD1 was found to be mutating more in the low-m6Ascore group, which resulted in a similar gene mutation status in the two groups. These findings suggested that TMB and m6Ascore are independent biomarkers/indicators for predicting ICB response. Interestingly, the combination of TMB and the risk characteristic we constructed can more clearly and accurately stratify SARC patients. Compared with other subgroups, patients with high TMB/high m6Ascore have the worst prognosis, which also shows that the TMB status does not affect the prognostic value of m6Ascore. m6Ascore has the potential to predict immunotherapy responsiveness, which may be independent of TMB (Zhang et al., 2020; Guan et al., 2021; Lin et al., 2021; Sun et al., 2021; Wu et al., 2021). In addition, the TIDE results suggested that the patients in the low-risk group may respond better to immunotherapy. Based on estimated IC50 values, the patients in the low-risk group showed sensitive chemotherapy responses to most drugs. Taken together, these findings suggest that this risk signature may play a role in the ICB treatment of SARC. Subsequent studies further confirmed that risk characteristic is independent prognostic factors in patients with sarcoma. Based on risk score and other clinically independent predictors, a nomogram for personalized clinical outcome prediction was established in the study, which was certified to perform well for predicting the 1-, 3-, and 5-year survival rates of SARC patients, showing a C-index of 0.744 (95% CI: 0.707–0.784).
Although potential biomarkers involved in tumorigenesis in a large number of samples were identified by the bioinformatics approach, it should be noted that this study also has some limitations as follows: 1) due to the lack of RNA-seq or microarray data in SARC patients, only TCGA data were included. Also, the SARC samples lacked some additional clinical follow-up information; therefore, factors such as the presence of other health conditions in patients to differentiate prognostic biomarkers were not included. 2) We internally verified the nomogram prediction model based on m6Ascore, and the findings of this study would be more meaningful if this model could be well validated externally with another real-world, independent, large-quantity, and high-quality cohort, and thus, a more diverse patient population could be extrapolated. However, the application of the prognostic prediction model based on m6Ascore required four types of data, namely, clinical data, RNA-seq (mRNAs and lncRNAs), and miRNA-seq, which involves high costs and is not easily feasible in practice. 3) Most importantly, experimental validation is needed to confirm these results and further explore the potential mechanism and role of these potential biomarkers in SARC. However, our findings showed that the nomogram model based on m6Ascore may be promising for clinical prediction of prognosis and might contain potential biomarkers for treatment response prediction for SARC patients, which remains an instructive and efficient way for predicting the accurate individual clinical outcomes of SARC patients.
Conclusion
In conclusion, in our study, a ceRNA network based on m6A-related genes was successfully constructed through bioinformatics analysis of TCGA database. Candidate biomarkers in the ceRNA network were used to establish a risk profile of m6A-related RNAs, which is significantly associated with the prognosis and immune microenvironment of SARC, and could effectively predict the prognosis and treatment efficacy of STSs. The results of this study suggest that these markers may play an important role in the therapeutic target and prognostic analysis of sarcoma patients.
Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.
Author contributions
HL and KW designed the research. HL, DL, XW, and JZ conducted the data analysis. HL and DL wrote the manuscript. KW and ZF revised and completed the final manuscript. All authors have read and agreed to the published version of the manuscript.
Funding
This research was funded by the Tianshan Innovative Research Team of Xinjiang Uygur Autonomous Region, China (2020D14020).
Acknowledgments
The authors acknowledge the data support from databases such as TCGA, TIMER, cBioPortal, GDC, miRTarBase, and LncBase.
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.
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.894080/full#supplementary-material
References
Alavi, S. N., Florou, V., Tinoco, G., Trent, J. C., and Wilky, B. A. (2018). A precision medicine approach in sarcoma: Identification of patients who may benefit from early use of pazopanib. Discov. Med. 25, 131–144.
Aung, T., Asam, C., and Haerteis, S. (2019). Ion channels in sarcoma: Pathophysiology and treatment options. Pflugers Arch. 471, 1163–1171. doi:10.1007/s00424-019-02299-8
Bacci, G., Ferrari, S., Mercuri, M., Bertoni, F., Picci, P., Manfrini, M., et al. (1998). Predictive factors for local recurrence in osteosarcoma: 540 patients with extremity tumors followed for minimum 2.5 years after neoadjuvant chemotherapy. Acta Orthop. Scand. 69, 230–236. doi:10.3109/17453679809000921
Barbieri, I., and Kouzarides, T. (2020). Role of RNA modifications in cancer. Nat. Rev. Cancer 20, 303–322. doi:10.1038/s41568-020-0253-2
Bone Cancer (2021). Bone cancer: Statistics, Available at:. https://www.cancer.net/cancer-types/bone-cancer-sarcoma-bone/statistics (Accessed September 25, 2021).
Cava, C., Bertoli, G., and Castiglioni, I. (2019). Portrait of tissue-specific coexpression networks of noncoding RNAs (miRNA and lncRNA) and mRNAs in normal tissues. Comput. Math. Methods Med. 2019, 9029351. doi:10.1155/2019/9029351
Chang, Z., Huang, R., Fu, W., Li, J., Ji, G., Huang, J., et al. (2020). The construction and analysis of ceRNA network and patterns of immune infiltration in Colon adenocarcinoma metastasis. Front. Cell. Dev. Biol. 8, 688. doi:10.3389/fcell.2020.00688
Chen, K., Wei, Z., Zhang, Q., Wu, X., Rong, R., Lu, Z., et al. (2019). Whistle: A high-accuracy map of the human N6-methyladenosine (m6A) epitranscriptome predicted using a machine learning approach. Nucleic Acids Res. 47, e41. doi:10.1093/nar/gkz074
Chen, S., Li, Y., Zhi, S., Ding, Z., Wang, W., Peng, Y., et al. (2020). WTAP promotes osteosarcoma tumorigenesis by repressing HMBOX1 expression in an m6A-dependent manner. Cell. Death Dis. 11, 659. doi:10.1038/s41419-020-02847-6
Chen, X. Y., Zhang, J., and Zhu, J. S. (2019). The role of m6A RNA methylation in human cancer. Mol. Cancer 18, 103. doi:10.1186/s12943-019-1033-z
Das, S., Camphausen, K., and Shankavaram, U. (2020). Cancer specific immune prognostic signature in solid tumors and its relation to immune checkpoint therapies. Cancers 12, 2476. doi:10.3390/cancers12092476
Ding, H., Zhao, J., Zhang, Y., Wang, G., Cai, S., and Qiu, F. (2018). Tumor mutational burden and prognosis across pan-cancers. Ann. Oncol. 29, viii16–viii17. viii16-viii17. doi:10.1093/annonc/mdy269.055
Doyle, L. A. (2014). Sarcoma classification: An update based on the 2013 world health organization classification of tumors of soft tissue and bone. Cancer 120, 1763–1774. doi:10.1002/cncr.28657
Eaton, B. R., Schwarz, R., Vatner, R., Yeh, B., Claude, L., Indelicato, D. J., et al. (2021). Osteosarcoma. Pediatr. Blood Cancer 68, e28352. doi:10.1002/pbc.28352
Ferguson, J. L., and Turner, S. P. (2018). Bone cancer: Diagnosis and treatment principles. Am. Fam. Physician 98, 205–213.
Gao, L., Zhao, Y., Ma, X., and Zhang, L. (2021). Integrated analysis of lncRNA-miRNA-mRNA ceRNA network and the potential prognosis indicators in sarcomas. BMC Med. Genomics 14, 67. doi:10.1186/s12920-021-00918-x
Gao, X., Zhao, H., Diao, C., Wang, X., Xie, Y., Liu, Y., et al. (2018). miR-455-3p serves as prognostic factor and regulates the proliferation and migration of non-small cell lung cancer through targeting HOXB5. Biochem. Biophys. Res. Commun. 495, 1074–1080. doi:10.1016/j.bbrc.2017.11.123
Guan, X., Xu, Z. Y., Chen, R., Qin, J. J., and Cheng, X. D. (2021). Identification of an immune gene-associated prognostic signature and its association with a poor prognosis in gastric cancer patients. Front. Oncol. 10, 629909. doi:10.3389/fonc.2020.629909
Hatina, J., Kripnerova, M., Houfkova, K., Pesta, M., Kuncova, J., Sana, J., et al. (2019). Sarcoma stem cell heterogeneity. Adv. Exp. Med. Biol. 1123, 95–118. doi:10.1007/978-3-030-11096-3_7
Hellmann, M. D., Callahan, M. K., Awad, M. M., Calvo, E., Ascierto, P. A., Atmaca, A., et al. (2018). Tumor mutational burden and efficacy of nivolumab monotherapy and in combination with ipilimumab in small-cell lung cancer. Cancer Cell. 33, 853–861. doi:10.1016/j.ccell.2018.04.001
Hellmann, M. D., Ciuleanu, T. E., Pluzanski, A., Lee, J. S., Otterson, G. A., Audigier-Valette, C., et al. (2018). Nivolumab plus Ipilimumab in lung cancer with a high tumor mutational burden. N. Engl. J. Med. 378, 2093–2104. doi:10.1056/NEJMoa1801946
Hisaoka, M., Matsuyama, A., Nagao, Y., Luan, L., Kuroda, T., Akiyama, H., et al. (2011). Identification of altered MicroRNA expression patterns in synovial sarcoma. Genes. Chromosom. Cancer 50, 137–145. doi:10.1002/gcc.20837
Hou, M., Guo, X., Chen, Y., Cong, L., and Pan, C. (2020). A prognostic molecular signature of N⁶-Methyladenosine methylation regulators for soft-tissue sarcoma from the cancer Genome Atlas database. Med. Sci. Monit. 26, e928400. doi:10.12659/MSM.928400
Hu, H., Xu, H., Lu, F., Zhang, J., Xu, L., Xu, S., et al. (2020). Exploring the effect of differentially expressed long non-coding RNAs driven by copy number variation on competing endogenous RNA network by mining lung adenocarcinoma data. Front. Cell. Dev. Biol. 8, 627436. doi:10.3389/fcell.2020.627436
Hu, Y., Yang, Z., Bao, D., Ni, J. S., and Lou, J. (2019). miR-455-5p suppresses hepatocellular carcinoma cell growth and invasion via IGF-1R/AKT/GLUT1 pathway by targeting IGF-1R. Pathol. Res. Pract. 215, 152674. doi:10.1016/j.prp.2019.152674
Huang, R., Meng, T., Chen, R., Yan, P., Zhang, J., Hu, P., et al. (2019). The construction and analysis of tumor-infiltrating immune cell and ceRNA networks in recurrent soft tissue sarcoma. Aging (Albany NY) 11, 10116–10143. doi:10.18632/aging.102424
Hui, J. Y. (2016). Epidemiology and etiology of sarcomas. Surg. Clin. North Am. 96, 901–914. doi:10.1016/j.suc.2016.05.005
Hui, L., Zhang, S., Wudu, M., Ren, H., Xu, Y., Zhang, Q., et al. (2019). CBLL1 is highly expressed in non-small cell lung cancer and promotes cell proliferation and invasion. Thorac. Cancer 10, 1479–1488. doi:10.1111/1759-7714.13097
Italiano, A., Mathoulin-Pelissier, S., Cesne, A. L., Terrier, P., Bonvalot, S., Collin, F., et al. (2011). Trends in survival for patients with metastatic soft-tissue sarcoma. Cancer 117, 1049–1054. doi:10.1002/cncr.25538
Ji, L., Chen, S., Gu, L., and Zhang, X. (2020). Exploration of potential roles of m6A regulators in colorectal cancer prognosis. Front. Oncol. 10, 768. doi:10.3389/fonc.2020.00768
Jiang, W., Zhang, L., Guo, Q., Wang, H., Ma, M., Sun, J., et al. (2019). Identification of the pathogenic biomarkers for hepatocellular carcinoma based on RNA-seq analyses. Pathol. Oncol. Res. 25, 1207–1213. doi:10.1007/s12253-019-00596-2
Jin, Y., Wang, Z., He, D., Zhu, Y., Hu, X., Gong, L., et al. (2021). Analysis of m6A-related signatures in the tumor immune microenvironment and identification of clinical prognostic regulators in adrenocortical carcinoma. Front. Immunol. 12, 637933. doi:10.3389/fimmu.2021.637933
Li, L., Xie, R., and Wei, Q. (2021). Network analysis of miRNA targeting m6A-related genes in patients with esophageal cancer. PeerJ 9, e11893. doi:10.7717/peerj.11893
Li, Q., Ren, C. C., Chen, Y. N., Yang, L., Zhang, F., Wang, B. J., et al. (2021). A risk score model incorporating three m6A RNA methylation regulators and a related network of miRNAs-m6A regulators-m6A target genes to predict the prognosis of patients with ovarian cancer. Front. Cell. Dev. Biol. 9, 703969. doi:10.3389/fcell.2021.703969
Li, Y., Qi, D., Zhu, B., and Ye, X. (2021). Analysis of m6A RNA methylation-related genes in liver hepatocellular carcinoma and their correlation with survival. Int. J. Mol. Sci. 22, 1474. doi:10.3390/ijms22031474
Lin, Y., Pan, X., Zhao, L., Yang, C., Zhang, Z., Wang, B., et al. (2021). Immune cell infiltration signatures identified molecular subtypes and underlying mechanisms in gastric cancer. NPJ Genom. Med. 6, 83. doi:10.1038/s41525-021-00249-x
Makdissi, F. B., Machado, L. V., Oliveira, A. G., Benvenuti, T. T., Katayama, M. L., Brentani, M. M., et al. (2009). Expression of E-cadherin, snail and Hakai in epithelial cells isolated from the primary tumor and from peritumoral tissue of invasive ductal breast carcinomas. Braz J. Med. Biol. Res. 42, 1128–1137. doi:10.1590/s0100-879x2009001200002
Nie, K., Zheng, Z., Wen, Y., Pan, J., Liu, Y., Jiang, X., et al. (2020). A novel ceRNA axis involves in regulating immune infiltrates and macrophage polarization in gastric cancer. Int. Immunopharmacol. 87, 106845. doi:10.1016/j.intimp.2020.106845
Niu, Y., Zhao, X., Wu, Y. S., Li, M. M., Wang, X. J., and Yang, Y. G. (2013). N6-methyl-adenosine (m6A) in RNA: An old modification with a novel epigenetic function. Genomics Proteomics Bioinforma. 11, 8–17. doi:10.1016/j.gpb.2012.12.002
Qi, X., Zhang, D. H., Wu, N., Xiao, J. H., Wang, X., and Ma, W. (2015). ceRNA in cancer: possible functions and clinical implications. J. Med. Genet. 52, 710–718. doi:10.1136/jmedgenet-2015-103334
Salmena, L., Poliseno, L., Tay, Y., Kats, L., and Pandolfi, P. P. (2011). A ceRNA hypothesis: The rosetta stone of a hidden RNA language? Cell. 146, 353–358. doi:10.1016/j.cell.2011.07.014
Sand, M., Skrygan, M., Sand, D., Georgas, D., Hahn, S. A., Gambichler, T., et al. (2012). Expression of microRNAs in basal cell carcinoma. Br. J. Dermatol. 167, 847–855. doi:10.1111/j.1365-2133.2012.11022.x
Sarcomas Tissue, Soft (2021). Statistics. Available at: https://www.cancer.net/ cancer-types/sarcoma-soft-tissue/statistics (Accessed September 25, 2021).
SEER (2021). Cancer statistics review, 1975-2017-SEER statistics, Available at:. https://seer.cancer.gov/csr/1975_2017/ (Accessed September 25, 2021).
Shang, X., Shi, L. E., Taule, D., and Zhu, Z. Z. (2021). A novel miRNA-mRNA Axis involves in regulating transcriptional disorders in pancreatic adenocarcinoma. Cancer Manag. Res. 13, 5989–6004. doi:10.2147/CMAR.S316935
Siegel, R. L., Miller, K. D., and Jemal, A. (2019). Cancer statistics, 2019. Ca. Cancer J. Clin. 69, 7–34. doi:10.3322/caac.21551
Sun, F., Yu, Z., Wu, B., Zhang, H., and Ruan, J. (2020). LINC00319 promotes osteosarcoma progression by regulating the miR-455-3p/NFIB axis. J. Gene Med. 22, e3248. doi:10.1002/jgm.3248
Sun, J., Shi, R., Zhang, X., Fang, D., Rauch, J., Lu, S., et al. (2021). Characterization of immune landscape in papillary thyroid cancer reveals distinct tumor immunogenicity and implications for immunotherapy. Oncoimmunology 10, e1964189. doi:10.1080/2162402X.2021.1964189
Tang, Y., Chen, K., Song, B., Ma, J., Wu, X., Xu, Q., et al. (2021). m6A-Atlas: a comprehensive knowledgebase for unraveling the N6-methyladenosine (m6A) epitranscriptome. Nucleic Acids Res. 49, D134–D143. doi:10.1093/nar/gkaa692
Wang, W., Mu, S., Zhao, Q., Xue, L., and Wang, S. (2019). Identification of differentially expressed microRNAs and the potential of microRNA-455-3p as a novel prognostic biomarker in glioma. Oncol. Lett. 18, 6150–6156. doi:10.3892/ol.2019.10927
Wu, J., Li, L., Zhang, H., Zhao, Y., Zhang, H., Wu, S., et al. (2021). A risk model developed based on tumor microenvironment predicts overall survival and associates with tumor immunity of patients with lung adenocarcinoma. Oncogene 40, 4413–4424. doi:10.1038/s41388-021-01853-y
Xu, F., Huang, X., Li, Y., Chen, Y., and Lin, L. (2021). m6A-related lncRNAs are potential biomarkers for predicting prognoses and immune responses in patients with LUAD. Mol. Ther. Nucleic Acids 24, 780–791. doi:10.1016/j.omtn.2021.04.003
Xuan, J. J., Sun, W. J., Lin, P. H., Zhou, K. R., Liu, S., Zheng, L. L., et al. (2018). RMBase v2.0: Deciphering the map of RNA modifications from epitranscriptome sequencing data. Nucleic Acids Res. 46, D327–D334. doi:10.1093/nar/gkx934
Yang, C., Wang, Y., Xue, W., Xie, Y., Dong, Q., and Zhu, C. (2020). Competing endogenous RNA (ceRNA) network analysis of autophagy-related genes in hepatocellular carcinoma. Pharmgenomics. Pers. Med. 13, 445–462. doi:10.2147/PGPM.S267563
Yang, H., Wei, Y. N., Zhou, J., Hao, T. T., and Liu, X. L. (2017). MiR-455-3p acts as a prognostic marker and inhibits the proliferation and invasion of esophageal squamous cell carcinoma by targeting FAM83F. Eur. Rev. Med. Pharmacol. Sci. 21, 3200–3206.
Yang, K., Wang, F., Li, K., Peng, G., Yang, H., Xu, H., et al. (2022). N6-methyladenosine modification-related long non-coding RNAs are potential biomarkers for predicting the prognosis of patients with osteosarcoma. Technol. Cancer Res. Treat. 21, 15330338221085354. doi:10.1177/15330338221085354
Yang, Y., Hsu, P. J., Chen, Y. S., and Yang, Y. G. (2018). Dynamic transcriptomic m6A decoration: Writers, erasers, readers and functions in RNA metabolism. Cell. Res. 28, 616–624. doi:10.1038/s41422-018-0040-8
Yi, X., Wang, Y., and Xu, S. (2020). MiR-455-3p downregulation facilitates cell proliferation and invasion and predicts poor prognosis of osteosarcoma. J. Orthop. Surg. Res. 15, 454. doi:10.1186/s13018-020-01967-1
Zhang, F., Yu, X., Lin, Z., Wang, X., Gao, T., Teng, D., et al. (2021). Using tumor-infiltrating immune cells and a ceRNA network model to construct a prognostic analysis model of thyroid carcinoma. Front. Oncol. 11, 658165. doi:10.3389/fonc.2021.658165
Zhang, L., Li, B., Peng, Y., Wu, F., Li, Q., Lin, Z., et al. (2020). The prognostic value of TMB and the relationship between TMB and immune infiltration in head and neck squamous cell carcinoma: A gene expression-based study. Oral Oncol. 110, 104943. doi:10.1016/j.oraloncology.2020.104943
Zhang, L., Tang, X. Z., Wan, J., Zhang, X., Zheng, T., Lin, Z., et al. (2021). Construction of a novel signature and prediction of the immune landscape in soft tissue sarcomas based on N6-methylandenosine-related LncRNAs. Front. Mol. Biosci. 8, 715764. doi:10.3389/fmolb.2021.715764
Zheng, F., Du, F., Qian, H., Zhao, J., Wang, X., Yue, J., et al. (2021). Expression and clinical prognostic value of m6A RNA methylation modification in breast cancer. Biomark. Res. 9, 28. doi:10.1186/s40364-021-00285-w
Zhou, Y. J., Zhu, G. Q., Lu, X. F., Zheng, K. I., Wang, Q. W., Chen, J. N., et al. (2020). Identification and validation of tumour microenvironment-based immune molecular subgroups for gastric cancer: Immunotherapeutic implications. Cancer Immunol. Immunother. 69, 1057–1069. doi:10.1007/s00262-020-02525-8
Zhu, K. P., Zhang, C. L., Ma, X. L., Hu, J. P., Cai, T., and Zhang, L. (2019). Analyzing the interactions of mRNAs and ncRNAs to predict competing endogenous RNA networks in osteosarcoma chemo-resistance. Mol. Ther. 27, 518–530. doi:10.1016/j.ymthe.2019.01.001
Zhu, N., and Hou, J. (2020). Assessing immune infiltration and the tumor microenvironment for the diagnosis and prognosis of sarcoma. Cancer Cell. Int. 20, 577. doi:10.1186/s12935-020-01672-3
Keywords: sarcomas, m6A, risk score, prognosis, chemotherapy, immunotherapy
Citation: Li H, Lin D, Wang X, Feng Z, Zhang J and Wang K (2022) The development of a novel signature based on the m6A RNA methylation regulator-related ceRNA network to predict prognosis and therapy response in sarcomas. Front. Genet. 13:894080. doi: 10.3389/fgene.2022.894080
Received: 11 March 2022; Accepted: 12 September 2022;
Published: 12 October 2022.
Edited by:
Ruth Roberts, ApconiX, United KingdomCopyright © 2022 Li, Lin, Wang, Feng, Zhang and Wang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Kai Wang, wangkaimath@sina.com