The detailed workflow for m6A patterns and subsequent analyses are shown in Figure S1A . 21 m6A methylation regulators (“writers”: CBLL1, KIAA1429, METTL14, METTL3, RBM15, RBM15B, WTAP, ZC3H13; “readers”: ELAVL1, FMR1, HNRNPA2B1, HNRNPC, IGF2BP1, LRPPRC, YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3; and “erasers”: ALKBH5 and FTO) were analyzed in STS. Gene Ontology (GO) enrichment analyses of these regulators were conducted, and the significantly enriched biological processes are summarized in Figure 1A . The locations across the chromosomes of the copy number variation (CNV) of the regulators are shown in Figure 1B . CNV was very common and mostly involved amplification, though FMR1, ZC3H13, RBM15, FTO, LRPPRC, and RBM15B had a high frequency of deletion ( Figure 1C ). The interaction patterns among the 21 m6A regulators were also analyzed using the Search Tool for the Retrieval of Interacting Genes/​Proteins (STRING) database ( Figure 1D ). Among the 237 STS samples in the The Cancer Genome Atlas Program (TCGA) cohort, 10 mutations were identified, mutation frequency of the 21 m6A regulators was 4.22%. IGF2BP1 had the highest mutation frequency, followed by RBM15 and YTHDC2 ( Figure 1E ).

Next, we investigated the difference in the expression of the 21 m6A regulators between the normal and tumor tissues. Fat and muscle tissue samples from the The Genotype-Tissue Expression (GTEx) database were used as adjacent normal tissue samples. Of the 21 regulators, 20 showed significant differences between STS and normal tissues, while ZC3H13 did not ( Figure 1F ). Kaplan–Meier (KM) survival analysis showed significant differences in overall survival between patients with high or low expression of the 21 regulators ( Figure S1B ). The above findings suggested that changes in the expression of m6A regulators may play a crucial role in occurrence and progression of STS.

The crosstalk among the 21 m6A regulators and their prognostic value in STS is comprehensively illustrated in the m6A regulator network ( Figure 2A ). All 21 m6A regulators were positively correlated with each other. Next, consensus clustering was performed using the “ConsensusClusterPlus” R package both TCGA and Gene Expression Omnibus (GEO) cohorts. K = 2 was selected based on the empirical cumulative distribution function (CDF) plots ( Figures 2B, C ). Thus, two m6A modification patterns, designated m6Acluster-A (n=284) and m6Acluster-B (n=284), were identified. STS samples with distinct m6A modification could be completely distinguished ( Figure 2D ). Kaplan–Meier (K-M) survival analysis of the two patterns revealed the clear survival advantage in m6Acluster-A, both TCGA-SARC and GSE21050 cohorts ( Figure 2E ). The expression of the 21 regulators was significantly different between the two patterns ( Figure 2F ). To verify the large difference in survival, the TCGA and GEO cohorts (as the validation dataset) respectively were employed to compare the survival between the two patterns ( Figure S2A ). Importantly, similar results were observed, indicating support for the two patterns identified by the clustering algorithm. Additionally, in the TCGA and GEO cohorts, there were significant differences in the expression of m6A regulators between the two patterns ( Figure S2B ).

Compared to m6Acluster-B, m6Acluster-A had increased enrichment of immune cells (using CIBERSORT, MCP-counter, xCell, EPIC, TIMER, quanTIseq and IPS algorithms), especially regarding anti-tumor immune cells, in the TCGA and GEO cohorts ( Figure 2G ). The ImmuneScore and StromalScore (evaluated by the ESTIMATE method) were compared between the two patterns ( Figure 2H ). The differences in immune cell infiltration between the two patterns were also respectively validated in the TCGA and GEO cohorts, and similar results were observed ( Figure S2C ).

To further explore the biological behaviors in the two patterns, Gene Set Variation Analysis (GSVA) and the “limma” package were used, which led to the identification of 84 differential pathways ( Figure S2D ). Typical biological pathways and immune signatures were compared between the two patterns to explore the potential mechanisms. Pathways related to immunity, metabolism, and stemness (cell cycle, DNA damage repair, DNA replication, and mismatch repair) showed significant differences between the two patterns. In particular, the epithelial-to-mesenchymal transition (EMT) and pan-fibroblast TGF-β response signaling pathways were significantly upregulated in m6Acluster-A, which had strong enrichment of CD8+ T cells, effector antigen processing machinery, and immune checkpoints ( Figure 2I ). DNA damage repair, DNA replication, and Wnt signaling pathways were significantly upregulated in m6Acluster-B. These results demonstrated that m6Acluster-A predominantly featured immune and stromal activation, and m6Acluster-B mainly featured DNA repair. Based on above results, we revealed two m6A modification patterns with distinct characteristics of immunity, metabolism, and stemness, which suggested that m6A modification might regulate immune microenvironment, metabolism processes, and tumor cell stemness to contribute to different behaviors of STS.

