Integrative Bulk, Single-Cell, and Spatial Transcriptomic Analyses Nominate PPARG as a Candidate Senescence-Related Prognostic Gene in Osteosarcoma

Summary

Using bulk, single-cell, and spatial transcriptomic analyses, together with machine-learning-based survival modeling and experimental validation, this study nominates PPARG as a candidate senescence-related prognostic gene in osteosarcoma and associates its reduced expression with adverse survival in TARGET-OS and vascular and microenvironmental features.

Abstract

Osteosarcoma remains challenging in metastatic, recurrent, or therapy-resistant disease. This study aimed to identify senescence-related prognostic genes and characterize their spatial contexts. Paired differential expression analysis was performed in GSE99671 using DESeq2, followed by intersection with the CellAge senescence gene set. Transcriptomic data from the TARGET-OS cohort were obtained from UCSC Xena. Candidate genes were evaluated using univariate Cox regression, Kaplan-Meier analysis, time-dependent receiver operating characteristic analysis, LASSO Cox regression, repeated LASSO analysis, and random survival forest modeling, with clinical covariates incorporated into adjusted Cox models. Functional enrichment, immune microenvironment analysis, single-cell transcriptomics, SP_BS3 spatial transcriptomics, GSE36001 expression validation, and qRT-PCR and Western blot validation in osteosarcoma 143B cells and osteoblast cells were used to characterize it. In GSE99671, 2,248 genes were differentially expressed with adjusted P < 0.05, and their intersection with the 866 CellAge genes yielded 105 senescence-related differentially expressed genes. In TARGET-OS, lower PPARG expression was associated with a higher mortality risk (univariate HR = 0.603, 95% CI = 0.454–0.802, P = 0.000494; adjusted HR = 0.224, 95% CI = 0.085–0.589, P = 0.00241). Adding PPARG to the clinical model increased the C-index from 0.707 to 0.829. PPARG was downregulated in both GSE99671 and GSE36001, and qRT-PCR and Western blotting confirmed lower PPARG mRNA and protein expression in osteosarcoma 143B cells than in osteoblast cells. Single-cell analysis localized PPARG to endothelial cells, pericytes, macrophages/monocytes, and tumor-associated stromal cells. Spatial transcriptomic analysis showed weak but significant positive correlations between PPARG expression and the CellAge senescence, endothelial-related, and pericyte-related scores. These findings nominate PPARG as a candidate CellAge-derived prognostic biomarker associated with unfavorable survival and vascular microenvironmental features in osteosarcoma and support investigation of its potential relevance to risk stratification and the senescence-related tumor microenvironment. GSE36001 provided external expression validation only; independent survival validation was not performed.

Introduction

Osteosarcoma is the most common primary malignant bone tumor in children, adolescents, and young adults¹. Although multi-agent chemotherapy combined with surgery has improved outcomes for localized disease^2,3, patients with metastatic, recurrent, or therapy-resistant osteosarcoma continue to have poor long-term survival^3,4. Emerging evidence indicates that oxidative stress-induced epigenetic remodeling can facilitate metastatic adaptation and tumor progression, highlighting the complex molecular plasticity underlying aggressive cancer phenotypes⁵. Robust biomarkers that are both clinically interpretable and biologically informative remain limited. Therefore, identifying molecular features that capture osteosarcoma heterogeneity and prognostic risk across multiple layers of data remains important.

Cellular senescence is a stable cell-cycle arrest program induced by telomere dysfunction, DNA damage, oxidative stress, oncogene activation, and therapeutic pressure⁶. Senescence can restrict aberrant proliferation; however, senescent cells can also reshape the tumor microenvironment through inflammatory, chemokine, growth factor, and extracellular matrix remodeling programs^7,8. In osteosarcoma, senescence-related genes may reflect both tumor-cell-intrinsic stress states and non-malignant microenvironmental compartments, yet their prognostic relevance and spatial organization have not been systematically evaluated.

