Osteosarcoma is the most common primary malignant bone tumor, arising from mesenchymal tissue, and carries a distinctly bimodal age distribution: it peaks in adolescents and young adults, then appears again in older patients, often secondary to Paget's disease or prior radiation. Among patients aged 0-24, the incidence is approximately 4.4 per million. When treated with surgery alone, five-year survival hovers around 15-17%. The introduction of combination chemotherapy in the 1970s dramatically changed outcomes, enabling roughly 70% of patients with localized disease to be cured. However, for the 20-30% with metastatic or recurrent disease, five-year survival has remained stubbornly near 20% for the past three decades, reflecting a near-complete failure to identify effective second-line strategies.
The missing prognostic layer: For localized osteosarcoma, chemotherapy response assessed at surgery remains the best available predictor of prognosis. But no validated molecular biomarker exists to stratify risk at diagnosis, guide treatment intensity decisions, or identify patients who should be enrolled in experimental protocols before recurrence. This gap motivates the search for transcriptomic markers that can be read from routinely available RNA sequencing data.
Pseudogenes as an unexplored resource: Pseudogenes are genomic sequences homologous to known protein-coding genes but rendered non-protein-coding through accumulated mutations, such as frameshift mutations, insertions, deletions, and premature stop codons. For decades they were dismissed as "gene fossils" or "junk DNA." However, evidence has accumulated that pseudogenes are actively transcribed and regulate their parent genes and broader transcriptomes through several mechanisms: they generate small interfering RNAs (siRNAs), act as competitive endogenous RNAs (ceRNAs) that sequester shared microRNAs, produce antisense transcripts, and directly modulate protein-coding gene expression. In hepatocellular carcinoma, high expression of pseudogene RP11-564D11.3 correlates with poor prognosis. In high-grade ovarian cancer, SLC6A10P expression predicts recurrence. Osteosarcoma had not been systematically examined.
This study by Liu et al. (published in Genes, 2019) represents the first systematic analysis of pseudogene expression as a prognostic resource in osteosarcoma, using machine learning to distill a 1,333-pseudogene transcriptome down to a four-member signature with clinical-grade predictive performance.
All data were drawn from the Therapeutically Applicable Research to Generate Effective Treatments (TARGET) database, a publicly available NCI-funded resource that pairs RNA sequencing data with clinical follow-up information for pediatric cancers. The full TARGET osteosarcoma dataset contains clinical records for 274 patients and matched RNA-seq expression profiles for 101 patients. The significant gap between the clinical and molecular datasets reflects the reality that RNA sequencing was not available for all archival specimens.
Cohort construction: Samples with incomplete clinical follow-up were excluded. Pseudogenes with low average expression, defined as average transcripts per million (TPM) at or below 1.0 across all samples, were also removed, since low-expressed features contribute primarily noise to survival models. After both filters, the working cohort included 94 osteosarcoma patients with 1,333 pseudogene expression profiles. Expression values were log2-transformed as log2(TPM + 1) to stabilize variance before statistical modeling.
Clinical characteristics: The 94-patient cohort included 54 males and 40 females. Median age at diagnosis was 14.5 years (range 3-39 years), consistent with the adolescent peak of osteosarcoma incidence. The primary tumor site was the leg or foot in 83 of 94 patients (88.3%), with the remaining 11 at other sites. Metastatic disease at enrollment was present in 22 patients (23.4%), and 72 were non-metastatic. Among the 94 patients, 57 were alive and 37 had died at the time of data extraction, yielding an event rate sufficient for survival modeling despite the modest cohort size.
The reliance on a single publicly available dataset is both a practical necessity and a recognized limitation. The TARGET osteosarcoma cohort is one of the largest publicly available datasets with paired RNA-seq and survival data for this disease. Using TPM normalization, which accounts for both gene length and sequencing depth, ensures that expression comparisons across samples are valid. The log2 transformation brings the distribution of expression values closer to normality, satisfying assumptions of the downstream Cox regression framework.
The prognostic signature was constructed using a sequential, three-stage variable selection pipeline designed to move from a large candidate feature space down to a compact, clinically deployable model. This stepwise approach is well-established in biomarker discovery because it prevents the statistical artifacts that arise when thousands of features are directly entered into a multivariate survival model.
Stage 1 - Univariate Cox regression: Each of the 1,333 pseudogenes was individually tested for association with overall survival using univariate Cox proportional hazards regression. Pseudogenes achieving p below 0.05 after false discovery rate (FDR) adjustment were retained. This yielded 125 survival-related pseudogenes, an initial reduction to about 9% of the original feature space. Among these, the majority (91 of 125) were identified as risk factors with hazard ratios above 1, suggesting that elevated expression is associated with worse outcomes, while the remaining 34 were protective.
Stage 2 - LASSO regression: The 125 candidate pseudogenes were then subjected to least absolute shrinkage and selection operator (LASSO) Cox regression. LASSO is a penalized regression method that adds an L1 penalty term to the likelihood function, which forces the coefficients of weakly predictive features to shrink exactly to zero. This achieves sparse variable selection, reducing the risk of overfitting that is inherent when sample sizes (n = 94) are small relative to the number of candidate predictors (n = 125). Ten-fold cross-validation was used to select the optimal regularization parameter lambda (lambda.min = 0.13), which yielded 15 pseudogenes with non-zero coefficients.
Stage 3 - Multivariate Cox regression: The 15 LASSO-selected pseudogenes were entered into a multivariate Cox proportional hazards regression. Only pseudogenes with p below 0.05 in this final model were retained for the signature. This reduced the set to four pseudogenes: RPL11-551L14.1, RPL7AP28, RP4-706A16.3, and RP11-326A19.5. A continuous risk score was computed for each patient as a linear combination of the four pseudogene expression values weighted by their multivariate Cox regression coefficients. All analyses were conducted in R version 3.5.1.
The final prognostic signature integrates expression levels of four pseudogenes into a single continuous risk score. Each pseudogene contributes to the score through a beta coefficient derived from multivariate Cox regression, representing its independent association with overall survival after accounting for the effects of the other three members.
Individual pseudogene hazard ratios: In the multivariate model, the four pseudogenes had the following hazard ratios for overall survival: RPL11-551L14.1, HR 0.65 (95% CI: 0.44-0.95), a protective factor where high expression associates with lower death risk; RPL7AP28, HR 0.32 (95% CI: 0.14-0.76), the strongest protective factor in the model; RP4-706A16.3, HR 1.89 (95% CI: 1.35-2.65), the only risk factor where high expression associates with worse survival; and RP11-326A19.5, HR 0.52 (95% CI: 0.37-0.74), another protective factor. Three of the four pseudogenes are independently protective, while one (RP4-706A16.3) is a risk factor.
Risk score formula: The continuous risk score is computed as: Risk score = (expression of RPL11-551L14.1 x -0.4327) + (expression of RPL7AP28 x -1.1344) + (expression of RP4-706A16.3 x 0.6360) + (expression of RP11-326A19.5 x -0.6503). The sign of each coefficient matches the directionality of the hazard ratio: negative coefficients for the three protective pseudogenes (high expression lowers the score toward better prognosis) and a positive coefficient for the risk factor RP4-706A16.3 (high expression raises the score toward worse prognosis).
Risk stratification: Patients were divided into high-risk (n = 47) and low-risk (n = 47) groups using the median risk score of 1.23 as the cut-off. Visual inspection of the risk score distribution confirmed that as risk score increases, survival time decreases and death events accumulate, and expression of the three protective pseudogenes decreases while expression of RP4-706A16.3 increases. None of the four pseudogenes had been reported in prior osteosarcoma studies, making them novel candidate biomarkers warranting experimental follow-up.
The core validation of the four-pseudogene signature used two complementary approaches: Kaplan-Meier survival analysis to assess whether the risk groups have meaningfully different survival trajectories, and receiver operating characteristic (ROC) curve analysis to quantify the signature's accuracy at specific clinical time points (3-, 5-, and 8-year survival).
Kaplan-Meier analysis: Patients in the high-risk group had significantly worse overall survival compared to the low-risk group (p below 0.0001, two-sided log-rank test). The separation of the two curves was clinically meaningful and statistically robust, supporting the signature's ability to distinguish meaningfully different prognostic groups from a single blood or tissue RNA-seq sample taken at diagnosis.
ROC curve analysis: Time-dependent ROC analysis produced AUC values of 0.885 for 3-year survival prediction, 0.878 for 5-year survival prediction, and 0.796 for 8-year survival prediction. The higher AUC at earlier time points is expected: the model was trained on overall survival follow-up data, so its predictive accuracy is strongest in the window where most events occur, and diminishes slightly at the longest time horizon where fewer data points inform the estimate. An AUC of 0.878 at 5 years places this signature among the better-performing molecular prognostic tools reported for osteosarcoma, where prior gene-based and lncRNA-based models have typically reported 5-year AUCs in the 0.70-0.85 range.
Comparison with a gene-based comparator: The authors also constructed an analogous risk score using protein-coding genes (PCGs) through the same pipeline and compared it to the pseudogene signature by ROC curve. The AUC values were close, with the pseudogene signature slightly outperforming, though the authors note that the small sample size limits the interpretability of this comparison. This finding suggests that pseudogene expression carries at least equivalent prognostic information to conventional gene expression for osteosarcoma survival, which is notable given that pseudogenes are rarely included in clinical genomic panels.
For a prognostic biomarker to have clinical value, it must demonstrate that its association with survival is not simply a proxy for other known prognostic factors. If the signature only performed well in, say, female patients or non-metastatic patients, its utility would be limited. The authors therefore performed pre-specified subgroup analyses stratified by three clinical variables: gender, age at diagnosis, and metastatic status.
Gender subgroups: In the male subgroup (n = 54), high-risk patients had significantly shorter overall survival than low-risk patients (p below 0.01), with a 5-year AUC of 0.753. In the female subgroup (n = 40), the same direction of effect was seen with a similar p-value and a 5-year AUC of 0.805. The distribution of risk scores did not differ significantly between males and females (p above 0.05), indicating that the signature is not confounded by sex.
Age subgroups: Patients were divided at age 18, yielding 72 patients under 18 and 22 patients 18 or older. In both groups, high-risk patients had significantly shorter overall survival. Five-year AUC values were 0.888 for the under-18 group and 0.861 for the 18-and-older group, both higher than the overall cohort AUC of 0.878. This reflects the fact that the under-18 group is larger and contributes more events to estimate the model. Risk score distributions did not differ between age groups (p above 0.05).
Metastatic status subgroups: In the non-metastatic subgroup (n = 72), the signature separated high- and low-risk groups with a shorter survival in the high-risk group. In the metastatic subgroup (n = 22), the same directionality was maintained, though the smaller sample size limits the power of this comparison. The 5-year AUC in the non-metastatic group was 0.857, and in the metastatic group it was 0.750. These findings confirm that the signature adds prognostic information beyond what is captured by knowing whether the patient has metastatic disease at diagnosis, which itself is a strong adverse prognostic factor. Stage information was not available for all patients and was therefore not included in subgroup analyses.
Beyond demonstrating statistical performance, the authors developed two components that move the signature closer to clinical use: a nomogram for individualized survival prediction and a co-expression analysis to identify the biological processes through which the four pseudogenes may influence osteosarcoma biology.
Nomogram and calibration: The nomogram integrates the pseudogene risk score with available clinical variables to generate individualized probability estimates for 3-, 5-, and 8-year survival. This tool converts the continuous risk score into a point-based system that a clinician can use directly, without understanding the underlying mathematical model. Calibration plots, which compare the nomogram's predicted survival probabilities with actual observed outcomes, confirmed that the nomogram's predictions align well with reality across the range of predicted risks. The R package rms was used for nomogram construction and calibration.
Co-expression analysis and functional annotation: Pearson correlation coefficients were calculated between each of the four pseudogenes and all protein-coding genes (PCGs) in the dataset. PCGs significantly correlated (either positively or negatively) with each pseudogene were collected and submitted to the Database for Annotation, Visualization, and Integrated Discovery (DAVID) for gene ontology (GO) biological process enrichment and KEGG pathway analysis (p below 0.05 threshold).
The co-expression results revealed that the four pseudogenes are associated with three broad biological themes: regulation of malignant phenotype (processes involved in tumor growth, invasion, and cell cycle control), immune regulation (immune effector processes, T cell and NK cell activity pathways), and DNA/RNA editing (RNA processing, splicing, and modification). These functional annotations are biologically coherent: osteosarcoma progression is driven by dysregulated proliferation and evasion of immune surveillance, and pseudogene-mediated regulation of miRNA activity and ceRNA networks could plausibly affect all three categories. The association with immune regulation is particularly noteworthy given the growing interest in immunotherapy for osteosarcoma, where responses to checkpoint blockade have been modest but mechanistically informative.
Single dataset, no external validation: The most significant limitation is that the entire analysis, including model construction, performance evaluation, and subgroup validation, was conducted on the same 94-patient cohort from the TARGET database. There is no independent external test set from a separate institution. Internal cross-validation, used during the LASSO step, helps prevent extreme overfitting at the feature selection stage but does not substitute for external validation. All reported AUC values should be interpreted as optimistic estimates of real-world performance. A reduction of 5-10 percentage points in AUC is plausible when the model is applied to a truly independent dataset with different patient demographics, RNA-seq protocols, or clinical workflows.
Pure data mining without experimental validation: No wet-lab experiments were performed to confirm the expression findings or to establish the mechanisms by which the four pseudogenes influence osteosarcoma biology. The hazard ratios and AUC values are derived entirely from bioinformatic analysis of existing public data. Experimental confirmation, including in vitro knockdown or overexpression of the four pseudogenes in osteosarcoma cell lines and in vivo validation in animal models, would be necessary before any mechanistic claims can be made.
Missing stage information: Disease stage is one of the most clinically important prognostic variables in osteosarcoma. The TARGET dataset did not have complete stage information for all 94 patients, preventing inclusion of stage in subgroup validation or nomogram construction. This is a meaningful gap: the four-pseudogene signature's independence from stage cannot be formally demonstrated, which limits its positioning relative to stage-based risk assessment.
Future directions: The authors identify three priorities: first, prospective validation of the four-pseudogene signature in independent, ideally multi-center osteosarcoma cohorts with complete clinical annotation; second, experimental characterization of the regulatory mechanisms through which RPL11-551L14.1, RPL7AP28, RP4-706A16.3, and RP11-326A19.5 influence tumor biology, focusing on their ceRNA and siRNA activity; and third, integration of the pseudogene signature with other molecular data layers (somatic mutations, copy number alterations, methylation) to build multimodal prognostic models. The strong association of co-expressed PCGs with immune regulation also suggests a role for these pseudogenes in predicting immunotherapy responsiveness, an area that warrants prospective investigation given osteosarcoma's clinical need for new treatment options.