m⁶A regulator-mediated methylation modification patterns and tumor immune microenvironment in sarcoma

Landscape of m⁶A regulators in sarcoma

The workflow of our study is shown in Figure 1. Through literature search and analysis, we identified a total of 23 m⁶A regulators, including 13 readers, eight writers, and two erasers (Supplementary Table 1). We summarized the process by which m⁶A regulators could add or remove m⁶A modification sites and change the potential biological functions of RNA splicing, translation, and degradation (Figure 2A). The investigation of copy number variation (CNV) alteration frequency showed a prevalent CNV alteration in 23 m⁶A regulators, and ALKB homolog 5 (ALKBH5), METTL3, heterogeneous nuclear ribonucleoprotein (HNRNP)-C, and ELAV-like 1 (ELAVL1) showed widespread CNV frequency gains. In contrast, METTL15, zinc finger CCCH-type containing 13 (ZC3H13), fat mass and obesity-associated (FTO), leucine-rich pentatricopeptide repeat containing (LRPPRC), and RNA-binding motif protein 15B (RBM15B) showed widespread CNV frequency loss (Figure 2B). We further searched for the position of CNV alteration in the m⁶A regulator on the chromosome (Figure 2C). We can completely distinguish sarcoma samples (259 sarcoma samples from The Cancer Genome Atlas [TCGA] database) from normal samples (911 muscle and adipose tissue samples from the University of California Santa Cruz [UCSC] Xena database) based on the expression levels of these 23 m⁶A regulators (Figure 2D). The median overall survival (OS) of patients with sarcoma was found to be 2.57 years (IQR = 1.3–4.34 years). Among them, Vir-Like m⁶A methyltransferase associated (VIRMA), METTL14, METTL3, Wilms tumor 1 (WT1)-associated protein (WTAP), fragile-X mental retardation 1 (FMR1), HNRNPA2B1, HNRNPC, IGF2BP3, LRPPRC, YTH domain containing (YTHDC)-1, YTHDC2, and FTO were highly expressed in normal tissues, while RBM15, RBM15B, ZC3H13, ELAVL1, IGF2BP1, YTH m⁶A RNA-binding protein (YTHDF)-1, YTHDF3, and ALKBH5 were highly expressed in sarcoma tissues. The expression levels of Cbl-like 1 (CBLL1) and IGF2BP2 in normal and sarcoma tissues were not statistically significant. Gene Ontology (GO) enrichment results are shown in Supplementary Figure 1A. We then used Cox regression analysis and Kaplan-Meier (K-M) analysis to determine the relationship between m⁶A regulators and the prognosis of patients with sarcoma. Univariate Cox regression analysis showed that the expression levels of IGFBP2 (P-value = 0.002), VIRMA (P-value = 0.039), IGFBP1 (P-value <0.001), HNRNPC (P-value = 0.009), HNRNPA2B1 (P-value = 0.004), YTHDF2 (P-value <0.001), and IGF2BP3 (P-value = 0.005) were protective factors for patients with sarcoma (Supplementary Figure 1B and Supplementary Table 2). At the same time, we calculated the K-M value of m⁶A regulators in patients with sarcoma by each best cut-off point (Supplementary Table 3) and the survival curve of m⁶A regulators (Supplementary Figure 2A–2W).

Figure 1. Design and workflow of the study.

Figure 2. Landscape of N⁶-methyladenosine (m⁶A) regulators in patients with sarcoma. (A) Regulation of m⁶A regulators and their biological functions in RNA metabolism. A total of 23 known m⁶A regulators, including 13 readers, 8 writers, and 2 erasers. m⁶A methylation involves multiple stages of the RNA life cycle, including pre-mRNA shearing, pre-mRNA splicing, pre-mRNA transport, RNA translation, RNA degradation, etc. (B) Copy number variation (CNV) mutations in 23 m⁶A regulators. (C) The location of CNV alterations of m⁶A regulators on chromosomes. (D) Differential expression levels of 23 m⁶A regulators between normal tissues (Adipose and muscle tissues) and tumor tissues (sarcoma samples). Asterisks represent the statistical P-values (^*P < 0.05; ^**P < 0.01; ^***P < 0.001).

