Bayesian Unsupervised Clustering Identifies Clinically Relevant Osteosarcoma Subtypes

Briefings in Bioinformatics 2024 AI 8 Explanations View Original
Original Paper (PDF)

Unable to display PDF. Download it here or view on PMC.

Plain-English Explanations
Pages 1-2
Why Osteosarcoma Subtyping Has Been So Difficult

Osteosarcoma is the most common primary bone cancer in children and adolescents, affecting roughly 3-4 per million individuals globally each year. It typically originates in the metaphyseal intramedullary cavity of long bones and shows a bimodal age distribution, with the first peak occurring during adolescent growth spurts and a second, smaller peak in adults over 60. Despite being the dominant bone sarcoma in the pediatric and adolescent and young adult (AYA) population, it remains rare in absolute terms, which has historically made building large clinical trial cohorts exceptionally difficult.

The stagnant survival problem: The five-year survival rate for osteosarcoma has been stuck at approximately 50% for the last 45 years. The standard chemotherapy regimen, known as MAP (methotrexate, adriamycin/doxorubicin, and cisplatin), has not improved substantially since the 1970s. Multiple Phase I and II clinical trials investigating new agents have consistently failed to advance to Phase III, in part because patients with osteosarcoma are recruited as a single undifferentiated group rather than being stratified by molecular subtype. Without subtyping, it is impossible to know which drug might work for which patient.

The heterogeneity challenge: Osteosarcoma tumors are biologically heterogeneous, meaning individual tumors contain mixtures of different cell states and transcriptional programs. Classical unsupervised clustering algorithms, such as hierarchical cluster analysis and k-means, assume that each tumor sample belongs exclusively to one cluster. This assumption is fundamentally incompatible with a disease where a single tumor may simultaneously exhibit features of multiple biological subtypes. The result is that prior subtyping studies, while suggestive, have systematically understated intratumoral heterogeneity.

This paper, from researchers at the University of East Anglia, addresses that limitation directly. By applying a Bayesian probabilistic model called latent process decomposition (LPD) to RNA sequencing data from 165 primary osteosarcomas, the authors identify three clinically meaningful molecular subtypes, one of which is associated with significantly worse outcomes under MAP chemotherapy.

TL;DR: Osteosarcoma's 5-year survival has not improved past 50% in 45 years. Classical clustering fails because individual tumors are heterogeneous mixtures of cell states. This paper uses Bayesian latent process decomposition (LPD) on RNA-seq data from 165 primary osteosarcomas across four datasets to identify three clinically relevant subtypes, one of which shows poor response to MAP chemotherapy.
Pages 2-3
Latent Process Decomposition: A Bayesian Mixed-Membership Model for Heterogeneous Tumors

Latent process decomposition (LPD) is a Gaussian mixed-membership Bayesian model that was originally developed for analyzing cDNA microarray data. Unlike hard-assignment clustering methods, which force each sample into a single predefined cluster, LPD treats each sample as a probabilistic mixture of multiple underlying "processes." The gene expression profile of any individual tumor is modeled as a weighted combination of these latent signatures, with weights drawn from a Dirichlet distribution. A sample can therefore be 70% Process 1 and 30% Process 2, reflecting its true biological heterogeneity rather than artificially assigning it to just one group.

Input data and model optimization: The authors used the TARGET osteosarcoma RNA-seq dataset (phs000468, n = 88 primary tumors from the Therapeutically Applicable Research to Generate Effective Treatments initiative) as the discovery cohort. The expression data were reduced to the top approximately 500 transcripts showing the greatest variance across samples. Two model hyperparameters were optimized simultaneously: the number of latent processes and a sigma parameter controlling process spread. Optimization was performed by assessing hold-out validation log-likelihood across a grid of hyperparameter combinations, with the optimal point selected at the plateau before overfitting began. For the TARGET dataset, the optimal configuration was three processes with a sigma value that prevented over-splitting.

