The outcomes of adolescents/young adults with osteosarcoma have not improved in decades. The chaotic karyotype of this rare tumor has precluded the identification of prognostic biomarkers and patient stratification. We reasoned that transcriptomic studies should overcome this genetic complexity. RNA sequencing (RNA-seq) of 79 osteosarcoma diagnostic biopsies identified stable independent components that recapitulate the tumor and microenvironment cell composition. Unsupervised classification of the independent components stratified this cohort into favorable (G1) and unfavorable (G2) prognostic tumors in terms of overall survival. Multivariate survival analysis ranked this stratification as the most influential variable. Functional characterization associated G1 tumors with innate immunity and G2 tumors with angiogenic, osteoclastic, and adipogenic activities as well as PPARγ pathway upregulation. A focused gene signature that predicted G1/G2 tumors from RNA-seq data was developed and validated within an independent cohort of 82 osteosarcomas. This signature was further validated with a custom NanoString panel in 96 additional osteosarcomas. This study thus proposes new biomarkers to detect high-risk patients and new therapeutic options for osteosarcoma.

Significance:

These findings indicate that the osteosarcoma microenvironment composition is a major feature to identify hard-to-treat patient tumors at diagnosis and define the biological pathways and potential actionable targets associated with these tumors.

Osteosarcoma, the most common primary bone cancer in adolescents and young adults, presents highly heterogeneous genomic and transcriptomic landscapes as a result of multiple chromosomal rearrangements (1). Such heterogeneity has complicated exploratory studies aiming to determine the key oncogenic drivers underlying disease emergence or progression, and no routinely used prognostic biomarkers or robust stratification have been established thus far. In the meantime, there have been no major changes to treatments for 40 years, and one-third of patients with osteosarcoma still experience treatment failure, primarily with metastatic relapses (2). Over the last 20 years, the focus has been on the genetic exploration of the disease by different DNA analysis techniques from comparative genomic hybridization (CGH) array (3) to whole-exome sequencing (WES; refs. 1, 4–8) or whole-genome sequencing (9). These genomic studies have reported numerous genetic events (mainly copy-number variations, CNV) related to slight cancer cell fitness increases without any clear tumor stratification. Recent studies on clonal dynamics suggest that genetic heterogeneity in osteosarcoma is due to neutral selection leading to high polyclonality despite common phenotypic features (10). Altogether, these observations raise the possibility that major phenotypic osteosarcoma traits could emerge from the epigenetic or/and transcriptional reprogramming of a polyclonal community. In osteosarcoma, this plasticity is probably partially driven by the cross-talk of the tumor cells with its specialized tumor microenvironment (TME) through ligand/receptor interactions and specific physicochemical stresses (acidosis, hypoxia; ref. 11). Therefore, deciphering the tumor cell interaction with the TME composition (12) and immune infiltrates (13) should guide the discovery of new biomarkers associated with tumor fate.

In this context, we reasoned that in osteosarcoma, gene expression at the RNA level, in addition to providing information on the TME composition, should better describe osteosarcoma tumors as a final read-out of the epigenetic and transcriptional modulation of a highly rearranged genetic landscape.

Population and materials