Identification of m⁶A methylation modification patterns

The m⁶A regulator network describes the association between the 23 m⁶A regulators, their interactions, and their prognostic significance in patients with sarcoma (Figure 3A). The results of the m⁶A regulatory network indicate that the interaction between writers, readers, and erasers may play unique roles in the formation of different m⁶A regulator expression patterns and affect the prognosis of patients with sarcoma [21, 22]. We used non-negative matrix factorization (NMF) consistent clustering to classify the sarcoma samples and identified four m⁶A modification clusters (Figure 3B and Supplementary Figure 3A–3J), namely, m⁶A-cluster-A, m⁶A-cluster-B, m⁶A-cluster-C, and m⁶A-cluster-D. K-M analysis showed that m⁶A-cluster-C had a significant survival advantage, while m⁶A-cluster-D showed a poor prognosis (Figure 3C). A heat map showed the differential expression patterns of 23 m⁶A regulators in four different m⁶A modification clusters (Figure 3D).

Figure 3. Establishment of m⁶A methylation modification patterns. (A) Regulatory functions and interactions of the 23 m⁶A regulators in sarcoma. The size of the circle represents the influence of different regulators on prognosis, while the Log-rank test was used to calculate values ranging from P < 1e⁻⁴, P < 0.001, P < 0.01, P < 0.05, and P > 0.05. Purple dots in circles, prognostic risk factors; green dots in circles, prognostic favorable factors. Erasers, readers, and writers are indicated by red, orange, and gray dots, respectively. The regulators' lines show their interactions, with negative correlations marked in blue and positive correlations in red. (B) The consensus matrix heatmap defined four m⁶A methylation modification clusters from 259 patients with sarcoma. (C) Kaplan-Meier curves of the overall survival rates of 259 patients with sarcoma with four m⁶A methylation modification clusters. (D) The heat map shows the expression levels of 23 m⁶A regulators in different m⁶A modification clusters and various clinicopathological features.

TIME landscape characteristics in four m⁶A regulator expression patterns

To explore the biological behavior of the four different m⁶A modification types, we performed pairwise gene set variation analysis (GSVA) enrichment analysis (Figure 4A–4F). GSVA results showed that m⁶A-cluster-A was significantly rich in material metabolism (amino acid and lipid metabolism) and human diseases (immune-related, such as systemic lupus erythematosus and primary immunodeficiency). It is worth noting that material metabolism-related processes, especially glucose metabolism, are enriched in m⁶A-cluster-B. The pathways enriched by m⁶A-cluster-C are phosphoinositide transduction and metabolism and human diseases (circulatory system related). The pathways enriched by m⁶A-cluster-D are mainly metabolism-related pathways (such as glycan metabolism and lipid metabolism) and human disease-related signaling pathways (basal cell carcinoma and diabetes). Surprisingly, GSVA results showed that the signaling pathways with high expression between different clusters were mostly material metabolism- and immune-related, further confirming a direct or indirect link between m⁶A regulator expression patterns and TIME. In addition, we constructed a visible box-line plot by the “ssGSEA” package to compare the differences in the relative abundance of 28 immune infiltrating cell subpopulations in 4 different m⁶A clusters (Figure 4G). Statistically significant differences in immune cell infiltration results indicate the reliability of the m⁶A cluster in elucidating the immunophenotyping of TIME. Three-dimensional PCA (3D PCA) analysis was performed in four different m⁶A clusters based on the expression of 23 m⁶A regulators. 3D PCA and its projection on the three planes (Figure 4H–4K) showed that the distribution of four different m⁶A modification clusters was disorderly.