Validation datasets: The same 500 high-variance transcripts and optimized LPD configuration were applied to three independent external datasets: GREEN (n = 14, including primary and metastatic samples plus circulating tumor cells), PERRY (n = 35 primary tumors sequenced on Illumina HiSeq 2000), and SCOTT (n = 52 samples including primaries, metastatic tissue, cell lines, and normal bone). All four datasets used Illumina HiSeq 2000 platforms except TARGET, which used an Affymetrix Human Exon ST Array, providing a cross-platform validation test.

Differential expression and pathway analysis: Samples were assigned to the process in which they had the highest membership weight. Differentially expressed (DE) transcripts were identified using DESeq2 with pairwise comparisons between the three subtypes. Inclusion criteria were a log2 fold change greater than 1.5, an adjusted P-value below 0.01 (Benjamini-Hochberg correction), and detection in at least 80% of model runs. Gene Ontology (GO) enrichment analysis used clusterProfiler 4.0. Statistical comparisons of clinical variables used ANOVA with Bonferroni correction for continuous variables and chi-square tests for categorical variables.

TL;DR: LPD assigns each tumor a probabilistic mixture of latent processes rather than forcing hard cluster assignment. Discovery used TOP 500 variance transcripts from 88 TARGET RNA-seq primaries; hyperparameters optimized by hold-out log-likelihood. Validated in three independent datasets (GREEN n=14, PERRY n=35, SCOTT n=52). DE transcripts required log2FC greater than 1.5, adjusted P less than 0.01, present in at least 80% of runs. Analysis performed in R v4.3.2 using DESeq2 and clusterProfiler 4.0.
Pages 3-5
Three Osteosarcoma Subtypes with Distinct Survival and Clinical Profiles

LPD identified three osteosarcoma processes in the TARGET discovery dataset. Samples were classified by dominant process membership into TARGET LPD-1 (n = 26), TARGET LPD-2 (n = 39), and TARGET LPD-3 (n = 23). The three subtypes differed significantly in sex distribution (chi-square P = 0.01) and vital status (P = 0.004), but not in age at diagnosis (median ages: 12.5, 15, and 15 years respectively, P = 0.184) or ethnicity (P = 1.00).

The poor prognosis subtype: TARGET LPD-1 emerged as the most clinically distinct group. It showed a marked female predominance (17 of 26 patients female) and, critically, the worst vital status distribution: 16 of 26 patients in LPD-1 were dead at the time of analysis, compared to only 7 of 39 in LPD-2 and 6 of 23 in LPD-3. Kaplan-Meier survival analysis confirmed that LPD-1 membership was associated with significantly inferior survival probability. Cox regression incorporating process assignment as a covariate further substantiated the independent prognostic value of subtype membership.

Favorable and intermediate subtypes: TARGET LPD-2 represented the largest subgroup and showed a predominantly male composition (30 of 39 patients) with the best survival outcomes across the cohort. TARGET LPD-3 showed intermediate clinical characteristics. The biological distinctions between these three groups were not simply a reflection of known clinical covariates such as stage or metastatic status, suggesting that the LPD-detected transcriptional signatures capture meaningful underlying biology beyond what is currently assessed in routine clinical practice.

The identification of a subtype with significantly worse outcomes under the MAP chemotherapy backbone has direct clinical implications. Patients classified as LPD-1 at diagnosis may be candidates for treatment intensification, alternative regimens, or early enrollment in clinical trials targeting the molecular vulnerabilities specific to that subtype, rather than receiving the same standard therapy applied uniformly to all osteosarcoma patients.

TL;DR: LPD identified three subtypes in 88 TARGET primaries: LPD-1 (n=26, poor prognosis, 16/26 dead, female-predominant), LPD-2 (n=39, best outcomes, male-predominant), and LPD-3 (n=23, intermediate). Sex distribution differed significantly (P=0.01) and vital status differed significantly (P=0.004). Age and ethnicity did not differ. LPD-1 patients may warrant alternative treatment strategies beyond standard MAP chemotherapy.
Pages 5-7
679 Differentially Expressed Transcripts and Eight Cross-Dataset Biomarker Candidates

