Soft tissue sarcomas are rare malignancies arising in the connective tissues that surround and support other body structures, including fat, muscle, blood vessels, and nerves. They collectively account for less than 1% of all malignant tumors, with an incidence of roughly 1 per 50,000 population per year. Despite this rarity, the WHO classification recognizes more than 50 distinct histological subtypes, each with its own prognosis and preferred treatment regimen. The clinical challenge is that many of these subtypes look nearly identical under the microscope, making differential diagnosis both difficult and consequential.
The morphological overlap problem: Synovial sarcoma (SS) and malignant peripheral nerve sheath tumor (MPNST) are a classic example. They can be morphologically indistinguishable, and their immunohistochemical profiles frequently overlap. Confirming a synovial sarcoma currently requires molecular testing for the SS18-SSX translocation, a process that is laborious and time-consuming. Similarly, myxofibrosarcoma (MFS) and undifferentiated pleomorphic sarcoma (UPS) share significant molecular and morphological features, and accurately grading leiomyosarcoma carries major prognostic implications that standard histological assessment alone does not fully resolve.
The opportunity in genomics data: Over the past decade, large-scale open-access transcriptome sequencing datasets have become available for sarcomas. The Cancer Genome Atlas (TCGA) soft tissue sarcoma cohort, with 206 revised cases carrying complete clinical data, is the largest publicly available resource of this kind. The authors recognized that machine learning algorithms, which can detect non-obvious patterns across thousands of gene expression variables simultaneously, were ideally suited to extract diagnostic, prognostic, and therapeutic insights from these datasets in ways that conventional statistical approaches could not.
The study applies four distinct computational strategies: a deep neural network to compare sarcoma transcriptomes to normal tissue, a random forest classifier to identify subtype-distinguishing diagnostic genes, a Kaplan-Meier-plus-k-nearest neighbor approach to identify prognostic genes, and a regulatory network reconstruction paired with Connectivity Map (CMAP) drug data to identify repurposable therapies. The results are then validated in independent cohorts and in laboratory cell line experiments.
Three gene expression datasets were integrated. The primary dataset was TCGA RNA-seq count data (downloaded February 2018 from the NIH GDC portal): 206 revised soft tissue sarcoma cases comprising soft tissue leiomyosarcoma (STLMS, n=53), uterine leiomyosarcoma (ULMS, n=27), dedifferentiated liposarcoma (DDLPS, n=50), undifferentiated pleomorphic sarcoma (UPS, n=44), myxofibrosarcoma (MFS, n=17), synovial sarcoma (SS, n=10), and MPNST (n=5). STLMS was the most common subtype and included FNCLCC grade 1 (n=11), grade 2 (n=35), and grade 3 (n=7) cases. The second dataset was Genotype-Tissue Expression (GTEx) project version 7, containing 9,662 normal tissue samples across 31 tissue types, used as the reference for the deep neural network. The third was a French Sarcoma Group expression array from GEO (accession GSE21050) with DDLPS (n=62) and LMS (n=84) samples, used for independent prognostic validation.
Normalization and filtering: Across all datasets, genes with low expression (RNA-seq counts below 2 counts per million in all samples, or expression array values below 2) were removed. Remaining data were normalized using the weighted trimmed mean of M-values (TMM) via Limma R package (v3.6), then log2-transformed. Where datasets were merged, joint normalization was applied to minimize batch effects.
Machine learning implementations: The deep neural network used TensorFlow (v1.6) and Keras (v2.1.4) with one converging hidden layer. t-SNE clustering used the Rtsne package at perplexity=60 and theta=0.5. Random forest analysis used the randomForest R package (v4.6) with Breiman's algorithm, with feature importance quantified by the Gini index. K-nearest neighbor analysis used the Caret package (v6.0) with repeated cross-validation, testing k values from 1 to 30. Connectivity Map analysis used the expression2kinase pipeline to infer regulatory transcription factors (ChEA database) and kinases from differentially expressed genes, then matched them to drugs in the CMAP database.
Validation approaches: Diagnostic markers from the random forest were validated in an independent cohort of 9 frozen tissue samples (5 SS and 4 MPNST, all molecularly confirmed) using quantitative RT-PCR normalized to housekeeping gene HPRT1. Prognostic genes were validated in the French Sarcoma Group cohort using k-NN trained on one half and tested on the other. HMMR protein was validated by immunohistochemistry on a tissue microarray (TMA) of 70 leiomyosarcoma cases scored by automated ImageJ color deconvolution and by a blinded pathologist. HDAC inhibitor sensitivity was tested in cell viability assays using three LMS lines (JA192, LMS04, LMS05) and one SS line (SYO-1).
To understand the differentiation state and likely tissue of origin of each sarcoma subtype, the authors trained a deep neural network on the GTEx normal tissue data. Principal components were computed across all 9,662 GTEx samples, and those values were used to train the neural network, which achieved a tissue-type prediction accuracy of 98% on the held-out normal tissue test data. The trained network was then applied to principal components from the 206 TCGA sarcoma samples, asking which normal tissue each tumor's expression profile most resembled.
Key neural network findings: Uterine leiomyosarcoma showed the expected overlap with uterine and cervical tissue. Soft tissue leiomyosarcoma overlapped with blood vessel expression patterns, consistent with the hypothesis that a subset of STLMS arise from small to medium-sized veins. The most biologically informative finding was that both MPNST and SS showed expression profiles most similar to neural tissue (brain and nerve). MPNST's neural identity was expected given its origin from peripheral nerve sheaths, but SS's neural signature was a new molecular confirmation of prior speculation about neural differentiation in this tumor, whose cell of origin has historically been unclear. MFS and UPS showed unexpectedly large overlap with adipose tissue, more so than DDLPS itself, possibly because TCGA DDLPS samples were enriched for the dedifferentiated (less adipocytic) component.
A parallel differentially expressed gene analysis identified 331 to 7,784 DEGs per subtype (Benjamini-Hochberg adjusted p < 0.05, logFC > 0), with DDLPS having the most (7,784) and STLMS the fewest (331). Gene Ontology enrichment via EnrichR confirmed biological coherence: STLMS DEGs enriched for muscle development GO terms, MPNST DEGs for neuronal development, and ULMS for cell cycle processes rather than muscle differentiation, suggesting proliferation dominates over differentiation in the uterine variant. DDLPS, UPS, and MFS shared many GO terms including immune system pathways, likely reflecting shared inflammatory infiltration in these tumors.
The t-SNE analysis on all 206 TCGA samples identified three natural molecular clusters in unsupervised fashion: (1) STLMS and ULMS, (2) SS and MPNST, (3) DDLPS, UPS, and MFS. MFS and UPS had the largest DEG overlap at 1,201 shared genes, followed by STLMS and ULMS at 210. The unsupervised clustering required no clinical labels and yet recapitulated exactly the three groupings that clinicians already recognize as the most morphologically challenging diagnostic boundaries in sarcoma pathology.
Once the three molecular groups were established, the authors trained separate random forest classifiers to distinguish the subtypes within each group. For each comparison, differentially expressed genes (adjusted p < 0.05) served as input features; samples were divided randomly into training and test sets, and classifiers were evaluated on held-out test data. Variable importance for each gene was quantified by the Gini index, which measures how much a given variable reduces impurity across all decision tree splits.
Prediction accuracy: Random forests distinguished STLMS from ULMS, SS from MPNST, and DDLPS from UPS/MFS with accuracy exceeding 95% in each comparison. The single exception was the pairwise comparison between UPS and MFS, where the forest achieved 88% accuracy, consistent with the 1,201 shared DEGs between these entities. This mirrors clinical experience: pathologists often find the UPS-MFS boundary among the most difficult calls in soft tissue pathology.
Top diagnostic genes by group: For group 1 (STLMS vs. ULMS), the highest-importance genes were HOXA11 and its antisense RNA HOXA11-AS, both established regulators of uterine development and homeostasis, making their higher expression in ULMS biologically coherent. For group 2 (SS vs. MPNST), the top genes were NEURL1 and NPAS1 (upregulated in SS), reflecting neural differentiation, and SCD (stearoyl-CoA desaturase, upregulated in MPNST), an enzyme in fatty acid biosynthesis whose differential expression may relate to the 48% rate of chromosome 10 deletions in MPNST, which removes the NEURL1 locus. For group 3, MDM2 and CDK4 appeared as expected DDLPS markers (reflecting 12q13-15 amplification), along with neighboring amplified genes FRS2, TSPAN31, and CTDSP2. JADE2 had the highest importance score for UPS vs. MFS but with overlapping expression distributions, consistent with the limited separability between these two entities.
qRT-PCR validation: The top three SS-vs-MPNST markers (NEURL1, SCD, NPAS1) were validated in an entirely independent cohort of 9 frozen tissue samples (5 SS, all SS18-SSX translocation-positive; 4 MPNST, all confirmed by nerve association, NF1 status, or H3K27me3 loss) using quantitative RT-PCR normalized to HPRT1. Expression patterns matched the TCGA-derived predictions in both direction and magnitude, confirming the diagnostic utility of these three-gene panel in independent tissue.
For all annotated sarcoma subtypes except MPNST (only 5 samples available), the authors screened all 24,168 expressed genes for prognostic significance. For each gene, the optimal expression cutoff was calculated using the maxstat R package (v0.7) applying the Hothorn-Lausen method, which finds the cutoff maximizing the log-rank statistic for disease-free interval (DFI), with p-values derived from the Streitberg algorithm. DFI, defined as time to local recurrence or distant metastasis, was the primary outcome because it had been validated previously as one of the most informative endpoints in soft tissue sarcoma studies.
Identified prognostic genes: In total, 429 genes met the threshold of p < 0.001 for DFI. The distribution was heavily subtype-specific: SS had the most with 166 prognostic genes, followed by STLMS (74 genes) and ULMS (34 genes), with smaller counts in DDLPS, UPS, and MFS. There was almost no cross-subtype overlap. Only two genes (KLF6 and MT1F) were shared between UPS and SS, and one gene (NPM2) between ULMS and MFS. No genes were shared between STLMS and ULMS, despite both being leiomyosarcomas, reinforcing that uterine and soft tissue LMS are biologically distinct diseases with different clinical behavior. Of the 429 genes, 228 had previously been reported as prognostic in non-sarcoma tumors in the Protein Atlas, and 201 were entirely novel.
Comparison to CINSARC: The CINSARC prognostic signature uses 67 genes associated with genome complexity and cell proliferation, designed as a pan-sarcoma prognosticator. Only one gene, CDCA3, overlapped between the 429 subtype-specific genes identified here and the CINSARC set. The authors attribute this divergence to design differences: CINSARC targets general proliferation, while this study's genes are subtype-specific and incorporate time-to-event DFI data. In direct comparison in the French Sarcoma Group validation cohort, the subtype-specific k-NN approach produced substantially better p-values: LMS k-NN p=0.045 vs. CINSARC p=0.24; DDLPS k-NN p=0.02 vs. CINSARC p=0.14; UPS k-NN p=0.012 vs. CINSARC p=0.038.
For the French Sarcoma Group cross-validation, genes significant in both the TCGA (DFI, p < 0.05) and the first French Sarcoma cohort (metastasis-free interval, p < 0.05) were carried forward. The k-NN algorithm was trained on the first half of the French cohort and tested on the second. Key cross-validated prognostic genes were: for LMS, HMMR, MXD4, and BRCA2; for DDLPS, KLF6; for UPS, PCMTD2, TNXA, TMEM65, and SNRNP48. These represent the strongest candidates for clinical prognostic biomarker development.
Among all computationally identified prognostic genes, HMMR (hyaluronan-mediated motility receptor, also known as RHAMM) was selected for independent protein-level validation in leiomyosarcoma. HMMR was significant for DFI in the TCGA data (p < 0.05) and also significant for metastasis-free interval in the French Sarcoma Group cohort. HMMR had previously been reported to associate with poor survival in liver, pancreatic, and lung cancer, and forms a protein complex with BRCA1 or BRCA2 and other DNA repair factors involved in homologous recombination.
Tissue microarray design and scoring: IHC for HMMR protein was performed on two TMA cohorts totaling 70 leiomyosarcoma cases: 32 from a previously described cohort and 38 newly constructed. The 70 cases included 43 females and 27 males with a mean diagnosis age of 62 years; 5 were uterine LMS and the rest soft tissue LMS, graded FNCLCC grade 1 (n=10), grade 2 (n=23), grade 3 (n=31), with 1 ungraded. All cases were stained simultaneously with a polyclonal rabbit anti-HMMR antibody (Sigma-Aldrich HPA040025) at 1:1000 dilution to ensure comparability. Automated scoring used ImageJ color deconvolution to separate hematoxylin and DAB channels; a cut-off DAB intensity score of 20 defined high vs. low expression. The second cohort was additionally scored manually by a blinded pathologist using intensity categories (1=weak, 2=moderate, 3=strong), and manual and automated scores showed concordance.
Outcome results: High HMMR protein expression was associated with significantly shorter disease-free interval in the full 70-case cohort (p = 0.0061). In multivariate Cox regression incorporating FNCLCC histological grade, HMMR staining remained an independent predictor (p = 0.0039) while FNCLCC grade did not reach significance (p = 0.285). This is clinically important: FNCLCC grading is the established standard for sarcoma prognosis, and HMMR's superior performance in multivariate analysis indicates it captures biological behavior that histological grade alone misses.
The biological mechanism likely relates to the recently recognized "BRCAness" signature in LMS, which involves defects in homologous recombination repair. Two of the top LMS prognostic genes from the genome-wide screen, HMMR and BRCA2, are both involved in homologous recombination. Overexpression of HMMR may represent a compensatory upregulation in tumors with defective repair pathway components, making it both a marker of repair pathway dysfunction and a potential therapeutic vulnerability. The authors suggest that the BRCAness connection warrants further study for LMS-specific treatment strategies.
To identify repurposable drugs, the authors applied a regulatory network reconstruction pipeline to each sarcoma subtype's signature genes. Using expression2kinase (maayanlab.net/X2K), differentially expressed genes were input into the ChIP-seq Enrichment Analysis (ChEA) database to infer which transcription factors most likely regulate those genes. Kinase Enrichment Analysis then identified upstream kinases. These transcription factors and kinases were matched to compounds in the Connectivity Map (CMAP) database, which links small molecule perturbation signatures to gene expression changes, to identify drugs predicted to reverse or recapitulate the sarcoma expression state.
Identified compounds: Doxorubicin, the standard first-line systemic therapy for soft tissue sarcomas, was identified as a potential effective therapy for most subtypes, serving as an internal positive control that validates the methodology. Trichostatin A (TSA), a pan-HDAC inhibitor, was predicted to be effective across all soft tissue sarcoma subtypes tested. A second HDAC inhibitor, Vorinostat, was specifically predicted for UPS and ULMS. Tanespimycin, an HSP90 inhibitor in clinical trials for solid tumors, was predicted for UPS, ULMS, and MPNST. The chord diagram visualization showed TSA and doxorubicin had the broadest subtype connections, while tanespimycin and vorinostat were more selective.
Cell viability assays: While HDAC inhibitor sensitivity in SS was established (HDAC inhibition disrupts the SS18-SSX oncoprotein leading to apoptosis), LMS sensitivity had not been systematically studied. The authors tested TSA and quisinostat (a more potent second-generation HDAC inhibitor) on three LMS lines (JA192, an STLMS line; LMS04, a ULMS line; LMS05, an STLMS line) and SYO-1 (SS, positive control). IC50 values for TSA ranged from 39 nM (SYO-1) to 474 nM (JA192), with LMS lines intermediate (LMS04: 229 nM; LMS05: 178 nM), reflecting greater TSA variability across lines.
Quisinostat showed uniformly low IC50 values across all four lines: SYO-1 at 15 nM, JA192 at 41 nM, LMS04 at 34 nM, and LMS05 at 39 nM. The narrow 15-41 nM IC50 range, consistent across ULMS and both STLMS lines, indicates that quisinostat's activity does not depend on subtype anatomy within LMS. This consistency, combined with the CMAP prediction from gene expression analysis, makes quisinostat a strong candidate for further preclinical investigation in LMS and for potential early-phase clinical trials in this historically chemoresistant disease subtype.
The most significant structural limitation is the extreme imbalance in sample sizes across subtypes within the TCGA cohort. MPNST had only 5 samples, which precluded prognostic gene identification entirely. SS had only 10 samples for the random forest training, making those results less statistically stable. LMS with 80 combined cases and DDLPS with 50 cases produced the most reliable findings. This imbalance reflects the underlying disease epidemiology, since these are genuinely rare tumors, but it means that conclusions about the smallest subtypes require independent replication before clinical application.
Cross-validation gaps for prognostic genes: Prognostic gene cross-validation in the French Sarcoma Group cohort was only possible for LMS, DDLPS, and UPS, because these were the only subtypes with both sufficient case numbers and available follow-up data in the external dataset. SS, despite having the most identified prognostic genes (166), had no available independent expression dataset with adequate clinical follow-up. There was also a systematic mismatch between TCGA outcome data (DFI, including loco-regional recurrence) and French Sarcoma Group data (metastasis-free interval only), meaning genes with prognostic value specifically for local recurrence risk may have been filtered out at the cross-dataset significance step.
Bulk RNA-seq and tumor microenvironment: All transcriptome data used in this study came from bulk RNA-seq, which averages expression across all cells in a tumor sample, including immune cells, stromal fibroblasts, endothelial cells, and normal adipocytes. Some tissue-type similarity signals detected by the deep neural network may therefore reflect microenvironmental composition rather than intrinsic tumor biology. The authors specifically note this limitation for the MFS and UPS overlap with adipose tissue, and for SS's salivary gland similarity (the latter potentially explained by the 2 biphasic SS cases in the cohort). Single-cell RNA sequencing and data from the Human Cell Atlas will be needed to definitively map tissue of origin in morphologically ambiguous sarcoma subtypes.
Translational gaps for clinical use: The random forest-identified diagnostic markers (NEURL1, SCD, NPAS1 for SS vs. MPNST) were validated only at the mRNA level by qRT-PCR. Clinical pathology requires IHC-ready protein markers applicable to standard formalin-fixed paraffin-embedded tissue. The HMMR result demonstrates that this translation is achievable from computational discovery to validated protein marker, but NEURL1 and NPAS1 require antibody development and FFPE-validation studies. The HDAC inhibitor findings, though compelling at the cell line level with quisinostat IC50 values of 15-41 nM across LMS lines, require in vivo testing in xenograft or patient-derived xenograft models and ultimately prospective clinical trials before they can change practice. The authors propose quisinostat's consistent sub-50 nM IC50 across three independent LMS lines justifies advancement to mouse models and early-phase trials.