Figure 4. Gene set variation analysis (GSVA), single sample gene set enrichment analysis, and principal component analysis (PCA). The heat map shows the GSVA scores of representative hallmark passages by comparing two by two m⁶A regulator expression patterns. (A) m⁶A cluster-A vs. m⁶A cluster-B; (B) m⁶A cluster-A vs. m⁶A cluster-C; (C) m⁶A cluster-A vs. m⁶A cluster-D; (D) m⁶A cluster-B vs. m⁶A cluster-C; (E) m⁶A cluster-B vs. m⁶A cluster-D; (F) m⁶A cluster-C vs. m⁶A cluster-D; (G) Abundance of infiltrating cells in each of the four m⁶A modification patterns. Asterisks represent the statistical P-values (^*P < 0.05; ^**P < 0.01; ^***P < 0.001). (H) Three-dimensional PCA (3D PCA) results of four m⁶A modification clusters, showing significant differences in the transcriptome between different modification clusters. (I) Projection of 3D PCA on the first principal component (PC1) and the second principal component (PC2). (J) Projection of 3D PCA on PC1 and the third principal component (PC3). (K) Projection of 3D PCA on PC2 and PC3.

Generation of m⁶A gene signatures and functional annotation

To deeply explore the potential biological functions of the four m⁶A regulator expression patterns, we analyzed the differentially expressed genes (DEGs) of the two m⁶A modification clusters using the “limma” package and took the intersection (Figure 5A and Supplementary Table 4). Finally, we identified 22 m⁶A-related DEGs (Supplementary Table 5). Surprisingly, based on the results of GO (Figure 5B) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (Figure 5C), these 22 m⁶A DEGs showed that tumor metabolism was closely related to m⁶A modification, which reinforces that m⁶A modification plays a non-negligible role in TME. We further performed univariate Cox analysis on 22 DEGs and obtained 19 prognosis-related DEGs (Supplementary Table 6). To further validate the mechanism between m⁶A and TIME, we performed consensus clustering analysis based on the 19 m⁶A prognosis-related DEGs. We classified the samples into three different genomic subtypes and named these three clusters as m⁶A gene clusters A, B, and C (Figure 5D and Supplementary Figure 4A–4J). Heat map showing differences in the expression of 19 m⁶A prognosis-related DEGs in the three m⁶A gene clusters (Figure 5E). K-M analysis results showed that m⁶A gene cluster B had a significant survival advantage, while m⁶A gene cluster C had a poor prognosis (Figure 5F). Immediately afterward, we found significant differences in the expression of m⁶A regulators among the three m⁶A gene clusters (Figure 5G).

Figure 5. Development of the m⁶A score and exploration of functional annotation. (A) Venn diagram showing 22 m⁶A-associated differentially expressed genes (DEGs) between the six clusters compared. (B) Functional annotation of m⁶A-related DEGs using Gene Ontology (GO) enrichment analysis. (C) Functional annotation of m⁶A-related DEGs using Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis. (D) The consensus matrix heatmap defined three m⁶A gene clusters. (E) The heat map shows the expression levels of 22 m⁶A-related DEGs in different m⁶A gene clusters and various clinicopathological features. (F) Kaplan-Meier curves of the overall survival rates of 259 patients with sarcoma with three m⁶A gene clusters. (G) Expression levels of 23 m⁶A regulators in each of the three m⁶A gene clusters. Asterisks represent the statistical P-values (^*P < 0.05; ^**P < 0.01; ^***P < 0.001).

Development of the m⁶A score and its clinical significance

Although our studies have further confirmed the relationship between m⁶A modification in TIME and prognosis, these analyses cannot predict the pattern of m⁶A modification in patients with sarcoma and cannot assess the prognosis of these patients. Therefore, we developed a scoring scheme called the m⁶A score to quantify the m⁶A regulator expression pattern in each sarcoma patient (Supplementary Table 7). We used a Sankey plot to illustrate the workflow (Figure 6A). The Kruskal-Wallis test showed that there were significant differences in m⁶A scores between the m⁶A clusters. The median of m⁶A cluster C was the lowest, and the median of m⁶A cluster D was the highest, indicating that a high m⁶A score may be closely related to the characteristics of material metabolism, while a low m⁶A score may be related to immune cell infiltration (Figure 6B). At the same time, m⁶A gene cluster C showed the highest median score of m⁶A, and m⁶A cluster B showed the lowest median score (Figure 6C). The results of the K-M curve showed that patients with high m⁶A scores had a shorter survival time than those with a low score group (Figure 6D). In different clinicopathological characteristics (age, race, sex, site, margin status, and metastasis status), it was further confirmed that the prognosis of patients with high m⁶A scores was worse than that of patients with low m⁶A scores (Supplementary Figure 5A–5L). We drew the receiver operating characteristic (ROC) curve to further explore the prediction function of the m⁶A score for patients with sarcoma. We found that the AUC of m⁶A score was 0.748, which is much larger than the AUC of other clinicopathological features, and the m⁶A score can predict the prognosis of patients with sarcoma (Figure 6E). The m⁶A score was validated in GSE63157, and the AUC value of the validation cohort was 0.675 (Figure 6F). Additionally, the validation cohort GSE63157, survival curves showed that low-risk patients had a favorable prognosis than high-risk patients (P < 0.001, Figure 6G).