PPARG encodes peroxisome proliferator-activated receptor gamma, a ligand-activated nuclear receptor involved in lipid metabolism, inflammatory regulation, cellular differentiation, and immune modulation⁹. The role of PPARG in cancer is context-dependent¹⁰. In some settings, PPARG is associated with differentiation and anti-inflammatory states, whereas in others, it may support adaptive tumor or stromal programs. However, its expression pattern, prognostic value, and cellular and spatial localization in osteosarcoma remain incompletely characterized.

In the present study, differentially expressed genes in GSE99671 were identified and intersected with the CellAge senescence gene set, yielding 105 senescence-related differentially expressed genes. TARGET-OS survival data, multiple machine-learning and survival-modeling approaches, and clinical adjustment were then used to nominate PPARG as the core gene. PPARG was further characterized using bulk functional and immune microenvironment analyses, single-cell transcriptomics, spatial transcriptomics, external expression validation from GSE36001, and qRT-PCR and Western blot validation in the osteosarcoma cell line 143B and human osteoblasts.

Protocol

This study used public datasets and cell lines for analysis and validation and did not involve human participants or clinical tissue samples; therefore, ethics approval was not required.

Differential expression analysis in GSE99671 and CellAge intersection
Raw count data and sample grouping information for GSE99671 were obtained from GEO^11,12. Eighteen pairs of osteosarcoma and matched normal samples were analyzed. Paired differential expression analysis was performed using DESeq2 with the design formula ~ pair_id + condition, where pair_id accounted for paired-individual effects and condition compared tumor tissue with normal tissue¹³. Genes with read counts of at least 10 in at least three samples were retained. Differential expression was defined as an adjusted P < 0.05, and a stricter threshold of adjusted P < 0.05 and |log2FC| ≥ 1 was applied for visualization. Differentially expressed genes were intersected with the 866 CellAge senescence genes after gene symbols were converted to uppercase¹⁴. Differential expression analysis was performed in R (version 4.3.2) using DESeq2 (version 1.40.2), and adjusted P values were calculated using the Benjamini-Hochberg method.

TARGET-OS cohort and prognostic modeling
TARGET-OS transcriptomic and clinical data were obtained from UCSC Xena¹⁵. Expression values for the 105 senescence-related differentially expressed genes were extracted. Eighty-five patients with complete survival time, survival status, and candidate gene expression data were included, with 27 death events. Standardized gene expression values were analyzed using univariate Cox regression¹⁶, Kaplan-Meier survival analysis, and time-dependent receiver operating characteristic (ROC) analysis¹⁷. LASSO Cox regression¹⁸, repeated LASSO stability analysis, and random survival forest modeling were used to assess selection stability and variable importance¹⁹. The integrated hub score and clinical-integrated ranking were calculated using the explicit binary criteria detailed below. Survival analyses were performed in R using the survival (version 3.5-7), timeROC (version 0.4), glmnet (version 4.1-8), and randomForestSRC (version 3.2.2) packages. For the univariate Cox screening of the 105 candidate genes, Benjamini-Hochberg false discovery rate (FDR) correction was applied, and genes with FDR < 0.05 were considered statistically significant.

Model preprocessing, PPARG grouping, and time-dependent ROC analysis
Among 85 patients with 27 deaths, genes with zero variance were excluded, missing candidate-gene expression values were median-imputed, and each candidate gene was z-score standardized. For Kaplan-Meier analysis, expression was dichotomized at the cohort median: values strictly above the median were assigned to the high-expression group and values at or below the median to the low-expression group (PPARG: n = 42 high and n = 43 low). Log-rank tests were two-sided. Time-dependent ROC analyses used timeROC with event cause = 1, marginal inverse-probability-of-censoring weighting, evaluation times of 365, 1,095, and 1,825 days, and iid = FALSE. To ensure larger marker values consistently indicate higher risk, standardized expression values were used for genes with positive Cox coefficients and multiplied by -1 for protective genes with negative coefficients.

