Leiomyosarcoma (LMS) is a rare and highly aggressive soft tissue sarcoma that arises from cells with smooth muscle differentiation. It can appear in the uterus, retroperitoneum, gastrointestinal tract, or extremities, and each anatomical site carries its own clinical behavior and prognosis. Despite multimodal treatment strategies combining surgery, radiation, and chemotherapy, advanced LMS carries a poor prognosis due to high recurrence rates and substantial resistance to conventional agents. The broader context is sobering: in 2024 alone, an estimated 2,001,140 new cancer cases and 611,720 cancer-related deaths were expected in the United States, placing the disease as the second leading cause of death after heart disease, and LMS represents one of the more aggressive subsets within this landscape.
The precision medicine gap: Innovations such as immune checkpoint blockade and cell-based therapies have not yet achieved meaningful efficacy in LMS. Objective response rates to PD-1/PD-L1 inhibitors in LMS fall below 10%, and PD-L1 expression is low across most tumors. This limited immunotherapy response reflects the significant heterogeneity and immune evasion mechanisms that characterize LMS. A major bottleneck is the scarcity of LMS-specific clinical cohorts and multi-omics datasets, which has historically limited understanding of the molecular drivers of disease progression.
The stemness angle: This study takes a distinct approach by focusing on tumor stemness, the degree to which malignant cells resemble undifferentiated stem cells. The stemness index (mRNAsi) is a quantitative metric derived using a one-class logistic regression (OCLR) machine learning model trained on human pluripotent stem cell data. A high mRNAsi suggests that tumor cells have reverted to a dedifferentiated state associated with greater cellular plasticity, therapeutic resistance, and metastatic potential. By integrating scRNA-seq data from the GEO database with bulk RNA-seq and clinical data from the TCGA-SARC cohort, the authors sought to identify which LMS tumor cells are most stem-like, map how those cells interact with the immune microenvironment, and build a prognostic gene signature that can stratify patient outcomes.
The study is published in Frontiers in Oncology in 2025 and represents one of the most comprehensive single-cell analyses of stemness-related heterogeneity in LMS to date, identifying six candidate prognostic genes with experimental validation in the SK-LMS-1 cell line.
The study draws on two major data sources. Bulk RNA-seq data came from the TCGA-SARC cohort, which contains normalized gene expression values (TPM) for 263 sarcoma patients and 2 normal samples alongside clinical outcome data. Single-cell RNA-seq data were retrieved from the GSE212527 dataset in the Gene Expression Omnibus (GEO) database. This dataset includes three sarcoma and four leiomyosarcoma samples; the authors selected the four LMS samples for downstream analysis, yielding a final integrated dataset of 52,781 cells after quality filtering.
scRNA-seq preprocessing: All single-cell data were processed using the Seurat R package (v4.4.0). Cells with fewer than 200 detected features or present in fewer than 3 cells were excluded. Cells with mitochondrial gene expression exceeding 20% were filtered out as low-quality or dying cells. Mitochondrial genes (prefixed "MT") and ribosomal protein genes ("RPS" and "RPL") were removed from downstream analysis. Batch effects across the four LMS samples were corrected using the Harmony algorithm. Dimensionality reduction was performed using the top 50 principal components from PCA, followed by UMAP for unsupervised clustering, yielding 37 distinct cell clusters. Cell type annotation used curated marker genes from the CellMarker 2.0 database, and marker genes for each cluster were identified with the Wilcoxon rank-sum test (min.pct = 0.1, |log2FC| > 0.25, adjusted p-value < 0.05).
Identifying malignant cells with CopyKAT: A critical methodological step was distinguishing malignant tumor cells from non-malignant stromal and immune cells within the single-cell data. The authors used CopyKAT (Copy Number Karyotyping of Tumors), a Bayesian segmentation algorithm that infers genome-wide DNA copy number variations (CNVs) from scRNA-seq data at approximately 5 megabase resolution. Since aneuploidy occurs in approximately 88-90% of human cancers but is absent in normal diploid stromal cells, CNV-based classification reliably separates malignant from non-malignant populations. This step identified 16,284 aneuploid (malignant) cells, representing 30.9% of the total dataset.
Stemness index calculation: For the 16,284 aneuploid cells, the mRNAsi was computed using a one-class logistic regression (OCLR) model trained on 11,774 human pluripotent stem cell samples from the Progenitor Cell Biology Consortium (PCBC). The gelnet R package implemented the OCLR with elastic net regularization (alpha = 0.5) and lambda optimized by 5-fold cross-validation. The top 500 genes by coefficient weight were retained to define the stemness signature. Each tumor cell's mRNAsi was then calculated as the Spearman correlation between the stemness score vector and the cell's gene expression profile. Cells were classified into high and low mRNAsi groups based on the distribution of scores visualized on UMAP.
Once high and low mRNAsi groups were defined among the aneuploid tumor cells, the authors performed differential gene expression analysis using FindMarkers in Seurat (Wilcoxon test, |log2FC| > 0.5, adjusted p-value < 0.05) to identify genes distinguishing the two groups. These differentially expressed genes (DEGs) were then projected onto the TCGA-SARC bulk cohort (n = 263) for prognostic modeling, creating a bridge between single-cell biology and clinical outcome data.
Univariate Cox regression screen: All DEGs from the high vs. low mRNAsi comparison were subjected to univariate Cox proportional hazards regression in the TCGA-SARC cohort. This step screened for genes whose expression was significantly associated with overall survival. Importantly, genes upregulated in the high mRNAsi group were predominantly classified as risk-associated (hazard ratio > 1), while genes upregulated in the low mRNAsi group showed no consistent prognostic association. This asymmetry reinforces the hypothesis that high-stemness gene programs are specifically linked to worse outcomes.
Lasso regression for feature selection: Genes passing the univariate Cox screen entered a Lasso (Least Absolute Shrinkage and Selection Operator) regression with alpha = 1 (pure L1 penalty), implemented using the glmnet package. The optimal regularization parameter lambda.min was determined by 10-fold cross-validation and identified as 0.056, at which 24 genes with non-zero coefficients were retained for multivariate modeling. Lasso shrinks less informative coefficients toward zero, effectively selecting a sparse, interpretable gene set from a high-dimensional pool of candidates.
Multivariate Cox regression and risk score construction: The 24 Lasso-selected genes entered multivariate Cox regression, which identified six genes that independently contributed to prognosis: SRPK1, DSE, CTBP1, PSMD10, BOP1, and HACD4. A stemness risk score was constructed as a weighted linear combination of each gene's expression value multiplied by its multivariate Cox regression coefficient. Five of the six genes (SRPK1, DSE, CTBP1, PSMD10, BOP1) were classified as risk genes with positive coefficients, meaning higher expression correlates with worse survival. HACD4 was identified as a protective gene with a negative coefficient, meaning higher expression correlates with better survival outcomes.
Model validation: The prognostic model's discriminative performance was assessed using time-dependent ROC curves at 1, 3, and 5 years in the TCGA-SARC cohort, yielding AUC values of 82%, 77%, and 79%, respectively. Patients were stratified into high-risk and low-risk groups using the median stemness risk score as a cutoff. Kaplan-Meier survival analysis confirmed significantly poorer survival in the high-risk group. It is important to note that no independent external validation cohort was used, which represents a stated limitation of the study.
After classifying the 16,284 aneuploid cells into high and low mRNAsi groups, reclustering identified seven distinct cell clusters, with clusters 2, 3, 4, 5, and 6 exhibiting high mRNAsi levels. Differential expression analysis between the two groups revealed biologically coherent patterns that align with known features of aggressive tumor behavior.
High mRNAsi signature genes: Genes specifically upregulated in the high mRNAsi group included MKI67, MGST1, and ALOX5AP. MKI67 encodes the Ki-67 antigen, a canonical marker of cellular proliferation widely used in clinical pathology to assess mitotic activity. MGST1 is associated with endoplasmic reticulum stress responses, and ALOX5AP is linked to lipid metabolism and inflammatory signaling. Together, these markers suggest that high-stemness LMS cells are actively proliferating and engaging metabolic stress pathways, behaviors consistent with aggressive, therapy-resistant phenotypes.
Low mRNAsi signature genes: Genes preferentially expressed in the low mRNAsi group included TRDN and FXYD1, both involved in sodium-calcium channel homeostasis. These markers reflect a more differentiated smooth muscle cell identity, consistent with the hypothesis that low-stemness cells retain more of the normal tissue phenotype from which LMS arises. The functional divergence between high and low mRNAsi cells suggests that the two populations occupy opposite ends of a differentiation spectrum within the same tumor.
Gene Set Variation Analysis (GSVA): Pathway-level analysis using GSVA revealed striking enrichment of the KEGG_AXON_GUIDANCE pathway in the high mRNAsi group. This pathway integrates Rho GTPase, PI3K/AKT, and MAPK/ERK signaling cascades, all of which are well-established drivers of tumor invasion and migration. In contrast, the KEGG_DNA_REPLICATION pathway was significantly enriched in the low mRNAsi group. DNA replication pathway enrichment in low-stemness cells may reflect a more ordered, terminally dividing cell population, while the invasion-associated signaling in high-stemness cells directly links stemness to metastatic potential. These GSVA results provide mechanistic grounding for the clinical risk associations observed in the prognostic model.
To understand how stemness-related heterogeneity influences the broader tumor ecosystem, the authors applied CellChat, an algorithm that quantifies intercellular communication by modeling known ligand-receptor interactions. Cell communication strength was quantified using interaction probability scores, and statistical significance was set at p < 0.05. Differential interaction analysis compared the high and low mRNAsi aneuploid cell groups in the context of their communication with other cell types present in the LMS tumor microenvironment, including epithelial cells, fibroblasts, cancer-associated fibroblasts (CAFs), macrophages, NK cells, and other immune populations.
Overall communication complexity: The high mRNAsi group showed a significantly greater number of cell-cell interactions with epithelial cells, fibroblasts, and immune cells compared to the low mRNAsi group. Fibroblast-immune cell interactions were markedly stronger in the high mRNAsi context, while epithelial-immune cell interactions were more prominent in the low mRNAsi setting. This divergence suggests that high-stemness LMS cells remodel the stromal and immune landscape in ways that favor tumor progression, specifically by enhancing fibroblast-mediated immune suppression.
Key ligand-receptor pairs: Among the most notable interactions in the high mRNAsi group were FN1-(ITGA8+ITGB1), FN1-CD44, COL4A2-CD44, and ANGPTL1-(ITGA8+ITGB1) operating between fibroblasts and immune cells. Fibronectin (FN1) interactions with integrin complexes and CD44 promote tumor angiogenesis, cell adhesion, migration, and invasion. COL4A2 (type IV collagen) interactions with CD44 further support extracellular matrix remodeling associated with metastasis. These ligand-receptor pairs represent concrete molecular links between stromal cell communication and tumor aggressiveness in the high-stemness compartment.
TGFbeta and EGF signaling pathways: Comparative analysis of signaling flow identified significant upregulation of TGFbeta, EGF, and IL2 signaling pathways in the high mRNAsi group relative to the low mRNAsi group. Validation of ligand-receptor gene expression confirmed cell type-specific patterns: AREG and HBEGF (EGF ligands) were highly expressed in NK cells and macrophages, respectively, while EGFR showed broader expression across stromal and immune populations. For the TGFbeta pathway, TGFBR2, TGFBR1, and TGFbeta1 were highly expressed in CAFs, macrophages, and myeloid cells, with TGFBR2 particularly abundant in CAFs. Both the frequency and intensity of TGFbeta and EGF signaling were higher in the high mRNAsi group, confirming the functional activity of these pathways in the high-stemness tumor compartment.
The six prognostic genes identified through the regression pipeline are not arbitrary; each has a documented biological role in tumor progression, and their specific functions help explain why high expression correlates with poor LMS outcomes. The authors also used TIMER2.0 with the CIBERSORT algorithm (LM22 signature matrix, covering 22 human hematopoietic cell phenotypes) to assess immune infiltration associations, and GeneMANIA to predict regulatory interaction networks.
BOP1 (Block of Proliferation 1): BOP1 is a recognized oncogene across multiple cancer types, including hepatocellular carcinoma, where it promotes epithelial-to-mesenchymal transition (EMT) and correlates with poor prognosis in pan-cancer analyses. In the TCGA-SARC cohort, BOP1 was significantly upregulated in tumor samples versus normal tissue. Its knockdown in SK-LMS-1 cells achieved 83.09% silencing efficiency by qRT-PCR and significantly reduced cell migration, Transwell invasion, and colony formation (p < 0.001).
CTBP1 (C-Terminal Binding Protein 1): CTBP1 is a transcriptional co-repressor involved in EMT and tumor plasticity. It has been linked to sarcomatoid transformation in hepatocellular carcinoma, making its elevated expression in LMS biologically coherent. Knockdown efficiency was 78.72%, with significant reductions in all three functional assays. SRPK1 (Serine/Arginine Protein Kinase 1) facilitates tumor progression in gastric and breast cancers by modulating RNA splicing and has been documented at elevated levels in synovial sarcoma; knockdown efficiency was 79.59%.
DSE (Dermatan Sulfate Epimerase): DSE has been reported as a tumor suppressor in melanoma through regulation of extracellular matrix interactions and immune cell infiltration, but in the high mRNAsi LMS context it is upregulated as a risk gene (knockdown efficiency 78.36%). PSMD10 contributes to tumorigenesis through proteasomal regulation and has been implicated in liver and thyroid cancers; its knockdown achieved 72.43% efficiency. HACD4, by contrast, is the sole protective gene in the signature, highly expressed in normal tissue and associated with better survival outcomes.
Immune infiltration and regulatory network: TIMER2.0 analysis showed that the five risk genes were positively correlated with M0 macrophages and naive CD8+ T cells, while exhibiting negative correlation with other macrophage subtypes. This immune infiltration pattern suggests that these genes may promote a pro-tumorigenic immune state characterized by macrophage and T cell dysfunction. Functional enrichment (GO analysis) confirmed significant involvement of these genes in cell invasion and cell cycle pathways. GeneMANIA network prediction identified TGIF1, PLD1, and PMSC4 as key interaction partners, genes with established roles in TGFbeta signaling, phospholipid metabolism, and proteasomal function respectively.
Small scRNA-seq sample size: The single-cell component of this study relies on only four LMS samples drawn from one GEO dataset (GSE212527). While 52,781 cells were analyzed, those cells all come from a very limited number of patients. This constrains the generalizability of cell population frequencies, stemness score distributions, and cell-cell communication patterns to the broader LMS patient population. Rare cell subpopulations, treatment history differences, and site-specific biology (uterine vs. retroperitoneal LMS, for example) cannot be reliably captured with four samples.
Batch effects and residual technical noise: Single-cell RNA-seq data generated from different tissue processing protocols, sequencing depths, and platforms carry inherent batch effects. The Harmony algorithm was applied to correct inter-sample variation, but residual batch-specific noise may still influence cell clustering, stemness score calculation, and downstream differential expression results. The authors acknowledge that this technical limitation could affect the reliability of cell type annotations and the stability of the stemness-based stratification.
Lack of external validation for the prognostic model: The six-gene stemness risk score was constructed and tested entirely within the TCGA-SARC cohort using internal cross-validation. No independent external cohort was used to validate the AUC values of 82% at 1 year, 77% at 3 years, and 79% at 5 years. Internal validation consistently overestimates real-world predictive performance, and the absence of external validation is a significant limitation for a prognostic biomarker intended for clinical translation. The model also relied on FPKM normalization without platform standardization across different RNA-seq datasets, which could affect portability.
In vitro validation in a single cell line: Functional experiments were conducted exclusively in SK-LMS-1, a commercially available LMS cell line. While this provides proof-of-concept evidence that the five risk genes influence migration, invasion, and colony formation, single cell-line studies do not capture the genetic and phenotypic diversity of primary LMS tumors. No in vivo validation (xenograft or patient-derived xenograft models) was performed, and HACD4, the protective gene, was not included in the knockdown experiments. Causal mechanistic evidence linking these genes to the stemness-immune microenvironment axis in vivo remains to be established.
Expanding scRNA-seq cohorts and multi-omics integration: The most immediate priority is validating the stemness-based findings in larger, multi-institutional single-cell datasets that capture diverse LMS subtypes, anatomical sites, and treatment histories. Integrating genomic, epigenomic, and proteomic layers with single-cell transcriptomics would provide a more complete picture of the regulatory mechanisms driving high-stemness cell identity. Longitudinal sampling from the same patients before and during treatment would clarify whether the mRNAsi shifts in response to chemotherapy or targeted therapy, a dynamic question this cross-sectional dataset cannot address.
In vivo functional validation and therapeutic targeting: The causal roles of BOP1, CTBP1, DSE, PSMD10, and SRPK1 in LMS progression need to be tested in mouse xenograft models and, ideally, patient-derived xenograft (PDX) models. Overexpression studies, not just knockdown, would help confirm oncogenic function. These genes also represent candidate therapeutic targets: SRPK1 inhibitors, for instance, have been explored in preclinical oncology settings, and PSMD10's proteasomal regulatory function places it within a pathway already targeted clinically. Drug screening against the six-gene panel in LMS cell lines and organoids could identify exploitable pharmacological vulnerabilities.
Immune checkpoint strategies for high-stemness LMS: The finding that high mRNAsi tumor cells exhibit elevated TGFbeta and EGF signaling and a more complex immune microenvironment opens specific therapeutic hypotheses. Combination strategies that target TGFbeta or EGFR signaling alongside immune checkpoint inhibition may overcome the immune evasion that limits single-agent PD-1/PD-L1 blockade in LMS. The observation that high-stemness cells correlate with M0 macrophage infiltration suggests that macrophage repolarization strategies (e.g., CSF1R inhibition or CD47 blockade) could also be explored to shift the immunosuppressive tumor microenvironment.
Clinical validation and risk stratification: Longitudinal cohort studies with larger patient numbers are needed to determine whether the stemness risk score can be used for early detection, monitoring treatment response, or predicting relapse in clinical settings. The authors suggest that cross-platform validation using RNA-seq data from independent sarcoma datasets, combined with prospective collection of serial biospecimens, would provide the evidence base required to advance the six-gene model toward clinical utility. Integration with imaging biomarkers (MRI-based radiomics, PET metabolic tumor volume) could further refine risk stratification beyond what genomic signatures alone provide.