Figure 6. Development of the m⁶A score and exploration of m⁶A-related clinical features. (A) Sankey diagrams of different m⁶A modification clusters, m⁶A gene clusters, m⁶A scores, and different prognosis status. (B) Differences in the m⁶A scores among four m⁶A modification clusters (P-value <0.001, Kruskal-Wallis test). (C) Differences in the m⁶A scores among three m⁶A gene clusters (P-value <0.001, Kruskal-Wallis test). (D) Survival analyses for low-and high-m⁶A score groups using Kaplan-Meier curves (P-value < 0.001, Log-rank test). (E) Comparison of the area under the receiver operating characteristic curve between the m⁶A score and clinicopathological characteristics. (F) The receiver operating characteristic curve of m⁶A score in the validation set GSE63157 (Area under curve = 0.675). (G) Survival analyses for low- and high-m⁶A score groups using Kaplan-Meier curves in validation set GSE63157 (P-value < 0.001, Log-rank test).

Role of the m⁶A score in immunotherapy

Using the Tumor Immune Estimation Resource (TIMER), we evaluated the association between the prognostic features and immune cell infiltration. The results indicated that DC (Cor = –0.175, P-value = 0.005, Figure 7A), macrophages (Cor = –0.320, P-value <0.001, Figure 7B), CD4⁺ T cells (Cor = –0.177, P-value = 0.005, Figure 7C), and neutrophils (Cor = −0.279, P-value <0.001, Figure 7D) were significantly negatively correlated with the m⁶A score. The other two types of immune cells, CD8⁺ T cells (Figure 7E) and B cells (Figure 7F) were not significantly different from the m⁶A score. Immune cell infiltration analysis results of TIMER further confirmed that a low m⁶A score might be related to immune cell infiltration. In addition, using CIBERSORT, we assessed the difference in immune cell infiltration between the high-m⁶A score group and the low-m⁶A score group. The results of CIBERSORT showed that the infiltration degrees of activated B cells, memory T cells, CD4⁺ T cells, M0 and M2 macrophages, and mast cells were higher in the high-m⁶A score group than in the low-m⁶A score group; in contrast, the infiltration degrees of naïve B cells, M1 macrophages, and resting mast cells were higher in the low-m⁶A score group than in the high-m⁶A score group (Figure 7G). Immune cell infiltration analysis results of CIBERSORT showed that different m⁶A scores had different immune cell types infiltrated in TIME. Based on the m⁶A score, we divided six patients from Tianjin Medical University Cancer Institute and Hospital into a low-m⁶A score group (four samples) and a high-m⁶A score group (two samples) and plotted their immune cell infiltration heat map (Figure 7H). By analyzing the degree of immune cell infiltration in the external validation set of the six patients, we further confirmed the conclusion that macrophage M2 and neutrophils were higher in the high-m⁶A score group than in the low-m⁶A score group (Figure 7I). The results of ICI-related biomarker analysis showed that the expression of TIM-3 (p = 0.005) was higher in the high-m⁶A score group than in the low-m⁶A score group (Figure 7J). The validation set GSE63157 high-m⁶A score group also had a higher level of TIM-3 expression than the low-m⁶A score group (Figure 7K). However, there was no statistical difference in the expression of PD1 (Figure 7L), PD-L1 (Figure 7M), CTLA4 (Figure 7N), T cell immunoglobulin and ITIM domain (TIGIT) (Figure 7O), and lymphocyte activation gene-3 (LAG3) (Figure 7P) between the high- and low-m⁶A score groups.