LASSO and repeated LASSO
Cox LASSO was fitted with glmnet using family = "cox", alpha = 1, prior z-score standardization (therefore standardize = FALSE), five-fold cross-validation, type.measure = "deviance" and random seed 123. The primary coefficient solution used lambda.min. Stability analysis repeated the same five-fold cross-validation 300 times; repeat b used seed 1000 + b (b = 1,...,300). For each gene, selection frequency was the proportion of repeats with a nonzero coefficient at lambda.min; nonzero selection at lambda.1se was also recorded.

Random survival forest
A survival forest was fitted to all 105 standardized candidate genes with randomForestSRC (version 3.2.2), using seed 123, ntree = 1,000, importance = TRUE, and na.action = "na.impute". Package defaults were retained for survival data: log-rank splitting, mtry = 11 (the ceiling of the square root of 105 predictors), minimum terminal node size = 15, nsplit = 10 random split points, sampling without replacement with a 0.632 sampling fraction, and anti-split variable importance.

Integrated hub score
Each of the 105 CellAge-related differentially expressed candidate genes received one point for each of six binary criteria: (1) differential-expression/CellAge intersection membership (all candidates received this point because adjusted P < 0.05 was required before scoring); (2) nominal univariate Cox P < 0.05; (3) Kaplan-Meier log-rank P < 0.05; (4) mean 3- and 5-year time-dependent AUC ≥ 0.65; (5) lambda.min repeated-LASSO selection frequency at or above the 70th percentile across candidates and > 0; and (6) random-survival-forest importance at or above the 70th percentile across candidates and > 0. All criteria had equal unit weight, giving a hub score from 0 to 6; genes scoring ≥ 4 advanced to clinical adjustment (18 genes). Univariate Cox FDR values were also reported, and FDR < 0.05 was used to denote multiple-testing significance, but the prespecified scoring indicator used nominal P < 0.05.

Clinical-integrated ranking
After merging expression and clinical records, the adjusted analyses included 40 patients with complete covariate records and 13 deaths. The final score was the initial hub score plus one point for each of seven criteria: adjusted Cox P < 0.05, adjusted Cox P < 0.10, sensitivity Cox P < 0.05 after removing definitive surgery, sensitivity Cox P < 0.10, mean 3- and 5-year AUC ≥ 0.65, likelihood-ratio-test P < 0.10 for the clinical-plus-gene versus clinical-only model, and delta AIC < 0. Because the 0.05 and 0.10 thresholds were nested, a P < 0.05 contributed two points, thereby giving greater weight to clearly significant adjusted and sensitivity Cox evidence. The total range was 0 to 13; ties were resolved by the smaller adjusted Cox P value, followed by the larger mean 3-/5-year AUC. PPARG received all six initial points and all seven clinical integration points (13/13), ranking first. The C-index improvement was reported descriptively and was not assigned a score.

Clinical adjustment
Candidate hub genes were integrated with TARGET-OS clinical variables, including sex, age, disease status at diagnosis, primary tumor site, specific tumor region and definitive surgery. Forty patients with complete expression and clinical records, including 13 deaths, entered the clinical-adjusted analyses. Cox models containing only clinical variables were compared with models including clinical variables plus gene expression. C-index, Akaike information criterion (AIC), and likelihood ratio test P values were used to evaluate model improvement. Sensitivity analysis was performed after removing the surgery variable.