To further explore the potential significance of each of the 21 regulators, their correlations with immunity, metabolic pathways, and stemness was analyzed. Regarding immunity, we analyzed the correlations between the expression of the 21 m6A regulators and the infiltration of 28 immune cells ( Figure S3A ).

Regarding metabolic pathways, 40 differential metabolic pathways were identified by the “limma” R package. Subsequently, 18 prognosis-related metabolic pathways were selected by univariate Cox analysis and the randomSurvivalForest algorithm ( Table S2 ). The relationship between the error rate and the number of classification trees is shown in Figure S3B . After ranking these metabolic pathways by importance according to the out-of-bag error, five metabolic pathways with relative importance >0.5 were considered in the subsequent analysis ( Figure S3C ). Most metabolic pathways were negatively correlated with the 21 m6A regulators ( Figure S3D ).

Regarding stemness, there were significant correlations between the expression of the 21 regulators and the six stemness indices ( Figure S3E ). The highest correlation coefficient was between RBM15 and mRNAsi. Additionally, HNRNPC, YTHDF2, and HNRNPA2B1 were significantly positively correlated with mDNAsi and mRNAsi. K-M survival analysis for the six stemness indices showed a survival advantage in the lower level of stemness index group ( Figure S3F ). Significant correlations between the expression of the 21 regulators with immune microenvironment, metabolic pathways, and stemness levels indicated that these regulators might play important roles in regulation of m6A modification in terms of immunity, metabolism, and stemness for STS.

To identify immune/metabolism/stemness phenotype-related genes related to the m6A modification patterns, WGCNA was used to identify biologically meaningful modules corresponding to phenotype-related genes. The 12 phenotypes investigated were ImmuneScore (calculated using the ESTIMATE method), five metabolic pathways (selected in the random survival forest analysis), and the six stemness indices.

First, by comparing the two m6A modification patterns, 2183 differentially expressed genes (DEGs) (|log₂FC|>0.5 and FDR<0.05) were identified to be used in WGCNA ( Table S3 ). In the subsequent WGCNA, there were five modules (merged dynamic) ( Figure 3A ). Of the 12 phenotypes, the three that were most correlated with module genes were selected for further analysis. A heatmap revealed the three key modules (MEblue, MEbrown, and MEyellow for ImmuneScore, Retinoic Acid metabolism pathway, and mRNAsi, respectively) ( Figure 3B ). We further analyzed the correlations of the hub genes in the three modules ( Figure S4A ). Ultimately, we identified 579 immune phenotype (ImmuneScore)-related genes, 326 metabolism phenotype-related (retinoic acid metabolism) genes, and 286 stemness phenotype (mRNAsi)-related genes ( Table S4 ).

Three unsupervised consensus clustering analyses were performed based on the immunity/metabolism/stemness phenotype-related genes in the three modules, with the optimal number of clusters being selected based on the corresponding CDF curve ( Figure S4B ). As a result, Immunity groups A–D, Metabolism groups A–C, and Stemness groups A–C, respectively, were defined ( Figure 3C ). Regarding immunity, the K-M survival analysis showed that STS patients in the Immunity A and D groups had better prognoses than those in the Immunity B and C groups ( Figure 3D ). The Immunity A and D groups had massive infiltration of anti-tumor immune cells (CD8+ T cells, macrophages, cytotoxic cells, dendritic cells, and Th1 cells) ( Figure S4C ) and lower expression of m6A regulators ( Figure 3E ). Regarding metabolism, the K-M survival analysis showed that the Metabolism A group had a better prognosis than the Metabolism B and C groups ( Figure 3D ). The three subgroups exhibited different metabolism processes ( Figure S4C ). The Metabolism A group had lower expression of the m6A regulators than the other two groups ( Figure 3E ). Regarding stemness, the K-M survival analysis showed that the Stemness B group had a better prognosis than the Stemness A and C groups ( Figure 3D ). Nevertheless, 15 typical tumor stemness-related biological processes in the three subgroups were comparable ( Figure S4C ). The Stemness B group had lower expression of the m6A regulators than the other two groups ( Figure 3E ).

