Bayesian unsupervised clustering identifies clinically relevant osteosarcoma subtypes

Abstract

Identification of cancer subtypes is a critical step for developing precision medicine. Most cancer subtyping is based on the analysis of RNA sequencing (RNA-seq) data from patient cohorts using unsupervised machine learning methods such as hierarchical cluster analysis, but these computational approaches disregard the heterogeneous composition of individual cancer samples. Here, we used a more sophisticated unsupervised Bayesian model termed latent process decomposition (LPD), which handles individual cancer sample heterogeneity and deconvolutes the structure of transcriptome data to provide clinically relevant information. The work was performed on the pediatric tumor osteosarcoma, which is a prototypical model for a rare and heterogeneous cancer. The LPD model detected three osteosarcoma subtypes. The subtype with the poorest prognosis was validated using independent patient datasets. This new stratification framework will be important for more accurate diagnostic labeling, expediting precision medicine, and improving clinical trial success. Our results emphasize the importance of using more sophisticated machine learning approaches (and for teaching deep learning and artificial intelligence) for RNA-seq data analysis, which may assist drug targeting and clinical management.

Keywords: RNA-seq; heterogeneity; latent process decomposition; osteosarcoma; precision medicine.

Figure 1

Latent process decomposition model optimization, subtype assignment and clinical outcome. (a) Hyperparameter optimization for the TARGET dataset. LPD assesses the explanatory power of different combinations of sigma values (process spread) and the number of processes. The optimal combination is determined as the point of maximum log-likelihood before the onset of overfitting, visually identified as a plateau in the curves. For the TARGET dataset, the optimal parameters were three processes and a sigma value of −0.0001. (b) Sample assignment to subtypes. Bar plot illustrates sample assignment to the three identified subtypes based on their degree of membership (gamma value). Higher gamma values indicate stronger membership in a specific subtype reflecting the extent to which each subtype captures sample-specific transcriptomic variability. (c) Kaplan–Meier curves illustrate the survival probability over time for each subtype. Pairwise comparisons between subtypes are shown with log-rank p-values and sample sizes provided for each comparison.

Figure 2

Correlation of gene expression profiles between poor prognosis TARGET LPD-1 and corresponding subtypes. Scatter plots comparing the expression levels of the top 500 most variable transcripts across the entire TARGET dataset between TARGET LPD-1 and the corresponding most similar subtypes from the GREEN (GREEN LPD-1), PERRY (PERRY LPD-2), and SCOTT (SCOTT LPD-1) datasets. Trend lines and Pearson correlation coefficients (r) with corresponding P-values are displayed for each comparison.

Figure 3

Overlap of DE transcripts. Venn diagram illustrating the overlap of DE transcripts between TARGET LPD-1 and the most closely correlated subtypes from the GREEN, PERRY, and SCOTT datasets. The diagram quantifies the number of DE transcripts in each dataset and identifies eight transcripts shared across all four poor prognoses datasets.

Figure 4

Comparative evaluation of traditional clustering methods. (a) Silhouette analysis to determine the optimal number of clusters for hierarchical and k-means clustering in the TARGET dataset. Three clusters were identified as optimal, with six clusters showing similar performance. (b) Kaplan–Meier survival curves comparing patient survival based on hierarchical and k-means clustering groups using both three and six clusters as suggested by the silhouette analysis. Log-rank test was used to assess statistical significance.