Figure 7. Role of the m⁶A score in immunotherapy. Infiltration abundances of six types of immune cells (Pearson correlation analysis). (A) Dendritic Cells (Cor = –0.175, P-value = 0.005). (B) Macrophages (Cor = –0.320, P-value <0.001). (C) CD4⁺ T cells (Cor = –0.177, P-value = 0.005). (D) Neutrophil (Cor = –0.279, P-value <0.001). (E) CD8⁺ T cells (Cor = –0.020, P-value = 0.755). and (F) B cells (Cor = –0.047, P-value = 0.454). (G) Box plots visualizing significantly different immune cells between high- and low-m⁶A score groups. Infiltration degrees of naive B cells, M1 macrophages, and resting mast cells were higher in the low-m⁶A score group than the high-m⁶A score group. Meanwhile, the infiltration degrees of memory B cells, activated memory CD4 T cells, M0 and M2 macrophages, and activated mast cells were higher in the high-m⁶A score group than the low-m⁶A score group. (H) Ten immune cells infiltration heat map of the six patients from Tianjin Medical University Cancer Institute and Hospital. (I) Infiltration degrees of M2 macrophages and neutrophils were higher in the high-m⁶A score group than the low-m⁶A score group in the six patients validation set. The box plot shows the differences in (J) hepatitis A virus cellular receptor 2 (HAVCR2) between high- and low-m⁶A score groups in TCGA database, (K) HAVCR2 between high- and low-m⁶A score groups in the validation set, GSE63157, and (L) programmed cell death 1 (PD-1), (M) CD274, (N) cytotoxic T lymphocyte-associated antigen-4 (CTLA-4), (O) T cell immunoglobulin and ITIM domain (TIGIT), and (P) lymphocyte activation gene-3 (LAG3) between high- and low-m⁶A score groups in The Cancer Genome Atlas (TCGA) database. Asterisks represent the statistical P-values (^*P-value <0.05; ^**P-value <0.01; ^***P-value <0.001).

Data collection and processing

Gene expression datasets and corresponding clinical datasets (including age, race, gender, site, margin status, metastasis status, status, and survival months) for patients with sarcoma were downloaded from TCGA database (https://www.cancer.gov/). Patients with sarcomas without related overall survival were excluded. Gene expression datasets of normal muscle and adipose tissue samples were downloaded from the UCSC Xena database (https://xenabrowser.net/). To achieve comparability between gene expression data, we converted gene expression datasets of UCSC from Fragments Per Kilobase Million (FPKM) format to transcripts per kilobase million (TPM) format. Genome mutation data (CNV) of TCGA-SCAR were also downloaded from the UCSC database. By searching the literature related to m⁶A methylation modification, 23 confirmed m⁶A regulators were included in our study [7, 33, 34]. We downloaded the gene file of “c2.cp.kegg.v7.4.symbols.gmt” from the Molecular Signatures Database (MSigDB) for GSVA analysis. The 23 m⁶A regulators included 13 readers, 8 writers, and 2 erasers (Supplementary Table 1).

Data collection for validation

We obtained gene expression profiles and clinical data for an independent cohort, GSE63157, from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). GSE63157 was used as an external validation cohort. We collected STS patients admitted to the Tianjin Medical University Cancer Institute and Hospital between 2016 and 2019. Tissue specimens from all patients were reviewed by pathologists and diagnosed with STS. All patients underwent whole-exome sequencing (WES). We quantified the scores of 10 immune cell types from the RNA-seq data of the six samples. The retrospective investigation was conducted in accordance with the Declaration of Helsinki and approved by the Ethics Committee of Tianjin Medical University Cancer Institute and Hospital (Approval No. E2019144). All patients provided signed informed consent. Genetic sequencing data were obtained from the tissue samples of six patients (the trial registration was NCT04126993) [35]. The Yuce Bio Company conducted gene sequencing procedures. Genomic DNA from formalin-fixed, paraffin-embedded (FFPE) sections from biopsy samples or whole blood control samples were extracted using the Gene Read DNA FFPE Kit (Qiagen, Germantown, MD, USA) and the Mag-Bind Blood and Tissue DNA HDQ 96 Kit (Qiagen), respectively. Library preparations were performed using the KAPA Library Quantification Kit (Roche, Indianapolis, IN, USA), and target enrichment was performed using the Target Seq Enrichment Kit (iGene Tech, Beijing, China), and sequencing was performed on a NovaSeq (Illumina, San Diego, CA, USA). The raw reads of WES-seq were processed using SOAPnuke (version 1.5.6, parameters: -l 20 -q 0.1 -n 0.1) to remove ambiguous reads and/or low-quality reads. These qualified sequence reads were then aligned to the human reference genome (UCSC hg38) using BWA-mem (BWA, version 0.7.12). We quantified the scores of 10 immune cell types from the RNA-seq data of the six samples.