Using the TCGA-SARC cohort, 531 CpG sites in the 21 m6A regulator genes were identified ( Table S5 ). Subsequently, 41 prognosis-related CpG sites were identified by univariate Cox regression. These sites were used to identify DNA methylation subgroups. K = 2 was selected as the most suitable choice based on the consistency of each cluster and the CDF curve ( Figure S5A ). Thus, the DNA methylation site clustering analysis identified two distinct subgroups, designated DNAmethy-Cluster-A and -B ( Figure S5B ). The heatmap shows the differences in the methylation sites in each subgroup ( Figure S5C ), with higher DNA methylation levels in DNAmethy-Cluster-A. The annotated distribution of clinical traits in the heatmap shows that the two subgroups had unique characteristics. K-M survival analysis showed that DNAmethy-Cluster-A had a higher survival rate ( Figure S5D ). Furthermore, six of the m6A regulators (RBM15B, KIAA1429, YTHDF2, HNRNPA2B1, HNRNPC, and ELAVL1) exhibited lower expression in DNAmethy-Cluster-A ( Figure S5E ).

We further explored the differences in immune cell infiltration, metabolic pathways, and stemness between the two DNA methylation subgroups. As expected, DNAmethy-Cluster-A had an immune-activated phenotype characterized by abundant immune cell infiltration ( Figure S5F ). This subgroup was significantly enriched in metabolic pathways, including phenylalanine metabolism, tryptophan metabolism, and nicotinate and nicotinamide metabolic pathways ( Figure S5G ). Furthermore, this group had lower DNA-related stemness index, as shown in the boxplot in Figure S5H . These results indicated the key roles of the m6A modification in DNA methylation.

STS patients were classified into two m6A modification patterns by consensus clustering based on the expression of 21 m6A regulators. We then examined the potential m6A modification pattern-related gene expression changes between the two patterns, identifying 204 DEGs (log₂FC>1 and FDR<0.05) ( Table S6 ). GO enrichment analysis of these DEGs revealed significant enrichment of T cell mediated immunity, negative regulation of immune response, and positive regulation of cell cycle ( Figure S6A ). The results further demonstrated that the DEGs were characterized by m6A modification, immunity, metabolism pathways, and stemness. The results also confirmed that m6A modification played a key role in the TME. Among the 204 m6A modification pattern-related DEGs, 141 prognosis-related genes were identified by univariate Cox regression ( Table S7 ). Subsequently, these 141 genes were subjected to unsupervised consensus clustering analysis and two stable phenotypes were obtained ( Figures S6B, C ). Ultimately, the samples were divided into two distinct m6A gene signature subgroups, designated geneCluster-A and geneCluster-B ( Figure S6D ). The geneCluster-A was associated with better prognosis ( Figure S6E ). There were significant differences in the expression of the 19 m6A regulators between the two subgroups ( Figure S6F ).

Considering the individual heterogeneity and complexity of m6A modification, we quantified the m6A modification pattern of individual STS patients using principal component analysis based on the 141 abovementioned genes. Thus, m6Ascore was defined for each STS patient. We visualized the changes in the attributes of individual patients in different clusters using an alluvial diagram ( Figure S6G ). To assess the prognostic value of m6Ascore, samples were divided into high- and low-m6Ascore using the optimal cutoff (1.46) determined by the “survminer” R package. Survival was higher in the low-m6Ascore group ( Figure S6H ), as verified in the TCGA and GEO cohorts ( Figure S6I ), and the expression levels of the 17 regulators were also significantly different between the two subgroups ( Figure S6J ). Both m6ACluster-A and geneCluster-A had a lower m6Ascore ( Figure S6K ). These results indicate m6Ascore could be used to predict prognosis in STS.

The histological subtypes of STS in the TCGA-SARC and GSE21050 cohorts mainly included undifferentiated pleomorphic sarcoma (UPS; 35.48%), dedifferentiated liposarcoma (DDLPS; 20.11%), desmoid tumor (DT; 0.4%), leiomyosarcoma (LMS; 35.67%), myxofibrosarcoma (MFS; 4.74%), malignant peripheral nerve sheath tumors (MPNST; 1.7%), and synovial sarcoma (SS; 1.9%). K-M survival analysis of these seven subtypes showed that UPS (with a lower m6Ascore) had improved survival, while LMS (with a higher m6Ascore) had poorer survival ( Figure S7A ). The stacked column chart shows the distribution of the histological subtypes in the high- and low-m6Ascore groups ( Figure S7B ). UPS, DDLPS, and MFS patients were mainly in the low-m6Ascore group, while LMS patients were mainly in the high-m6Ascore group ( Figure S7C ). K-M survival analysis of high- and low-m6Ascore subgroups in each histological subtype showed that low-m6Ascore subgroups had a better prognosis, but the difference was only significant for UPS ( Figure S7D ).