Functional enrichment and immune microenvironment analysis
TARGET-OS samples were stratified according to PPARG expression. Differential expression between the PPARG-high and PPARG-low groups was used to generate ranked gene lists for gene set enrichment analysis (GSEA)²⁰. The pathways displayed were Nemeth Inflammatory Response LPS Up, Burton Adipogenesis 5, Burton Adipogenesis 6, Krieg KDM3A Targets Not Hypoxia, Reactome: Transcriptional Regulation By TP53, Fulcher Inflammatory Response Lectin Vs LPS Dn, Hollmann Apoptosis Via CD40 Dn, Zhou Inflammatory Response Live Dn, WP: Fatty Acids And Lipoproteins Transport In Hepatocytes, Sweet Lung Cancer KRAS Up, KEGG Medicus Pathogen HIV Tat To TLR2/4 NF-kB Signaling Pathway and Reactome: Fatty Acids. A bulk CellAge senescence score was calculated, and Spearman correlation analysis was used to assess associations between PPARG and senescence genes or immune microenvironment signatures²¹. Differences in microenvironment scores between PPARG-high and PPARG-low groups were evaluated using non-parametric tests with multiple-testing correction.

Single-cell transcriptomic analysis
A published human osteosarcoma single-cell transcriptomic dataset was analyzed using a preprocessed object for which quality control, dimensionality reduction, clustering, and manual annotation had already been completed^22,23,24. A total of 68,336 cells and 32,297 genes were included. For main-text interpretation, annotations were simplified into 13 major cell types: B cells, CAFs, cycling cells, endothelial cells, erythroid cells, macrophages/monocytes, malignant osteosarcoma cells, myogenic cells, neutrophils, osteoclast-like cells, pericytes, T/NK cells, and tumor-associated stromal cells. Dimensionality reduction, feature expression, dot, and violin plots were used to visualize PPARG localization. PPARG-positive cells were defined as cells with expression greater than zero. Differences among cell types were assessed using Kruskal-Wallis and Wilcoxon rank-sum tests with Benjamini-Hochberg correction. Single-cell analyses were performed in R using Seurat (version 5.0.1).

Spatial transcriptomic analysis
The SP_BS3 spatial transcriptomic sample was used to construct a spatial expression object^25,26. Quality control thresholds were nFeature_Spatial ≥ 200 and percent.mt ≤ 30, leaving 4,572 spots for analysis. Data were normalized, 3,000 highly variable genes were selected, and data scaling, principal component analysis, neighborhood graph construction, spatial spot clustering, and dimensionality reduction were performed. A spatial CellAge senescence score was calculated after removing PPARG from the gene set to avoid circular correlation. Endothelial, pericyte, macrophage/monocyte, tumor-associated stromal, malignant osteosarcoma, and osteoclast-like signatures were constructed and scored. Spearman correlation analysis was used to assess associations between PPARG expression and spatial scores. Label transfer was performed using the single-cell dataset as the reference and the spatial dataset as the query to infer predicted cell-type scores for each spot^23,24. Spatial transcriptomic analyses were performed in R using Seurat (version 5.0.1), and Benjamini-Hochberg FDR correction was applied to all spatial correlation P values.

External expression validation in GSE36001
The GEO dataset GSE36001 was used exclusively as an independent expression validation cohort; because survival outcomes were unavailable, it was not used for prognostic validation^11,27. The dataset included 19 osteosarcoma samples and six normal control samples. GPL6102 platform annotation was used to convert probe identifiers to gene symbols. When multiple probes mapped to the same gene, the probe with the highest average expression was retained. Differential expression between tumor and normal groups was assessed using limma²⁸. Analyses were performed in R using limma (version 3.56.2), and adjusted P values were calculated using the Benjamini-Hochberg method.

