Osteosarcoma (OS) is the most common primary bone malignancy in adolescents and young adults. While neoadjuvant chemotherapy has improved outcomes, the 5-year survival rate for patients with metastases at diagnosis remains less than 30%. Identifying therapeutic targets at the molecular level requires moving beyond individual gene expression analysis to understanding how sets of interacting genes (modules) are dysregulated in each patient.
Traditional gene expression studies compare group averages between disease and control cohorts, but this population-level approach cannot capture the molecular heterogeneity that is characteristic of osteosarcoma and other rare tumors. A personalized approach, comparing a single patient's expression profile against a large reference dataset of normal samples, enables the identification of patient-specific aberrations that would be missed in group-level analyses.
The individualized module aberrance score (iMAS) method addresses this need by integrating protein-protein interaction (PPI) network structure with expression data to score the dysregulation of gene modules in individual patients. This approach was previously validated in colorectal and lung cancer using the Average Z equation method, and this study applies it to osteosarcoma for the first time.
Gene expression data were obtained from the EMBL-EBI ArrayExpress database. Two datasets were combined: E-GEOD-33382 and E-GEOD-28974, providing a total of 99 samples. The study used 15 normal bone tissue samples as the control group and 84 osteosarcoma patient samples as the disease group. The expression platform was the Illumina Human-6 v2.0 expression BeadChip using nuIDs as identifiers.
Data preprocessing used the Linear Models for Microarray Data (LIMMA) package. Background adjustment and quantile normalization were performed using the robust multi-array average (RMA) method. An outlier-robust median polish procedure was used to estimate model parameters, and where multiple probes mapped to the same gene symbol, expression values were averaged, yielding 23,214 unique genes for analysis.
The human protein-protein interaction (PPI) map was obtained from the STRING database version 9.0, which contained 1,048,576 PPI pairs. Filtering to high-confidence interactions (combined score greater than or equal to 0.8) and intersecting with the 23,214 expressed genes yielded a PPI subnet of 37,381 interactions across 6,665 nodes. Edge weights were assigned as the absolute Spearman correlation coefficient between connected gene pairs based on expression data.
Modules were identified from the PPI network using a clique-merging algorithm. All maximal cliques of size at least k were enumerated via a fast depth-first search (CLIQUES algorithm). Each clique was scored by a weighted interaction density score: the sum of all edge weights within the clique divided by the number of possible edges, reflecting both connectivity and expression co-regulation. Cliques were sorted by score in descending order, and those with pairwise overlap above a threshold of 0.5 were merged into modules, yielding a final set of discrete, non-redundant gene modules.
The iMAS (individualized module aberrance score) was computed for each patient-module combination. For each individual OS patient, their expression data were quantile-normalized against the combined dataset, and a Z-score was computed for each gene relative to the accumulated normal samples (ANS) cohort of 15 control tissues. The module-level iMAS was defined as the average Z-score across all genes in the module: iMAS = (sum of Z-scores) / n, where n is the number of genes in the module. This creates a patient-specific signature of module-level dysregulation.
Module significance was assessed by Wilcoxon test comparing module scores between OS patients and normal samples, with p-values corrected for multiple testing using the false discovery rate (FDR). The significance thresholds were FDR less than 0.05 and log fold change greater than 1. Modules reaching adjusted p less than 0.01 were designated as significantly differentially expressed. GO biological process enrichment was performed using BiNGO in Cytoscape version 3.2.0 at p less than 0.01.
The clique-merging algorithm identified 83 modules containing 2,084 genes from the osteosarcoma PPI network. Of these, 61 modules were significantly different between OS patients and normal controls by Wilcoxon test at the adjusted p less than 0.01 threshold. These 61 modules represent a compressed representation of the osteosarcoma transcriptomic landscape: rather than analyzing thousands of individual genes, clinical information is captured by a manageable set of functional gene clusters.
GO biological process analysis of the 61 significant modules identified 1,568 enriched biological processes, with the most significant terms reaching p values below 1.0 x 10^-36. The relation network of biological processes at p less than 1.0 x 10^-31 revealed a densely connected cluster of core cellular functions, with circle size representing the degree of interconnection. Modules M1, M2, and M5 showed the clearest separation between OS samples and the accumulated normal sample cohort, while M3 and M4 showed relatively weaker discrimination.
Cluster analysis of the 84 OS samples using iMAS technology on Beer's data distinguished 5 module clusters (labeled 1 through 5). The 15 normal samples clustered together with consistent module scores (shown in red), while the OS samples (shown in blue) distributed across the 5 clusters based on their individual module aberrance patterns. This heterogeneity in module dysregulation across patients reflects the known clinical heterogeneity of osteosarcoma and supports the potential value of personalized treatment assignment based on cluster membership.
A Support Vector Machine (SVM) classifier (C-classification kernel) was used to evaluate whether the 61 significant modules could replace the full gene expression dataset for distinguishing OS patients from normal controls. Expression data were randomly split into training and test groups in a 6:4 ratio, and 5-fold cross-validation was used to estimate performance.
The SVM applied to the screened 61 modules achieved an accuracy of 100%, a specificity of 1.00, and a sensitivity of 1.00: every OS sample and every normal sample was correctly classified without error. This perfect classification indicates that the module-level representation retains the full discriminatory information present in the original 23,214-gene dataset while dramatically reducing feature dimensionality.
The high SVM accuracy demonstrates that modules can serve as surrogate features for individual gene panels in OS diagnosis. Because modules represent biologically coherent functional units rather than arbitrary gene combinations, they provide interpretable biological context for classification decisions. The SVM's ability to handle high-dimensional data with relatively small sample sizes (15 normal, 84 OS) is well-suited to the rare tumor setting characteristic of osteosarcoma research.
Traditional differentially expressed gene (DEG) analysis compares group averages and cannot capture patient-specific molecular aberrations. The iMAS method's key innovation is computing module scores for each individual patient relative to accumulated normal samples (ANS), rather than comparing group averages. This enables the identification of dysregulated modules that are specific to individual patients even when those aberrations are not shared across the entire OS cohort.
Ahn et al. previously demonstrated that the Average Z equation used in iMAS can reveal noticeable expression aberrations with clinical significance and distinguish survival-relevant pathways statistically in colorectal and lung cancer. The current study extends this framework to osteosarcoma, where tumor heterogeneity is high and individual-level pathway analysis could guide personalized therapy selection that standard group-level DEG studies would miss.
The module-centric approach also provides analytical advantages in rare tumors with limited sample sizes. Instead of attempting to power statistical tests across thousands of individual genes with only 84 OS samples, iMAS aggregates signal across gene sets whose co-regulation is supported by PPI network topology. This aggregation improves statistical stability and enables meaningful analysis even when individual gene-level tests lack power.
The study combines two public datasets (E-GEOD-33382 and E-GEOD-28974) to reach 84 OS samples and 15 normal controls. The control group of only 15 accumulated normal samples is small for a reference cohort, and may not capture the full range of normal bone tissue expression variability. The authors themselves note that the distribution outcomes of the 61 modules have not been verified in an independent abundance of samples and recommend that screened modules be tested in new studies.
Despite perfect SVM accuracy, the model was tested on the same dataset from which modules were derived. No truly independent external validation cohort was used to confirm that the 61 modules generalize to osteosarcoma patients from other institutions, sequenced on different platforms, or treated in different clinical contexts. This creates a risk of overfitting to the specific biological and technical characteristics of the two combined datasets.
The iMAS framework does not integrate clinical outcome data such as survival, metastasis, or chemotherapy response. While the 5 OS patient clusters identified by iMAS likely represent biologically distinct subgroups, whether those subgroups have different clinical outcomes or respond differently to treatment remains unknown. Linking iMAS module clusters to clinical endpoints would be the critical validation step for therapeutic decision-making applications.
The immediate next step is validating the 61 identified modules in an independent osteosarcoma cohort with clinical outcome annotation. If module clusters correlate with survival, metastatic potential, or chemotherapy response, they could serve as biomarkers for personalized treatment assignment: patients in high-risk clusters receiving intensified therapy or novel agents, while low-risk clusters receive standard care. The framework is scalable as large national OS genomic initiatives (e.g., NCI TARGET) provide more samples.
Integration with other data modalities could deepen the personalized medicine application. Adding DNA methylation, copy number variation, or somatic mutation data alongside gene expression would allow iMAS to capture multi-omic sources of module dysregulation, which is particularly important in osteosarcoma given its complex and highly rearranged genome. Methylation-driven silencing of module hub genes may be as clinically important as transcriptional dysregulation.
The SVM classification framework could be extended from binary OS vs. normal discrimination to multi-class prediction tasks such as identifying patients likely to develop metastasis, predicting chemotherapy response, or stratifying patients by molecular subtype. Given the rarity of osteosarcoma (4 to 5 cases per million per year), consortium-based collection of samples across institutions using the iMAS framework could accelerate the identification of clinically actionable personalized biomarkers.