Molecular characterization and prognostic modeling of liquid-liquid phase separation-related genes in osteosarcoma based on single-cell sequencing and weighted gene co-expression network analysis

Introduction

Intracellular biological processes (BPs) are tightly regulated in space and time to ensure proper physiological function (1). Membrane-bound organelles have traditionally been considered the primary means of spatially compartmentalizing biochemical reactions (2-4). However, an alternative mode of subcellular organization, liquid-liquid phase separation (LLPS), has emerged as a pivotal mechanism for the formation of membrane-less organelles (5). These structures, driven by the physicochemical forces underlying LLPS, play essential roles in cellular function, gene regulation, and environmental adaptation (6,7).

LLPS depends on multivalent interactions among biomacromolecules, enabling proteins and nucleic acids to aggregate into dynamic biomolecular condensates (8-11). These condensates can form and dissolve rapidly, facilitating the fine-tuning of diverse cellular and biochemical processes in response to internal and external cues (12,13). Accumulating evidence highlights the critical involvement of LLPS in numerous physiological and pathological contexts, including neurodegenerative diseases, viral infections, and cancer (14,15).

In cancer biology, LLPS-derived biomolecular condensates modulate cell signaling and gene expression, thereby influencing tumor progression. A notable example is the tumor suppressor p53, the transcriptional activity of which is sustained by LLPS-mediated nuclear condensate formation. While wild-type p53 maintains its functionality in the liquid-phase state, oncogenic mutations can lead to protein misfolding and a transition from a liquid to a solid-like state, impairing DNA-binding capacity and tumor-suppressive function (16,17). Another hallmark of LLPS in tumors is the formation of stress granules (SGs) (18). Although the mechanisms by which SGs contribute to tumorigenesis remain to be fully elucidated, cancer cells exhibit an enhanced capacity to form SGs compared to normal cells. SGs support tumor survival under adverse conditions and contribute to phenotypes such as chemoresistance (19).

As the predominant malignant bone tumor, osteosarcoma represents a leading cause of cancer-linked deaths in pediatric and adolescent populations (20,21). While the combination of chemotherapy and surgical resection has improved 5-year survival rates, the emergence of chemoresistance limits long-term outcomes (22,23). In light of the lack of targeted therapies, elucidating novel pathogenic mechanisms is essential for the development of effective interventions. Recent studies suggest that LLPS may play previously unrecognized roles in cancer, including osteosarcoma.

Therefore, this study elucidated the role of LLPS in osteosarcoma utilizing transcriptomic and single-cell sequencing data. Through a combination of bioinformatics algorithms, our study identified several LLPS-related prognostic markers and developed a robust prognostic model. These findings present fresh insights into the molecular pathogenesis of osteosarcoma and highlight potential therapeutic targets. The results were further validated by reverse transcription quantitative polymerase chain reaction (RT-qPCR), immunohistochemistry (IHC) and proliferation assays. We present this article in accordance with the TRIPOD+AI and MDAR reporting checklists (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1366/rc).

Methods

Identification of LLPS-related genes

The LLPS database (http://llps.biocuckoo.cn/), a dedicated online resource for LLPS, was utilized to identify LLPS-related genes. In this study, LLPS-related genes involving three distinct mechanisms (client, regulator, and scaffold) that had been reviewed were included, and only those associated with Homo Sapiens were extracted. A total of 3,611 candidate LLPS-related genes were obtained after compilation (table online: https://cdn.amegroups.cn/static/public/tcr-2025-1366-1.xlsx) and defined as the LLPS gene set.

Data acquisition and processing

Single-cell RNA sequencing data of osteosarcoma (GSE162454), comprising six osteosarcoma samples, were from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/). Downstream single-cell analysis was performed using the Seurat package. Cells with fewer than 200 or more than 8,000 unique features, as well as those with mitochondrial gene content exceeding 25%, were filtered out. Quality-controlled data were subjected to principal component analysis (PCA) for dimensionality reduction, followed by clustering via uniform manifold approximation and projection (UMAP). Canonical marker genes were utilized for cell type classification. To evaluate the expression patterns of LLPS-related genes, the PercentageFeatureSet function was used to calculate the relative abundance of these specific genes in each single cell, expressed as the percentage of their total expression among all detected genes in that cell. Cells were subsequently divided into high- and low-expression cohorts based on median percentage values. Inter-group differentially expressed genes (DEGs) were detected via the FindMarkers function (|log₂FC| =0.5; P<0.05).

Two additional osteosarcoma gene expression datasets, GSE21257 (GPL10295; n=53) and GSE39058 (GPL14951; n=47), both of which contain gene expression matrices and survival information, were downloaded via the getGEO function in the GEOquery package. These datasets were used for prognostic model construction and validation.