CNV analysis and differential analysis of m⁶A regulators and establishment of m⁶A regulator network

To explore the extent of mutation of 23 m⁶A regulators in sarcoma, we performed CNV analysis and labeled the positions of 23 m⁶A regulators in human chromosomes. Gene Ontology (GO) enrichment analysis based on 23 m⁶A regulators was performed. Combined with the datasets of sarcoma samples in TCGA and the datasets of muscle and adipose tissues in UCSC, we performed differential analysis. Immediately afterward, we combined the survival data of patients with sarcoma with K-M analysis for each m⁶A regulator and explored the interrelationships between the m⁶A regulators. Based on the results of the above analysis, an m⁶A regulator network was established.

Consistent cluster analysis of 23 m⁶A regulators

NMF consensus clustering adds non-negative constraints to the decomposed matrix based on matrix factorization. For example, matrix A (M_A) is decomposed into two non-negative matrices B (M_B) and C (M_C) (i.e., M_A ≈ M_B × M_C, M_B ≥ 0, M_C ≥ 0) [36]. NMF consensus clustering was applied to identify different m⁶A regulator expression patterns by the “ConsensusClusterPlus” package, and we named the m⁶A regulator expression pattern the m⁶A cluster. The selection of the optimal number of clusters was determined by the cophenetic, consensus index, and silhouette coefficients. According to the consensus clustering results, we performed K-M analysis and plotted a heat map of patients in different m⁶A clusters.

GSVA

GSVA is a non-parametric unsupervised analysis method that assesses whether different metabolic pathways are enriched among different clusters by combining the analysis of differences in gene expression between clusters [37]. We analyze the differences in biological processes between four different m⁶A regulator expression patterns by the “GSVA” package. Adjusted P-values less than 0.05, and log fold change (log FC) less than 0.1 were considered statistically significant.

Estimation of immune cell infiltration and Principal Component Analysis (PCA)

The enrichment score calculated by single-sample gene set enrichment analysis (ssGSEA) represents the relative abundance of each immune cell in each m⁶A cluster [38]. PCA uses the idea of dimensionality reduction to transform K dimensions into M dimensions (K > M) with a small loss of information and uses M dimensions to explain the variance-covariance structure of multiple variables, thereby simplifying the system structure [39]. In our study, the N dimensions represent 23 m⁶A modulators, and the K dimensions represent two dimensions, which can be visualized. At the same time, we used different colors to distinguish the different m⁶A clusters.

Identification of m⁶A related DEGs

To identify m⁶A-related genes, we divided patients into four different m⁶A clusters (Clusters A, B, C, and D) based on the expression of 23 m⁶A regulators. Next, we performed a two-by-two comparison of different clusters (Cluster AvsB, AvsC AvsD, BvsC, BvsD, and CvsD) to identify DEGs. We combined the DEGs intersection of the 6 groups (Cluster AvsB, AvsC AvsD, BvsC, BvsD, and CvsD) by jvenn [40]. Moreover, we obtained prognosis-related DEGs by performing univariate Cox regression analysis. The log FC and adjusted P-values were used to evaluate the significance of m⁶A DEGs. The filter criterion was set to log FC less than 2.00, and the adjusted P-value was less than 0.05.