Biological samples at diagnosis were prospectively collected from patients (up to 50 years) registered in the therapeutic approved French OS2006/sarcoma09 trial (NCT00470223; https://www.clinicaltrials.gov/ct2/show/NCT00470223; ref. 14). This study was carried out in accordance with the ethical principles of the Declaration of Helsinki and with Good Clinical Practice guidelines. The study was done in 40 French centers in accordance with the ethical principles of the Declaration of Helsinki and with Good Clinical Practice guidelines. Written informed consent for blood and tumor samples was obtained from the patients or their parents/guardians if they were under 18 years of age upon enrollment. The information given to the patients, the written informed consent used, the collection of samples, and the research project were approved by an independent ethics committee and Institutional Review Boards. As part of the ancillary biological studies, RNA sequencing (RNA-seq), CGH array, and WES were performed at the Gustave Roussy Cancer Campus.

Between 2009 and 2014, 79 frozen osteosarcoma biopsy samples at diagnosis were collected and analyzed by RNA-seq. The clinical and survival characteristics of the 79 patients were comparable with those of the whole osteosarcoma population registered in the OS2006 trial. The M/F sex ratio was 0.55, and the median age was 15.43 years (4.71–36.83). Most primary tumors were located in the limb (94.9%), and 20.2% of patients had metastases at diagnosis. Chemotherapy was administered as per the OS2006 trial: the MTX-etoposide/ifosfamide (M-EI) regimen was used in 86% of patients, the adriamycin/platinum/ifosfamide (API-AI) regimen was used in 12.6% of patients, and 34.1% of patients received zoledronic acid according to the randomization arm. Seventy-seven patients underwent surgery for the primary tumor, and poor histologic response was observed in 24% of the patients. Twenty-nine patients relapsed, with a median latency of 1.55 years (range, 0.09–3.62). Twenty-two patients died, with a median latency of 2.25 years (range, 0.07–6.9). The median follow-up was 4.8 years. Overall, the 3-year progression-free survival (PFS) and overall survival were 65.82% and 82.27%, respectively.

The median age of the 32 patients from the MAPPYACTS trial was 15.25 years (9.7–21.7) with a M/F sex ratio of 0.68. Most of the samples come from progressive relapse (n = 30) and were biopsied from lung metastases (n = 17). No survival and progression data were available at the time of this study.

Nucleic acid extraction

DNA and RNA were isolated using an AllPrep DNA/RNA Mini Kit (Qiagen) according to the manufacturer's instructions. Quantification/qualification were performed using a Nanodrop 2000 spectrophotometer (Thermo Fisher Scientific) and Bioanalyzer DNA 7500 (Agilent Technologies).

RNA-seq

RNA-seq libraries were prepared with the TruSeq Stranded mRNA Kit following the recommendations. The key steps consist of poly(A) mRNA capture with oligo dT beads using 1 μg of total RNA, fragmentation to approximately 400 bp, DNA double-strand synthesis, and ligation of Illumina adaptors followed by amplification of the library by PCR for sequencing. Library sequencing was performed using Illumina sequencers (NextSeq 500 or HiSeq 2000/2500/4000) in 75 bp paired-end mode in both techniques, and data sequencing was processed by bioinformatics analyses. For the optimized detection of potential fusion transcripts by RNA-seq, an in-house designed metacaller approach was used.

Gene expression analysis of RNA-seq data

The quality of stranded paired-ended RNA-seq libraries was evaluated with FastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Reads were aligned with STAR, and quantification was performed with RSEM using the GRCh37 ENSEMBL mRNA dataset as reference sequences.

ICA was performed using BIODICA, one of the fastest and reproducible implementations of fastICA, including ICASSO stability analysis (15). Genes represented by less than 100 reads were filtered out from the independent component analysis (ICA), and the gene expression matrix was log scaled and scaled by genes before decomposition in 50 components. Pearson correlations between independent components (IC) to assess inter-IC relationships were calculated from the metagene matrix. For each IC, we defined contributory genes as genes within 3 SDs of the mean of the corresponding IC metagene vector. For each component, gene set enrichment analysis was performed using MSigDB (16) for positional enrichment and g:Profiler (17) for functional enrichment. For each component, only the most significant enrichment detected among the gene ontology (GO), REACTOME or tissue factor categories is reported. For positional enrichment, only cytobands with a minimum FDR of 1e−4 are reported. When several cytobands from the same chromosomal arm were enriched, only the less significant P values are reported. To build the gene network from the gene expression matrix, we relaxed the criteria of contributory gene selection to increase the sensitivity and capture IC interactions. We selected the contributory genes within 2.5 SDs of the mean of the corresponding IC metagene vector and infer the network using Pearson distance to calculate the edge lengths. Edges with a distance above 0.25 were discarded. Network and subnetwork illustrations and analyses were produced with Cytoscape v3.6.1 (18). We used the REACTOME FI viz plugin (19) to detect clusters in the network with the spectral partition-based network clustering algorithm and estimate their functional enrichments in GO biological process terms.

Hierarchical classification of the metagene matrix and corresponding heatmap was generated with the pheatmap R package using Pearson distance and Ward construction method. To project the ICs in a two-dimensional orthogonal space, we used a principal component analysis of metagene vectors. The Pls-da function from the mixOmics R package (20) was used to select the most contributory ICs to the G1/G2 patient classification.

Kaplan–Meier overall survival and PFS curves were generated using the survminer and survival R packages. The log-rank test was used to compare survival distributions between groups and to provide corresponding P values.

Differential mRNA expression was estimated with the DESeq2 R package (21) from the raw read count table.

A Gaussian regression model using the glmnet R package (parameters: type.measure = “mse,” alpha = 0.8, family = “Gaussian”) was used to select copy-number alterations (CNA; inferred from the CGH array in pair with the RNA-seq samples) modeling the best IC. Likewise, we detected dependencies between ICs corresponding to large chromosomal transcriptional modulation and CNAs from the same chromosomal region.

A logistic regression model using the glmnet R package (22) was used to define a focused signature of 15 genes able to discriminate G1 from G2 RNA-seq libraries (parameters: type.measure = “mse,” alpha = 0.35, family = “binomial”). A cross-validation strategy was used to select the best lambda value. An extended signature of 37 genes has also been generated to design the custom NanoString panel. Our purpose was to maximize the probability to detect the expression of enough genes, with another technology, to predict the G1/G2 stratification. In both signatures, two gage genes (GAGE12D and GAGE12G) and two histone genes (HIST2H2AA3 and HIST2H2AA4) had high sequence similarity. To reduce redundancies in the signature, we only considered one of the homologous genes. We then predicted G1 and G2 tumors with the focused signature from the 82 osteosarcoma RNA-seq data points generated by the TARGET consortium and compared the signature in 42 relapsed osteosarcoma samples in the MAPPYACTS trial (NCT02613962). For TARGET and MAPPYACTS, gene expression was scaled by the mean and the variance of the gene expression in our 79 samples. For TARGET samples, Kaplan–Meier overall survival and PFS curves were generated using the survminer and survival R packages. The log-rank test was used to compare survival distributions between groups and to provide corresponding P values.

Oligonucleotide array CGH assay

In all experiments, sex-matched normal DNA from a pooled human female or male (Promega) was used as a reference. Oligonucleotide array CGH (aCGH) processing was performed as detailed in the manufacturer's protocol (version 7.5; http://www.agilent.com). Equal amounts (500 ng) of tumor and normal DNAs were fragmented with AluI and RsaI (Fermentas). The fragmented DNAs were labeled with cyanine Cy3-deoxyuridine triphosphate (dUTP) or Cy5-dUTP. Hybridization was carried out on SurePrint G3 Human CGH Microarray 4 × 180K (Agilent Technologies) arrays for 24 hours at 65°C in a rotating oven (Robbins Scientific) at 20 rpm. Hybridization was followed by appropriate washing steps.

Oligonucleotide aCGH mixed (joint + individual) preprocessing

Scanning of glass microarrays was performed with an Agilent G2505C DNA Microarray scanner at 100% PMT with 3 μm resolution at 20°C in a low ozone concentration environment. Data were extracted from scanned TIFF images using Feature Extraction software (v11.5.1.1, Agilent), along with protocol CGH_1105_Oct12. All further data processing steps were performed in the R statistical environment v3.4 (http://cran.r-project.org). Acquired raw intensities were transformed to log2(test/ref). Joint normalization of the whole cohort of raw CGH profiles was performed using the cghseg package (v1.0.2.1; ref. 23) using default parameters; a common “wave-effect” track was computed and then subtracted from all individual profiles through a LOWESS regression. As a second step, an individual normalization step was performed for each profile by subtracting precomputed GC-content tracks through a LOWESS regression. Joint segmentation of the whole cohort of normalized CGH profiles was performed using a penalized least squares regression implemented in the copynumber package (v1.16.0; ref. 24). To set a unique value to the gamma (penalty) parameter for the whole cohort, the optimal gamma value was computed for each individual profile using 100 cross-validation folds, and the median of the optimal gamma values was chosen (median = 18, SD = 4.2, range = 11:41). Joint segmentation resulted in 1,604 segments. Individual centering of profiles was performed using an in-house method selecting the most centered mode in the distribution density of probes' log2(ratio) values. No calling of aberrations was performed.

Oligonucleotide aCGH analysis

All genomic coordinates were established on the UCSC human genome build hg19 (25). Hierarchical clustering of samples was performed in R with the segmented data using Euclidean distances and Ward construction method. Clustering of samples based on nonnegative matrix factorization and spherical k-means were performed using the NMF (v0.20.6; ref. 26) and spherical k-means (v0.2.10; ref. 27) packages, respectively. Minimum common region analysis was performed using GISTIC2 (v2.0.22; ref. 28). Differential analyses were performed using nonparametric statistical tests (Wilcoxon for two classes, Kruskal–Wallis for N classes), and all P values were FDR-adjusted using the Benjamini–Hochberg method.

Multivariate proportional hazards model

To model the associations between overall survival and the major clinical factors (sex, histologic response, metastatic status, tumor size, treatment, chemotherapy, and pubertal status) and the G1/G2 signature, we used stability selection with boosting using the Cox model (29). The selection frequencies for these predictors were computed from 1,000 bootstrap samples. We fixed the number q of selected variables per boosting run to 2 and the threshold to define stable variable |{\pi _{th}}$| to 0.9 (which can be relaxed, increasing the risk of false-positive predictor selection). The per-family error rate, corresponding to the expectation of the maximum number of false-positive selected predictors, and the significance level were computed from the two previously defined quantities (q and |{\pi _{th}}$|⁠) from the definition provided by Meinshausen and Bühlmann (30). These analyses were performed using the mboost R package.

The impact of the G1/G2 stratification on the PFS (defined as the time between the inclusion and the occurrence of relapse or death) was assessed using multivariable Cox model adjusted for the standard risk factors (age, sex, metastasis at diagnosis, histologic response, tumor height, Zometa treatment, chemotherapy, puberty status). The proportional hazards assumption was checked using Schoenfeld residuals. This analysis was performed using the coxph function of the R package survival.

Targeted RNA expression analysis

Direct gene expression analysis was conducted with a NanoString custom approach and following the supplier's recommendations. A specific panel of 41 probes of interest including five housekeeping genes (CNOT1, EIF4G2, SF1, SLC39A1, and SURF4) was designed according to NanoString. RNA integrity was measured with a fragment analyzer system (RNA concentration, RNA quality number, percentage of RNA fragments > 300 nt), and RNA concentrations and purity were measured with a Nanodrop ND8000.

Because of the number of targets (<400), 100 ng of total RNA was hybridized to NanoString probes. Our cohort composed of 176 samples was split into two batches of hybridizations, including a no template control (NTC, water) sample and a universal RNA sample (Agilent Technologies, P/N: 740,000) per batch. Before hybridization, NanoString positive and negative controls were added to samples as spike-in controls. Probes and mRNA were hybridized for 16 hours at 65°C prior to processing on the NanoString nCounter preparation station to remove excess probes and immobilize biotinylated hybrids on a cartridge coated with streptavidin. Cartridges were scanned at a maximum scan resolution [555 fields of view (FOV)] on the nCounter Digital Analyzer (NanoString Technologies) to count individual fluorescent barcodes and quantify RNA molecules. NanoString nSolver 4.0 software was used to control raw data and to normalize data. The imaging QC criterion was higher than 93% (threshold: 75% of FOV). The binding density was monitored as positive and negative control signals. Counts obtained for each target in NTC samples (1–9 counts) were deduced to the corresponding target in the samples of interest. The geometric average of housekeeping gene signals was used to normalize the RNA content. The comparison of universal RNA was in platform agreement (r2: 0.98792). Of 176 samples, three samples with low performance were rejected from the primary analysis: one sample had a low binding density (lower threshold: 0.1; EX148_ARN010: 0.09), and two samples had a high mRNA content normalization factor (threshold: 20; EX148_ARN071: 25.7 and EX148_ARN165: 870.6).

G1/G2 strata prediction from the targeted RNA expression panel

To detect suspicious discrepancies, we first compared the targeted and nontargeted RNA expression of the 35 genes as well as the five housekeeping genes. Likewise, we filtered out seven genes with a P value of correlation (Pearson method) higher than 10−5. From the 166 samples, 70 samples, already analyzed by RNA-seq and therefore classified G1 or G2, were chosen as the training set to retrain our logistic model with elastic-net regularization using glmnet on the 28 targeted RNA expression panel. The alpha value was set as 0.8, and the lambda value was chosen by a cross-validation strategy. Then, we predicted the G1/G2 strata of the 96 remaining samples, and Kaplan–Meier overall survival and PFS curves were generated using the survminer and survival R packages. The log-rank test was used to compare survival distributions between groups and to provide corresponding P values.

Data and materials availability

The data are available in the European Genome-Phenome Archive as the EGAS00001005326 study. The code used to compute the analyses and the figures is available at https://github.com/gustaveroussy/OSTEOMICS_pub.

We performed the RNA-seq of 79 primary osteosarcoma tumors sampled at diagnosis from patients accrued in and treated with MTX (n = 69) or API-AI (n = 10) by the first-line OS2006 trial (NCT00470223). These samples recapitulated the main clinical feature distribution of the whole cohort (Supplementary Table S1). To study the transcriptome landscape of both osteosarcoma and TME cells, we first dissociated their respective transcriptome programs. To accomplish this, we decomposed the RNA gene expression matrix by ICA (BIODICA implementation; ref. 15) with the ICASSO stabilization procedure (15) into 50 ICs, also called gene modules.

ICs recapitulate biological functions and their interactions

We first tested whether specific ICs were separately correlated with clinical variables. With the exception of IC2, which was associated with sex and obvious gene enrichment from the Y chromosome (Fig. 1A), no other ICs were significantly associated with clinical items after P-value adjustment. We next characterized ICs by functional enrichment analysis using the most contributory genes of each IC (Supplementary Table S2). We observed that 90% of the ICs were significantly associated with either large chromosomal regions under transcriptional modulation with the MSigDB database (16) or specific cellular/molecular functions with g:Profiler (Fig. 1A; ref. 17). Several ICs appeared specific to altered gene expression involving statistically significant large chromosomal regions. Some of these regions recapitulate known CNVs observed in osteosarcomas (4), as confirmed by a comparison with 70 CGH arrays produced in pair with our RNA-seq samples (Fig. 1A; Supplementary Table S3; Supplementary Fig. S1). Contributory genes from the other ICs recapitulate the transcriptional program of tumor and TME cells such as the bone microenvironment (osteoclasts/bone resorption IC25; osteoblast/ossification IC44), angiogenesis (IC13), the immune response (IC41), neuron projection (IC48), and muscle (IC1). Pearson correlation analysis between ICs revealed a proximity between some of them with two main groups depicted by red and blue boxes (Fig. 1A). Such IC clustering shows that osteosarcoma tumors share two specific IC patterns and suggests a strong relationship among the TME, cancer cells, and large chromosomal region modulations. To decipher these potential functional interactions, we inferred a coexpression network using only the genes contributing significantly to ICs (>3 SDs from the mean). Colored by ICs, genes belonging to the same ICs demonstrate strong interconnection in the network, confirming the existence of a continuum between IC latent variables and real measured gene expression (Fig. 1B). In most cases, we identified matches between the ICs and functional subnetwork annotations performed with REACTOME FI viz (Fig. 1C), confirming the potency of ICA to decompose the gene expression matrix from heterogeneous tumors into functionally relevant gene modules. In addition, annotation of the coexpression network helped us classify genes from ICs into more precise subnetworks. IC41, initially flagged as a general immune response, was subdivided into “adaptive immune response/neutrophil degranulation,” “type 1 IFN response,” and “adaptive immune response.” IC9 was subdivided into “neutrophil degranulation/defensin” and “adipogenesis/peroxisome proliferator-activated receptor (PPAR) signaling pathways.” At a larger scale, the network is structured around several hubs of highly connected ICs/subnetworks, probably mirroring the active cross-talk between the TME and tumor cells.

Figure 1.

A, Correlation matrix of the metagene vectors corresponding to the 50 ICs. The table reports the biological function (g:Profiler) or chromosomal regions (MSigDB) significantly enriched in each component. Chromosomal regions highlighted in green have also been detected as CNAs in paired samples analyzed by CGH array using glmnet regression. B and C, Show the same Inferred gene coexpression network generated with genes contributing the most to each IC. B, Nodes are colored according to their IC membership. C, Nodes are colored according to the clusters found with the spectral partition-based network clustering algorithm of the REACTOME FI viz plugin in Cytoscape. The legend shows the results of the functional enrichment analysis of the major clusters for GO biological processes generated with the REACTOME FI viz plugin.

Figure 1.

A, Correlation matrix of the metagene vectors corresponding to the 50 ICs. The table reports the biological function (g:Profiler) or chromosomal regions (MSigDB) significantly enriched in each component. Chromosomal regions highlighted in green have also been detected as CNAs in paired samples analyzed by CGH array using glmnet regression. B and C, Show the same Inferred gene coexpression network generated with genes contributing the most to each IC. B, Nodes are colored according to their IC membership. C, Nodes are colored according to the clusters found with the spectral partition-based network clustering algorithm of the REACTOME FI viz plugin in Cytoscape. The legend shows the results of the functional enrichment analysis of the major clusters for GO biological processes generated with the REACTOME FI viz plugin.

Close modal

Stable ICs stratify patients by prognosis

We questioned whether IC/subnetwork interconnections reflecting the TME composition were linked to clinical variables. We categorized tumors by hierarchical classification using stable ICs (stability index > 0.5; Fig. 2A). This unsupervised analysis defined two groups of tumors associated with different average living statuses and discriminated 35 low-risk tumors, coined G1 (8.5% death rate, DR), from a higher-risk group, coined G2 (43.1% DR; Fig. 2A). We confirmed the significantly different overall survival, estimated by the Kaplan–Meier method, between the G1 and G2 groups using the log-rank test (P = 0.00042; Fig. 2B). With a median follow-up of 4.8 years, the 3-year overall survival rates were 100% and 67.8% for G1 and G2, respectively. A difference was also observed for PFS, although the difference was less significant (P = 0.042; Fig. 2C). We next challenged the prognostic value of our stratification with a multivariate proportional hazards model including known osteosarcoma prognostic factors. The model tested in this stability analysis identified the G1/G2 stratification at diagnosis as the most contributory factor to overall survival (with 68.5% inclusion frequency; Fig. 2D). A multivariate analysis confirmed the prognostic value of this stratification (HR = 6.3) and its independency from clinical features described as prognostic factors at baseline (metastasis status at diagnosis, HR = 2.8) or after first-line chemotherapy (histologic response, HR = 7.1).

Figure 2.

A, Hierarchical classification of metasample vectors for the 34 most stable ICs (stability index > 0.5). The x-axis corresponds to the 79 RNA-seq samples. The sample and IC dendrograms are colored according to patient groups G1 (pink) and G2 (blue). Clinical annotations on the top of the heatmap show, from top to bottom, for each patient, sex, puberty status, age, histopathologic subtype, tumor height class, histologic response after first line of chemotherapy, presence of metastasis at diagnosis, presence of metastasis at any timepoint, patient that relapsed, living patients. B, Overall survival curves for G1 (pink) and G2 (blue) patients. C, PFS curves for G1 (pink) and G2 (blue) patients. D, Barplot illustrating the proportion of models including clinical variables and G1/G2 stratification for the prediction of the overall survival of the patients in the cohort. E, Univariate and multivariate survival analyses of the G1/G2 stratification and relevant clinical variables. F, Principal component analysis of IC metagene vectors. The names of the ICs with the higher contribution to the principal components are indicated on the plot. Bold and transparent dots show deceased and living patients, respectively. Dots are colored according to G1 (pink) or G2 (blue) membership. G, Loadings of the most contributitory ICs as estimated by partial least squares discriminant analysis to stratify patients as G1 or G2.

Figure 2.

A, Hierarchical classification of metasample vectors for the 34 most stable ICs (stability index > 0.5). The x-axis corresponds to the 79 RNA-seq samples. The sample and IC dendrograms are colored according to patient groups G1 (pink) and G2 (blue). Clinical annotations on the top of the heatmap show, from top to bottom, for each patient, sex, puberty status, age, histopathologic subtype, tumor height class, histologic response after first line of chemotherapy, presence of metastasis at diagnosis, presence of metastasis at any timepoint, patient that relapsed, living patients. B, Overall survival curves for G1 (pink) and G2 (blue) patients. C, PFS curves for G1 (pink) and G2 (blue) patients. D, Barplot illustrating the proportion of models including clinical variables and G1/G2 stratification for the prediction of the overall survival of the patients in the cohort. E, Univariate and multivariate survival analyses of the G1/G2 stratification and relevant clinical variables. F, Principal component analysis of IC metagene vectors. The names of the ICs with the higher contribution to the principal components are indicated on the plot. Bold and transparent dots show deceased and living patients, respectively. Dots are colored according to G1 (pink) or G2 (blue) membership. G, Loadings of the most contributitory ICs as estimated by partial least squares discriminant analysis to stratify patients as G1 or G2.

Close modal

To refine our analysis, we selected ICs contributing the most to the prognostic G1/G2 stratification by their projection in the first two principal component spaces (Fig. 2F). As confirmed by a supervised partial least squares discriminant analysis (Fig. 2G), IC39 best represents the G1/G2 stratification, but several other ICs participate in the distinction between the two groups (Fig. 2F and 2G). The first principal component (PCA1) confirmed the high association of several ICs with prognosis, while the second component (PCA2) was not significantly associated with any clinical item. Therefore, although PCA2 might be associated with an interesting state of tumor cells (IC48: neuroprojection, IC49: EGR1-regulated genes, IC50: KIT and YES1 expression, IC42: E2F1 regulated genes), we further focused our work on the PCA1 functions related to survival. Altogether, our observations suggest that the complex traits related to our stratification and patient prognosis emerged from the interaction of several ICs/subnetworks.

Stratified ICs are related to specific biological functions

To understand at the biological level this hidden complex tumor trait detectable at diagnosis, we functionally characterized the two prognostic groups at the network scale. We overlaid the log2 fold change gene expression between G1 and G2 tumors on the inferred networks (Fig. 3A; Supplementary Fig. S1B).

Figure 3.

A, Selection of subnetworks showing the greatest differential gene expression between G1 and G2 tumors. Nodes are colored according to their log2-fold change. The color scale from red to blue shows the genes more highly expressed in G1 or G2. B, Gene expression distribution for each of the subnetworks and their predicted biological function. Adjusted P value computed by gene set enrichment analysis using the log2-fold change of genes from the network (G1 vs. G2) to detect significant enrichment in subnetwork gene sets. C, Volcano plot illustrating the log2-fold change (x-axis) by the inverse of the log10 of the adjusted P value generated by the differential gene expression analysis between G1 and G2 tumors. Gene names are indicated for genes with adjusted P values lower than 0.001. The color scale from red to blue shows the genes more highly expressed in G1 or G2. D, GISTIC analysis of the CNAs found to be significantly associated with G1 or G2 tumors.

Figure 3.

A, Selection of subnetworks showing the greatest differential gene expression between G1 and G2 tumors. Nodes are colored according to their log2-fold change. The color scale from red to blue shows the genes more highly expressed in G1 or G2. B, Gene expression distribution for each of the subnetworks and their predicted biological function. Adjusted P value computed by gene set enrichment analysis using the log2-fold change of genes from the network (G1 vs. G2) to detect significant enrichment in subnetwork gene sets. C, Volcano plot illustrating the log2-fold change (x-axis) by the inverse of the log10 of the adjusted P value generated by the differential gene expression analysis between G1 and G2 tumors. Gene names are indicated for genes with adjusted P values lower than 0.001. The color scale from red to blue shows the genes more highly expressed in G1 or G2. D, GISTIC analysis of the CNAs found to be significantly associated with G1 or G2 tumors.

Close modal

Two IC41 subnetworks, significantly enriched in genes involved in the “innate immune response/IFN 1 response” and “inflammatory response,” were found upregulated in the favorable G1 tumor (Fig. 3A and B). Interestingly, IC39, a complex component in terms of function and the most influential to the G1 group, appears to be related to the epigenetic reprogramming of the immune TME (e.g., HDAC11/HDAC6 negatively regulates IL10, KDM6A positively regulates IL6; refs. 31, 32). In addition, some G1 tumors also express specific cancer testis antigens (CTA; IC26), known as a source of neoantigens (33). In contrast, the unfavorable G2 group exhibited connected subnetworks summarizing other aspects of the TME previously described as related to poor prognosis and a prometastatic phenotype, such as osteoclast differentiation (IC25; ref. 34), tumor angiogenesis/VEGFR (IC13; ref. 35) and neutrophil degranulation (IC9; ref. 36). Two other isolated subnetworks were strongly upregulated in the G2 group (Fig. 3A and B): (i) the fibroblastic network (IC1), which could reflect either the myofibroblastic reprogramming of osteosarcoma stem cells that is crucial for lung metastasis (37), or the cancer-associated fibroblasts contributing to the mesenchymal-like phenotype and metastasis (38); and (ii) the adipocyte/PPAR signaling pathway (IC9), which is involved in inflammatory states and metabolic reprogramming favoring cancer progression (39).

To complement this functional approach and test the validity of our findings, we performed a differential gene expression analysis between G1 and G2 (Fig. 3C; Supplementary Table S4). Consistent with the above results, in G1 tumors, we observed the reexpression of specific CTAs, including the GAGE family clustered on Xp11.23, in addition to genes implicated in spermatogenesis (e.g., RHOXF2B, GTSF1, RHOXF2, and DAZ3). This reexpression evokes some major epigenomic changes at the methylation level (40) as required by the control of pre-meiosis or meiosis (41). In accordance, we observed a significant upregulation of some genes (Padjusted < 1.3 × 10−4) involved in the RNA gene silencing pathway (e.g., MAEL, DDX4, and TDRD1) and known to suppress transposable element expression during meiosis (42) through the piRNA pathway. Similarly, the functional characteristics of G2 group tumors were confirmed by differential gene expression analysis with the upregulation of genes implicated in adipocyte differentiation, osteoclastogenesis, and angiogenesis. Some of these factors have been associated with unfavorable prognosis in osteosarcoma [e.g., PLA2G2A (43), PAQR3 (44), and HP (45)].

Thus, the findings of all three complementary functional analysis methods were consistent, underlying the major contribution of the TME to osteosarcoma progression, with the innate immune response associated to the favorable G1 tumors while angiogenesis, osteoclastogenesis, and adipogenesis were upregulated in unfavorable G2 tumors.

Stratified ICs are associated with distinct large chromosomal regions

In addition to these distinct TME compositions, we observed altered gene expression involving large chromosomal regions, putatively related to CNA dosage effects (Fig. 1A) and associated with the G1 and G2 stratification groups. From the most influential ICs to the G1 group, several ICs, associated to this large chromosomal region, reflected osteosarcoma tumorigenesis and response to chemotherapy (46) such as the modulation of cytoband 12q14.1 (IC34:CDK4, OS9) or cytobands 4qter and 6p21.1–22.1 [IC6: RUNX2, CDC5L (47), UBR2]. Other G1 group contributory ICs were characterized by the dysregulated expression of telomeric regions (IC24: 15qter, 21qter, 12qter; IC34: 4qter; IC35: 13qter). This feature was not detected in the unfavorable G2 group, although the most positively contributory gene of IC19 was DAXX, which, as part of the ATRX/DAXX/HistoneH3.3 complex, regulates telomere maintenance through the alternative lengthening of telomeres (48). For the unfavorable G2 group, three main chromosomal regions were dysregulated on chromosomes 6p (IC19), 8q (IC23), and 22q (IC33). The first two are well-known high copy gain or amplification events associated with osteosarcoma oncogenesis (4, 49), more frequently observed in recurrent/metastatic osteosarcoma than in primary osteosarcoma (50), and previously linked to unfavorable prognosis. Cytoband 6p is a recurrent amplified region in osteosarcoma (51) that dysregulates the expression level of the oncogene CDCL5 (47). Cytoband 8q copy-number gain is strongly suspected to participate in tumorigenesis through MYC-driven superenhancer signaling (52) and to be a prognostic factor in osteosarcoma (53).

Altogether, those ICs with dysregulated chromosomal regions could emphasize a major contribution of CNVs to G1/G2 stratification. To estimate this potential influence, we integrated our results with 70 CNV profiles (CGH array) paired with our RNA-seq samples. The differential analysis of the aCGH profile using GISTIC2.0 (28) characterized four chromosomes with regions differentially and significantly altered between prognostic groups (Padjusted < 0.1; Supplementary Table S5; Fig. 3D). Hence, gains of 2p21–23 and 22q11.22 cytobands and losses of 11p12-p15 cytobands were detected in the G2 group, while 13q12 region loss was more pronounced in the G1 group than in the G2 group. Surprisingly, among these regions, only the chr22:21859431–22318144 (8.27E-03) region was also identified at the expression level (IC33) as participating in the G1/G2 stratification. We also noticed that the chr8q cytoband, including MYC, was close to significance in the aCGH profile analysis, supporting again the multiple reports about MYC amplification involvement in osteosarcoma (52).

Validation of the osteosarcoma RNA-seq prognostic signature

To validate the robustness of our stratification in an independent cohort, we identified a 15-gene signature predicting G1/G2 tumors by machine learning from the initial gene expression count matrix, with logistic regression regularized by elastic net (Supplementary Table S6). We tested the ability of this signature to predict G1/G2 tumors in an independent cohort of 82 pediatric patients with osteosarcoma whose gene expression count table and paired clinical data are available via open access at https://ocg.cancer.gov/programs/target/projects/osteosarcoma. To check the prediction validity, we compared the overall survival of predicted G1 and G2 tumors and generated related log-rank P values (Fig. 4B). Consistent with observations in the OS2006 cohort, the results in this independent cohort demonstrated that predicted G2 tumors were significantly associated with worse prognosis than predicted G1 tumors (P: 0.00039), supporting the predictive power of this gene signature, regardless of the first-line treatment used. Indeed, in this independent cohort, the first-line chemotherapy used, primarily based on the MAP regimen, was different from that used in the OS2006 cohort (14). To challenge the stratification with a different technology than RNA-seq, we designed a custom NanoString panel based on an extended signature of 28 genes (Supplementary Table S7). We analyzed with this panel 166 samples from the OS2006 cohort. Seventy samples, already analyzed by RNA-seq, were used to retrain our elastic-net model. The NanoString retrained models were used to predict G1/G2 stratification of the 96 remaining tumors. Strikingly, the G1/G2 stratification of these 96 osteosarcomas was again associated to significantly different overall survival (P: 0.023) and PFS (P: 0.036) rates (Fig. 4C; Supplementary Fig. S3C).

Figure 4.

A, Overall survival curves of predicted G1 (pink) and G2 (blue) patients from an independent cohort of 82 patients with osteosarcoma. B, Overall survival curves of predicted G1 (pink) and G2 (blue) patients using a custom NanoString panel on 96 samples from the OS2006 cohort. D, The barplot on the left shows the proportion of G1 and G2 tumors from the 29 of 79 patients who relapsed with tumors submitted to RNA-seq and enrolled in the OS2006 trial. The bar plot on the right shows the proportion of G1 and G2 tumors sampled at relapse among 32 patients with osteosarcoma enrolled in the MAPPYACTS trial.

Figure 4.

A, Overall survival curves of predicted G1 (pink) and G2 (blue) patients from an independent cohort of 82 patients with osteosarcoma. B, Overall survival curves of predicted G1 (pink) and G2 (blue) patients using a custom NanoString panel on 96 samples from the OS2006 cohort. D, The barplot on the left shows the proportion of G1 and G2 tumors from the 29 of 79 patients who relapsed with tumors submitted to RNA-seq and enrolled in the OS2006 trial. The bar plot on the right shows the proportion of G1 and G2 tumors sampled at relapse among 32 patients with osteosarcoma enrolled in the MAPPYACTS trial.

Close modal

Finally, to investigate the reversibility of this prognostic signature throughout disease progression, we used the focused RNA-seq–based signature to stratify 32 patients sampled at relapse and sequenced during the MAPPYACTS trial (Supplementary Table S8). We observed a higher proportion of G2 tumors at relapse (65.6%) than at diagnosis in the OS2006 (55.6%) and TARGET (47.5%) cohorts (Fig. 4D). Interestingly, a similar proportion of unfavorable G2 prognosis (68.9%) was observed for the 29 patients with metastasis from the OS2006 RNA-seq cohort at diagnosis. Therefore, our result suggests the persistence of this stratification despite the disease progression but also an overrepresentation in metastatic disease.

Our work, using an unsupervised machine learning strategy, defines a repertoire of ICs describing the transcriptional program of osteosarcoma tumors and TMEs at diagnosis. Tumor relative composition in ICs stratifies the cohort into favorable and unfavorable prognosis tumors, defining a transcriptomic osteosarcoma classification linked to biological functions. Thus, the functional characterization of the ICs confirmed the observations that favorable prognosis tumors are associated with specific innate immune expression (54) and that unfavorable prognosis tumors is depicted by a TME combining angiogenic, osteoclastic, and adipogenic activities, prone to induce metastases (11, 55,56). At the clinical level, this stratification is congruous with the observed efficacy of multityrosine kinase inhibitors with antiangiogenic activity in relapsed osteosarcoma (57), while the observed inefficacy of zoledronate in front-line osteosarcoma treatment (14) is thought to be partially linked to its action on the immune system (54). At the genetic level, although CNVs are considered major drivers in OS, we found a relatively low number of CNVs associated to the differential disease progression in both groups. From three regions, only the chromosome 22q alteration was not previously described as prognostic in osteosarcoma and may need further investigation (4). This region contains the PI4KAP2, RIMBP3B, RIMBP3C, UBE2L3, YDJC, CCDC116, SDF2L1, MIR301B, MIR130B, PPIL2, YPEL1, MAPK1, PPM1F, and TOP3B genes. The MAPK signaling pathway has been previously proposed as an important driver of the metastatic stage, whereas the PI3K-Akt signaling pathway may participate in the early and late stages of osteosarcoma evolution (58).

Gene expression in oncology often shows high dispersion between samples, which leads to overfitted models or classifications based on noise rather than signal, questioning result reproducibility in other cohorts as well as the introduction of such models for patient use. To confirm this stratification and its relation to prognosis, we used a focused signature of 15 genes within an independent cohort of 82 patients with osteosarcoma primary tumors. Interestingly, despite the different therapy regimens used in the two cohorts, the signature was able to well discriminate favorable and unfavorable prognosis patients, suggesting the existence of two tumor phenotypes, already present at diagnosis, that will respond differentially to treatments. Together with the major influence of the TME in the stratification, our observations suggest an early cross-talk, preceding diagnosis, of the TME with specific tumor cells that could lock the pathology into distinct treatment sensitivities and metastatic potentials. The persistence of the stratification through the disease progression, supported by the overrepresentation of the unfavorable phenotype in tumors at relapse, suggests a perpetuation of this unfavorable cross-talk in distant metastatic sites. Zhou and colleagues reported recently similar observations with a comparable TME composition between primitive tumors and lung metastasis at single-cell level (13).

A valid criticism of RNA-seq–based signatures concerns the reproducibility in laboratories with different technical setups and its deployment as routine clinical practice. To overcome this limitation and confirm the statistical significance of our signature, we designed a custom NanoString panel to validate the association of this stratification with the overall survival and the PFS on additional 96 tumors. With this panel, we propose a quick and reproducible assay to stratify large cohort.

In this context, our results might guide new preclinical explorations and enable treatment stratification in osteosarcomas with drugs already tested in patients but so far without the assistance of companion biomarkers. For instance, immune modulation-based therapies (e.g., mifamurtide) might be beneficial to G1 group, when anti-osteoclastic and anti-angiogenic therapies could be more efficient to treat the high risk G2 group. In addition, this study has shed light on potential new actionable targets from underlying biological pathways rarely reported in osteosarcoma. Among them, PPARγ, master regulator of the adipogenesis, plays also a central role in the bone remodeling by regulating oppositely the osteoclast and osteoblast differentiation (59). Together with its activity on the anti-inflammatory M2 polarization (60), PPARγ is therefore known to upregulate pathways that recapitulate the unfavorable G2 TME and likewise, constitutes a promising therapeutic target to switch the TME from unfavorable to favorable. On the other side, in the favorable G1, we observed an upregulation of the piRNA pathway and CTAs. The CTAs, as potential antigens, are promising target for immunotherapies. The confirmation of their expressions at the protein level in osteosarcoma could, in the near future, offers the opportunity to include patients in new trials. The observed partial expression of the piRNA pathway require further investigations to determine whether, like during the gamete development, it could play a role in the silencing of the transposons and contributes to the osteosarcoma cell physiology.

Overall, our results suggest that personalized therapies in osteosarcoma should be based not only on somatic genetic alterations in the tumor (e.g., mutations/CNVs) but also on RNA expression profiling to consider both the DNA abnormalities transcribed at the RNA level and the TME landscape (23). This work should lead to the development of prognostic biomarkers in the near future to help clinicians identify high risk patients at diagnosis and drive the setup of associated personalized therapy in osteosarcoma.

A. Marchais reports a patent for EP21305080 pending. G. Vassal reports grants from French National Cancer Institute, Fondation ARC, and Imagine for Margo during the conduct of the study. L. Brugières reports grants from French National Cancer Institue (InCa), Ligue Nationale Contre le Cancer, Federation Enfants et Santé, Société Française des Cancers de l'Enfant (SFCE), NOVARTIS, and CHUGAI during the conduct of the study. N. Gaspar reports grants from SFCE 2013, ligue ADO 2015, Etoile de Martin 2015, and MVE 2017 during the conduct of the study; in addition, N. Gaspar has a patent for EP21305080 pending. No disclosures were reported by the other authors.

A. Marchais: Conceptualization, data curation, formal analysis, supervision, validation, investigation, visualization, methodology, writing–original draft, writing–review and editing. M.E. Marques Da Costa: Data curation, formal analysis, validation, investigation, methodology, writing–original draft. B. Job: Data curation, formal analysis, visualization, writing–review and editing. R. Abbas: Data curation, formal analysis, visualization, writing–review and editing. D. Drubay: Data curation, formal analysis, visualization, writing–review and editing. S. Piperno-Neumann: Resources, investigation, writing–review and editing. O. Fromigué: Data curation, investigation, writing–review and editing. A. Gomez-Brouchet: Resources, investigation, writing–review and editing. F. Redini: Resources, investigation, writing–review and editing. R. Droit: Data curation, formal analysis, methodology. C. Lervat: Resources, investigation, writing–review and editing. N. Entz-Werle: Resources, investigation, writing–review and editing. H. Pacquement: Resources, investigation, writing–review and editing. C. Devoldere: Resources, investigation, writing–review and editing. D. Cupissol: Resources, investigation, writing–review and editing. D. Bodet: Resources, investigation, writing–review and editing. V. Gandemer: Resources, investigation, writing–review and editing. M. Berger: Resources, investigation, writing–review and editing. P. Marec-Berard: Resources, investigation, writing–review and editing. M. Jimenez: Resources, data curation, funding acquisition, investigation, writing–review and editing. G. Vassal: Resources, funding acquisition, investigation, writing–review and editing. B. Geoerger: Resources, funding acquisition, investigation, writing–review and editing. L. Brugières: Conceptualization, funding acquisition, investigation, project administration, writing–review and editing. N. Gaspar: Conceptualization, resources, supervision, funding acquisition, investigation, methodology, writing–original draft, project administration, writing–review and editing.

Association Etoile de Martin, la Ligue contre le cancer, the Société Française des Cancers de l'Enfant et l'adolescent (SFCE) with the Fédération Enfants et Santé, Imagine for Margo, Une main vers l'espoir, Fondation ARC, and "Parrainage médecin-chercheur" Gustave Roussy.

For the omics studies conducted by the U1015 Gustave Roussy team, the authors acknowledge:

The funders: Ligue Nationale contre le Cancer, Société Française des Cancers et Leucémies de l'Enfant (SFCE) and Enfants et Santé, Association L'Etoile de Martin Etoile de Martin, Association Imagine For Margo.

The MAPPYACTS protocol is funded by the French national cancer institute, Imagine for Margo, and the foundation ARC.

For the OS2006 trial, the authors are indebted:

- To all the children, young people, and parents who accepted to participate in this study.

- To the data management and data analysis team: Nadir CHEURFA, Christine LARUE, Camille DUMONT, Saïd MAALEM, Guillaume DANTON, Violeta GAVEIKAITE, Muriel BEAUBRUN, Pascale JAN, Muriel WARTELLE, Anne-Sophie VEILLARD, Audrey MAUGUEN, Caroline DOMERG, and Bob-Valéry OCCEAN.

- To the sponsor coordinating team who were in charge of the on-site source data verification and data collection and vouch for the accuracy of the data: Marta JIMENEZ, Karine BUFFARD, Jessy DELAYE, Céline MAHIER, Naïma BONNET, Patricia NEZAN, Khadija CHERIF, Bochra HACHANI, Julie GARRABEY, Flavie LE BOURGEOIS, Julien KAKOU, Dilan USTUNER, Sabrina CIESLIK, Raquel CAPITAO, Rim EL ABED, Charlotte GUYADER, Nicolas DE SOUSA, Mélanie PARAS, Valérie PONS, UNICANCER, Paris.

- To the principal investigators of the OS2006 trial: Laurence BRUGIERES, Gustave Roussy, Villejuif; Sophie PIPERNO-NEUMANN, Institut Curie, Paris.

- To the statistician of the trial: Marie-Cécile LE DELEY, Centre Oscard Lambret, Lille; Aurélie CHEVANCE Gustave Roussy, Villejuif; Rachid ABBAS, Gustave Roussy, villejuif.

- To all investigators who participated in the trial: Jacques-Olivier BAY, Centre Jean Perrin, Clermont Ferrand; Claire BERGER, CHU St Etienne, St Etienne; Jean-Pierre BERGERAT, Hôpital de Hautepierre, Strasbourg; François BERTUCCI, Institut Paoli Calmettes, Marseille; Binh BUI NGUYEN, Institut Bergonié, Bordeaux; Jean-Yves BLAY, Centre Léon Bérard, and Claude Bernard University, Lyon; Emmanuelle BOMPAS, Institut de Cancérologie de l'Ouest, Saint Herblain; Liana-Stefania CARAUSU, CHU Morvan, Brest; Loïc CHAIGNEAU, CHU Besançon, Besançon; Christine CHEVREAU, Institut Claudius Regaud, Toulouse; Olivier COLLARD, Institut de Cancérologie de la Loire, St Etienne; Nadège CORRADINI, Hôpital Mère-enfant, Nantes; Gérard COUILLAULT, CHU de Dijon, Dijon; Didier CUPISSOL, Institut Régional du Cancer de Montpellier/Val d'Aurelle, Montpellier; Corinne DELCAMBRE, Centre François Baclesse, Caen; Anne DEVILLE, Hôpital de l'Archet 2, Nice; Catherine DEVOLDERE, CHU Amiens Picardie, Amiens; Florence DUFFAUD, Hôpital La Timone, Marseille; Jean-Christophe EYMARD, Institut Jean Godinot, Reims; Natacha ENTZ-WERLé, CHU Hautepierre, Strasbourg; Virginie GANDEMER, CHU de Rennes, Rennes; Jean-Claude GENTET, CHU La Timone, Marseille; Stéphanie GORDE GROSJEAN, Hopital Américain, Reims; Cécile GUILLEMET, Centre Henri Becquerel, Rouen; Nicolas ISAMBERT, Centre G.F Leclerc, Dijon; Antoine ITALIANO, Institut Bergonié, Bordeaux; Justyna KANOLD, CHU Estaing, Clermont Ferrand; Pierre KERBRAT, Centre Eugene Marquis, Rennes; Véronique LAITHIER, CHU Jean Minjoz, Besançon; Valérie LAURENCE, Institut Curie, Paris; Axel LE CESNE, Institut Gustave Roussy, Villejuif; Odile LEJARS, Hôpital de la Clocheville, Tours; Cyril LERVAT, Centre Oscar Lambret, Lille; Claude LINASSIER, CHRU Bretonneau, Tours; Patrick LUTZ, Hôpital de Hautepierre, Strasbourg; Perrine MAREC-BéRARD, Centre Leon-Berrard, Lyon; Frédéric MILLOT, Hôpital Jean Bernard, Poitiers; Odile MINCKES, CHU Côte de Nacre, Caen; Hélène PACQUEMENT, Insitut Curie, Paris; Isabelle PELLIER, CHU d'Angers, Angers; Nicolas PENEL, Centre Oscar Lambret, Lille; Christophe PIGUET, CHU Dupuytren, Limoges; Dominique PLANTAZ, CHU Grenoble, Grenoble; Graziella RAIMONDO, Croix Rouge Française, Margency; Isabelle RAY-COQUARD, Centre Léon Bérard, Lyon; Maria RIOS, Centre Alexis Vautrin, Nancy; Hervé RUBIE, Hôpital des enfants, Toulouse; Laure SAUMET, Hôpital Arnaud de Villeneuve, Montpellier; Claudine SCHMITT, Hôpital d'enfants, Nancy; Patrick SOULIE, Institut de Cancérologie de l'Ouest, Nantes; Marie-Dominique TABONE, Hôpital Trousseau, Paris; Antoine THYSS, Centre Antoine Lacassagne, Nice; Jean-Pierre VANNIER, CHUHôpital Charles Nicolle, Rouen; Cécile VERITÉ-GOULARD, Hôpital des enfants, Bordeaux.

- To the clinical research assistants and research nurses who collected the data: Frédéric DUEE, Centre Jean Perrin, Clermont Ferrand; Catherine SOLER, CHU St Etienne, St Etienne; Christine KRAVANJA and Hélene ETIENNE, Hôpital de Hautepierre, Strasbourg; Nathalie NICOLAS, Institut Paoli Calmettes, Marseille; Nicolas ROUGE and Simon DUPUIS, Institut de Cancérologie de l'Ouest, Nantes; Imene HEZAM, Institut Gustave Roussy, Villejuif; Sabrina SELLAN-ALBERT, Institut Bergonié, Bordeaux; Typhaine THOMAS, CHU Morvan, Brest; Esma KARAYIGIT, CHU Besançon, Besançon; Caroline TRESSENS, Institut Claudius Regaud, Toulouse; Séverine MARCHAND, Institut de Cancérologie de la Loire, St Etienne; Jean-Pierre RANARIVELO, CHU, Nantes; Frédérique DEBOMY, CHU de Dijon, Dijon; Angélique CABOCO, Institut Régional du Cancer de Montpellier/Val d'Aurelle, Montpellier; Sophie DANET, Centre François Baclesse, Caen; Isabelle CHAMPENOIS, Hôpital de l'Archet 2, Nice; Alicia PLAYE, CHU Amiens Picardie, Amiens; Josiane ONEZIME, Hôpital La Timone, Marseille; Véronique DANGMANN, Institut Jean Godinot, Reims; Isabelle DEBROISE, CHU de Rennes, Rennes; Sylvie ABED, Hôpital pour enfants la Timone, Marseille; Eric MARQUIS, Hopital Américain, Reims; Sabine LORIVEL, Centre Henri Becquerel, Rouen; Sophie DEBARD, Centre G.F Leclerc, Dijon; Nadège ROUEL-RABIAU, CHU Estaing, Clermont Ferrand; Euriell FORTIN, Centre Eugene Marquis, Rennes; Stéphanie CLERC-GUAY, CHU Jean Minjoz, Besançon; Pauline SMIS, Centre Oscar Lambret, Lille; Christine COLOMBAT, CHRU Bretonneau, Tours; Florence VINCENT, Hôpital de Hautepierre, Strasbourg; Cécile GIRAUD and Sandrine GUILLEMIN, Centre Léon Bérard, Lyon; Sandrine RAFERT and Violaine GOYEAU, Hôpital Jean Bernard, Poitiers; Esther LE BRETON, CHU Côte de Nacre, Caen; Anne MACLOU, Institut Curie, Paris; Josué RAKOTONJANAHARY, CHU d'Angers, Angers; Julie GRENIER, CHU Dupuytren, Limoges; Sandrine BILLET, CHU Grenoble, Grenoble; Mireille HIRSCH, Centre Alexis Vautrin, Nancy; Nathalie COUTEAU, Hôpital des enfants, Toulouse; Karine PERRRIN-LOPEZ, Hôpital Arnaud de Villeneuve, Montpellier; Sandrine PALL-KONDOLFF, Hôpital d'enfants, Nancy; Christine RAGU, Hôpital Trousseau, Paris; Colin DEBAIGT, Centre Antoine Lacassagne, Nice; Cécile GIROD, CHU-Hôpital Charles Nicolle, Rouen; Agathe LARRIEU, Hôpital des enfants, Bordeaux.

- To the pathology review board: Sébastien AUBERT, CHRU Lille; Corinne BOUVIER, CHU La Timone, Marseille; Gonzague DE PINIEUX, CHU Tours; Frédérique DIJOUD, CHU Dijon; Jean Marc GUINEBRETIERE, Institut Curie Saint Cloud; Anne GOMEZ-BROUCHET, CHU Toulouse; Frédérique LAROUSSERIE, Hôpital Cochin, Paris; Luc MARCELIN, CHU Strasbourg.

- To all the other pathologists involved in the project: Monique AUGRAND, Hôpital Robert Boulin, Libourne; Florence BURTIN, CHU Rennes; Elisabeth CASSAGNAU, CHU Nantes; Jean Michel COINDRE, Institut Bergonié, Bordeaux; Carole CORDONNIER, CHU Amiens; Aurore COULOMB, Hôpital Trousseau, Paris; Louise GALMICHE, Hôpital Necker, Paris; Bernadette KANTELIP, CHU Besançon; Jerzy KLIJANIENKO, Institut Curie, Paris; Patrice JOSSET, Hôpital Trousseau, Paris; Sébastien LEPREUX, CHU Bordeaux; Aurélien MARAN, CHU Montpellier; Dominique MICHIELS, CHU Nice; Karine MONTAGNE, CHU Brabois, Nancy; Michel MOREAU, CHU Toulouse; Corinne PASQUIER, CHU Grenoble; Michel PEUCHMAUR, Hôpital Robert Debré, Paris; Dominique RANCHERE, Centre Léon Bérard Lyon; Jean Michel ROMEO, Paris; Pierre ROUSSELOT, CHU Caen; Jean Christophe SABOURIN, CHU Rouen; Christine SCHAEFFER, CHU Dijon; Virginie VERKARRE, Hôpital Necker, Paris; Jean Michel VIGNAUD, Hôpital Central, Nancy; Dominique ZACHAR, CHU Reims; Elie ZAFRANI, Hôpital Henri Mondor, Creteil.

- To the Independent Data Monitoring Committee: Mark BERNSTEIN, Serge LEYVRAZ, Agnès LINGLART, Monia OUALI, and Jeremy WHELAN.

- To the funders of the OS2006 trial: INCa, NOVARTIS, CHUGAI, Ligue Nationale contre le Cancer, Fédération Enfants et Santé, Société Française des Cancers et Leucémies de l'Enfant.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

Note: Supplementary data for this article are available at Cancer Research Online (http://cancerres.aacrjournals.org/).

1.
Lorenz
S
,
Barøy
T
,
Sun
J
,
Nome
T
,
Vodák
D
,
Bryne
J-C
, et al
.
Unscrambling the genomic chaos of osteosarcoma reveals extensive transcript fusion, recurrent rearrangements and frequent novel TP53 aberrations
.
Oncotarget
2016
;
7
:
5273
88
.
2.
Trama
A
,
Botta
L
,
Foschi
R
,
Ferrari
A
,
Stiller
C
,
Desandes
E
, et al
.
Survival of European adolescents and young adults diagnosed with cancer in 2000–07: population-based data from EUROCARE-5
.
Lancet Oncol
2016
;
17
:
896
906
.
3.
Atiye
J
,
Wolf
M
,
Kaur
S
,
Monni
O
,
Böhling
T
,
Kivioja
A
, et al
.
Gene amplifications in osteosarcoma-CGH microarray analysis: CGH microarray analysis of osteosarcoma
.
Genes Chromosomes Cancer
2005
;
42
:
158
63
.
4.
Perry
JA
,
Kiezun
A
,
Tonzi
P
,
Van Allen
EM
,
Carter
SL
,
Baca
SC
, et al
.
Complementary genomic approaches highlight the PI3K/mTOR pathway as a common vulnerability in osteosarcoma
.
Proc Natl Acad Sci U S A
2014
;
111
:
E5564
73
.
5.
Kovac
M
,
Blattmann
C
,
Ribi
S
,
Smida
J
,
Mueller
NS
,
Engert
F
, et al
.
Exome sequencing of osteosarcoma reveals mutation signatures reminiscent of BRCA deficiency
.
Nat Commun
2015
;
6
:
8940
.
6.
Bousquet
M
,
Noirot
C
,
Accadbled
F
,
Sales de Gauzy
J
,
Castex
MP
,
Brousset
P
, et al
.
Whole-exome sequencing in osteosarcoma reveals important heterogeneity of genetic alterations
.
Ann Oncol
2016
;
27
:
738
44
.
7.
Chen
KS
,
Kwon
WS
,
Kim
J
,
Heo
SJ
,
Kim
HS
,
Kim
HK
, et al
.
A novel TP53-KPNA3 translocation defines a de novo treatment-resistant clone in osteosarcoma
.
Cold Spring Harb Mol Case Stud
2016
;
2
:
a000992
.
8.
Chiappetta
C
,
Mancini
M
,
Lessi
F
,
Aretini
P
,
De Gregorio
V
,
Puggioni
C
, et al
.
Whole-exome analysis in osteosarcoma to identify a personalized therapy
.
Oncotarget
2017
;
8
:
80416
28
.
9.
Behjati
S
,
Tarpey
PS
,
Haase
K
,
Ye
H
,
Young
MD
,
Alexandrov
LB
, et al
.
Recurrent mutation of IGF signalling genes and distinct patterns of genomic rearrangement in osteosarcoma
.
Nat Commun
2017
;
8
:
15936
.
10.
Gambera
S
,
Abarrategi
A
,
González-Camacho
F
,
Morales-Molina
Á
,
Roma
J
,
Alfranca
A
, et al
.
Clonal dynamics in osteosarcoma defined by RGB marking
.
Nat Commun
2018
;
9
:
3994
.
11.
Corre
I
,
Verrecchia
F
,
Crenn
V
,
Redini
F
,
Trichet
V
.
The osteosarcoma microenvironment: a complex but targetable ecosystem
.
Cells
2020
;
9
:
976
.
12.
Petitprez
F
,
Sun
C-M
,
Lacroix
L
,
Sautès-Fridman
C
,
de Reyniès
A
,
Fridman
WH
.
Quantitative analyses of the tumor microenvironment composition and orientation in the era of precision medicine
.
Front Oncol
2018
;
8
:
390
.
13.
Zhou
Y
,
Yang
D
,
Yang
Q
,
Lv
X
,
Huang
W
,
Zhou
Z
, et al
.
Single-cell RNA landscape of intratumoral heterogeneity and immunosuppressive microenvironment in advanced osteosarcoma
.
Nat Commun
2020
;
11
:
6322
.
14.
Piperno-Neumann
S
,
Le Deley
M-C
,
Rédini
F
,
Pacquement
H
,
Marec-Bérard
P
,
Petit
P
, et al
.
Zoledronate in combination with chemotherapy and surgery to treat osteosarcoma (OS2006): a randomised, multicentre, open-label, phase 3 trial
.
Lancet Oncol
2016
;
17
:
1070
80
.
15.
Kairov
U
,
Cantini
L
,
Greco
A
,
Molkenov
A
,
Czerwinska
U
,
Barillot
E
, et al
.
Determining the optimal number of independent components for reproducible transcriptomic data analysis
.
BMC Genomics
2017
;
18
:
712
.
16.
Liberzon
A
,
Subramanian
A
,
Pinchback
R
,
Thorvaldsdottir
H
,
Tamayo
P
,
Mesirov
JP
.
Molecular signatures database (MSigDB) 3.0
.
Bioinformatics
2011
;
27
:
1739
40
.
17.
Raudvere
U
,
Kolberg
L
,
Kuzmin
I
,
Arak
T
,
Adler
P
,
Peterson
H
, et al
.
g:Profiler: a web server for functional enrichment analysis and conversions of gene lists (2019 update)
.
Nucleic Acids Res
2019
;
47
:
W191
8
.
18.
Shannon
P
.
Cytoscape: a software environment for integrated models of biomolecular interaction networks
.
Genome Res
2003
;
13
:
2498
504
.
19.
Haw
R
,
Loney
F
,
Ong
E
,
He
Y
,
Wu
G
.
Perform pathway enrichment analysis using ReactomeFIViz
.
Methods Mol Biol
2020
;
2074
:
165
79
.
20.
Le Cao
K-A
,
Rohart
F
,
Gonzalez
I
,
Dejean
S
.
mixOmics: Omics Data Integration Project
.
R package version 6.1.1. [Internet]
;
2016
.
Available from
: http://www.bioconductor.org/packages/release/bioc/html/mixOmics.html.
21.
Love
MI
,
Huber
W
,
Anders
S
.
Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2
.
Genome Biol
2014
;
15
:
550
.
22.
Friedman
J
,
Hastie
T
,
Tibshirani
R
.
Regularization paths for generalized linear models via coordinate descent
.
J Stat Softw
2010
;
33
:
1
22
.
23.
Picard
F
,
Lebarbier
E
,
Hoebeke
M
,
Rigaill
G
,
Thiam
B
,
Robin
S
.
Joint segmentation, calling, and normalization of multiple CGH profiles
.
Biostatistics
2011
;
12
:
413
28
.
24.
Nilsen
G
,
Liestøl
K
,
Loo
PV
,
Moen Vollan
HK
,
Eide
MB
,
Rueda
OM
, et al
.
Copynumber: efficient algorithms for single- and multi-track copy number segmentation
.
BMC Genomics
2012
;
13
:
591
.
25.
Karolchik
D
.
The UCSC genome browser database
.
Nucleic Acids Res
2003
;
31
:
51
4
.
26.
Gaujoux
R
,
Seoighe
C
.
A flexible R package for nonnegative matrix factorization
.
BMC Bioinformatics
2010
;
11
:
367
.
27.
Hornik
K
,
Feinerer
I
,
Kober
M
,
Buchta
C
.
Spherical k-means clustering
.
J Stat Softw
2012
;
50
:
1
22
.
28.
Mermel
CH
,
Schumacher
SE
,
Hill
B
,
Meyerson
ML
,
Beroukhim
R
,
Getz
G
.
GISTIC2.0 facilitates sensitive and confident localization of the targets of focal somatic copy-number alteration in human cancers
.
Genome Biol
2011
;
12
:
R41
.
29.
Hofner
B
,
Boccuto
L
,
Göker
M
.
Controlling false discoveries in high-dimensional situations: boosting with stability selection
.
BMC Bioinformatics
2015
;
16
:
144
.
30.
Meinshausen
N
,
Bühlmann
P
.
Stability selection: stability Selection
.
J R Stat Soc Ser B Stat Methodol
2010
;
72
:
417
73
.
31.
Cheng
F
,
Lienlaf
M
,
Perez-Villarroel
P
,
Wang
H-W
,
Lee
C
,
Woan
K
, et al
.
Divergent roles of histone deacetylase 6 (HDAC6) and histone deacetylase 11 (HDAC11) on the transcriptional regulation of IL10 in antigen presenting cells
.
Mol Immunol
2014
;
60
:
44
53
.
32.
de Groot
AE
,
Pienta
KJ
.
Epigenetic control of macrophage polarization: implications for targeting tumor-associated macrophages
.
Oncotarget
2018
;
9
:
20908
27
.
33.
Salmaninejad
A
,
Zamani
MR
,
Pourvahedi
M
,
Golchehre
Z
,
Hosseini Bereshneh
A
,
Rezaei
N
.
Cancer/testis antigens: expression, regulation, tumor invasion, and use in immunotherapy of cancers
.
Immunol Invest
2016
;
45
:
619
40
.
34.
Kelleher
FC
,
O'Sullivan
H
.
Monocytes, macrophages, and osteoclasts in osteosarcoma
.
J Adolesc Young Adult Oncol
2017
;
6
:
396
405
.
35.
Yang
J
,
Yang
D
,
Sun
Y
,
Sun
B
,
Wang
G
,
Trent
JC
, et al
.
Genetic amplification of the vascular endothelial growth factor (VEGF) pathway genes, including VEGFA, in human osteosarcoma: VEGFA amplification in osteosarcoma
.
Cancer
2011
;
117
:
4925
38
.
36.
Mollinedo
F
.
Neutrophil degranulation, plasticity, and cancer metastasis
.
Trends Immunol
2019
;
40
:
228
42
.
37.
Zhang
W
,
Zhao
J-M
,
Lin
J
,
Hu
C-Z
,
Zhang
W-B
,
Yang
W-L
, et al
.
Adaptive fibrogenic reprogramming of osteosarcoma stem cells promotes metastatic growth
.
Cell Rep
2018
;
24
:
1266
77
.
38.
Quail
DF
,
Joyce
JA
.
Microenvironmental regulation of tumor progression and metastasis
.
Nat Med
2013
;
19
:
1423
37
.
39.
Harmon
GS
,
Lam
MT
,
Glass
CK
.
PPARs and lipid ligands in inflammation and metabolism
.
Chem Rev
2011
;
111
:
6321
40
.
40.
Fratta
E
,
Coral
S
,
Covre
A
,
Parisi
G
,
Colizzi
F
,
Danielli
R
, et al
.
The biology of cancer testis antigens: putative function, regulation and therapeutic potential
.
Mol Oncol
2011
;
5
:
164
82
.
41.
McFarlane
RJ
,
Wakeman
JA
.
Meiosis-like functions in oncogenesis: a new view of cancer
.
Cancer Res
2017
;
77
:
5712
6
.
42.
Ozata
DM
,
Gainetdinov
I
,
Zoch
A
,
O'Carroll
D
,
Zamore
PD
.
PIWI-interacting RNAs: small RNAs with big functions
.
Nat Rev Genet
2019
;
20
:
89
108
.
43.
Mintz
MB
,
Sowers
R
,
Brown
KM
,
Hilmer
SC
,
Mazza
B
,
Huvos
AG
, et al
.
An expression signature classifies chemotherapy-resistant pediatric osteosarcoma
.
Cancer Res
2005
;
65
:
1748
54
.
44.
Ma
Z
,
Wang
Y
,
Piao
T
,
Li
Z
,
Zhang
H
,
Liu
Z
, et al
.
The tumor suppressor role of PAQR3 in osteosarcoma
.
Tumor Biol
2015
;
36
:
3319
24
.
45.
Kwon
J-O
,
Jin
WJ
,
Kim
B
,
Ha
H
,
Kim
H-H
,
Lee
ZH
.
Haptoglobin acts as a TLR4 ligand to suppress osteoclastogenesis via the TLR4–IFN-β axis
.
J Immunol
2019
;
202
:
3359
69
.
46.
Martin
JW
,
Chilton-MacNeill
S
,
Koti
M
,
van Wijnen
AJ
,
Squire
JA
,
Zielenska
M
.
Digital expression profiling identifies RUNX2, CDC5L, MDM2, RECQL4, and CDK4 as potential predictive biomarkers for neo-adjuvant chemotherapy response in paediatric osteosarcoma
.
PLoS One
2014
;
9
:
e95843
.
47.
Lu
X-Y
,
Lu
Y
,
Zhao
Y-J
,
Jaeweon
K
,
Kang
J
,
Xiao-Nan
L
, et al
.
Cell cycle regulator gene CDC5L, a potential target for 6p12-p21 amplicon in osteosarcoma
.
Mol Cancer Res
2008
;
6
:
937
46
.
48.
Hoelper
D
,
Huang
H
,
Jain
AY
,
Patel
DJ
,
Lewis
PW
.
Structural and mechanistic insights into ATRX-dependent and -independent functions of the histone chaperone DAXX
.
Nat Commun
2017
;
8
:
1193
.
49.
Selvarajah
S
,
Yoshimoto
M
,
Ludkovski
O
,
Park
PC
,
Bayani
J
,
Thorner
P
, et al
.
Genomic signatures of chromosomal instability and osteosarcoma progression detected by high resolution array CGH and interphase FISH
.
Cytogenet Genome Res
2008
;
122
:
5
15
.
50.
Yen
C-C
,
Chen
W-M
,
Chen
T-H
,
Chen
W Y-K
,
Chen
P C-H
,
Chiou
H-J
.
Identification of chromosomal aberrations associated with disease progression and a novel 3q13.31 deletion involving LSAMP gene in osteosarcoma
.
Int J Oncol
2009
;
35
:
775
88
.
51.
Lau
CC
,
Harris
CP
,
Lu
X-Y
,
Perlaky
L
,
Gogineni
S
,
Chintagumpala
M
, et al
.
Frequent amplification and rearrangement of chromosomal bands 6p12-p21 and 17p11.2 in osteosarcoma
.
Genes Chromosomes Cancer
2004
;
39
:
11
21
.
52.
Chen
D
,
Zhao
Z
,
Huang
Z
,
Chen
D-C
,
Zhu
X-X
,
Wang
Y-Z
, et al
.
Super enhancer inhibitors suppress MYC driven transcriptional amplification and tumor progression in osteosarcoma
.
Bone Res
2018
;
6
:
11
.
53.
Tarkkanen
M
,
Elomaa
I
,
Blomqvist
C
,
Kivioja
AH
,
Kellokumpu-Lehtinen
P
,
Böhling
T
, et al
.
DNA sequence copy number increase at 8q: a potential new prognostic marker in high-grade osteosarcoma
.
Int J Cancer
1999
;
84
:
114
21
.
54.
Gomez-Brouchet
A
,
Illac
C
,
Gilhodes
J
,
Bouvier
C
,
Aubert
S
,
Guinebretiere
JM
, et al
.
CD163-positive tumor-associated macrophages and CD8-positive cytotoxic lymphocytes are powerful diagnostic markers for the therapeutic stratification of osteosarcoma patients: an immunohistochemical analysis of the biopsies fromthe French OS2006 phase 3 trial
.
Oncoimmunology
2017
;
6
:
e1331193
.
55.
Heymann
M-F
,
Lézot
F
,
Heymann
D
.
The contribution of immune infiltrates and the local microenvironment in the pathogenesis of osteosarcoma
.
Cell Immunol
2019
;
343
:
103711
.
56.
Zhu
L
,
McManus
MM
,
Hughes
DPM
.
Understanding the biology of bone sarcoma from early initiating events through late events in metastasis and disease progression
.
Front Oncol
2013
;
3
:
230
.
57.
Duffaud
F
,
Mir
O
,
Boudou-Rouquette
P
,
Piperno-Neumann
S
,
Penel
N
,
Bompas
E
, et al
.
Efficacy and safety of regorafenib in adult patients with metastatic osteosarcoma: a non-comparative, randomised, double-blind, placebo-controlled, phase 2 study
.
Lancet Oncol
2019
;
20
:
120
33
.
58.
Wang
D
,
Niu
X
,
Wang
Z
,
Song
C-L
,
Huang
Z
,
Chen
K-N
, et al
.
Multiregion sequencing reveals the genetic heterogeneity and evolutionary history of osteosarcoma and matched pulmonary metastases
.
Cancer Res
2019
;
79
:
7
20
.
59.
Ahmadian
M
,
Suh
JM
,
Hah
N
,
Liddle
C
,
Atkins
AR
,
Downes
M
, et al
.
PPARγ signaling and metabolism: the good, the bad and the future
.
Nat Med
2013
;
19
:
557
66
.
60.
Bouhlel
MA
,
Derudas
B
,
Rigamonti
E
,
Dièvart
R
,
Brozek
J
,
Haulon
S
, et al
.
PPARgamma activation primes human monocytes into alternative M2 macrophages with anti-inflammatory properties
.
Cell Metab
2007
;
6
:
137
43
.

Supplementary data