To characterize the molecular underpinnings of the three LPD subtypes, the authors performed pairwise DE analysis across all three group comparisons. They identified 679 significantly DE transcripts: 336 upregulated and 343 downregulated, all exceeding a log2 fold change of 1.5, with adjusted P-values below 0.01, and detected in at least 80% of LPD model runs. Gene Ontology enrichment analysis identified 69 significantly altered biological processes associated with these transcripts.

Extracellular matrix dominance: The most overrepresented biological processes involved extracellular matrix (ECM) structure and organization pathways. ECM remodeling is well-established as a driver of invasion, metastasis, and treatment resistance in sarcomas. The enrichment of ECM-related transcripts in the poor prognosis LPD-1 subtype points toward a tumor microenvironment characterized by stromal remodeling that may facilitate disease progression and limit drug penetration.

The core eight-gene signature: By applying the LPD model to the three independent validation datasets and cross-referencing the DE transcripts of poor prognosis subtypes across all four cohorts (TARGET LPD-1, GREEN LPD-1, PERRY LPD-2, SCOTT LPD-1), the authors identified eight transcripts consistently dysregulated across all poor prognosis groups. This cross-dataset core signature represents promising candidates for biomarker development and therapeutic targeting. Several of the eight genes have previously been implicated in osteosarcoma biology: ANGPT1 upregulation has been associated with non-metastatic disease, CGREF1 overexpression is linked to poor prognosis, KAZALD1 is overexpressed in osteoblastic osteosarcoma and cancer-associated fibroblasts, and CILP is downregulated in cartilage tumors.

Targetability potential: The identification of consistently dysregulated transcripts across independent datasets transforms the LPD subtypes from purely descriptive classifications into actionable biological categories. Genes appearing in the eight-gene core signature are candidates for follow-up functional studies, drug target validation, and ultimately incorporation into precision medicine frameworks that match specific molecular vulnerabilities to approved or experimental agents.

TL;DR: 679 DE transcripts (336 up, 343 down, log2FC greater than 1.5, adjusted P less than 0.01) characterized the three subtypes. GO enrichment identified 69 altered biological processes, with ECM pathways most overrepresented. Cross-dataset comparison of poor prognosis subtypes across all four cohorts identified 8 consistently dysregulated transcripts including ANGPT1, CGREF1, KAZALD1, and CILP as biomarker and drug target candidates.
Pages 7-8
LPD Subtypes Replicate Across Three Independent Datasets

External validation used the GREEN, PERRY, and SCOTT datasets applying the same 500-transcript input and the same hyperparameter optimization approach as in the TARGET discovery analysis. Because LPD is applied separately to each dataset, the model identifies the optimal number of processes independently for each cohort rather than forcing a predetermined three-subtype solution. This tests whether the biological signals identified in TARGET are real and reproducible, not artifacts of a particular dataset or the forced imposition of a specific cluster number.

GREEN dataset results: LPD identified two subtypes in GREEN (n = 14). GREEN LPD-1 (n = 10) was predominantly composed of primary and metastatic tissue samples (chi-square P = 0.006). GREEN LPD-2 (n = 4) exclusively contained circulating tumor cell samples. The gene expression profile of GREEN LPD-1 correlated strongly with TARGET LPD-1 (Pearson correlation P less than 0.001), indicating that the poor prognosis transcriptional signature is detectable even in a small mixed-composition cohort. No survival data were available for the GREEN dataset.

PERRY and SCOTT dataset results: LPD identified three subtypes in PERRY (n = 35) and four in SCOTT (n = 52). In PERRY, Kaplan-Meier analysis showed a non-significant trend toward poorer survival for PERRY LPD-1 and LPD-2 compared to LPD-3 (P = 0.15), likely reflecting limited statistical power with only 35 patients. In SCOTT, significant age differences were found across subtypes (ANOVA P = 0.005), with SCOTT LPD-1 and LPD-3 being significantly younger than LPD-4. SCOTT LPD-2 was enriched for osteosarcoma cell lines (chi-square P less than 0.001), which LPD correctly separated from primary tissue samples. The poor prognosis subtypes PERRY LPD-2 and SCOTT LPD-1 both showed high gene expression similarity to TARGET LPD-1 (Pearson correlation P less than 0.001 for both).