Establishment of m⁶A gene clusters and m⁶A score

We then established different m⁶A gene clusters to analyze the m⁶A DEGs. The process is as follows: (I) Analysis of differential genes using unsupervised clustering; (II) The number and stability of gene clusters were also defined by the consensus clustering algorithm. Based on the results of the m⁶A gene clusters, we performed a K-M analysis. Furthermore, to further quantify the m⁶A methylation modification in patients with sarcoma, we constructed a novel scoring system using prognosis-related m⁶A DEGs and named it the m⁶A score. The m⁶A score was constructed as follows: (I) We performed PCA and extracted the principal components 1 and 2 (PC1 and PC2); (II) Calculate the m⁶A score using the following equation:

$${\text{m}}^6\text{A score}=\sum ({\text{PC1}}_i+{\text{PC2}}_i)$$

where i is the prognosis-related m⁶A DEGs.

Based on optimal cut-off points, patients with sarcoma were separated into high- and low m⁶A score groups. To further verify the m⁶A score, we used the Kruskal-Wallis test to evaluate whether the m⁶A score was different in the m⁶A cluster and the m⁶A gene cluster. To further demonstrate the characteristics of the m⁶A score, we assessed the survival status between high- and low-m⁶A score groups by K-M analysis. We also assessed survival status in different clinicopathological characteristics, including age (≤65, >65 years), race (white and others [Asian, black, and African American]), sex (female and male), site (soft tissues and parts other than soft tissues [uterus and retroperitoneum, etc.]), margin status (R0 and R1/2), and metastasis status (yes and no/unknown), according to high- and low-m⁶A score groups. To further explore the predictive function of the m⁶A score for patients with sarcoma, we drew the receiver operating characteristic curve (ROC) and calculated the area under the curve (AUC) for the m⁶A score and the different clinicopathological characteristics. To further test the accuracy of the m⁶A score, we calculated the m⁶A score values of the independent cohort (GSE63157) and plotted the ROC curve and K-M analysis.

Correlation between m⁶A score and tumor immune microenvironment

First, we analyzed the differences in ICI-related biomarkers (CTLA-4, PD-1/PD-L1, hepatitis A virus cellular receptor 2/T cell immunoglobulin domain and mucin domain-3 [HAVCR2/TIM-3], lymphocyte activation gene 3 [LAG3], and T cell immunoreceptor with immunoglobulin and ITIM domain [TIGIT]) between the high- and low-m⁶A score groups. Second, we used the CIBERSORT database to identify a complex association between 22 different immune cells and different m⁶A score groups [41]. Third, TIMER [42] (https://cistrome.shinyapps.io/timer/) was used to systematically analyze and estimate the abundance of the six immune cell immune infiltrates in the m⁶A score. Finally, we used the set GSE63157 from GEO and six samples from Tianjin Medical University Cancer Institute and Hospital to verify the relationship between the m⁶A score and TIME.

Statistical analyses

Statistical analyses in this study were performed using Ri386-4.0.3. Differences in quantitative data and normally distributed variables were compared using the t-test, and differences in non-normally distributed variables were compared using the Wilcoxon rank-sum test. Differences were compared for more than two groups of variables using one-way analysis of variance and the Kruskal-Wallis test. Prognostic analysis was performed using the Kaplan-Meier survival analysis and Cox proportional hazards model. Pearson’s analysis was used for the correlation analysis. A P-value <0.05 (two-tailed) was considered to indicate statistical significance. The Benjamini-Hochberg method was used to control for FDR for multiple hypothesis testing.

Ethics approval and consent to participate

TCGA, UCSC, and GEO data did not involve animal experiments, human specimens, or ethics-related issues. Ethical approval for the six patients was obtained from the research ethics committee of the Cancer Institute and Hospital of Tianjin Medical University prior to the study. Written informed consent was obtained from all patients and/or their families. All specimens were handled and stored anonymously according to ethical and legal standards.