Osteosarcoma is the most common primary bone malignancy, with a bimodal age distribution: one peak in adolescence and a second in older adulthood. In the United States, bone cancers account for 5% of cancer mortality in children under 15 and rise to 15% in those aged 15-19, contributing to more than 1,900 cancer-related deaths annually in children and adolescents. Early, accurate diagnosis is critical because it directly affects survival and the likelihood of limb-sparing surgery versus amputation.
The diagnostic challenge: Plain x-ray is the primary frontline tool for assessing bone tumors, and osteosarcoma presents with recognizable features including a sun-burst pattern, periosteal elevation forming Codman's triangle, and new bone formation in soft tissue. However, these features overlap with other aggressive but benign lesions, and distinguishing malignant from benign tumors by imaging alone requires significant radiological expertise. An incorrect call can delay curative therapy or expose a patient to unnecessary surgery.
The multimodal opportunity: Combining imaging data with molecular data offers a path toward more objective and reproducible classification. This paper proposes fusing x-ray image features with RNA sequencing (RNA-seq) transcriptomic data in a machine learning model. The hypothesis is that the radiological presentation of a tumor and its underlying gene expression profile together carry more discriminative information than either data source alone, and that a random forest (RF) classifier trained on their combined features can reliably separate osteosarcoma from benign bone lesions.
Prior work by Shen et al. demonstrated the viability of this multimodal concept by combining x-ray features with metabolomic data, achieving an AUC of 0.94. This paper extends that approach by substituting RNA-seq for metabolomics, a more clinically scalable molecular assay available from routine tumor biopsy specimens.
The dataset was assembled in collaboration with Shanghai Jiao Tong University (SJTU) and collected under a Shanghai General Hospital IRB-approved protocol with written informed consent from all participants. Every sample in the dataset contains both a plain x-ray image in DICOM format and a matched RNA-seq library, making this a fully paired multimodal cohort. This pairing is essential: the classifier relies on combining features from both modalities for each individual patient, so missing either data type for any patient would exclude that case entirely.
Cohort composition: The dataset comprises 17 patients: 6 with benign bone tumors and 11 with osteosarcoma. The benign group had a mean age of 40.33 years (standard deviation 12.48) and consisted of 5 males and 1 female. The osteosarcoma group was substantially younger, with a mean age of 22.73 years (standard deviation 16.14), reflecting the known adolescent peak, and comprised 4 males and 7 females. The class imbalance (roughly 2:1 osteosarcoma-to-benign ratio) is a known limitation that the authors explicitly acknowledge.
RNA-seq library preparation: RNA sequencing libraries were generated in 200 bp paired-end mode using the Illumina HiSeq platform. The resulting FASTQ files were subjected to quality assessment before alignment. All 17 samples passed quality control; no adapter trimming was applied and no samples were excluded based on quality metrics. This indicates the raw sequencing data were of sufficient quality to proceed directly to alignment without preprocessing adjustments.
The deliberate decision to include only samples with both imaging and molecular data, rather than supplementing with imaging-only or RNA-seq-only cases, reflects a commitment to a true multimodal design. However, it also keeps the dataset very small, which has direct consequences for the statistical reliability of the reported performance metrics and the generalizability of the findings.
Before any features could be extracted from the x-ray images, the tumor and surrounding bone had to be delineated from the background. The authors used a semi-automated segmentation approach called Graph Cut with Lazy Snapping. In this method, the user places foreground seeds (marking the tumor region) and background seeds (marking the surrounding non-tumor tissue). These user-defined seeds serve as hard constraints for a graph-based energy minimization algorithm, which then assigns foreground or background labels to every pixel in the image to minimize the Gibbs Energy E(X).
Lazy Snapping and superpixel preprocessing: Before Graph Cut labeling, the watershed algorithm pre-segments the image into superpixels, producing a compressed graph where each node represents a small uniform region rather than an individual pixel. This dramatically reduces the computational cost of the energy minimization while preserving the boundaries relevant for tumor delineation. The energy function E(X) has two components: a likelihood energy E1(x) penalizing unlikely label assignments for each node, and a prior energy E2(xi, xj) penalizing discontinuous label transitions between neighboring nodes.
Feature extraction from segmented regions: Once segmentation was complete, two categories of features were extracted: features from the tumor segmentation alone (morphological and texture features) and features computed using both the tumor and bone segmentations (tumor border clarity and distance to joint). Morphological features included mean pixel intensity, standard deviation, entropy, kurtosis, skewness, convex area, eccentricity, perimeter, major axis length, minor axis length, and solidity. Texture features were computed from the gray-level co-occurrence matrix (GLCM), including contrast, correlation, energy, and homogeneity, as well as wavelet decomposition coefficients from approximation, horizontal, vertical, and diagonal detail sub-bands.
The combination of GLCM texture descriptors and wavelet-based statistics captures both local pixel-level heterogeneity and multi-scale spatial patterns in the tumor appearance, reflecting known radiological correlates of osteosarcoma aggressiveness such as cortical destruction and periosteal reaction.
The RNA-seq processing pipeline followed a four-step workflow: quality assessment, genome alignment and count generation, normalization, and feature reduction. Raw FASTQ files were first evaluated using FastQC version 0.11.8 with default parameters, which produces per-base sequence quality scores, GC content distributions, and over-representation metrics. All samples passed with acceptable overall quality; no adapter trimming was needed.
Alignment and count generation: Reads were aligned to the human reference genome (hg38) using STAR with a genome index built from the Gencode V29 GTF annotation file. The maximum number of multiple alignments permitted per read was set to 20. All 17 samples achieved a ratio of uniquely mapped reads greater than 76%, indicating high alignment quality. The aligned BAM files were coordinate-sorted and processed with the QoRTs Java utility to generate transcript-level count files, which were then read into R using the QoRTs library.
Normalization with DESeq2: Raw counts were normalized using DESeq2's median ratio method. For each sample, a size factor is calculated as the geometric mean of per-transcript count ratios across all samples, and raw counts are divided by that size factor. A critical data-leakage prevention measure was applied: for cross-validation folds, the test set size factors were derived from the training set size factors rather than computed from the full dataset. This prevents the test samples from influencing their own normalization, which would artificially inflate performance metrics.
Literature-guided feature reduction: The raw transcript count matrix contains thousands of features, far exceeding the 17-sample dataset size and creating severe overfitting risk. To reduce this, a prior literature search identified 211 genes whose upregulation or downregulation is associated with osteosarcoma tumorigenesis and metastasis. Only transcript counts mapping to these 211 genes were retained as RNA-seq features, filtering out the vast majority of the transcriptome and limiting the molecular feature space to biologically plausible candidates.
After image feature extraction and RNA-seq feature reduction, the two modalities were merged into a single feature vector per patient. This combined representation means each of the 17 observations contains both the x-ray-derived morphological and texture features and the counts from transcripts associated with the 211 osteosarcoma-relevant genes. The combined feature space was then fed into a random forest (RF) classifier, an ensemble method that builds multiple decision trees on bootstrapped subsets of the training data and aggregates their predictions.
Random forest hyperparameters: The RF was configured with 50 trees. At each decision split, the number of features considered was the square root of the total feature count, the standard approach for classification tasks in random forests. The minimum leaf size was set to 1, allowing trees to grow until each leaf contains a single observation. These are largely default hyperparameter settings rather than tuned values, reflecting the practical constraint that the tiny dataset size makes robust hyperparameter optimization infeasible.
Dimensionality reduction with PCA: Before training each RF fold, principal component analysis (PCA) was applied to the combined feature matrix. PCA compresses the high-dimensional feature space into a smaller set of orthogonal principal components (PCs) that capture most of the variance in the data, directly addressing the risk of overfitting when features vastly outnumber samples. To prevent data leakage, PCA was fit only on training samples, and the resulting transformation coefficients were applied to the test samples. Test sets were also scaled by training-set statistics before PCA projection. The number of PCs used as model input was varied to assess performance sensitivity.
Two cross-validation strategies: Performance was evaluated under both three-fold cross-validation and leave-one-out cross-validation (LOOCV). For three-fold CV, observations were divided into three folds containing 2 benign and 3-4 osteosarcoma cases each, with no overlap. In LOOCV, all observations except one form the training set and the held-out observation is tested; this is repeated for every individual. LOOCV is generally considered more appropriate than k-fold CV for very small datasets because it maximizes the training set size at each iteration, though it can produce optimistic estimates of generalization performance.
Under three-fold cross-validation with one principal component, the RF model achieved an AUC of 0.7272, an F1 score of 0.8462, sensitivity of 1.0 (all osteosarcoma cases correctly identified), specificity of 0.3333, and overall accuracy of 0.9647. The high sensitivity with low specificity reflects the class imbalance: with only 6 benign cases, incorrectly classifying 4 of them as osteosarcoma still yields acceptable overall accuracy because osteosarcoma dominates the dataset. Performance did not significantly improve when additional PCs were added beyond one under this cross-validation scheme, suggesting that the variance captured by the first PC is the most informative, with additional PCs likely adding noise.
Leave-one-out cross-validation results: Under LOOCV with 2 PCs, the model achieved an AUC of 0.9015, representing a substantial improvement over the three-fold result. The increase in performance with LOOCV is expected because LOOCV trains on 16 of 17 samples at each iteration, providing the model with almost the entire dataset, whereas three-fold CV trains on only two-thirds. The authors note that performance improved monotonically as the number of PCs increased from 1 to 2 but did not improve meaningfully beyond 2, suggesting that variance beyond the second PC reflects noise rather than discriminative signal.
Contextualizing the numbers: The AUC of 0.9015 under LOOCV is a promising result, but it must be interpreted cautiously. LOOCV estimates can be optimistic for very small datasets, and with only 17 samples, a single misclassification in the benign group shifts specificity by 16.7 percentage points. The authors are transparent about this: they note that a single incorrect call among the 6 benign cases produces large swings in specificity. The high sensitivity (1.0 in three-fold CV) is clinically meaningful because missing an osteosarcoma diagnosis is the most dangerous error in this setting, but the low specificity means many benign cases would be flagged incorrectly in a real screening scenario.
Compared to the earlier Shen et al. work using x-ray plus metabolomics (AUC 0.94), the current RNA-seq approach performs somewhat lower under three-fold CV but comparably under LOOCV. Given the different dataset sizes and cross-validation strategies involved, direct numeric comparison between the two studies is difficult, but both support the general principle that multimodal fusion outperforms imaging alone in this classification task.
Sample size and class imbalance: The most significant limitation of this study is the dataset of 17 patients with a 2:1 class imbalance (11 osteosarcoma, 6 benign). At this scale, performance metrics are highly sensitive to individual misclassifications. Specificity, in particular, is unreliable when calculated over only 6 benign cases: a single wrong call changes specificity from 1.0 to 0.833, and two wrong calls drop it to 0.667. The authors acknowledge this directly, noting that "a difference of 1 incorrect classification could result in large changes in specificity." This limitation means the reported AUC values, while promising, cannot be taken as reliable estimates of the classifier's generalization performance on a broader population.
Dimensionality collapse via PCA: The use of PCA to collapse thousands of RNA-seq transcript features down to 1-2 principal components is a pragmatic solution to the curse of dimensionality in a small dataset, but it comes at a cost. PCA is an unsupervised technique that maximizes explained variance without any regard for class separability. The most informative PCs for distinguishing osteosarcoma from benign tumors may not correspond to the PCs that explain the most variance in the overall feature matrix. Supervised dimensionality reduction techniques (linear discriminant analysis, UMAP with supervision, or LASSO-penalized logistic regression for feature selection) might better identify the transcript features most relevant to classification.
Semi-automated segmentation dependency: The image processing pipeline requires a human operator to place foreground and background seeds for the Graph Cut segmentation. This introduces variability across operators and centers, and makes the pipeline non-trivially deployable in settings without radiological expertise. Fully automated segmentation using convolutional neural networks (U-Net architectures) trained on annotated bone tumor images would be a necessary step toward a deployable clinical tool.
Path to clinical translation: Despite these limitations, the results suggest that the multimodal framework combining plain x-ray features with RNA-seq transcriptomics can contribute to computer-aided diagnosis of osteosarcoma. The natural next steps are external validation on a larger independent cohort, exploration of alternative molecular assays (such as targeted gene panels that are cheaper and faster than whole-transcriptome RNA-seq), integration of deep learning-based image features rather than hand-crafted morphological and texture descriptors, and replacement of semi-automated segmentation with fully automated approaches. The collaboration between SJTU and University of Michigan groups provides a foundation for such larger prospective studies.