Based on the aforementioned strong associations of the m6A regulators with immunity, metabolism, and stemness in STS patients, we further investigated the correlations of m6Ascore with immunity, metabolism, and stemness. As expected, m6Ascore was significantly correlated with the ImmuneScore and immune cells, including anti-tumor cells (CD8+ T cells, macrophages, Th1 cells, natural killer cells, dendritic cells, and TH17 cells) and pro-tumor immune cells (Th2 cells) ( Figure 4A ). The heatmap of immune cell infiltration (based on CIBERSORT, MCP-counter, xCell, EPIC, TIMER, quanTIseq and IPS algorithms) indicated that the low-m6Ascore group had higher immune cell infiltration, especially regarding anti-tumor-related cells (CD8+ T, dendritic, natural killer, and Th1 cells) ( Figure 4B ). The low-m6Ascore group also had a higher ImmuneScore and StromalScore (calculated using the ESTIMATE method) ( Figure 4C ).

Next, we used the GSVA and “limma” R package to analyze the differences in 114 metabolic pathways between the high- and low-m6Ascore groups, which identified 37 metabolic pathways ( Figure 4D ). Compared to the low-m6Ascore group, the high-m6Ascore group was significantly enriched in the Propanoate metabolism, Lysine degradation, and Glycogen degradation pathways. In addition, there were differences in stemness-related pathways between the two groups ( Figure 4E ).

We compared the somatic mutations in the TCGA-SARC cohort and found that low-m6Ascore group had a higher tumor mutation rate than the high-m6Ascore group (68.61% versus 65.31%) ( Figure 4F ). Differentially mutated genes between the two groups are displayed as a forest plot ( Figure S7E ). The low m6Ascore group had a markedly higher TMB ( Figure 4G ). Higher TMB was associated with increased survival ( Figure 5A ), and low-m6Ascore combined with higher TMB was also associated with better survival ( Figure 5B ). In addition to TMB, we further studied other immunogenic biomarkers and found that intratumor heterogeneity (LOH), DNA damage including homologous recombination deficiency (HRD), tumor neoantigen burden (TNB), intratumor heterogeneity (ITH), and aneuploidy were significantly negatively correlated with m6Ascore ( Figure 4G ). In summary, the differences in tumor immunogenicity between the high- and low-m6Ascore groups were significant ( Figure 5C ).

We further explored CNV between the two m6Ascore groups. The low-m6Ascore group had a higher focal-level gain (p<0.01) and loss (p=0.02) burden and a higher arm-level gain burden (p=0.04) compared to the high-m6Ascore group ( Figure 5D ). Figure 5E shows the distribution of the G-score (based on the frequency and amplitude of the gains and losses) across all chromosomes in the high- and low-m6Ascore groups. Focal amplifications and deletions in various chromosomal regions were detected in both the low- and high-m6Ascore groups ( Figures 5F and S7F ). These results show that the low-m6Ascore group had relatively high immunogenicity, while the high-m6Ascore group had relatively low immunogenicity. Moreover, Gene Set Enrichment Analysis (GSEA) showed that the Hedgehog signaling, Myogenesis, and Spermatogenesis pathways were substantially enriched in the high-m6Ascore group, while the Hypoxia, IL6-JAK-STAT3 signaling, Inflammatory response, KRAS signaling, and P53 pathways were enriched in the low-m6Ascore group ( Figure 5G ).

We further assessed the ability of m6Ascore to predict the clinical benefit of immunotherapy. The low-m6Ascore group had higher expression of immune checkpoint-related genes (especially regarding PDCD1, PD1, and CTLA4) than the high-m6Ascore group ( Figure 6A ). This suggested that patients with different m6Ascores may have different responses to immune checkpoint inhibitors.

We then used the Tumor Immune Dysfunction and Exclusion (TIDE) algorithm to predict the likelihood of response to immunotherapy, and it demonstrated that the low-m6Ascore group had a lower TIDE score and may therefore be more likely to respond to immunotherapy than the high-m6Ascore group (p<0.01) ( Figure 6B ). m6Ascore was significantly positively correlated with TIDE score (cor=0.24, p<0.01) ( Figure 6C ). We also used Subclass Mapping (SubMap) algorithm to compare the expression profile of the two m6Ascore groups with an independent cohort of 47 melanoma patients treated with immunotherapy. The low-m6Ascore group was more likely to respond to anti–PD-1 antibody treatment (nominal p<0.01, Bonferroni-corrected p<0.01) ( Figure 6D ).