The cross-platform replication is particularly notable: TARGET used an Affymetrix array platform while the validation datasets used Illumina sequencing, yet the core transcriptional signatures were reproduced. This suggests that the osteosarcoma subtypes detected by LPD reflect robust biological differences that transcend technical variation in measurement platform.

TL;DR: LPD replicated poor prognosis subtypes across all three validation datasets. GREEN LPD-1 (n=10) correlated with TARGET LPD-1 (P less than 0.001). PERRY showed non-significant survival trend (P=0.15, underpowered). SCOTT identified 4 subtypes with significant age differences (P=0.005); poor prognosis PERRY LPD-2 and SCOTT LPD-1 both correlated with TARGET LPD-1 (P less than 0.001). Replication held across different sequencing platforms (Affymetrix vs. Illumina).
Pages 8-9
LPD Outperforms Hierarchical and k-Means Clustering for Prognostic Stratification

A core claim of the paper is that LPD offers meaningful advantages over the classical unsupervised methods most commonly used for cancer subtyping. To test this directly, the authors applied both hierarchical clustering (with Ward's D2 linkage) and k-means clustering (with the Hartigan-Wong algorithm and 50 restarts) to the same TARGET dataset used for LPD discovery. The optimal number of clusters for each method was determined using silhouette score analysis.

Silhouette analysis and cluster selection: For both hierarchical and k-means methods, silhouette analysis identified three clusters as optimal, with six clusters as a secondary suggestion with similar performance. This is an important methodological detail: even when applying the same target cluster count as LPD (three), the traditional methods failed to recapitulate the clinically meaningful subtypes identified by the Bayesian approach.

Survival curve comparison: Kaplan-Meier survival curves generated from the hierarchical and k-means cluster assignments showed no significant differences in survival between any of the resulting groups (P greater than 0.05 for all comparisons across both three-cluster and six-cluster solutions). This is a striking contrast to LPD, which identified a subtype with significantly inferior survival. The failure of traditional methods is attributable precisely to their hard-assignment assumption: by forcing heterogeneous samples into single clusters, they average out the biological signal that LPD preserves through mixed membership.

This head-to-head comparison provides the clearest evidence for the paper's central argument. When the same input data are processed by competing algorithms, only LPD identifies prognostically relevant osteosarcoma subtypes. The implication is that the field's previous reliance on hierarchical or k-means clustering may have masked real biological heterogeneity and contributed to the failure to identify clinically actionable subgroups in prior studies.

TL;DR: Hierarchical clustering (Ward's D2) and k-means (Hartigan-Wong, 50 restarts) both identified 3 clusters as optimal by silhouette score on the same TARGET data. Neither method produced significant survival differences between clusters (P greater than 0.05). LPD identified prognostically significant subtypes from the same input data. The hard-assignment assumption of classical methods averages out the biological signal that LPD's mixed-membership model preserves.
Pages 9-10
Small Datasets, Incomplete Clinical Annotation, and the Rarity Problem

Limited discovery cohort size: The LPD model was developed on the TARGET dataset containing 88 primary osteosarcoma samples. While this is one of the larger publicly available osteosarcoma RNA-seq datasets, it remains small by the standards of AI and machine learning studies in common cancers such as breast or colorectal cancer, where training sets of thousands of samples are routine. Small discovery cohorts increase the risk that the identified subtypes reflect dataset-specific noise rather than universal osteosarcoma biology, and they limit the precision of DE analysis and GO enrichment findings.

Incomplete clinical annotation in validation cohorts: The external validation datasets had incomplete clinical data. The GREEN dataset provided no survival information at all. The PERRY dataset had only 35 patients with available survival data, which was insufficient to achieve statistical significance even when survival trends aligned with the expected direction (P = 0.15). The SCOTT dataset lacked survival data entirely, limiting validation to gene expression and sample composition comparisons. The inability to independently confirm the survival association in adequately powered validation cohorts is the most important clinical limitation of the study.

Cross-dataset technical challenges: The TARGET dataset used Affymetrix arrays while the validation datasets used Illumina HiSeq 2000, introducing potential batch effects and platform-specific biases. Although the authors found that core transcriptional signatures replicated across platforms, the differing library preparation methods and sequencing parameters required careful curation of input features. The requirement to use only transcripts detected across all platforms constrains the expressiveness of the model and may exclude biologically relevant genes that are platform-specifically detectable.

Access to tissue and biological material: Osteosarcoma is a rare disease, and the authors explicitly note that obtaining paired tissue and comprehensive clinical follow-up data is particularly challenging. Biopsy material is often limited in quantity, and posttreatment samples commonly exhibit extensive chemotherapy-related damage that degrades RNA quality. These practical constraints on data collection are unlikely to be resolved quickly and will continue to limit the scale of osteosarcoma computational studies in the near term.

TL;DR: Key limitations: discovery cohort limited to 88 primary tumors (TARGET); validation cohorts had no survival data (GREEN, SCOTT) or underpowered survival analysis (PERRY n=35, P=0.15). Cross-platform differences between Affymetrix (TARGET) and Illumina (validation datasets) required careful feature curation. Disease rarity and limited tissue access are structural barriers to larger validation studies.
Pages 10-11
From Subtype Discovery to Clinical Stratification and Precision Treatment

Clinical trial stratification: The authors' central translational argument is that LPD-based subtyping should be used to stratify patients in future osteosarcoma clinical trials rather than enrolling them as a uniform group. The 45-year stagnation in survival outcomes is at least partly attributable to testing new drugs against an insufficiently characterized patient population. A trial that enrolls patients based on LPD subtype membership, treating LPD-1 patients differently from LPD-2 or LPD-3 patients, would have a fundamentally better chance of detecting subtype-specific treatment effects. This mirrors the transformation that molecular subtyping brought to breast cancer trials, where separating HER2-positive from triple-negative from luminal subtypes dramatically improved the specificity and success rate of therapeutic development.

Expanding multi-omics integration: RNA-seq transcriptomics captures gene expression but not somatic mutations, copy number alterations, DNA methylation, or protein-level changes. Emerging computational platforms, including iCluster, Bayesian consensus clustering, and copula mixed models, can integrate multiple genomic data types simultaneously. Future extensions of LPD to osteosarcoma could incorporate mutation profiles, methylation arrays, and proteomic data to produce richer and more stable subtype definitions. The Cancer Genome Atlas and similar platforms already provide multi-omics data for several cancer types, and a similarly comprehensive osteosarcoma multi-omics resource would substantially advance this agenda.

Drug response prediction and the eight-gene signature: The eight consistently dysregulated transcripts identified across all four poor prognosis subtypes are high-priority candidates for functional follow-up studies. Validating these genes as druggable targets using cell lines and patient-derived organoids, then testing targeted agents against LPD-stratified tumor models, represents a concrete experimental pipeline. The authors note that LPD-based subtyping could serve as a companion diagnostic framework for predicting response to new and experimental drugs, transforming osteosarcoma drug development from empirical chemotherapy trials into hypothesis-driven precision oncology.

Rare cancer computational infrastructure: The broader lesson from this work is that rare cancers require computational methods specifically adapted for small, heterogeneous datasets. Standard machine learning approaches developed and optimized for common cancers with thousands of samples do not transfer straightforwardly to diseases where the largest single-institution cohorts contain tens to low hundreds of cases. Bayesian approaches that encode prior biological knowledge and handle uncertainty explicitly are better suited to the rare cancer context. The authors advocate for wider adoption of such methods across the sarcoma field and beyond, particularly for pediatric and AYA cancers where rarity, heterogeneity, and funding constraints compound the analytical challenges.

TL;DR: Key future directions: use LPD subtype as a stratification variable in osteosarcoma clinical trials rather than enrolling all patients uniformly; integrate multi-omics data using iCluster or similar frameworks; pursue functional validation of the 8-gene cross-dataset signature as drug targets; develop Bayesian and probabilistic methods designed specifically for rare, small-cohort cancers where standard ML approaches underperform.