Our overall study design is illustrated in Figure 1 .

We collected a training cohort coupled with two validation cohorts. All datasets are available from public data sources.

The training cohort comprised 98 clinically annotated patient cases produced by the Therapeutically Applicable Research to Generate Effective Treatments (TARGET) Osteosarcoma project (https://ocg.cancer.gov/programs/target/projects/osteosarcoma). Transcripts per million (TPM) values and read counts were downloaded from the TARGET data matrix. Furthermore, TPM values were scaled with Min-Max normalization.

The dataset GSE21257 containing clinical information of 53 osteosarcoma patients (25), was obtained from the Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo/). It was treated as a validation cohort to test the molecular subtyping and drug targets. The other validation dataset provided by the European Genome-phenome Archive (EGA, https://ega-archive.org) contains RNA-seq data from 48 specimens of osteosarcoma patients labeled with primary, relapsed, or metastatic (17). The accession number is EGAS00001003887.

The gene expression matrix of the 98 osteosarcoma samples was used to identify the molecular subtypes using the non-negative matrix factorization (NMF) consensus cluster method (26). NMF is a practical machine-learning approach for identifying molecular classification and has been widely used in many subtyping cases (27). Firstly, we filtered the genes expressed (read counts≥10) in less than 24 samples, and 26824 genes remained. Secondly, after min-max normalization, all TPM data were rescaled from 0 to 1. Thirdly, we performed principal component analysis (PCA) to select the top 1600 genes based on the sum of PC1 and PC2 rotation values. These 1600 genes and expression profiles were further subjected to NMF v.0.22.0 in R for unsupervised consensus-clustering. Non-smooth NMF (nsNMF) algorithm was employed with 100 iterations for the rank survey between 2 and 6 clusters.

The R package Survminer plotted the Kaplan-Meier survival curves. This tool employed a log-rank test to compare subtypes’ overall and event-free survival rates.

Using the single sample gene set enrichment analysis (ssGSEA) method, xCell and ESTIMATE calculated the tumor-infiltrating immune cell content of cancer samples. xCell could compute the abundance of 64 different cell types and scores of immune, stroma, and microenvironment (28). ESTIMATE could take advantage of the unique properties of the transcriptional profiles of cancer samples to infer tumor cellularity as well as the different infiltrating normal cells (29). They were two efficient publicly available algorithmic tools.

The differentially expressed genes (DEGs) were identified as highly expressed in one subtype compared to all three other subtypes. The limma package in R performed this computation. Adjusted P-values were applied for multiple testing corrections through the default Benjamini-Hochberg false discovery rate (FDR) method (30). Adjusted P-value < 0.05 and |fold change (FC)| > 1.5 were set as the cutoff values to determine DEGs.

Gene set enrichment analysis (GSEA) was performed on all genes. GSEA was run with Molecular Signatures Database (MSigDB) set V7.4, C2: curated gene sets (browse 6290 gene sets) (31). In addition, enrichment analyses of gene ontology (GO), Reactome, and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway, were conducted on the DEG by using the R package clusterProfiler (version 4.2.2) (32).

The PPI (protein-protein interactions) network was constructed via an open-source platform Cytoscape (33). Using the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database (https://string-db.org), the PPI network illustrated how proteins interrelate functionally and physically with each other by encoding the gene list as input (34).

Druggable genes for osteosarcoma were mined using the Drug–Gene Interaction Database (DGIdb 4.0, https://www.dgidb.org). This database contained drug targets, indications, MOA, and drug status from over thirty reputable sources, including DrugBank, FDA, and NCI (35). The hazard ratios (HR) were computed with log-rank tests using univariate analysis based on a cox proportional hazards regression model. In addition, we used the Survminer R package to calculate the log-rank p values and the optimal cut-off values to stratify patients into high/low gene expression groups and further visualize the survival differences between groups. High-risk genes were identified with HR>1 and log-rank p<0.05.

The U2OS and Saos-2, human osteosarcoma cell lines, were obtained from the American Type Culture Collection (ATCC, VA, USA). The U2OS was cultured in high glucose DMEM supplemented with 10% fetal bovine serum (FBS), and the Saos-2 was cultured in McCoy’s 5A supplemented with 15% FBS. They were incubated at 37°C in 5% CO_2. 1% penicillin and streptomycin (Gibco, USA) were added to the base media as supplements.

Knockdown of SQLE in cells was achieved by using Turbofect transfection reagent (Thermo Fisher Scientific, OR, USA). 2 μL Turbofect reagent was mixed with 500 μL Opti-MEM (Thermo Fisher Scientific) and combined with 1 μg siRNA (GenePharma, Shanghai, China). The mixture was incubated at room temperature for 20 minutes and added dropwise to the cells. The targeted oligos were as follows (5’-3’):

siNegative control: UUCUCCGAACGUGUCACGUTT.

Total RNA was extracted by Trizol reagents (Thermo Fisher Scientific). After reverse transcription by HiScript III All-in-one RT SuperMix (Vazyme, Nanjing, China), the mRNA levels were measured by SYBR master mix dye (Vazyme). The relative gene expression level was analyzed using the 2-ΔΔ Ct method, and GAPDH was used as the internal reference gene. The Q-PCR primer sequences were as follows (5’-3’): SQLE, forward, GGCATTGCCACTTTCACCTAT; reverse, GGCCTGAGAGAATATCCGAGAAG. GAPDH, forward, TGCACCACCAACTGCTTAGC; reverse, GGCATGGACTGTGGTCATGAG. QPCR was conducted after 24h siRNA transfection.

Total protein was extracted by RIPA lysis buffer (Thermo Fisher Scientific). After BCA quantification, equal volumes of protein samples were separated by 10% SDS-PAGE and transferred to a 0.45 nitrocellulose filter membrane (Millipore, MA, USA). The immunoblots were sequentially probed with the SQLE primary antibody (Cat No. 12544-1-AP, Proteintech, Wuhan, China) and secondary antibody (Cat No. SA00001-2, Wuhan, China). Finally, the detection was performed using an ECL chemical luminescent detection kit (Thermo Fisher Scientific), and the bands were further analyzed using ImageJ software. The levels of the target protein were normalized to β-actin expression. Western blot analyses were conducted after 48h siRNA transfection.

Terbinafine was purchased from Selleckchem (Houston, TX, USA). Cells were plated in 96-well plates at a density of 3,000 cells; the indicated siRNAs or terbinafine were introduced to cells after 12 h. When it came to the endpoint of the assay, 10% CCK8 reagent (Dojindo, Kyushu Island, Japan) was added into the wells; after 2 h incubation at 37°C, the OD values were detected at 450 nm.

Cells treated by siRNAs or terbinafine were plated in 6-well plates at a density of 2,000 cells and then cultured for two weeks. The cells were fixed and stained with 0.1% crystal violet (Beyotime, Beijing, China).

15,000 cells were plated in the upper chamber for the silencing cells, while DEME containing 10% FBS was added to the lower chambers. For the inhibitor-treated cells, cells in the upper chamber were treated with terbinafine (25 μM, 50 μM) for 24h in advance, while DEME containing 10% FBS was added to the lower chambers. After 24h incubation, the crossed cells were fixed and stained with 0.1% crystal violet.

Treated cells were washed with PBS twice and resuspended in isopropanol containing 1% Triton X-100 for 30 min at room temperature. After centrifugation at 12,000 rpm for 15 min, the supernatants were dried under nitrogen. The powders were dissolved with assay buffers. Cholesterol measurements were performed according to the instruction guide of the Amplex™ Red cholesterol assay kit (Thermo Fisher Scientific).

The subtype diagnostic model was constructed using a supervised learning method SVM. Permutation-based feature importance test (PermFIT) implemented in the R package “deepTL” (https://github.com/SkadiEye/deepTL) was proposed for assisting the interpretation of individual features in complex frameworks including SVM (36). Herein, we conducted PermFIT-SVM on expression profiles of DEG with a threshold of p-value<0.05 and thereby obtained 16, 25, 32, and 69 feature genes for S-I, S-II, S-III, and S-IV subtypes, respectively.

ProMS (https://github.com/bzhanglab/proms) is a unified and effective computational framework containing an SVM classifier for feature selection with the help of an omics view (e.g., RNA-seq) (37). We applied ProMS-SVM on feature genes of each subtype for training and testing the subtype diagnostic model. Preliminary runs were repeated 20 times, from K=1 to K=5, with the percentiles of fifth, 10th, 15th, 20th, and 25th. Then, runs were repeated 100 times to confirm the model’s accuracy. Hyperparameters were tuned using grid search with 3-fold cross-validation to select the best model via measuring the area under the receiver operating curve (AUC) value.

Utilizing patient information and expression data of 13 genes from SDM, the model’s genes were rigorously screened out using the R package “glmnet” (38) to repeat 5-fold cross-validation 100 times. The optima λ (λ.min) was selected to minimize the penalization in each iteration of the LASSO regression analysis. The number of model genes was determined by the most frequent gene numbers when setting λ.min and coefficient ≠ 0. Then, genes with the most frequencies were retained and combined as significant predictors of mortality. Finally, another 100 iterations were conducted to determine the coefficient of each signature gene.

The NMF algorithm, an effective consensus-clustering method, was employed for osteosarcoma subtype identification. The cophenetic score and silhouette width of the rank survey profiles indicated a potential option of four-subtype ( Supplementary Figures 1A, B ). Furthermore, we conducted 400 iterations for the clustering runs in the case of rank=4. The consensus clustering heat maps ( Supplementary Figure 2 ) illustrated that four-subtype is an appropriate solution for osteosarcoma ( Figure 2A ; Table S1 ). Hereinafter referred to as S-I, S-II, S-III, and S-IV. The TARGET database provided clinical information and genomic alteration ( Table S2 ). According to immune scores computed by xCell ( Table S3 ) and expression level of immune-related genes ( Supplementary Figures 3A, B ), S-II had the highest tumor-infiltrating immune cell abundance. In addition, overall survival and disease-free survival varied among subtypes ( Figure 2B ). Patients of S-IV showed the most unfavorable outcome.

To obtain the biological signatures of each subtype, we performed functional enrichment analyses on all genes and DEG ( Tables S4 , S5 ). GSEA conducted on all genes ( Figure 2C ) demonstrated that patients of S-I tended to exhibit more osteoclastic bone resorption (e.g., ACP5, SIGLEC15, CTSK), which was consistent with the clinical manifestations of osteosarcoma at the early stage (39). Multiple immune pathways were enriched in S-II (e.g., IGKC, S100A9, CD14), while S-III was prone to manifest salient cell proliferation features (e.g., TUBB2B, MYC). Patients of S-IV showed markedly enhanced lipid metabolism (e.g., SQLE, CD36, LPL), including cholesterol biosynthesis, which was believed to be a critical driven factor in the development of some solid tumors (40, 41). Apart from GSEA analysis, functional enrichment of DEG also revealed that S-I mainly correlated with bone resorption and osteoclast differentiation. Remarkable enrichment in immune-related functions was found in S-II. Various processes involved primarily in cancer cell proliferation were exhibited in S-III. Several functional categories, including the metabolism of lipids, were enriched in S-IV ( Figure 2D ).

Cholesterol metabolism produces metabolites with various biological functions and plays complex roles in tumorigenesis (42). Clinical trials manipulating cholesterol metabolism to inhibit tumorigenesis and reinvigorate anti-tumor immunity are ongoing (43). Overall, cholesterol homeostasis is maintained by the balance of de novo biosynthesis, cholesterol uptake, bile acid metabolism, esterification, and efflux (44). To decipher whether cholesterol metabolism was reprogrammed in patients of different osteosarcoma subtypes, we identified expression features of critical regulators that contribute to various modules ( Figure 3A ; Table S6 ). Cholesterol accumulation was facilitated by cholesterol biosynthesis and cellular uptake of LDL-cholesterol, while cholesterol consumption was mainly regulated by cholesterol esterification and efflux (45). When cholesterol accumulates to a high level, cells relieve cellular stress through esterification and efflux (46). In S-IV patients, mediators involved in cholesterol biosynthesis (HMGCS, HMGCR, SQLE) and cholesterol consumption (ACAT1, ACAT2, ABCA1, ABCG1) were highly expressed, while regulators promoting cellular uptake of LDL-cholesterol (LDLR, Niemann-Pick C proteins including NPC1 and NPC2) were not active. Especially, NPC2 was significantly depleted in S-IV. Besides, massive intracellular cholesterol limited the translocation of the SREBP2 to the nucleus, reducing its transcriptional target LDLR (47). The LDLR depletion suppressed the cholesterol uptake pathway and its downstream regulators NPC1/NPC2. These findings suggested cholesterol biosynthesis was particularly active in S-IV, leading to cholesterol accumulation. As a result of negative feedback regulation, cellular uptake of cholesterol was inhibited.

In addition to the TARGET cohort, GSEA performed on the EGA cohort ( Table S7 ) indicated that the cholesterol biosynthesis pathway was more enriched in recurrent or pulmonary metastases than in primary lesions ( Figure 3B ). Furthermore, according to gene expression level comparison, almost all the genes involved in this pathway were observably highly expressed in relapsed or metastatic tumors ( Figure 3C ; Table S8 ). Since recurrence or metastasis lead to treatment failure in most cases (48), and S-IV patients were predicted with poor prognosis, a hint was uncovered that active cholesterol biosynthesis could probably confer an unfavorable outcome.

As patients of S-IV tended to exhibit active cholesterol biosynthesis with a poor prognosis, and this pathway was strongly associated with a poor prognosis in the EGA cohort, we speculated that S-IV patients could benefit from targeting the cholesterol biosynthesis process.

Genes highly expressed in S-IV, showing high-risk scores for mortality (HR>1, log-rank p<0.05), and druggable in the DGIdb database, were identified as potential candidates for targeted therapy. 24 candidates fit the criteria above ( Figure 4A ; Table S9 ). They were ordered by decreasing HR values and top 15 were presented ( Figure 4B ). Among these candidates, CD36, SQLE, and LPL participated in cholesterol metabolism. Each over-representation led to an apparent inferior overall and disease survival ( Figure 4C ). SQLE is one of the rate-limiting enzymes in cholesterol biosynthesis (49). Meanwhile, Food and Drug Administration (FDA)-approved drugs targeting SQLE have been used in trials exploring the treatment of various cancers (50, 51). The impacts of SQLE on the prognosis of the GSE21257 dataset and the EGA cohort were also examined; high expression portended a poor prognosis ( Figure 4D ) and a significant trend of recurrence or metastasis ( Figure 4E ). Hence, we proposed that SQLE was a suitable drug target for osteosarcoma patients of S-IV.

To evaluate the contribution of SQLE to the phenotype associated with tumor growth, we knocked down SQLE via specific siRNAs in osteosarcoma cells, including U2OS and Saos-2 cell lines. To confirm the knockdown effect of SQLE, Q-PCR and western blot were used to detect SQLE expression at the mRNA and protein levels, respectively ( Figures 5A, B ; Supplementary Figures 4A, B ). Experimental results indicated that the silencing of SQLE could substantially inhibit cell proliferation, colony formation, and migration ( Figures 5C–G , Supplementary Figures 4C–E ). Furthermore, we examined the effect of terbinafine, a specific inhibitor of SQLE, on the cell phenotype. According to a conventional pharmacodynamic measure IC50 (drug concentration causing 50% inhibition), terbinafine was found to inhibit cell viability in a time- and dose-dependent manner ( Figure 6A and Supplementary Figure 5A ). Next, we set the concentration gradient of terbinafine based on the value of IC50 at 48h. The colony formation ability decreased dose-dependently ( Figures 6B, C ), and the migrated cells were reduced with the gradient of terbinafine challenge ( Figures 6D, E ; Supplementary Figures 5B, C ). Meanwhile, both SQLE knockdown and terbinafine treatment decreased intracellular cholesterol ( Figure 5H ; Supplementary Figure 4F , Figure 6F , and Supplementary Figure 5D ). Due to the negative feedback regulation of intracellular cholesterol, the protein abundance of SQLE was dose-dependently enhanced by terbinafine challenge in U2OS and Saos-2 cells. ( Supplementary Figures 6A, B ). These results signified that targeting SQLE could reduce proliferating and migrating ability of osteosarcoma via reducing intracellular cholesterol. Thus, SQLE could be a potential therapeutic target for osteosarcoma patients.

Identifying subtypes of osteosarcoma patients at initial biopsy could help make more precise therapeutic strategies early and optimize treatment efficacy (52). We used the PermFIT to determine 142 feature genes for subtypes ( Table S10 ) and employed the machine learning approach SVM by ProMS to build the subtype diagnostic model (SDM). This model was composed of 13 genes with high or low expression levels ( Table S11 ). Patients were subtyped in sequence of S-IV, S-III, S-II, and S-I ( Figure 7A ). To assess the accuracy of the SDM, we applied this model to distinguish the discovery cohort. Most patients were grouped into identical subtypes in this model compared with the NMF approach ( Figure 7B ; Table S12 ). We also applied SDM to the GSE21257 dataset; survival curves suggested that S-IV generated poorer outcomes than the rest ( Figure 7C ; Table S13 ). These results reflected that our SDM achieved excellent prediction accuracy.

To recognize grim prognosis using fewer genes, we established a prognostic classifier using the LASSO Cox regression analysis. 13 genes measured at the SDM were included in the LASSO regression selection. Finally, only 4 of the 13 genes with the most significant appearances were retained as the signatures of our LASSO model ( Tables S14 , 15 ), including PTPRZ1 (HR, 3.51; 95% CI, 1.76-7), SLC8A3 (HR, 4.1; 95% CI, 1.9-8.83), SIGLEC15 (HR, 0.31; 95% CI, 0.15-0.67), and CCR1 (HR, 0.12; 95% CI, 0.02-0.89) ( Figure 7D ; Table S14 ). After coefficient determination, we used the following formulas for this 4-gene model to calculate the risk score:

Risk Score = 1.486×SLC8A3+0.541×PTPRZ1−0.508×CCR1−0.552×SIGLEC15

The model’s high/low scores divided the TARGET cohort into two groups with sharply different prognoses (HR, 5.93; 95% CI, 3.01-11.68; P < 0.001). The 5-year survival rate tended to zero in the high-risk group, while it exceeded 60% in the low-risk group. This 4-gene model also separated GSE21257 patients with favorable or unfavorable prognosis based on its risk scores; high score generated statistically significant mortality (HR, 5.4; 95% CI, 2.06-14.13; P < 0.001) than the low ( Figure 7E ). Thus, the 4-gene model could act as a stable prognostic predictor across various datasets.