We next explored the prognostic value of the m6Ascore in immune checkpoint inhibitor therapy by classifying patients receiving immune checkpoint inhibitor therapy in the TCGA-SKCM cohort to high or low m6Ascore groups. Patients with high m6Ascores had significantly worse survival than those with lower m6Ascores in TCGA-SKCM cohort (p=0.03) ( Figure S7G ). However, response event outcomes were missing in the clinical information from the TCGA-SKCM cohort, so we further validated the predictive performance of m6Ascore in immunotherapy using an external cohort. The IMvigor210 cohort of 348 anti–PD-L1 antibody (atezolizumab)-treated muscle-invasive bladder cancer patients was used to further validate the value of m6Ascore for predicting the clinical benefit of immunotherapy. Based on our scoring strategy, the m6Ascore of each patient in the IMvigor210 cohort was calculated. The low-m6Ascore group had a significant survival advantage, implying that the m6Ascore reflects sensitivity to immunotherapy ( Figure 6E ). The low-m6Ascore group also mainly included patients who responded to immunotherapy ( Figure 6F ). In addition, the m6Ascore was significantly lower in the complete/partial immunotherapy response group than the non-response group (stable/progressive disease) ( Figure 6G ).

Chemotherapy regimens generally involve a combination of several anti-cancer drugs. We assessed the response of the low and high-m6Ascore groups to three common chemotherapeutic drugs for STS: docetaxel, doxorubicin, and gemcitabine. We trained a predictive model on a GDSC cell line dataset using ridge regression, with a satisfactory predictive accuracy evaluated by 10-fold cross-validation. The low-m6Ascore group was predicted to be more sensitive to docetaxel (P <0.01) and gemcitabine (P <0.01) ( Figure 6H ).

To further determine the performance of the m6Ascoring system in various cancer types. The m6Ascores of 10327 samples of 32 cancer types was calculated. Univariate Cox regression indicated that m6Ascore was a favorable factor in ACC, PRAD, MESO, LAML, SKCM, and STAD, and a risk factor in LUAD, PAAD, LGG, KICH, and KIRP ( Figure S8A ). The K-M survival analyses showed that there was a significant difference in overall survival between the high and low m6Ascore groups in the 32 cancer types (p<0.05) ( Figure S9 ). Lower m6Ascore improved prognosis in ESCA, COAD, KIRP, BLCA, READ, PAAD, THYM, UCEC, BRCA, KIRC, LUAD, LGG, KICH, CESC, LIHC, and UVM.

The pan-cancer analysis showed that all 32 cancers exhibited a significant correlation between m6Ascore and ImmuneScore (calculated using the ESTIMATE method) ( Figure S8B ). Next, the correlations between m6Ascore and the proportions of 28 immune cells (calculated using the xCell method) were analyzed in the 32 cancer types. The correlation trends in 32 cancer types differed ( Figure S8C ). The proportions of regulatory T cells, M2 macrophages, and Th2 cells (which are all pro-tumor, immunity-suppressing cells) were correlated with m6Ascore in TGCT, PRAD, OV, ACC, GBM, KIRP, LAML, LUAD and LUSC. The proportions of CD8+ T, dendritic, natural killer, and Th1 cells were correlated with m6Ascore in ACC, CESC, LUAD, TGCT and THYM. We further investigated the correlations between m6Ascore and five important metabolic pathways (selected in the random survival forest analysis) in the 32 cancer types and found that 21 cancers were significantly associated with the Retinoic acid metabolism pathway ( Figure S8D ). Additionally, 31 cancers (all except OV) exhibited significant associations between m6Ascore and the six stemness indices ( Figure S8E ), with negative correlations in DLBC and GBM and positive correlations in BRCA, CESC, ESCA, STAD, and LUSC.

TMB, microsatellite instability (MSI), and expression levels of immune checkpoint-related genes can be used to predict the response to immune checkpoint blockade immunotherapy. Of the 32 cancer types, 18 exhibited a significant correlation between m6Ascore and TMB, as shown in radar charts ( Figure S8F ), 13 exhibited a significant correlation between m6Ascore and MSI ( Figure S8G ), and nine (ACC, LIHC, LUAD, LUSC, MESO, OV, PAAD, STAD, and TGCT) exhibited significant correlations between m6Ascores and both TMB and MSI. We further investigated the correlations between m6Ascore and 15 immune checkpoint-related genes in the 32 cancer types, and there was a significant correlation between PD-L1 (CD274) expression and m6Ascore in 29 cancer types ( Figure S8H ), which again confirmed the ability of m6Ascore to predict the clinical benefit of immunotherapy.