qRT-PCR and Western blot validation
Experimental validation was performed using the human osteosarcoma cell line 143B and human osteoblast cells. Osteosarcoma cells were cultured in Dulbecco's modified Eagle's medium supplemented with 10% fetal bovine serum and 1% penicillin-streptomycin at 37 °C in a humidified atmosphere containing 5% CO₂, and passaged with 0.25% trypsin-EDTA at 80%–90% confluence. Human osteoblast cells were maintained under the recommended culture conditions. All cell lines were confirmed to be free of mycoplasma contamination. For qRT-PCR, total RNA was extracted using a phenol-guanidinium-based RNA extraction reagent, and RNA concentration and purity were assessed spectrophotometrically. One microgram of total RNA was reverse-transcribed using a reverse-transcription reagent according to the recommended protocol. qRT-PCR was performed using fluorescent DNA-binding dye-based chemistry under the following cycling conditions: initial denaturation at 95 °C for 30 s, followed by 40 cycles at 95 °C for 5 s and 60 °C for 30 s, with melt-curve analysis to confirm amplification specificity. Each reaction was performed in technical triplicate, and three independent biological experiments were conducted. GAPDH was used as the internal control, and relative PPARG expression was calculated using the 2^-ΔΔCt method²⁹. The PPARG forward primer was 5'-CGAAGACATTCCATTCACAAGAACAG-3', and the reverse primer was 5'-AGATGCAGGCTCCACTTTGATTG-3'.

Western blotting was performed to examine PPARG protein expression. Cells were lysed in radioimmunoprecipitation assay buffer supplemented with protease inhibitors, and protein concentrations were determined using a bicinchoninic acid assay. Equal amounts of protein (30 µg per lane) were separated by 10% sodium dodecyl sulfate-polyacrylamide gel electrophoresis and transferred onto polyvinylidene difluoride membranes. After blocking with 5% non-fat milk for 1 h at room temperature, membranes were incubated overnight at 4 °C with primary antibodies against PPARG (1:1,000) and GAPDH (1:5,000), followed by incubation with a horseradish peroxidase-conjugated secondary antibody (1:5,000) for 1 h at room temperature. Protein bands were visualized using chemiluminescence detection, and three independent experiments were performed. Band intensities were quantified using image-analysis software³⁰. Between-group differences were analyzed using two-tailed unpaired Student's t tests. Data were presented as the mean ± standard deviation (SD) from three independent experiments, and P < 0.05 was considered statistically significant. Statistical analyses of experimental data were performed using statistical analysis software (version 9.0).

Statistical analysis
Unless otherwise specified, all bioinformatic analyses were performed in R (version 4.3.2). Two-sided P values < 0.05 were considered statistically significant. Correlations were assessed using Spearman's rank correlation coefficient (ρ). Multiple-testing correction was performed using the Benjamini-Hochberg FDR method where applicable. Experimental data were presented as the mean ± SD and were compared using two-tailed unpaired Student's t tests. Experimental statistical analyses were performed using statistical analysis software (version 9.0).

Results

GSE99671 identified 105 CellAge senescence-related differentially expressed genes
GSE99671 included 36 samples from 18 paired tissues. After low-count filtering, 16,683 genes were retained. At an adjusted P < 0.05, 2,248 genes were differentially expressed. Under the stricter threshold of adjusted P < 0.05 and |log2FC| ≥ 1, 594 genes were significant, including 102 upregulated and 492 downregulated genes in tumors (Figure 1A,B). I...

Discussion

