The sequencing and corresponding clinical data for STS were acquired via the GEO and TCGA databases. The transcriptome data within the GTEx database for normal adipose and muscle tissues was employed as the control group. The UCSC bioinformatic pipeline (TOIL RNA-seq) was utilized to perform co-analysis of the datasets (Wang et al., 2018). The somatic mutation and CNVs were obtained through TCGA-SARC cohort as well. TARGET Pan-Cancer (PANCAN) cohort was utilized for pan-cancer analysis. We identified two bulk RNA-sequencing datasets (GSE39055 and GSE176307) and one scRNA-seq dataset through GEO database. For examining relationships among SRGs with drug responses, the study included a patient population who had undergone immunotherapeutic therapy using concomitant inhibition of both PD-1 and CTLA-4.

A total of 279 gene candidates associated with cellular senescence were acquired through the use of the CellAge database, a comprehensive cellular senescence repository (https://genomics.senescence.info/cells/). In order to perform Least Absolute Shrinkage and Selection Operator (LASSO) regression analysis, we implemented the “glmnet” package. Additionally, a tenfold cross-validation approach was employed for determining the optimal penalty regularization parameter λ. To minimize empirical error, we utilized the support vector machine recursive feature elimination (SVM-RFE) method, which follows the principle of structural risk minimization, another machine learning approach. Genomic locations of SRGs were virtualized by using “Rcircos” package. Unsupervised clustering was performed to identify clusters and subtypes associated with SRGs. MaxK = 9 and repetitions = 1,000 was set as key parameters for R package “ConsensusClusterPlus” (Wilkerson and Hayes, 2010).

For comparing the signatures of pathways among different clusters or subtypes, GSVA analysis was performed by utilizing the R package “GSVA” (Hänzelmann et al., 2013) with predefined datasets from the MSigDB. We analyzed the data on the basis of package “limma”, while modified t-statistics were employed for visualization. We utilized the “clusterProfiler” package for performing GO annotation and set the false discovery rate (FDR) threshold as 0.05 (Yu et al., 2012). For visualizing the correlation among SRGs in STS, we utilized the “corrplot” package. For the purpose of visualizing the relationships among SRGs and prognosis, we utilized the R package “igraph”.

We determined DEGs among the senescence-related clusters or subtypes using the R package “limma”. The p-values were adjusted based on Benjamini-Hochberg method for addressing the issue of multiple comparisons. We defined a statistically significant difference as the adjusted p-value <0.05. Specific threshold was adopted to assess the expression difference of genes.

We conducted an analysis for assessing the infiltration levels of immunocytes in STS based on single-sample gene set enrichment analysis (ssGSEA), as previously described (Bindea et al., 2013). The degree of immune infiltrations was scaled to range from 0 to 1. Relationships of TME with other biological process was evaluated by using signatures of tumor mutation burden (TMB) (Mariathasan et al., 2018). We also utilized the R package “ESTIMATE” to compute ESTIMATE scores based on gene signatures (Yoshihara et al., 2013), to assess the levels stromal and immune infiltrations within tumor. Potential signatures involved in immunotherapy response and cancer-immunity cycles were retrieved from previous research (Chen and Mellman, 2013; Qi et al., 2023). The association between senescence scores and GSVA scores was evaluated by using the R package “ggcor”.

For the quantification of senescence, we introduced the senescence scoring system, following the establishment of senescence-related clusters. We picked up the DEGs among different senescence-related clusters and conducted the PCA to obtain senescence score. On this setting, genes related to most set factors will display a high score and vice versa, as previously reported (Zhang et al., 2020; Chong et al., 2021). Equation used to calculate the senescence socre was Σ (PC1i + PC2i) (i indicates the expression levels for the selected genes based on the PCA).

The single-cell sequencing dataset GSE131309 (Jerby-Arnon et al., 2021) was analyzed following the standard pipelines, using the Seurat package. We applied LogNormalize (scale factor = 10,000) for normalizing the gene expression. Subsequently, we used FindVariableGenes to recognize the 2,000 highly variable genes. Cell sub-population was labeled using the annotation methods described in previous publication (Jerby-Arnon et al., 2021).

The Genomics of Drug Sensitivity in Cancer (GDSC) is the publicly available source that contains molecular features of cancers predicting response to anti-cancer drugs. The database records thousands of tumor cell lines and 518 compounds. For determining corresponding IC50 and drug sensitivity score, we used the “pRRophetic” and “oncoPredict” (Iorio et al., 2016; Maeser et al., 2021).

Four STS cell lines were used in this study, including SW-872, hSS-005R, SW-982 and HSF. SW-872, SW-982 and HSF were obtained from ATCC while hSS-005R was established by our laboratory. In the 37°C environment with 5% CO₂, the cell lines were grown in Dulbecco’s modified Eagle medium (DMEM) supplemented with 10% fetal bovine serum (FBS).

We collected and further sequenced four pairs of STS samples as well as their corresponding normal tissues using Oxford Nanopore Technologies (Oxford, United Kingdom). The sequencing matrix is available under accession number GSE198568. We further verified corresponding expression for SRGs in this dataset.

We used the RNA Express Total RNA Kit (M050, NCM Biotech, China) to extract total RNA from the cells. Following the removal of genomic DNA, we used The RevertAid First Strand cDNA Synthesis kit (K1622, Thermo Fisher Scientific, United States) for synthesizing cDNA from the total RNA. Detailed steps were similar to previous study (Qi et al., 2022; Qi et al., 2023). As shown in Supplementary Table S1, this study utilized the following primers.

The Hanbio (Shanghai, China) designed and synthesized the siRNA for E2F1 knockdown, as well as the corresponding negative control (siNC). The cells were seeded into a 6-well plate with the concentration of 1.5 × 10⁵ cells per well. The 50 nmol of siRNA or siNC was further transfected into the cells using 5 μL of Lipofectamine 2000 reagent. Following a 6-h transfection period, fresh medium was further added to the culture medium. The siRNA sequences were listed in Supplementary Table S2. We collected or analyzed the cells 48 h after transfection and used them for subsequent experiments.

We plated cells into a 96-well plate with the concentration of 2000 cells per well and cultured overnight. The cell counting kit-8 (CCK-8) was utilized for testing the cells that underwent transfection after being cultured after the indicated periods of time. For each well, we added 100 μL of 10% CCK-8 solution was added, followed by 1.5 h of incubation. The absorbance at 450 nm for each well was detected by using the Tecan Spark^® multimode microplate reader.

After transfection with siRNA, cells were seeded in 6-well plates with the concentration of 1,000 cells per well and subsequently cultivated for the indicated durations. At the end of experiment, cellular clones were fixed utilizing 4% paraformaldehyde (PFA) and subsequently subjected to staining by using 0.2% crystal violet solution with an additional period of 15 min.

In order to assess the migratory capabilities of the cells, we utilized the wound healing assay. Cells were seeded into 6-well plates with the concentration of 4 × 10⁵ cells per well. Upon achieving approximately 80% cellular confluence, we created a scratch employing the 100 μL pipette tip. Subsequently, cells were cultivated in DMEM with 2% FBS. The area covered by the migrated cells was then measured after indicated duration, by the light microscope.

This study employed the R statistical analysis software (version 4.1.0) as part of its methodological approach. Spearman’s correlation test was employed for determining the associations among the studied SRGs. Parametric comparisons were calculated by utilizing Student’s t-tests while nonparametric comparisons were determined through the Wilcoxon signed-rank test. Comparisons across multiple groups were ascertained through the implementation of one-way ANOVA or the Kruskal-Wallis test. Patient survival outcomes were comparatively analyzed utilizing the Log-rank test. Utilizing both univariate and multivariate Cox regression analyses, the predictive determinants were identified. In order to establish the threshold value for the senescence score, the “surv_cutpoint” function within the “survminer” package was employed. After determining optimal senescence score cutoff value, The patients were stratified into two groups characterized by high and low senescence levels. Following this, Chi-square tests or Fisher’s exact tests were applied to compare their clinical attributes.

Collectively, 279 SRGs were identified based on the database. Two machine learning-based algorithms, namely LASSO regression analysis and SVM-RFE, were employed to analyze the input data, in which 45 and 73 SRGs were identified respectively (Figures 1A–D). 33 SRGs overlapped by result of two algorithms were used for downstream analysis (Figure 1E).

For shedding lighting on illumination of the characteristics of SRGs in STS. We explored the SCNV in pan-cancer level. Results indicated the SCNV gain in ASPH, BLVRA, CEBPB, E2F1, IGFBP3, PMVK, SENP7, SPOP, SRC, and ZMAT3 (Figure 1F). Interestingly, E2F1 was highly expressed in most cancer types (Figure 1G). TP53 and SMARCA4 showed high frequency of mutation among the 33 analyzed SRGs (Figure 1H). ASPH, CDKN2A, E2F1, P3H1, and PTTG1 were identified as prognostic indicators of unfavorable outcomes within multiple cancer types (Figure 1I). Additionally, the CNV condition and expression level of ASF1A, CDKN2A, MAGOHB, NADK, SPOP, and ARC were remarkably correlated in most cancer types (Figure 1J).

We identified mutations related to SRGs in TCGA-SARC cohort. Results demonstrated that a large number of samples present SRGs-associated mutations (altered in 100 of 237 samples) (Figure 2A). The frequency of CNV of SRGs is displayed in Figure 2B. PMVK showing the highest CNV gain while TP53 acquiring the highest CNV loss. Most SRGs locate at chromosome 1 and 17 (Figure 2C). Interaction among somatic mutations of SRGs was measured by somatic interaction function. KL and CBX8 exhibited high co-occurrence with each other (p < 0.05) (Figure 2D). We found that expressing levels of 33 SRGs were effective in discriminating between tumors and normal tissues (Figures 2E, F).

Expression pattern of SRGs was further analyzed at single-cell level using data from GSE131309 (Figure 3A). We noticed that ASPH, BHLHE40 and ITSN2 showed wide expression in different cell types while IGFBP3, CDKN2A, HK3 and MAP4K1 expressed in specific cell types (Figures 3B, C; Supplementary Figure S1). Consistently, expression of MAGOHB and E2F1 were upregulated while CEBPB was downregulated in STS cell lines and our clinical samples (Figures 3D–I).

TME is a critical mediator of within tumor characteristics, as it involves interactions among immune cells and cancer cells and can influence their fate. We conducted an analysis to examine associations among immune signatures and expression of genetic factors. The results of the analysis showed that NINJ1, MAP4K1, HK3, and CEBPB exhibited a positive correlation with the majority of immune cell types (Figure 4A). A comprehensive correlation analysis of SRGs was also conducted, which showed the wide positive association between most SRGs. However, N3H1, NINJ1, HK3, and CEBPB were negatively correlated with most SRGs (Figure 4B).

Consensus unsupervised clustering was performed on TCGA-SARC cohorts, which resulted in the division of patients into distinct groups according to the expressing levels of 33 SRGs (Supplementary Figure S2A–F). Consequently, four clusters were identified according to the optimal clustering number, containing 117 cases within cluster C1, 70 cases within cluster C2, 41 cases within cluster C3 and 36 cases within cluster C4. They showed difference in survival rate, in which the cluster C4 displaying best prognosis (Figure 4C). The four clusters had distinct expression patterns of 33 SRGs (Figure 4D). Further, we applied the GSVA analysis to compare the clusters enriched in different clusters. As a result, cluster C2 was found to be exhibit negative enrichment within pathways related to cytokine-cytokine receptor interaction, complement and coagulation cascades, lysosome, sphingolipid metabolism, and amino sugar and nucleotide sugar metabolism (Figures 4E, F; Supplementary Figure S2G–J). Interestingly, cluster C2 showed significant downregulated adaptive immune cell infiltration (Figure 4G).

To further reveal the characteristics of different clusters, we calculated the senescence-related DEGs, only four overlapping DGEs were identified among different clusters (Figure 5A). Unsupervised consensus clustering was carried out by utilizing DEGs, for identifying subtypes associated with senescence (Supplementary Figure S3A–F). Consequently, the study identified three distinct subtypes, denoted as S1, S2, and S3, containing 72, 137, 48 patients respectively. They showed significant difference in outcomes and heterogeneity in clinical features (Figures 5B, C). Moreover, the GSVA analysis indicated that subtype S3 was positively enriched in sugar metabolism while negatively enriched in DNA repair pathways such as the DNA replication, mismatch repair, cell cycle, homologous recombination (Figure 5D; Supplementary Figure S3G).

To establish an individualized predictive model for senescence, we subsequently calculated the senescence score based on senescence related DEGs (Figure 6A). As expected, senescence-related clusters showed significant differences in senescence score (Figure 6B). A marked disparity in tumor mutation burden was observed between the high and low senescence scoring cohorts (Figure 6C). Notably, individual exhibiting high senescence scores exhibited unfavorable prognosis, as validated in three datasets (TCGA-SARC, GSE39055, GSE176307) (Figures 6D–F). The TME scores, encompassing stromal, immune, and ESTIMATE scores, demonstrated marked variation between the high and low senescence score cohorts (Figure 6G). Interestingly, we found a positive correlation among the senescence scoring and tumor mutation burden (Figure 6H). Distinct clinical characteristics were observed between the high and low senescence score groups, including age (p = 0.048) and histology (p < 0.001) (Figure 6I). Moreover, a positive correlation was discovered between senescence score and several immune cell subtypes, including activated CD4⁺ T cell, γδT cell, MDSC and macrophage (Figure 6J). In addition, the multivariate Cox regression analysis demonstrated that the senescence scoring exhibited significant prognostic value as a unfavorable prognostic indicator within STS (Figure 6K; Supplementary Figure S4A–I).

We subsequently depicted the differences of genomic and transcriptomic profiles between high and low senescence scores groups. As illustrated, a lower mutation frequency was noticed within the cohort exhibiting higher senescence score, wherein alterations within 129 (66.15%) of 195 individuals (Figures 7A, B). Notably, the arm level amplification was more prevalent in high senescence score cohort while depletion was more common in low senescence score group (Figure 7E). GSVA analysis indicated the positive enrichment of unfolded protein response, TNF signaling via NFKB, p53 pathway, TGFβ pathway and MYC targets in high senescence score group (Figure 7C), but the pathways of myogenesis was negatively enriched (Figure 7C). Furthermore, we analyzed senescence score’s association with immunotherapy-related pathway and cancer immunity cycles. Senescence score was positively correlated with several subtypes of immune cells including CD8⁺ T cell, dendritic cell, macrophage and NK cell (Figure 7D). Meanwhile, a positive correlation between senescence score and several immunotherapy-related pathways was identified, including IFN-γ signature, APM signal and Proteasome (Figure 7D).

To examine the association among senescence scoring system and response to immunotherapy, we analyzed a cohort of patients who had received immunotherapy. The findings revealed an unfavorable prognosis among patients exhibiting elevated senescence scores (p = 0.013) and a higher non-response rate when treated with immunotherapy (p < 0.001) (Figures 7F, G). We further screened the GDSC database and noticed a statistically significant increase in IC50 values for CCT007093, KIN001.135, and Lapatinib among the group with high senescence scores (Figure 7H).

As aforementioned, E2F1 showed high expression in STS, thus we further explored its role in STS. Our findings revealed that demonstrated that the downregulation of E2F1 significantly suppressed STS proliferation, as evidenced by CCK-8 assay and clone formation (Figures 8A–C). Moreover, scratch assay demonstrated that suppression of E2F1 significantly impeded the migration of STS cell lines (Figure 8D).