The UCSC Xena browser (https://xenabrowser.net/datapages/) was used to download gene expression profiles of The Cancer Genome Atlas (TCGA)—sarcoma (SARC) cohort and normal tissues in the Genotype-Tissue Expression (GTEx) dataset (Goldman et al., 2020). FPKM values of RNA sequencing (RNA-Seq) data from TCGA and GTEx were normalized through log₂(FPKM+1) transformation. RNA-Seq data from two database were then processed and unified following sufficiently rigorous procedures, including the uniform realignment, the quantification of gene expression and the correction of batch effect (Wang et al., 2018). Clinical information of TCGA-SARC cohort has been made available for download at cBioPortal (https://www.cbioportal.org/) (Cerami et al., 2012). Within TCGA-SARC cohort, a total of 259 patients with STS were screened, composed of 104 patients having leiomyosarcoma (LMS), 59 patients having dedifferentiated liposarcoma (DDLPS), 49 patients having undifferentiated pleomorphic sarcoma (UPS), 25 patients having myxofibrosarcoma (MFS) and 22 patients having other STS. In GTEx, gene expression profiles of 911 normal human adipose and muscle were integrated with that of TCGA-cohort owing to the lack of normal tissues in TCGA. Additionally, the RNA-Seq profiles with clinical characteristics of GSE30929 were available in the GEO data repository for further validation.

In total, 37 PRGs were identified based on previous studies (Kang et al., 2014; Zhu et al., 2017; Feng et al., 2018; Tsuchiya et al., 2019; Xue et al., 2019; Zhang et al., 2020; Zheng and Kanneganti, 2020), and were listed in Supplementary Table S1. Differential expression of PRGs between SARC and corresponding normal tissue was conducted, utilizing the “limma” R package (Version 3.48.3). In clusters based on consensus clustering analysis of PRGs, DEGs between cluster one and cluster two were identified if | log₂ (fold change) | > 2 in expression value and false discovery rate (FDR) < 0.05. We established correlation network of PRGs using “corr” R package (Version 0.4.3) with the correlation coefficient set of 0.4. Somatic mutation of PRGs were visualized utilizing “Maftools” R package (Version 2.8.0). Circos plot was present to demonstrate the location of PRGs by “Circos” R package (Version 1.2.1).

The PPI network of PRGs was constructed with the minimum required interaction score set of 0.9 to ensure high confidence, utilizing STRING database (https://string-db.org/) (Szklarczyk et al., 2021). Moreover, PPI network of PRGs was further analyzed using Cytoscape software (version 3.8.2). Hub genes were then screened based on MCODE, with the indicators set as degree cutoff = 2, node score cutoff = 0.2, K-core = 2 and maximum depth = 100.

The “ConsensusClusterPlus” R package (Version 1.56.0) was introduced to conduct consensus clustering analysis, so as to identify pyroptosis-related subtypes of STS. Because k-means clustering analysis was stochastic, repetitions was set to 1,000 to ensure stable clustering (Wilkerson and Hayes, 2010). Differences in survival among clusters were visualized with Kaplan-Meier (KM) plots based on the R packages of “survival” (Version 3.2–11) with “survminer” (Version 0.4.9).

DEGs (| log₂ (fold change) | > 2 in expression value and FDR <0.05) between cluster one and cluster two based on pyroptosis-related consensus clustering were collected. With the “clusterProfiler” R package (Version 4.0.4) introduced, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis were performed.

The prognostic significance of DEGs within different clusters was evaluated by utilizing univariate COX regression analysis. In order to narrow down gene selection, DEGs with significant impact on survival (p < 0.01) were subsequently incorporated into the least absolute shrinkage and selection operator (LASSO) Cox regression analysis based on the “glmnet” R package (Version 4.1–2). The risk scores were further established according to the formula: ∑ i n X i ∗ Y i (Xi: coefficients of the gene i, Yi: expression values of the gene i).

The risk scores were subsequently discretized to divide TCGA-SARC cohort into high- and low-risk groups. Survival between two risk groups was analyzed with KM plot in TCGA-SARC cohort. For further external validation of risk scores, gene expression profiles of GSE30929 were entered into the formula, and we sorted these patients into high- and low-risk groups according to the same cut-off point. The “timeROC” R package (Version 0.4) was introduced for evaluating the predictive accuracy. Furthermore, the risk scores integrated with clinical characteristics including age, sex, race, histology, tumor site, tumor multifocality, surgical margin, tumor depth and radiotherapy were analyzed using multivariate COX regression analysis. The nomogram was performed to illustrate the prognostic model, which was further evaluated by the calibration curve.

The 16 immune cells infiltration and 13 related functions were quantified through ssGSEA in different pyroptosis-related risk groups of TCGA-SARC cohort and GSE30929, by utilizing the R package of “GSVA” (Version 1.40.1) and “GSEABase” (Version 1.54.0). The correlation between expression of DEGs in the gene signature and immune infiltrates was analyzed by the Tumor Immune Estimation Resource 2.0 database (TIMER2.0) (Li et al., 2020).

The human synovial sarcoma cell line (SW-982) was purchased from the American Type Culture Collection (ATCC). The human skin fibroblast cell line (HSF) with related media were purchased from Fenghui Biotechnology Co., Ltd (Hunan, China). The primary human synovial sarcoma cells (hSS-005R) were also established for the validation of PRGs. The human synovial sarcoma cells SW-982 and hSS-005R were cultured in Dulbecco’s modified Eagle medium (DMEM) (Gibco, United States) supplemented with 10% fetal bovine serum (FBS) (Gibco, United States) and 1% penicillin-streptomycin (NCM Biotech, China). Cells were cultured at 37°C with 5% CO₂ in a humidified incubator (Thermo Fisher Scientific, United States).

Total cellular RNA was extracted using the RNA Express Total RNA Kit (M050, NCM Biotech, China). For cDNA synthesis, reverse transcription was conducted with the RevertAid First Strand cDNA Synthesis kit (K1622, Thermo Fisher Scientific, United States). Subsequently, RT-qPCR was performed on the StepOne Plus (Applied Biosystems, United States) by utilizing SYBR Green qPCR Master Mix (2×) (Bimake, United States). The primers used for the RT-qPCR were listed in Table 1.

Statistical Analysis were conducted by utilizing R (Version 4.1.0). Differential gene expression between two groups was identified using Wilcoxon rank sum test, with p value calculated for each gene. Spearman’s correlation test and matrix were conducted to compare gene expression with each other. Survival differences were compared utilizing log-rank test with KM curve. The χ2 test or Fisher’s exact test was introduced to evaluate clinical characteristics between high- and low-risk groups. COX regression analysis was conducted to identify prognostic factors, and hazard ratio (HR) with 95% confidence interval (CI) were also computed. Statistical difference of p < 0.05 was defined significant.

The study design was illustrated in Figure 1. Totally 37 PRGs were included to detect differential expression between STS and normal tissues from TCGA-SARC and GTEx dataset (Supplementary Table S1). We found that 34 PRGs were differentially expressed (p < 0.05), among which 14 genes (CASP3, CASP5, CASP6, DHX9, GSDMA, GZMA, GZMB, IL1B, NLRC4, NLRP3, NLRP7, NOD2, PYCARD, TNF) were upregulated and 20 genes (APIP, CASP4, CASP8, CASP9, ELANE, FOXO3, GPX4, GSDMB, GSDMC, GSDME, IL18, IL6, NLRP1, NLRP2, NLRP6, NOD1, PJVK, PLCG1, PRKACA, SCAF11) were downregulated in STS group (Figure 2A, Supplementary Figure S1A). To validate PRGs in related cell lines, we performed RT-qPCR analysis (Figure 2B). The expression levels of several key PRGs including CASP3, IL1B and DHX9 were significantly higher in the human synovial sarcoma cells SW-982 and hSS-005R, compared with those in the human skin fibroblast cell line (HSF). The correlation network of 34 PRGs was present in Supplementary Figure S1B and Supplementary Figure S1C. The chromosome location of 34 PRGs was illustrated in Figure 2C. We then analyzed copy number variations (CNVs) of PRGs (Figure 2D). It could be found that 15 of 237 (6.33%) SARC samples displayed pyroptosis-related mutations. And the majority of mutations were missense mutations. The PPI network was subsequently established and seven hub genes (NLRP1, IL18, NLRC4, NLRP3, PYCARD, CASP5, IL1B) were identified (Figure 2E). Furthermore, these 34 PRGs were efficient to discriminate STS and normal tissues on the expression level (Figure 2F).

To elucidate different STS subtypes and corresponding clinical characteristics and prognosis, the TCGA-SARC cohort was clustered into two distinct clusters based on PRGs through consensus clustering analysis (Figure 3A, Supplementary Figure S2A–F). Within pyroptosis-related cluster 1, There were 135 STS patients and 124 STS patients were in pyroptosis-related cluster 2. Remarkably, overall survival (OS) curves of these two clusters indicated significantly great survivorship difference (p = 0.019, Figure 3B). In Figure 3C, the gene expressing level of PRGs in two distinct clusters were displayed.

According to the stringent selecting criterion of | log₂ (fold change) | > 2 in expression value and FDR <0.05, a total of 577 DEGs were identified between pyroptosis-related cluster one and cluster 2 (Figure 4A). There were also significant differences in clinicopathological characteristics including age, histology, metastatic status and survival between two distinct clusters (p < 0.05).

Subsequently, these 577 DEGs were subjected to the GO and KEGG enrichment analysis, in order to reveal biological processes and mechanisms of pyroptosis-related clusters. GO enrichment analysis indicated that DEGs were predominantly enriched in immune response-activating cell surface receptor signaling pathway, immune response-activating signal transduction, external side of plasma membrane and antigen binding (Figure 4B). Moreover, these DEGs were also significantly enriched in cytokine-cytokine receptor interaction, hematopoietic cell lineage and cell adhesion molecules (Figure 4C).

The prognostic significance of 577 DEGs between cluster one and cluster two were analyzed by utilizing univariate COX regression analysis. Accordingly, 42 genes were preserved based on the strict criteria (p < 0.01) and processed for subsequent analysis (Figure 5A). The LASSO COX regression analysis was then performed, and five key genes were eventually identified with the pyroptosis-related gene signature constructed (Figures 5B,C). The risk score = (-0.05167*CTSG exp.) + (0.06184*DUSP9 exp.) + (-0.02483*CLEC10A exp.) + (-0.02979*CPA3 exp.) + (-0.06014*CD1C exp.). Based on the median of the risk scores in TCGA-SARC cohort, 259 STS patients were divided into the low-risk group (n = 130) and the high-risk group (n = 129) (Figure 5D, Supplementary Figure S3A). Principal component analysis (PCA) and t-distributed stochastic neighbor embedding (t-SNE) indicated that low-risk and high-risk group could be clearly distinguished (Supplementary Figure S3C, Supplementary Figure S3E). The KM plot of OS rate demonstrated significant difference between two risk scores groups of TCGA-SARC cohort (p < 0.001, Figure 5E). Time-dependent receiver operating characteristics (ROC) curves were introduced for assessing model performance, and the area under curve (AUC) of 1-year, 3-year and 5-year OS rate were 0.683, 0.668 and 0.690, accordingly (Figure 5F).

For further validating the accuracy of the pyroptosis-related gene signature, corresponding data of GSE30929 was retrieved and the risk score of was calculated respectively. According to the same cut-off point of TCGA-SARC cohort, GSE30929 cohort was divided into different risk groups (Figure 5G, Supplementary Figure S3B). PCA and t-SNE also illustrated optimal degree of discrimination between high- and low-risk groups of GSE30929 (Supplementary Figure S3D, Supplementary Figure S3F). Remarkably, the disease-free survival (DFS) of two distinct risk groups demonstrated significant discrepancy (p < 0.001, Figure 5H). AUC were 0.679 for 1-year, 0.668 for 3-year and 0.640 for 5-year (Figure 5I).

The potential clinical utility of pyroptosis-based risk score was further investigated. Clinical characteristics of gender, age, tumor histology and tumor site between high- and low-risk groups were visualized in Figure 6A. Alluvial diagram also illustrated the relationship of pyroptosis-based cluster distribution, clinical characteristics, different risk groups and survival outcomes (Figure 6B). Furthermore, multivariate COX regression analysis integrated with clinical characteristics and pyroptosis-based risk score were performed to establish the prognostic model (Figure 6C, Supplementary Table S2). Based on the established prognostic model, a novel nomogram was subsequently constructed for predicting the survival probability of STS patients (Figure 6D). The 3-year and 5-year OS rate have proven to be relatively well predicted by the calibration curve of the nomogram (Figure 6D).

In order to compare the immune activity, the ssGSEA was applied to analyze the immune infiltration and functions of high- and low-risk groups. The ssGSEA score calculated could indicate the infiltration degrees of immune cells and pathways within TCGA-SARC cohort and GSE30929 cohort. In TCGA-SARC cohort, the infiltration degrees of CD8⁺ T cell, dendritic cell (DC), immature DC (iDC), macrophage, mast cell, neutrophil, natural killer (NK) cell, plasmacytoid DC (pDC), T helper (Th) cell, T follicular helper (Tfh) cell, Th1 cell, Th2 cell, tumour-infiltrating lymphocyte (TIL) and regulatory T cell (Treg) were significantly lower in the high pyroptosis-related risk group (p < 0.05, Figure 7A). All 13 related immune functions were also significantly decreased in the high pyroptosis-related risk group (p < 0.05, Figure 7C). GSE30929 dataset was subsequently included to analyze the immune activity. The results were similar with the majority of immune infiltration and function significantly decreased in the high-risk group (Figure 7B, Figure 7D). The correlation of key DEGs in the gene signature and immune infiltration was also illustrated (Supplementary Figure S4A–E).