This study nominated PPARG as a candidate senescence-related prognostic gene in osteosarcoma by integrating differential expression analysis, CellAge gene intersection, TARGET-OS survival modeling, clinical adjustment, multi-omic localization, and experimental validation. This analytical framework was consistent with the current emphasis on molecularly informed osteosarcoma research and curated senescence-gene resources for interpreting senescence-related biology^14,31<...

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
Anti-GAPDH primary antibody	Proteintech Group, Wuhan, China	10494-1-AP	Primary antibody used to detect GAPDH as a loading control in western blotting; dilution 1:5,000.
Anti-PPARG primary antibody	Proteintech Group, Wuhan, China	16643-1-AP	Primary antibody used for detection of PPARG protein by western blotting; dilution 1:1,000.
BCA protein assay kit	Beijing Solarbio Science & Technology Co., Ltd., Beijing, China	PC0020	Colorimetric assay used to determine total protein concentration before electrophoresis.
DESeq2	Bioconductor	Version 1.40.2	R package used for differential gene-expression analysis of count-based transcriptomic data.
Dulbecco's modified Eagle's medium (DMEM)	Beijing Solarbio Science & Technology Co., Ltd., Beijing, China	11995	Basal culture medium used for maintenance of 143B osteosarcoma cells.
ECL detection reagent	Beijing Solarbio Science & Technology Co., Ltd., Beijing, China	PE0010	Chemiluminescent substrate used for detection of HRP-conjugated antibodies in western blotting.
Fetal bovine serum (FBS)	Zhejiang Tianhang Biotechnology Co., Ltd. (Sijiqing), Huzhou, Zhejiang, China	11011-8611	Serum supplement added to culture medium to support cell growth and viability.
glmnet	CRAN	Version 4.1-8	R package used for penalized regression analyses, including LASSO and elastic-net modeling.
GraphPad Prism	GraphPad Software, San Diego, CA, USA	Version 9.0	Software used for statistical analysis, graph generation, and visualization of experimental data.
HRP-conjugated secondary antibody	Proteintech Group, Wuhan, China	SA00001-2	Horseradish peroxidase-conjugated secondary antibody used for western blot detection; dilution 1:5,000.
Human osteoblast cells	Cell Applications, Inc., San Diego, CA, USA	406-05A	Primary human osteoblast cells used as the nonmalignant comparison/control cell type.
Human osteosarcoma 143B cell line	American Type Culture Collection (ATCC), Manassas, VA, USA	CRL-8303	Human osteosarcoma cell line used for in vitro validation experiments and molecular assays.
ImageJ	National Institutes of Health (NIH), USA	Version 1.53	Image-analysis software used for quantitative analysis of experimental images.
limma	Bioconductor	Version 3.56.2	R package used for differential expression analysis and linear-model-based statistical testing.
Penicillin-streptomycin	Beijing Solarbio Science & Technology Co., Ltd., Beijing, China	P1400	Antibiotic supplement used in cell culture medium to reduce bacterial contamination.
PPARG and GAPDH primers	Sangon Biotech (Shanghai) Co., Ltd., Shanghai, China	Custom synthesis; sequences provided in the Methods	Custom oligonucleotide primers used for qPCR analysis of PPARG and GAPDH expression.
PVDF membranes, 0.45 µm	Beijing Solarbio Science & Technology Co., Ltd., Beijing, China	YA1701	Membranes used for protein transfer during western blotting.
R statistical software	R Foundation for Statistical Computing, Vienna, Austria	Version 4.3.2	Statistical computing environment used for bioinformatic analyses, model construction, and visualization.
randomForestSRC	CRAN	Version 3.2.2	R package used for random survival forest modeling and feature importance analysis.
Reverse transcription kit	Beyotime Biotech Inc., Shanghai, China	D7168M	Used for synthesis of complementary DNA (cDNA) from isolated RNA before quantitative PCR.
RIPA lysis buffer	Beijing Solarbio Science & Technology Co., Ltd., Beijing, China	R0010	Protein-extraction buffer used to lyse cells for western blot analysis.
Seurat	Satija Laboratory	Version 5.0.1	R package used for single-cell RNA-sequencing data processing, integration, clustering, and visualization.
survival	CRAN	Version 3.5-7	R package used for survival analysis, including Cox proportional hazards modeling.
SYBR Green qPCR master mix	Beijing Solarbio Science & Technology Co., Ltd., Beijing, China	SR1110	Fluorescent master mix used for quantitative real-time PCR amplification.
timeROC	CRAN	Version 0.4	R package used to generate time-dependent receiver operating characteristic curves and calculate predictive performance over time.
TRIzol reagent	Invitrogen, Thermo Fisher Scientific, Waltham, MA, USA	15596026CN	Reagent used for extraction of total RNA from cultured cells.
Trypsin-EDTA, 0.25%	Beijing Solarbio Science & Technology Co., Ltd., Beijing, China	T1300	Cell-dissociation reagent used for passaging and harvesting adherent cells.