UCSC Xena (https://xenabrowser.net/datapages/), which aggregates large-scale genomics datasets including The Cancer Genome Atlas (TCGA) and Therapeutically Applicable Research to Generate Effective Treatments (TARGET), was used to download the GDC TARGET-OS dataset. After excluding samples without survival data, 30 cases remained for further model validation.

Single-sample enrichment scoring of LLPS-related genes

The single-sample gene set enrichment analysis (ssGSEA) method from the GSVA package was applied to compute the LLPS score for each sample. Genes in each sample were ranked in descending order according to their expression levels, followed by a gene set enrichment analysis (GSEA)-like running sum and statistical calculation for the LLPS gene set. The results were normalized according to the GSVA recommendations to ensure comparability among samples. The input data were derived from the microarray expression matrix of GSE21257, and no scoring was performed when the number of detectable genes within a gene set was less than 10. A higher LLPS score indicated stronger coordinated upregulation of LLPS-related genes in the corresponding sample.

Weighted gene co-expression network analysis (WGCNA)

As a systems biology approach, WGCNA helps to characterize gene correlation patterns, detect gene modules, and examine their links to phenotypic characteristics. It aids in elucidating biological mechanisms, identifying biomarkers, and exploring molecular pathways. Gene co-expression modules were constructed through the WGCNA package, and Pearson correlation coefficients were calculated between module eigengenes and LLPS scores to identify modules most strongly related to LLPS activity.

Functional enrichment analysis

To unveil the biological functions and pathways potentially involving LLPS-related genes, an intersection was taken between DEGs from the single-cell analysis and genes within the LLPS-associated WGCNA modules. The intersection was visualized via a Venn diagram. Gene Ontology (GO) analysis on these genes was carried out utilizing the clusterProfiler package. The aim was to identify related BP, cellular components (CCs), as well as molecular functions (MFs). The Kyoto Encyclopedia of Genes and Genomes (KEGG) helped to analyze possible signaling pathways. All enriched GO and KEGG terms were considered significant, while the P and Q-value were <0.05.

Construction and validation of the prognostic model

Univariate Cox proportional hazards regression analysis was conducted to identify prognosis-associated genes from the intersected gene set (P<0.05). To prevent model overfitting and identify the most robust prognostic genes, least absolute shrinkage and selection operator (LASSO) regression with 10-fold cross-validation was employed. The optimal penalty parameter (“lambda.min”) was selected for model construction. The candidate genes selected via LASSO were further analyzed through multivariate Cox regression to identify the final hub genes. Each sample’s risk score was computed by applying the predict function to the hub genes’ expression levels. Based on the median risk score, samples were split into high- and low-risk categories. Survival differences across the groups were assessed via Kaplan-Meier (KM) analysis. Moreover, heatmaps presented the expression patterns of the hub genes, while receiver operating characteristic (ROC) curves were derived to test the risk model’s predictive accuracy.

Nomogram construction

A prognostic nomogram was generated using hub genes identified after multivariate Cox regression analysis to help forecast individual survival outcomes in the osteosarcoma population. Calibration curves were utilized to check how closely the predicted results matched the actual outcomes. Decision curve analysis (DCA) was carried out to evaluate the nomogram’s usefulness in clinical practice.

Evaluation of the independent prognostic value

Based on the clinical information from the GSE21257 dataset (including age, sex, and Huvos grade of osteosarcoma) and the LLPS risk score, univariate and multivariate Cox regression analyses were performed to determine whether the risk score served as an independent prognostic factor. Furthermore, prognostic differences between high- and low-risk groups within each clinical subgroup were compared to demonstrate the stability and broad applicability of its prognostic discriminative capacity.

Functional characterization and spatial distribution of hub genes

GSEA on every hub gene was performed to further characterize their functional implications, with enrichment plots generated for visualization. In addition, gene-specific UMAP plots based on the GSE162454 dataset were constructed to illustrate the spatial distribution of hub genes across single cells.

Immune-related analyses

To better understand the immune differences between the high- and low-risk cohorts, immune-related analyses were conducted on hub genes. The compositional differences of 22 immune cell types across groups were examined through CIBERSORT. Correlation heatmaps were used to visualize the associations both among immune cell populations and between hub genes and immune cells. Additionally, leveraging the ssGSEA method, enrichment scores were computed for 16 categories of tumor-infiltrating immune cells along with 13 immune-related functions (24). Differential expression analyses of immune checkpoint molecules and human leukocyte antigen (HLA) genes across risk groups were also conducted, with results presented in boxplots.

Drug sensitivity analysis

To identify potentially beneficial therapeutic agents for osteosarcoma patients, drug sensitivity analysis was conducted using the “cgp2016ExprRma” dataset from the Genomics of Drug Sensitivity in Cancer Database (GDSC). The sensitivity of each drug across the risk groups was estimated via the pRRopheticPredict function. A more stringent filtering criterion was adopted, retaining only results with a P<0.01 to ensure robust findings.

The experiments of RT-qPCR

Among the four predicted oncogenes, EXOSC8 exhibited the highest coefficient in the prognostic model, suggesting its potentially significant impact on osteosarcoma. Thus, EXOSC8 was selected for experimental validation.

The RT-qPCR was performed to evaluate the mRNA expression levels of EXOSC8 in osteoblasts and five commonly used osteosarcoma cell lines (143B, HOS, U2OS, SAOS2, MG63). Total RNA was extracted utilizing the RNA Extraction Kit (Beyotime, China), followed by reverse transcription using a cDNA Synthesis Kit (Takara, Japan). Amplification was conducted through FastStart Universal SYBR Green Master ROX (Roche, Mannheim , Germany). GAPDH served as the internal control, and the 2^−ΔΔCt approach was employed for relative expression evaluation. The EXOSC8 primers were designed and synthesized by Sangon Biotech Co., Ltd. (Shanghai, China) with the sequences: left primer 1-ACCATCAACAGATGCCCCTG; right primer 1-CACTTGGGCCTCTTCTCCAG.

Patient-derived samples and IHC

Tissue samples from three patients, including both osteosarcoma and adjacent non-tumorous tissues, were collected for immunohistochemical analysis. They received surgical resection at The First Affiliated Hospital of Guangxi Medical University. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. The study was approved by the Medical Ethics Committee of The First Affiliated Hospital of Guangxi Medical University (approval No. 2025-E0235). All patients or their legal guardian(s) had given informed consents to the use of their samples for research before study. The tissue sections were prepared by the Guangxi Medical University Ruigu Medical Science Inspection Center.

For IHC, an anti-EXOSC8 antibody was obtained from Proteintech Group, Inc. (Wuhan, China). Tissue sections underwent sequential incubation with primary and secondary antibodies, 3,3’-diaminobenzidine (DAB) chromogenic reaction, hematoxylin counterstaining, dehydration, and mounting. The protein expression of EXOSC8 was subsequently evaluated under an electron microscope.

Cell transfection

According to the experimental design, cells were seeded into plate wells. When cell confluence reached 60–70%, following the manufacturer’s instructions, the silencing plasmid targeting EXOSC8 (Table S1) and its corresponding control plasmid were mixed with TransMate transfection reagent, respectively, to prepare transfection complexes, which were then added to the cell culture plate. After 6 hours of transfection, the medium was replaced with fresh complete medium. The cells were subsequently cultured for 48 hours under conditions of 37 ℃ and 5% CO₂ for further experimental analysis.

Cell viability assay

HOS cells were seeded at a density of 5×10³ cells per well in a 96-well plate. When cell confluence reached 60–70%, plasmid transfection was performed according to the experimental groups (CON, NC, si-EXOSC8). After 48 hours of transfection, 10 µL of Cell Counting Kit-8 (CCK-8) reagent was added to each well, followed by incubation for 1 hour. Absorbance was measured at a wavelength of 450 nm using a microplate reader.

EdU cell proliferation assay

The EdU-488 Cell Proliferation Detection Kit was used. HOS cells were seeded at a density of 20×10³ cells per well in a 24-well plate. After 48 hours of plasmid transfection according to the experimental groups, EdU reagent was added and incubated for 2 hours. The cells were then fixed, permeabilized, and the reaction solution was added. Finally, after counterstaining with 4’,6-diamidino-2-phenylindole (DAPI), cell images were captured using an upright fluorescence microscope (Zeiss, Oberkochen, Germany).

Clonogenic formation

HOS cells were seeded at a low density (500–1,000 cells per well) in a 6-well plate, and plasmid transfection was performed according to the experimental groups. After transfection, the cells were continuously cultured for 10–14 days until visible clonal colonies formed. The colonies were fixed with 4% paraformaldehyde and stained with 0.1% crystal violet. The colony numbers were quantified using ImageJ.

Statistical analysis

Our bioinformatics analyses were enabled by R 4.3.2. Quantification of IHC results was conducted using Image J win 64. Graphical visualizations of experimental results were generated with GraphPad Prism 8.0. The study employed Cox regression analyses (both univariate and multivariate) to pinpoint prognosis-related genes and develop the prognostic model. Survival outcomes were assessed via KM analysis, with drug sensitivity comparisons across groups performed through Wilcoxon rank-sum testing. A two-tailed P value <0.05 denoted statistical significance (P<0.05 denoted by *, P<0.01 by **, and P<0.001 by ***).

Results

Identification of LLPS-related genes

A total of 3,611 genes related to LLPS were obtained from the LLPS online database. These were collectively defined as the LLPS-related gene set.

Single-cell RNA sequencing analysis

Before quality control, 50,780 cells were obtained; following filtering, 45,787 high-quality cells were retained, with mitochondrial gene expression accounting for less than 25% (Figure 1A,1B). A significant positive correlation was observed between the total gene expression and the number of detected genes, with correlation coefficients of 0.88 before quality control and 0.90 after, while both metrics showed no significant correlation with mitochondrial gene expression (Figure 1C,1D). A volcano plot was used to display the top 2,000 highly variable genes, including the top 10 most variable genes (Figure 1E). PCA revealed a uniform distribution of cells across the six samples, indicating minimal batch effects (Figure 1F). At a resolution of 0.1, cells were grouped into 13 clusters (Figure 1G). Differential expression analysis (DEA) was enabled by the FindAllMarkers function with a log₂ fold change threshold of 1, and the top 10 DEGs in each cluster were visualized using a heatmap (Figure 1H). Classical marker genes were employed to annotate cell types, which were grouped into six major categories: osteoblast, fibroblast/stromal, myeloid, T, B, and endothelial cell (Figure 1I). The relative proportions of these cell types across different samples were visualized in a bar plot (Figure 1J). Subsequently, the proportion of LLPS-linked genes expressed in every cell was computed. Based on the median proportion (42.81396), cells were split into high- and low-LLPS expression cohorts (Figure 1K). 2,670 DEGs were identified.

Figure 1 Quality control and analysis of single-cell RNA sequencing data. (A,B) Quality control before and after cell filtering. (C,D) Correlation between mitochondrial gene expression and total gene count before and after filtering. (E) Volcano plot of the top 2,000 highly variable genes. (F) PCA plot showing cell distribution across six samples. (G) Clustering analysis identifying 13 distinct clusters. (H) Heatmap of the top 10 DEGs in each cluster. (I) Cell type annotation resulting in six major cell types. (J) Bar plot showing the relative proportions of the six cell types. (K) Division of cells into high- and low-LLPS expression groups based on LLPS gene expression. DEG, differentially expressed gene; LLPS, liquid-liquid phase separation; OS, osteosarcoma; PCA, principal component analysis; UMAP, uniform manifold approximation and projection.

Transcriptomic data analysis

The LLPS enrichment score was quantified for every osteosarcoma sample in the GSE21257 dataset via ssGSEA. Subsequently, WGCNA was performed. Genes expressed in more than half of the samples were retained after calculating the average expression values across all genes. Hierarchical clustering using the hclust function was applied to identify and exclude outlier samples (Figure 2A), resulting in 52 samples and 12,499 genes included in the final analysis. When the soft-thresholding power was set to 2, the scale-free topology index R²>0.8, and the average connectivity of the network stabilized (Figure 2B). After merging similar modules, 12 non-grey gene modules were identified (Figure 2C-2E), among which the turquoise module exhibited the strongest correlation with LLPS (correlation coefficient =0.84, P<0.001), 3,780 genes from the turquoise module were retained for subsequent analyses (Figure 2F).

Figure 2 Transcriptomic data and GSEA. (A) Outlier sample GSM531298 was excluded; 52 samples were retained for WGCNA. (B) A soft-thresholding power of 2 was suitable for network construction. (C-E) 12 non-grey gene modules were identified. (F) The turquoise module exhibited the strongest association with LLPS. (G) Venn diagram of 599 overlapping LLPS-related DEGs. (H) Results of GO enrichment analysis. (I) Results of KEGG pathway enrichment analysis. BP, biological process; CC, cellular component; DEG, differentially expressed gene; GO, Gene Ontology; GSEA, gene set enrichment analysis; KEGG, Kyoto Encyclopedia of Genes and Genomes; LLPS, liquid-liquid phase separation; MF, molecular function; TCA, tricarboxylic acid cycle; WGCNA, weighted gene co-expression network analysis.

GSEA

The intersection of the 3,780 WGCNA-derived genes and the 2,670 DEGs yielded 599 LLPS-related DEGs (Figure 2G). These 599 genes were subjected to GO and KEGG enrichment analyses, with the top 10 significantly enriched BPs or pathways presented (Figure 2H,2I).

GO-BP analysis indicated that these genes are primarily involved in protein translation, modification, folding, and stability regulation. They participate in cytoplasmic and mitochondrial translation, ribosome biogenesis, and the regulation of ubiquitin-mediated protein degradation, post-translational modification, and protein localization, underscoring their pivotal role in maintaining protein homeostasis and cellular function. GO-CC analysis revealed the predominant enrichment in ribosomes, mitochondrial protein complexes, and focal adhesions, suggesting important roles in protein synthesis, energy metabolism, and cell migration—processes critical to cellular adaptation and cancer progression. GO-MF analysis further demonstrated their involvement in ribosomal structure, transcriptional regulation, protein stability, and energy metabolism, supporting their influence on protein synthesis, gene regulation, and metabolic processes.

KEGG pathway analysis indicated significant enrichment in protein synthesis (ribosome), energy metabolism (tricarboxylic acid cycle), protein degradation (proteasome), and cancer-related signaling pathways. Additionally, these genes were implicated in oxidative stress responses and N-glycan biosynthesis. Notably, associations with neurodegenerative diseases and viral infections were also observed, suggesting that these genes may regulate protein homeostasis, metabolic reprogramming, and cellular stress responses. Collectively, the KEGG findings highlight their potential involvement in the pathogenesis, progression, and therapeutic resistance of osteosarcoma.

Construction and validation of the LLPS-associated prognostic model

Following the normalization of all samples in the GSE21257 dataset (Figure 3A,3B), univariate Cox regression analysis was carried out using the intersecting genes with a significance threshold of P<0.05. This yielded a preliminary list of 120 prognosis-related genes (table online: https://cdn.amegroups.cn/static/public/tcr-2025-1366-2.xlsx). Setting the random seed to 2024, LASSO regression was executed with 10-fold cross-validation. As shown in Figure 3C,3D, the optimal lambda value was determined to be 0.1038315, identifying nine candidate genes: CPE, CCT6A, EXOSC8, SMARCB1, CTSH, ACO2, AASDH, SEPHS1, and C15orf40. Subsequently, by applying backward stepwise multivariate Cox regression, CPE and CTSH were eliminated, yielding a final model comprising seven hub genes. Specifically, CCT6A, EXOSC8, ACO2, and SEPHS1 were identified as potential oncogenes, while SMARCB1, AASDH, and C15orf40 appeared to exhibit tumor-suppressive properties.

Figure 3 Construction of the LLPS-associated prognostic model. (A,B) Batch effects were successfully mitigated following normalization. (C,D) The prognostic model was constructed using LASSO regression. (E) High-risk group displayed notably poorer survival outcomes. (F) Heatmap illustrated differential expression of hub genes across risk groups. (G) ROC curves demonstrated excellent model performance. (H-J) Validation using TARGET-OS dataset: similar expression patterns and survival trends. (K-M) Validation using the GSE39058 dataset: consistent results with the training cohort. AUC, area under the curve; LASSO, least absolute shrinkage and selection operator; LLPS, liquid-liquid phase separation; ROC, receiver operating characteristic.

Patients were stratified into high- and low-risk groups based on the median risk score. A significant difference in survival was observed between the two groups, with the high-risk group exhibiting poorer prognoses (Figure 3E). The expression profiles of the hub genes between risk groups were visualized using a heatmap, highlighting their potential biological relevance in osteosarcoma (Figure 3F). ROC curves were subsequently used to evaluate model performance, with the areas under the curve (AUCs) for 1-, 3-, and 5-year survival reaching 0.893, 0.950, and 0.947, which suggest excellent predictive performance (Figure 3G).

To validate the prognostic model, two independent external datasets, TARGET-OS and GSE39058, were employed. Validation results were presented using KM survival curves, heatmaps, and ROC curves (Figure 3H-3J for TARGET-OS and Figure 3K-3M for GSE39058), both of which demonstrated similar trends and consistent model performance.

Construction of a nomogram

To visually represent the contribution of hub genes to survival prognosis, a nomogram was constructed (Figure 4A). The prediction accuracy of the nomogram was evaluated through calibration curves and DCA (Figure 4B,4C). As illustrated, the calibration curve demonstrated good agreement with actual outcomes, while DCA suggested that clinical decision-making based on the LLPS-associated model would provide the greatest benefit to patients.

Figure 4 Nomogram construction. (A) Nomogram based on hub gene expression profiles. (B) Calibration curve indicating strong concordance with observed survival rates. (C) DCA illustrating the clinical utility of the LLPS-based model. DCA, decision curve analysis; LLPS, liquid-liquid phase separation.

Evaluation of independent prognostic value

As shown in the figures, following univariate and multivariate Cox regression analyses, both the osteosarcoma Huvos grade and the LLPS risk score were retained (forest plots, Figure 5A,5B). Therefore, these two indicators are considered to have independent prognostic value, with higher Huvos grades serving as protective factors and higher LLPS risk scores as risk factors. In subgroup analyses based on clinical features, the risk score remained independently associated with prognosis within age and Huvos stratifications (P<0.05), indicating stability of the model in major clinical contexts. No significant differences were observed in the female subgroup, which may be attributable to the limited sample size (Figure 5C-5H).

Figure 5 Evaluation of independent prognostic value. (A,B) Indicators with independent prognostic significance were identified based on univariate and multivariate Cox regression analyses. (C-H) Survival assessment within clinical subgroups, including age, sex, and Huvos grade. CI, confidence interval; HR, hazard ratio.

Subcellular localization and GSEA of hub genes

To further understand the role and spatial distribution of hub genes, the cellular localization and GSEA analysis of hub genes were performed. The cellular localization of hub Genes was based on single-cell sequencing data. As shown in Figure 6A-6G, CCT6A, EXOSC8, ACO2, and SEPHS1 were all highly expressed in osteoblasts to varying degrees. CCT6A, EXOSC8 and ACO2 were also highly expressed in T cells, myeloid cells, and fibroblast/stromal. AASDH was specifically highly expressed in fibroblast/stromal. The expression of C15orf40 was low in all cells. Interestingly, although the prediction model constructed in this study and the generally accepted view both believe that SMARCB1 is an antioncogene, sublocalization analysis suggested that it was highly expressed in osteoblasts. Given that there are very few studies on the mechanism of SMARCB1 in osteoblasts of osteosarcoma, this finding needs to be further validated, and such phenomenon may just occur by chance. In addition, GSEA enrichment analysis was performed for each hub gene to evaluate its enriched pathways (Figures S1-S7).

Figure 6 Subcellular localization of hub genes. (A-G) Distribution of hub genes across different cell types. UMAP, uniform manifold approximation and projection.

Immune-related analyses

To investigate the immune landscape related to risk stratification, the CIBERSORT method was initially utilized to assess the infiltration levels of 22 distinct immune cell populations in samples from both high- and low-risk cohorts (Figure 7A). A heatmap presented possible correlations among the 22 immune cell sorts, and between immune cells and the hub genes (Figure 7B,7C). As illustrated, the composition of immune cells varied considerably between groups. Notably, T cells CD4 memory resting exhibited a significant positive correlation with the hub genes, whereas T cells gamma delta showed a significant negative correlation. Based on 22 immune cell types, their correlation with LLPS scores was further analyzed. The results indicated that Monocytes, CD8⁺ T cells, and M2 Macrophages were negatively correlated with LLPS, whereas resting memory CD4⁺ T cells and M0 Macrophages were positively correlated; other cell types showed no significant correlation (Figure 7D). Furthermore, LLPS scores for each sample in GSE21257 were calculated using the ssGSEA method. Between the high- and low-LLPS score groups, CD8⁺ T cells and Monocytes exhibited more pronounced immune infiltration in the low-risk group, while resting memory CD4⁺ T cells and M0 Macrophages showed the opposite pattern (Figure 7E). These findings suggest complex and potentially biologically meaningful interactions among immune cell types and between immune cells and hub genes that warrant further exploration.

Figure 7 Immune analyses related to hub genes. (A) Immune cell composition across individual samples. (B,C) Correlation among 22 immune cell types and between immune cells and hub genes. (D) Correlation between 22 immune cell types and LLPS scores. (E) Differences in immune cell infiltration between high- and low-LLPS score groups. (F) Enrichment scores of 16 tumor-infiltrating immune cells. (G) Enrichment scores of 13 functional pathways related to immunity. (H) Differentially expressed immune checkpoints. (I) Differentially expressed HLAs. *, P<0.05; **, P<0.01; ***, P<0.001. HLA, human leukocyte antigen; LLPS, liquid-liquid phase separation.

Subsequently, ssGSEA was performed to evaluate the enrichment of immune cell infiltration within tumors, covering 16 immune cell subtypes and 13 functional pathways linked to immunity (Figure 7F,7G). Differences in immune checkpoint expression and HLA expression between the high- and low-risk cohorts were also examined (Figure 7H,7I). Among the 16 tumor-infiltrating immune cell types, significant differences were observed for activated dendritic cells (aDCs), mast cells, neutrophils, Type II helper T cells (Th2), and tumor-infiltrating lymphocytes (TILs), with enrichment scores being markedly higher in the low-risk cohort. In addition, marked intergroup differences were noted in several immune functions like antigen-presenting cell (APC) co-inhibition, APC co-stimulation, chemokine receptor activity, HLA expression, inflammation promotion, as well as Type II interferon response. These results suggest that diminished infiltration of specific immune cells and impaired immune functions may contribute to greater malignancy in osteosarcoma.

Additionally, seven immune checkpoints and 14 HLAs were differentially expressed between risk groups. The immune checkpoints with altered expression were primarily involved in immune regulation mediated by T and natural killer (NK) cells. The differentially expressed HLAs could be broadly categorized into major histocompatibility complex class I and class II (MHC I and MHC II) molecules for immune recognition and antigen presentation. Most of them were upregulated in the low-risk cohort, suggesting their potential protective role in osteosarcoma pathogenesis.

LLPS-associated drug sensitivity analysis

Drug sensitivity analysis identified five candidate compounds demonstrating differential responses between the high- and low-risk groups (Figure 8A-8E). The low-risk group was potentially more sensitive to SB 505124, pyrimethamine, and KIN001-135, whereas patients in the high-risk group may benefit more from CGP-60474 and dimethyloxalylglycine (DMOG). These compounds exert diverse pharmacological effects, including inhibition of the transforming growth factor-β (TGF-β) receptor, modulation of autophagy, regulation of the MAPK/PI3K signaling pathways, inhibition of CDK1/2/5, and stabilization of hypoxia-inducible factor-1α (HIF-1α). These mechanisms may influence key pathophysiological processes in tumor progression.

Figure 8 Drug sensitivity analysis. (A-E) Differential drug responses for five candidate compounds. IC₅₀, half maximal inhibitory concentration.

Experiments to verify EXOSC8

To validate the expression of EXOSC8, RT-qPCR analysis demonstrated significantly elevated mRNA levels in four osteosarcoma cell lines (excluding SAOS2) (Figure 9A). IHC further demonstrated EXOSC8 upregulation at the protein level in osteosarcoma tissues (Figure 9B,9C). HOS was selected for subsequent experiments due to the highest expression level. CCK-8, EdU, and colony formation assays demonstrated that silencing EXOSC8 significantly inhibited osteosarcoma cell proliferation (Figure 10A-10F).

Figure 9 Validation of EXOSC8 expression. (A) Elevated EXOSC8 expression at the mRNA level in osteosarcoma cell lines. (B) Protein-level validation of EXOSC8 overexpression in osteosarcoma tissues (NT: adjacent non-tumorous tissues; PT: tumor tissue). (C) Statistical analysis of immunohistochemical results. *, P<0.05; ***, P<0.001; ****, P<0.0001; ns, not significant. IHC, immunohistochemistry.

Figure 10 Effect of EXOSC8 on proliferation. (A) Silencing of EXOSC8 in HOS. (B-F) Cell proliferation was significantly inhibited after EXOSC8 silencing (E, the colonies were fixed with 4% paraformaldehyde and stained with 0.1% crystal violet). ***, P<0.001; ns, not significant. Con, control; DAPI, 4',6-diamidino-2-phenylindole; EdU, 5-Ethynyl-2’-deoxyuridine; HOS, one of the cell lines of osteosarcoma; NC, negative control; OD, optical density.

Discussion

Osteosarcoma is the most common malignant bone tumor featuring high therapeutic resistance and poor clinical outcomes. One of the major obstacles to improving treatment efficacy lies in its complex molecular mechanisms, which remain incompletely understood. Therefore, it is crucial to keep exploring the molecular mechanisms behind osteosarcoma and to discover new targets for treatment. LLPS has recently garnered growing attention in oncology. Originally a physiological phenomenon within cells, aberrant LLPS has been implicated in the formation of oncogenic condensates, potentially inducing tumor-associated phenotypes in otherwise normal cells (25).

The present study employed multiple computational approaches, including WGCNA and LASSO regression, to identify LLPS-related genes in osteosarcoma. Multivariate Cox regression analysis subsequently identified seven hub genes as potential biomarkers related to osteosarcoma prognosis, which were then incorporated into a prognostic model and nomogram. ROC and DCA demonstrated that the model possesses substantial clinical utility, offering a convenient tool for evaluating patient survival and informing clinical decision-making. Based on the median risk score, patients were stratified into high- and low-risk cohorts. Further immune profiling was executed to compare tumor microenvironment (TME) characteristics and immune cell infiltration between risk groups, aiming to elucidate underlying immunological mechanisms and variations in tumor immune regulation. Moreover, drug sensitivity analysis was performed to predict differential responses to conventional chemotherapeutic agents, targeted therapies, and immune checkpoint inhibitors across risk groups, thereby assessing the feasibility of individualized treatment strategies and contributing to precision oncology.

Ultimately, seven hub genes were identified as regulators of osteosarcoma pathogenesis and progression. Chaperonin containing TCP1 subunit 6A (CCT6A), encoding a molecular chaperone, is frequently overexpressed in various malignancies and acts as an oncogene by modulating glycolytic activity (26-29). Exosome component 8 (EXOSC8), encoding a subunit of the RNA exosome complex, is critical in RNA degradation and processing. Cui et al. reported its oncogenic role in colorectal cancer (30). Aconitase 2 (ACO2), a key enzyme in the tricarboxylic acid cycle, is highly expressed in most cancers according to a pan-cancer analysis of TCGA data (31). Interestingly, in non-small cell lung cancer (NSCLC), downregulation of ACO2 expression has been shown to immortalize bronchial epithelial cells (32). Selenophosphate synthetase 1 (SEPHS1) catalyzes the synthesis of selenophosphate from selenide and ATP and is implicated in various diseases (33,34). Yang et al. demonstrated that SEPHS1 overexpression enhances hepatocellular carcinoma cell invasiveness and is significantly related to poor prognosis (35). SEPHS1 mRNA has abundant rapidly proliferating cells like embryonic and cancer cells, and gene disruption leads to impaired cell proliferation (33,35,36). Our findings suggest that CCT6A, EXOSC8, ACO2, and SEPHS1 may function as oncogenes in osteosarcoma by promoting tumor initiation and progression.

SMARCB1, also known as INI1, encodes a protein that is one of the core subunits of the complex of SWItch/Sucrose Non-Fermentable chromatin remodeling complex (SWI/SNF). SMARCB1 is widely considered as a classic tumor suppressor gene, which is consistent with the results of this study. However, its high expression in osteoblasts observed in this study needs further verification, because this phenomenon may occur by chance. SMARCB1 loss destabilizes SWI/SNF, making it more susceptible to degradation (37). Low SWI/SNF activity can alter genome-wide chromatin accessibility, directly affecting gene expression (38). SWI/SNF activity is continuously required to maintain normal chromatin structure, and its inactivation may trigger a series of oncogenic molecular events (39). In contrast, the biological functions of AASDH and C15orf40 remain largely unexplored. However, based on our results, they were preliminarily classified as tumor suppressor genes, which merit further investigation in future studies.

Beyond studies on tumor cells themselves, current cancer research increasingly emphasizes the TME (40-42). In this study, based on the marked differences in immune cell infiltration between high- and low-risk osteosarcoma groups, immune infiltration could be considered in future clinical management. This approach would allow rapid prognostic assessment at initial diagnosis and provide a reference for comparison during follow-up to evaluate disease progression. Grouping based on high/low LLPS scores can also utilize corresponding immune cell infiltration levels to distinguish high- and low-LLPS osteosarcoma subtypes, potentially guiding treatment strategies. Analyses of immune molecules mostly suggest a protective role. Therefore, future interventions targeting LLPS to modulate immune molecule expression or abundance and remodel the osteosarcoma TME may offer a novel strategy for improving prognosis, aligning with current trends in cancer research.

However, this study has also some limitations. Due to the lack of necessary clinical information in GSE21257, the predictive value of the model for treatment responses could not be assessed. Additionally, the direct interactions between hub genes or immune-related molecules and LLPS were not explored. These gaps, however, provide directions for future research, which may further investigate these relationships.

In conclusion, this study identified seven LLPS-associated genes that are significantly linked to osteosarcoma prognosis. These genes may contribute to osteosarcoma development and progression through LLPS-related molecular mechanisms and represent promising targets for future research and clinical intervention, offering novel insights into the diagnosis and treatment of osteosarcoma.

Conclusions

This study identified seven DEGs related to LLPS that are of potential significance in osteosarcoma. These genes can serve as new diagnostic biomarkers and prognostic indicators from the perspective of LLPS. Our study presents new perspectives on osteosarcoma diagnosis and treatment.

Acknowledgments

We are deeply grateful to Haijun Tang for his extensive expertise and technical guidance, which greatly contributed to the success of this study.

Footnote

Reporting Checklist: The authors have completed the TRIPOD+AI and MDAR reporting checklists. Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1366/rc

Data Sharing Statement: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1366/dss

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1366/prf

Funding: This research was supported by the National Natural Science Foundation of China under Grant (grant Nos. 81760485 and 82160536).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1366/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. The study was approved by the Medical Ethics Committee of The First Affiliated Hospital of Guangxi Medical University (approval No. 2025-E0235). All patients or their legal guardian(s) had given informed consents to the use of their samples for research before study.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Zhang J, Qiu Z, Zhang Y, et al. Intracellular spatiotemporal metabolism in connection to target engagement. Adv Drug Deliv Rev 2023;200:115024. [Crossref] [PubMed]
2. Chen B, Lyssiotis CA, Shah YM. Mitochondria-organelle crosstalk in establishing compartmentalized metabolic homeostasis. Mol Cell 2025;85:1487-508. [Crossref] [PubMed]
3. Sigrist SJ, Haucke V. Orchestrating vesicular and nonvesicular membrane dynamics by intrinsically disordered proteins. EMBO Rep 2023;24:e57758. [Crossref] [PubMed]
4. Jain A, Zoncu R. Organelle transporters and inter-organelle communication as drivers of metabolic regulation and cellular homeostasis. Mol Metab 2022;60:101481. [Crossref] [PubMed]
5. Hyman AA, Weber CA, Jülicher F. Liquid-liquid phase separation in biology. Annu Rev Cell Dev Biol 2014;30:39-58. [Crossref] [PubMed]
6. Hirose T, Ninomiya K, Nakagawa S, et al. A guide to membraneless organelles and their various roles in gene regulation. Nat Rev Mol Cell Biol 2023;24:288-304. [Crossref] [PubMed]
7. Li Y, Liu Y, Yu XY, et al. Membraneless organelles in health and disease: exploring the molecular basis, physiological roles and pathological implications. Signal Transduct Target Ther 2024;9:305. [Crossref] [PubMed]
8. Giudice J, Jiang H. Splicing regulation through biomolecular condensates and membraneless organelles. Nat Rev Mol Cell Biol 2024;25:683-700. [Crossref] [PubMed]
9. Mohanty P, Kapoor U, Sundaravadivelu Devarajan D, et al. Principles Governing the Phase Separation of Multidomain Proteins. Biochemistry 2022;61:2443-55. [Crossref] [PubMed]
10. Liu Z, Qin Z, Liu Y, et al. Liquid-liquid phase separation: roles and implications in future cancer treatment. Int J Biol Sci 2023;19:4139-56. [Crossref] [PubMed]
11. Elena HL, Claudia CB, Enrique C, et al. Unveiling the role of biomolecular condensates in cellular function and cancer. Adv Biol Regul 2025; Epub ahead of print. [Crossref]
12. Fefilova AS, Fonin AV, Vishnyakov IE, et al. Stress-Induced Membraneless Organelles in Eukaryotes and Prokaryotes: Bird's-Eye View. Int J Mol Sci 2022;23:5010. [Crossref] [PubMed]
13. Demmerle J, Hao S, Cai D. Transcriptional condensates and phase separation: condensing information across scales and mechanisms. Nucleus 2023;14:2213551. [Crossref] [PubMed]
14. Boyko S, Surewicz WK. Tau liquid-liquid phase separation in neurodegenerative diseases. Trends Cell Biol 2022;32:611-23. [Crossref] [PubMed]
15. Tong X, Tang R, Xu J, et al. Liquid-liquid phase separation in tumor biology. Signal Transduct Target Ther 2022;7:221. [Crossref] [PubMed]
16. Datta D, Navalkar A, Sakunthala A, et al. Nucleo-cytoplasmic environment modulates spatiotemporal p53 phase separation. Sci Adv 2024;10:eads0427. [Crossref] [PubMed]
17. Chen C, Fu G, Guo Q, et al. Phase separation of p53 induced by its unstructured basic region and prevented by oncogenic mutations in tetramerization domain. Int J Biol Macromol 2022;222:207-16. [Crossref] [PubMed]
18. Li T, Zeng Z, Fan C, et al. Role of stress granules in tumorigenesis and cancer therapy. Biochim Biophys Acta Rev Cancer 2023;1878:189006. [Crossref] [PubMed]
19. Redding A, Grabocka E. Stress granules and hormetic adaptation of cancer. Trends Cancer 2023;9:995-1005. [Crossref] [PubMed]
20. Panez-Toro I, Muñoz-García J, Vargas-Franco JW, et al. Advances in Osteosarcoma. Curr Osteoporos Rep 2023;21:330-43. [Crossref] [PubMed]
21. Beird HC, Bielack SS, Flanagan AM, et al. Osteosarcoma. Nat Rev Dis Primers 2022;8:77. [Crossref] [PubMed]
22. Zhra M, Akhund SA, Mohammad KS. Advancements in Osteosarcoma Therapy: Overcoming Chemotherapy Resistance and Exploring Novel Pharmacological Strategies. Pharmaceuticals (Basel) 2025;18:520. [Crossref] [PubMed]
23. Aoki Y, Kubota Y, Masaki N, et al. Targeting Methionine Addiction of Osteosarcoma with Methionine Restriction to Overcome Drug Resistance: A New Paradigm for a Recalcitrant Disease. Cancers (Basel) 2025;17:506. [Crossref] [PubMed]
24. Liu S, Zhang Q, Liu W, et al. Prediction of Prognosis in Patients With Endometrial Carcinoma and Immune Microenvironment Estimation Based on Ferroptosis-Related Genes. Front Mol Biosci 2022;9:916689. [Crossref] [PubMed]
25. Zhang S, Pei G, Li B, et al. Abnormal phase separation of biomacromolecules in human diseases. Acta Biochim Biophys Sin (Shanghai) 2023;55:1133-52. [Crossref] [PubMed]
26. Ma J, Wei Q, Zhang L, et al. CCT6A functions as promising diagnostic biomarker and promotes cell proliferation in colorectal cancer. J Cancer 2024;15:5897-909. [Crossref] [PubMed]
27. Wang Y, Wang H. CCT6A Regulates Cervical Cancer Cell Glycolysis and Proliferation under Hypoxic Conditions via the Telomerase Cajal Body Protein 1/Telomerase Reverse Tranase. Gynecol Obstet Invest 2025;90:30-41. [Crossref] [PubMed]
28. Zeng W, Wu M, Cheng Y, et al. CCT6A knockdown suppresses osteosarcoma cell growth and Akt pathway activation in vitro. PLoS One 2022;17:e0279851. [Crossref] [PubMed]
29. Yu SK, Yu T, Wang YM, et al. CCT6A facilitates lung adenocarcinoma progression and glycolysis via STAT1/HK2 axis. J Transl Med 2024;22:460. [Crossref] [PubMed]
30. Cui K, Gong L, Zhang H, et al. EXOSC8 promotes colorectal cancer tumorigenesis via regulating ribosome biogenesis-related processes. Oncogene 2022;41:5397-410. [Crossref] [PubMed]
31. Wang Z, Zheng W, Chen Z, et al. Pan-Cancer analysis shows that ACO2 is a potential prognostic and immunotherapeutic biomarker for multiple cancer types including hepatocellular carcinoma. Front Oncol 2022;12:1055376. [Crossref] [PubMed]
32. Mirhadi S, Zhang W, Pham NA, et al. Mitochondrial Aconitase ACO2 Links Iron Homeostasis with Tumorigenicity in Non-Small Cell Lung Cancer. Mol Cancer Res 2023;21:36-50. [Crossref] [PubMed]
33. Bang J, Kang D, Jung J, et al. SEPHS1: Its evolution, function and roles in development and diseases. Arch Biochem Biophys 2022;730:109426. [Crossref] [PubMed]
34. Kang D, Lee J, Jung J, et al. Selenophosphate synthetase 1 deficiency exacerbates osteoarthritis by dysregulating redox homeostasis. Nat Commun 2022;13:779. [Crossref] [PubMed]
35. Yang S, Zhang H, Yang H, et al. SEPHS1 promotes SMAD2/3/4 expression and hepatocellular carcinoma cells invasion. Exp Hematol Oncol 2021;10:17. [Crossref] [PubMed]
36. Hao J, He AY, Zhao X, et al. Pan-Cancer Study of the Prognosistic Value of Selenium Phosphate Synthase 1. Cancer Control 2023;30:10732748231170485. [Crossref] [PubMed]
37. Radko-Juettner S, Yue H, Myers JA, et al. Targeting DCAF5 suppresses SMARCB1-mutant cancer by stabilizing SWI/SNF. Nature 2024;628:442-9. [Crossref] [PubMed]
38. Dreier MR, Walia J, de la Serna IL. Targeting SWI/SNF Complexes in Cancer: Pharmacological Approaches and Implications. Epigenomes 2024;8:7. [Crossref] [PubMed]
39. Liu NQ, Paassen I, Custers L, et al. SMARCB1 loss activates patient-specific distal oncogenic enhancers in malignant rhabdoid tumors. Nat Commun 2023;14:7762. [Crossref] [PubMed]
40. Hanahan D. Hallmarks of Cancer: New Dimensions. Cancer Discov 2022;12:31-46. [Crossref] [PubMed]
41. Wang Q, Shao X, Zhang Y, et al. Role of tumor microenvironment in cancer progression and therapeutic strategy. Cancer Med 2023;12:11149-65. [Crossref] [PubMed]
42. Goenka A, Khan F, Verma B, et al. Tumor microenvironment signaling and therapeutics in cancer progression. Cancer Commun (Lond) 2023;43:525-61. [Crossref] [PubMed]