Method for identifying hard-to-treat osteosarcoma patients at diagnosis and improving their outcome by providing new therapy

FIELD OF THE INVENTION

The present invention relates to the field of anticancer treatment. In particular, the present invention concerns the prognosis of osteosarcoma in an individual, and provides an in vitro method for determining this prognosis, as well as diagnostic kits to perform this method. Also described herein is the use of an inhibitor of the PPARy pathway as an antineoplastic treatment to improve the overall survival of patients with a poor prognosis osteosarcoma.

BACKGROUND OF THE INVENTION

Osteosarcoma, the most common primary bone cancer in adolescents and young adults, presents a highly heterogeneous genomic and transcriptomic landscape as a result of multiple chromosomal rearrangements ^{1 2} . Such heterogeneity has complicated exploratory studies to determine the key oncogenic drivers underlying disease emergence or progression, and no routinely used prognostic biomarkers or robust stratification has been defined so far. As a result, there have been no major changes to treatments for 40 years, and a third of patients with osteosarcoma still experience treatment failure, primarily with metastatic relapses ^3-7 . Over the last 20 years, the focus has been on genetic exploration of the disease by different DNA analysis techniques from Comparative Genomic Hybridization array (CGH) ^8-12 to Whole Exome (WES) ^{1 ,2,13-16} or Genome Sequencing (WGS) ^17-19 . These have reported numerous genetic events (mainly Copy Number Variations, CNV) related to slight/medium tumor fitness increase without any clear tumor stratification.

It has recently been demonstrated that despite common phenotypic features, primary osteosarcoma is highly polyclonal under near neutral selection ²⁰ . These observations raise the possibility that the major phenotypic osteosarcoma traits are an emergent characteristic from the epigenetic or/and transcriptional reprogramming of heterogeneous genomic rearrangements in a polyclonal community. In other words, progression is a consequence of tumor plasticity.

In addition, in the highly constrained bone environment of osteosarcoma, it is challenging to evaluate tumor cell co-evolution with the tumor microenvironment (TME) composition and stromal cell proportions or activation states ^{21 22} . SUMMARY OF THE INVENTION

To better understand the mechanisms underlying osteosarcoma and thus enrich the armamentarium against this disease, the inventors reasoned that, the transcriptomic landscape could overcome the genetic complexity of this disease and provide more interpretable information. With unsupervised machine learning strategy, they defined the repertoire of gene components describing osteosarcoma tumor clones and TME. They observed that component interactions through co-expression stratify the cohort into good and poor prognosis tumors.

Functional characterization of the components associates good prognosis tumors with specific innate immune expression and poor prognosis tumors with angiogenic, osteoclastic, and adipogenic activities, with distinct CNVs specific to each group. These distinct functional characteristics can be used to stratify treatment in osteosarcomas, for instance immune modulation (e.g. mifamurtide) for G1 group tumors (good prognosis), and anti-osteoclastic and anti-angiogenic therapies for G2 (hard-to- treat) group.

The inventors also identified underlying biological pathways of osteosarcoma, involving PPARy, piRNA or CTAs which highlight new actionable targets.

Finally, they identified a panel of 37 genes involved in osteosarcoma and confirmed the predictive power of a minimal prognostic signature of 15 genes.

According to a first aspect, the present invention thus relates to an in vitro method for determining the prognosis of osteosarcoma in an individual, comprising measuring the expression levels of a collection of signature genes from a biological sample taken from said individual, applying the expression levels measured to a predictive model associating expression levels of said collection of signature genes with osteosarcoma outcome, and evaluating the output of said predictive model to determine prognosis of osteosarcoma in said individual.

The invention also pertains to a diagnostic kit for predicting the progression of osteosarcoma in a subject by measuring the level of expression of a collection of signature genes from a biological sample taken from the individual.

Another object of the present invention is the use of a PPARy inhibitor, especially a PPARy antagonist, for treating patients with a poor prognosis of osteosarcoma.

In particular, the invention provides an in vitro method for determining prognosis of osteosarcoma in an individual, comprising: (a) measuring expression levels of a collection of signature genes from a biological sample taken from said individual, wherein said collection of signature genes comprises at least two genes selected from the group consisting of: AMZ2P1 , ASXL2, BIRC6, C11orf58, C1orf53, C5orf28, CCDC34, DCUN1 D2, DNMT3A, ECI1 , EIF3M, ESCO1 , FAM35DP, GAGE12D, GAGE12G, HADHA, HIST2H2AA3, HIST2H2AA4, MAPK1 , MEAF6, MYO1 B, MYOM3, NDP, PAIP1 , PCID2, PEG10, PEX26, POLR2C, PRR14L, RALGAPA2, SLC2A6, SLC7A4, ST7L, THAP9-AS1 , TIMP3, TTPAL, WNT4; (b) applying the expression levels measured in step (a) to a predictive model relating expression levels of said collection of signature genes with osteosarcoma outcome; and (c) evaluating the output of said predictive model to determine prognosis of osteosarcoma in said individual. In some embodiments, the collection of signature genes comprises CCDC34 and MEAF6. In some embodiments, the collection of signature genes comprises at least 5, preferably at least 7, more preferably at least 8 genes selected from the group consisting of: AMZ2P1 , ASXL2, BIRC6, C11orf58, C1orf53, C5orf28, CCDC34, DCUN1 D2, DNMT3A, ECI1 , EIF3M, ESCO1 , FAM35DP, GAGE12D, GAGE12G, HADHA, HIST2H2AA3, HIST2H2AA4, MAPK1 , MEAF6, MYO1 B, MYOM3, NDP, PAIP1 , PCID2, PEG10, PEX26, POLR2C, PRR14L, RALGAPA2, SLC2A6, SLC7A4, ST7L, THAP9-AS1 , TIMP3, TTPAL, WNT4. In some embodiments, the collection of signature genes comprises AMZ2P1 , C5orf28, CCDC34, GAGE12D, MEAF6, SLC2A6, SLC7A4 and THAP9-AS1. In some embodiments, the collection of signature genes comprises at least 10, preferably at least 12, more preferably at least 15 genes selected from the group consisting of: AMZ2P1 , ASXL2, BIRC6, C11orf58, C1orf53, C5orf28, CCDC34, DCUN1 D2, DNMT3A, ECI1 , EIF3M, ESCO1 , FAM35DP, GAGE12D, GAGE12G, HADHA, HIST2H2AA3, HIST2H2AA4, MAPK1 , MEAF6, MYO1 B, MYOM3, NDP, PAIP1 , PCID2, PEG10, PEX26, POLR2C, PRR14L, RALGAPA2, SLC2A6, SLC7A4, ST7L, THAP9-AS1 , TIMP3, TTPAL, WNT4. In some embodiments, the collection of signature genes comprises AMZ2P1 , C5orf28, CCDC34, ESCO1 , GAGE12D, HADHA, MEAF6, RALGAPA2, SLC2A6, SLC7A4, THAP9-AS1 and TIMP3. In some embodiments, the collection of signature genes comprises AMZ2P1 , ASXL2, C11orf58, C5orf28, CCDC34, ESCO1 , GAGE12D, HADHA, MEAF6, PAIP1 , RALGAPA2, SLC2A6, SLC7A4, THAP9-AS1 and TIMP3. In some embodiments, the collection of signature genes comprises AMZ2P1 , C5orf28, CCDC34, ESCO1 , FAM35DP, GAGE12D, HADHA, HIST2H2AA3, HIST2H2AA4, MAPK1 , MEAF6, POLR2C, RALGAPA2, SLC2A6, SLC7A4, ST7L, THAP9-AS1 , TIMP3 and TTPAL In some embodiments, the expression levels of said collection of signature genes are measured at diagnosis. In some embodiments, the expression levels of said collection of signature genes are measured at relapse.

In some embodiments, the gene expression levels of said signature genes are combined with one or more other parameters to predict progression of osteosarcoma in said individual. In some embodiments, the one or more other parameters are selected from the group consisting of drugs administered to the patient, age, tumour height and stage at diagnosis, presence or absence of metastasis at diagnosis and any combinations thereof.

In some embodiments, the biological sample is an osteosarcoma biopsy sample from the individual.

In some embodiments, the method further comprises the development of said predictive model using stability selection. In some embodiments, the method further comprises developing said predictive model using logistic regression. In some embodiments, the method further comprises developing said predictive model by selecting genes using stability selection with elastic-net regularized logistic regression.

The invention also pertains to a diagnostic kit for predicting progression of osteosarcoma in a subject, wherein said kit comprises at least one nucleic acid probe or oligonucleotide, which can be used in a method for measuring the level of expression of a collection of signature genes from a biological sample taken from said individual, wherein said collection of signature genes comprises at least two genes selected from the group consisting of: AMZ2P1 , ASXL2, BIRC6, C11orf58, C1orf53, C5orf28, CCDC34, DCUN1 D2, DNMT3A, ECI1 , EIF3M, ESCO1 , FAM35DP, GAGE12D, HADHA, HIST2H2AA3, HIST2H2AA4, MAPK1 , MEAF6, MYO1 B, MYOM3, NDP, PAIP1 , PCID2, PEG10, PEX26, POLR2C, PRR14L, RALGAPA2, SLC2A6, SLC7A4, ST7L, THAP9- AS1 , TIMP3, TTPAL, WNT4.

Another aspect of the present invention is the use of a PPARy inhibitor, especially a PPARy antagonist, as an antineoplastic treatment in a subject with osteosarcoma. In some embodiments, said subject is an adolescent or a young adult. In some embodiments, said subject is identified as a poor responder to chemotherapy by the method described above. In some embodiments, said inhibitor is used in combination with another antineoplastic treatments. In some embodiments, said PPARy inhibitor is T0070907 (CAS N°313516-66-4). BRIEF DESCRIPTION OF THE DRAWINGS

Figure 1 : A. Overall survival curves for G1 (pink) and G2 (blue) patients.

B. Progression free survival curves for G1 (pink) and G2 (blue) patients. C. Barplot illustrating the proportion of models including clinical variables and G1/G2 stratification for the prediction of the overall survival of patients of the cohort. D. Principal component analysis of independent component metagene vectors. Name of the ICs for the most contributive loading variables 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. E. Loadings of the most contributive ICs as estimated by partial least square discriminant analysis to stratify patients between G1 or G2.

Figure 2: Gene expression distribution for each of the subnetworks and their predicted biological function. Adjusted P-value computed by gene set enrichment analysis using the Iog2 fold change of genes from the network (G1 versus G2) to detect significant enrichment in subnetwork gene sets.

Figure 3: A. GISTIC analysis of the copy number alterations found significantly associated to G1 or G2 tumors. B. Overall survival curves of predicted G1 (pink) and G2 (blue) patients from an independent cohort of 82 osteosarcoma patients.

C. Overall survival curves of predicted G1 (pink) and G2 (blue) patients using custom Nanostring panel on 96 samples from OS2006 cohort D. Barplot on the left shows proportion of G1 and G2 tumors among 39 out of 79 patients enrolled in the OS2006 trial, with tumor submitted to RNA-sequencing and which relapsed. Barplot on the right shows proportion of G1 and G2 tumors sampled at relapse among 42 patients with osteosarcoma enrolled in MAPPYACTS trial.

Figure 4: A. Cell viability assessment of osteosarcoma cell lines under the pressure of a PPAR antagonist (T0070907 - round points) and two PPAR agonists (TGZ - squares or RGZ - triangles) at 48 and 72h, with concentrations of [0; 0.01 ; 0.1 ; 1 ; 5; 10; 17.5; 20; 25; 30; 50 or 100 pmol/L] by MTS assay. At least three independent tests were performed, with three replicates each. B. The protein expression of PPARy was measured by Western Blot in osteosarcoma cell lines. Cells from normal lung tissue were used as a positive control for PPARy expression and p-actin was used as a housekeeping gene control.

Figure 5: Scatter plot between each IC sample contribution and the log fold change of the CNA showing the greatest coefficient in the corresponding regression model. Figure 6: A. Barplot illustrating gene selection based on their correlation between RNA-seq and Nanostring expression. Blue bars show genes correlating with p-values lower than 10' ⁵ and selected for the Nanostring signature. The five genes in italic/bold corresponds to the housekeeping genes, chosen based on their low variance in RNA-seq, to normalized Nanostring data B. Progression free survival curves of predicted G1 (pink) and G2 (blue) patients using custom Nanostring panel on 96 samples from OS2006 cohort.

Figure 7: Differentially expressed genes (p-adj<0.01) in HOS cells after treatment with antagonist or agonist of PPARG.

Figure 8: PPARy expression in osteosarcoma PDXs by immunohistochemistry (IHC)

Figure 9: In vivo effects of T0070907 in osteosarcoma paratibial PDX model MAP-217 alone (upper panel) and combined to methotrexate (lower panel), on primary tumor (left) and lung metastasis (right). A. Delay (days) to reach 2.5 time the initial tumor size in control and T0070907 treated mice. B. Metastases evaluation by IHC (Humain KI67 staining) in control and T0070907 treated mice. C. Tumor growth curve of the primary tumor for single agents and combination.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

Unless otherwise indicated, the practice of the method and system disclosed herein involves conventional techniques and apparatus commonly used in molecular biology, microbiology, protein purification, protein engineering, protein and DNA sequencing, and recombinant DNA fields, which are within the skill of the art. Such techniques and apparatus are known to those of skill in the art and are described in numerous texts and reference works.

Unless defined otherwise herein, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art. Various scientific dictionaries that include the terms included herein are well known and available to those in the art. Although any methods and materials similar or equivalent to those described herein find use in the practice or testing of the embodiments disclosed herein, some methods and materials are described. The terms defined immediately below are more fully described by reference to the Specification as a whole. It is to be understood that this disclosure is not limited to the particular methodology, protocols, and reagents described, as these may vary, depending upon the context they are used by those of skill in the art. In the present text, the following general definitions are used:

As used herein, the singular terms "a," "an," and "the" include the plural reference unless the context clearly indicates otherwise.

"Nucleic acid sequence," "expressed nucleic acid," or grammatical equivalents thereof used in the context of a corresponding signature gene means a nucleic acid sequence whose amount is measured as an indication of the level of expression of genes. The nucleic sequence can be a portion of a gene, a regulatory sequence, genomic DNA, cDNA, RNA including mRNA and rRNA, or others. A particular embodiment utilizes mRNA as the primary target sequence. As is outlined herein, the nucleic acid sequence can be a sequence from a sample, or a secondary target such as, for example, a product of a reaction such as a PCR amplification product (e.g., "amplicon"). A nucleic acid sequence corresponding to a signature gene can be any length, with the understanding that longer sequences are more specific. Probes are made to hybridize to nucleic acid sequences to determine the presence or absence of expression of a signature gene in a sample.

As used herein the term "comprising" means that the named elements are included, but other element (e.g., unnamed signature genes) may be added and still represent a composition or method within the scope of the claim.

As used herein, the term "signature gene" refers to a gene whose expression is correlated, either positively or negatively, with disease extent or outcome or with another predictor of disease extent or outcome. A "signature nucleic acid" is a nucleic acid comprising or corresponding to, in case of cDNA, the complete or partial sequence of a RNA transcript encoded by a signature gene, or the complement of such complete or partial sequence. A signature protein is encoded by or corresponding to a signature gene of the disclosure.

The term "relapse prediction" is used herein to refer to the prediction of the likelihood of osteosarcoma recurrence in patients with no apparent residual tumor tissue after treatment. The predictive methods of the present disclosure can be used clinically to make treatment decisions by choosing the most appropriate treatment modalities for any particular patient. The predictive methods of the present disclosure also can provide valuable tools in predicting if a patient is likely to respond favorably to a treatment regimen, such as surgical intervention, chemotherapy with a given drug or drug combination, and/or radiation therapy.

The terms "subject", “individual” or “patient” herein refer to a human subject.

“Osteosarcoma” herein designates a type of bone cancer that begins in the cells that form bones. Osteosarcoma is most often found in the long bones - more often the legs, but sometimes the arms - but it can start in any bone. In very rare instances, it occurs in soft tissue outside the bone.

The term “antineoplastic treatments” herein designate any treatment for cancer except surgery. They include chemotherapy, hormonal and biological therapies, and radiotherapy.

As used herein, “treat”, “treatment” and “treating” refer to any reduction or amelioration of the progression, severity, and/or duration of cancer, particularly a solid tumor; for example in an osteosarcoma, reduction of one or more symptoms thereof that results from the administration of one or more therapies.

Other definitions will be specified below, when necessary.

According to a first aspect, the present invention pertains to an in vitro method for determining prognosis of osteosarcoma in an subject, comprising the steps of: a. measuring expression levels of a collection of signature genes from a biological sample taken from said subject, wherein said collection of signature genes comprises at least two genes selected from the group consisting of: AMZ2P1 , ASXL2, BIRC6, C11orf58, C1orf53, C5orf28, CCDC34, DCUN1 D2, DNMT3A, ECI1 , EIF3M, ESCO1 , FAM35DP, GAGE12D, GAGE12G, HADHA, HIST2H2AA3, HIST2H2AA4, MAPK1 , MEAF6, MYO1 B, MYOM3, NDP, PAIP1 , PCID2, PEG10, PEX26, POLR2C, PRR14L, RALGAPA2, SLC2A6, SLC7A4, ST7L, THAP9-AS1 , TIMP3, TTPAL, WNT4; and b. applying the expression levels measured in step (a) to a predictive model relating expression levels of said collection of signature genes with osteosarcoma outcome; and c. evaluating the output of said predictive model to determine prognosis of osteosarcoma in said subject.

In the above method, the predictive model can have been developed by any appropriate method known by the skilled person. Machine learning (ML) approaches have been demonstrated to be useful in medicine and can advantageously be used to build the predictive method. An example of such a method is disclosed in the experimental part below. Other examples of computational protocols for identifying and refining prognostic signatures have been described and can be used by the skilled person in the context of the present invention ^66-69 (Vey et al., Cancers, 2019; Xia et al., Nature Communications, 2019; Jiang et al., Cell Systems, 2018; Liu et al, 2020; EBioMedicine). The predictive model is advantageously developed using data collected from patients known to have osteosarcoma. The skilled person can obtain different prediction models corresponding to different clinical situations. For example, a model can be developed using data collected only from patients at diagnosis; another model can be developed using data from patients who underwent surgery of the primary tumors, etc..

In some embodiments, the gene expression data may be preprocessed by normalization, background correction and/or batch effect correction. The preprocessed data may then be analyzed for differential expression of genes for stratification of patients in a good prognostic group (group 1) versus a poor prognostic group (group 2).

In some embodiments, to develop a predictive model for the prognosis of osteosarcoma, the probes included in the final model are selected from an entire set of probes using stability selection. In some embodiments, the model is developed using logistic regression, for example using elastic-net regularized logistic regression. Elastic- net regression is a high dimensional regression method that incorporates both a LASSO (L1) and a ridge regression (L2) regularization penalty. The exact mix of penalties (LASSO vs. ridge) is controlled by a parameter a (a = 0 is pure ridge regression and a = 1 is pure LASSO). The degree of regularization is controlled by the single penalty parameter. Both LASSO and ridge regression shrink the model coefficients toward zero relative to unpenalized regression but LASSO can shrink coefficients to exactly zero, thus effectively performing variable selection. LASSO alone however, tends to select randomly among correlated predictors, which the addition of the ridge penalty helps prevent. It is known that the ridge penalty shrinks the coefficients of correlated predictors towards each other while the lasso tends to pick one of them and discard the others. In some embodiments, one may use the implementation of elastic-net logistic regression in the R package 'glmnet. '

The idea behind stability selection is to find 'stable' probes that consistently show to be predictive of recurrence across multiple data sets obtained by 'perturbing' the original data. Specifically, perturbed versions of the data are obtained by subsampling m < n subjects (n is the total number of subjects) without replacement. Regularized regression (or elastic-net in some embodiments) is then performed on each subsample version of the data to obtain the complete regularized path (/.e., the model coefficients as a function of the regularization penalty). The effect of the LASSO penalty is to shrink the vast majority of the probe coefficients to exactly zero; the probes with non-zero coefficients (predictive) across a sizable proportion of the subsample versions of the data are deemed stable predictors.

In some embodiments, to implement stability selection with elastic net regression, one may calibrate the tuning parameter a using repeated cross-validation (e.g., using R package for a 10-fold cross-validation). The skilled person will choose the tuning parameter a to provide good prediction and to include as many possible features as necessary while maintaining good prediction. In some embodiments, stability selection may be implemented using different numbers of subsamples of the data, each having a portion of the total sample size (each with roughly the same proportion of cases and controls as the original), in order to identify robust predictors for the final model. In some embodiments, standardization of the gene expression levels by their standard deviation (the default in glmnet to place all gene features on the same scale) is not done, since differential variability of the gene expression levels may be biologically important. In some embodiments, such standardization may be performed. In some embodiments, clinical variables such as, for example, the presence of metastases at diagnosis are force included (i.e., not subject to the elastic net regularization penalty).

According to a particular embodiment illustrated in the experimental part, the predictive model leads to a stratification of patients into 2 groups: G1 of good prognosis (OS at 3 years = 100%) and G2 of poorer prognosis (OS at 3 years <70%), requiring a closer monitoring.

Patients classified in the G1 group may benefit from a treatment with an immunomodulator (e.g., mifamurtide), while patients in the G2 group may benefit from a treatment with antiangiogenesis agents like multityrosine kinase inhibitors or monoclonal antibodies, anti-osteoclasts like biphosphonates oranti-RANKL, as well as PPARgamma antagonists (as illustrated in Example 2).

Tables 4 to 8 of the experimental part below provide parameters for 5 different models which can be used to perform the above method.

In a preferred embodiment, the collection (or panel) of signature genes includes at least CCDC34 and MEAF6. Table 7 shown in the experimental part provides parameters for a model based on the expression levels of these 2 genes. Of course, these are provided as an example and not limiting. The skilled person can refine these parameters using, for example, a different cohort of patients, or adapt them to values obtained with a different technique used for measuring the expression levels of the genes. This of course holds true for all the models provided in this application. In some embodiments, 2, 3, 4, 5, 6, 7, 8, 9, 10, II, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21 , 22, 23, 24, 25, 26, 27, 28, 29, 30, 31 , 32, 33, 34, 35, 36 or all of the 37 genes shown in Table 4 may be used in a predictive model. For illustrative purpose only, Table 4 shown in the experimental part provides parameters for a model based on the expression levels of these 37 genes.

In some embodiments, at least 5, preferably at least 7 and more preferably at least 8 genes of the list of Table 4 are included. In some embodiments, at least 10, preferably at least 12 and more preferably at least 15 genes of the list of Table 4 are included. In some embodiments, these two or more genes may be selected by their correlation with recurrence in the training data set to develop the predictive models. In some embodiments, the one or more genes may be selected by their reliability ranks. In some embodiments, the two or more genes may be selected by their predictive power rankings.

In some embodiments, the panel of signature genes comprises at least AMZ2P1 , C5orf28, CCDC34, GAGE12D, MEAF6, SLC2A6, SLC7A4 and THAP9-AS1. For illustrative purpose only, Table 6 shown in the experimental part provides parameters for a model based on the expression levels of these 8 genes.

In some embodiments, the panel of signature genes comprises at least AMZ2P1 , ASXL2, C11orf58, C5orf28, CCDC34, ESCO1 , GAGE12D, HADHA, MEAF6, PAIP1 , RALGAPA2, SLC2A6, SLC7A4, THAP9-AS1 and TIMP3. For illustrative purpose only, Table 5 shown in the experimental part provides parameters for a model based on the expression levels of these 15 genes.

In some embodiments, the panel of signature genes comprises at least AMZ2P1 , C5orf28, CCDC34, ESCO1 , FAM35DP, GAGE12D, HADHA, HIST2H2AA3, HIST2H2AA4, MAPK1 , MEAF6, POLR2C, RALGAPA2, SLC2A6, SLC7A4, ST7L, THAP9-AS1 , TIMP3 and TTPAL. For illustrative purpose only, Table 8 shown in the experimental part provides parameters for a model based on the expression levels of these 20 genes.

In some embodiments, the panel of signature genes comprises at least AMZ2P1 , C5orf28, CCDC34, ESCO1 , GAGE12D, HADHA, MEAF6, RALGAPA2, SLC2A6, SLC7A4, THAP9-AS1 and TIMP3. These 12 genes are present in the models shown in both Table 5 and Table 8.

The method of the invention, in any of its embodiments described above, can further comprise combining the gene expression levels of said signature genes with one or more other parameters to predict progression of osteosarcoma in said subject. Non-limitative examples of such clinical parameters which can be combined to the expression levels of the selected genes include the presence and amount of metastases at diagnosis, the size, height and/or stage of the tumor at diagnosis, drugs administered to the patient, age, and any combinations thereof.

In some embodiments, based on the expression levels for the set of probes determined by stability selection and the clinical variables, the predictive model is obtained by fitting logistic model using elastic-net regularized logistic regression. Glmnet solves the following problem over a grid of values of A covering the entire range. The tuning parameter A controls the overall strength of the penalty.

In some embodiments, instead of selecting a subset of genes from the panel disclosed in Table 4, a model may weight the genes differently in the logistic regression. In some embodiments, determining the prognosis of osteosarcoma for an individual involves applying expression levels of the collection of signature genes to the predictive model, which involves weighting said expression levels according to stability rankings of the collection of signature genes. In some embodiments, the method involves weighting expression levels according to predictive power rankings of the collection of signature genes.

The logistic regression model above expresses the specific way the expression levels and the clinical variables are combined to obtain a score for each individual. In some embodiments, expression levels are weighted in the elastic-net regularized logistic regression. In some embodiments, expression levels are weighted in using LASSO. The weighting here does not refer to the model coefficients (which can be thought of as weights for the expression levels and clinical variables), but rather to an additional mechanism for differentially accounting for variable importance in the logistic regression procedure. In this regard, alternative embodiments consider unweighted logistic regression, /.e., treating all genes equally, and weighted logistic regression, weighting by the stability selection frequencies.

In some embodiments, various clinical variables (e.g., drugs administered to the patient, age, tumour height and stage at diagnosis, presence of absence of metastasis at diagnosis) will be included in the same logistic model along with the signature genes. Coefficients will be defined for each variable (gene expression and clinical values). This logistic regression model will provide a probability of having a clinical recurrence given the provided gene expression scores and clinical variables. This probability will be a number between 0-1 , and it will indicate for each given patient the prognosis of the disease.

In some embodiments, in addition to identifying the coefficients of the predictive model, the disclosure identifies the most useful specificity and sensitivity a user wishes to have for a specific risk probability. Based on the desired specificity and sensitivity levels, the method will report the risk status of each patient. For example, the skilled person may find that given the specificity and sensitivity of our model, a patient with 45% chance of being classified in the good prognosis responders group (group 1) might be better of being classified in the poor progbosis responders group (group 2) rather than group 1 or vice versa. In other words, more user-friendly criteria can be chosen based on more detailed analysis in further datasets to determine the most practical interpretation of the risk probability depending on how much clinicians want to risk having a false positive or a false negative.

Individuals suspected of having any of a variety of bone cancers, such as osteosarcoma, osteocarcinoma or cancer of the bones joints and soft tissue as well as pediatric sarcomas, can be evaluated using a method of the disclosure. Exemplary cancers that can be evaluated using a method of the disclosure include, but are not limited to osteosarcoma, pediatric sarcoma a well as adult and pediatric relapsing cancers in bone tissue.

Exemplary clinical outcomes that can be determined from a model of the disclosure include, for example, response to a particular course of therapy such as surgical removal of a tumor, radiation, or chemotherapy.

The skilled person will appreciate that patient tissue samples containing osteosarcoma cells may be used in the methods of the present disclosure including, but not limited to those aimed at determining the prognosis of the disease. In these embodiments, the level of expression of the signature gene can be assessed by assessing the amount, e.g. absolute amount or concentration, of a signature gene product, e.g., protein and RNA transcript encoded by the signature gene and fragments of the protein and RNA transcript) in a sample obtained from a patient. The sample can, of course, be subjected to a variety of well-known post-collection preparative and storage techniques (e.g. fixation, storage, freezing, lysis, homogenization, DNA or RNA extraction, ultrafiltration, concentration, evaporation, centrifugation, etc.) prior to assessing the amount of the signature gene product in the sample.

Tissue samples useful for preparing a model for determining the prognosis of osteosarcoma include, for example, paraffin and polymer embedded samples, ethanol embedded samples and/or formalin and formaldehyde embedded tissues, although any suitable sample may be used. In general, nucleic acids isolated from archived samples can be highly degraded and the quality of nucleic preparation can depend on several factors, including the sample shelf life, fixation technique and isolation method.

If required, a nucleic acid sample having the signature gene sequence(s) are prepared using known techniques. For example, the sample can be treated to lyse the cells, using known lysis buffers, sonication, electroporation, etc., with purification and amplification as outlined below occurring as needed, as will be appreciated by those in the art. In addition, the reactions can be accomplished in a variety of ways, as will be appreciated by those in the art. Components of the reaction may be added simultaneously, or sequentially, in any order, with preferred embodiments outlined below. In addition, the reaction can include a variety of other reagents which can be useful in the assays. These include reagents like salts, buffers, neutral proteins, e.g. albumin, detergents, etc., which may be used to facilitate optimal hybridization and detection, and/or reduce non-specific or background interactions. Also reagents that otherwise improve the efficiency of the assay, such as protease inhibitors, nuclease inhibitors, anti-microbial agents, etc., can be used, depending on the sample preparation methods and purity.

In some embodiments, biological samples in addition to or instead of osteosarcoma tissue may be used to determine the expression levels of the signature genes. In some embodiments, the suitable biological samples include, but are not limited to, circulating tumor cells isolated from the blood, urine of the patients or other body fluids, exosomes, and circulating tumor nucleic acids.

According to a particular embodiment, the expression levels of said collection of signature genes are measured at diagnosis.

According to another particular embodiment, the expression levels of said collection of signature genes are measured at relapse.

In some embodiments, the gene expression levels of the signature genes may be measured multiple times. In some embodiments, the dynamics of the expression levels may be used in combination of the signature genes expression levels to better predict the clinical outcome. One skilled in the art understands various approaches may be used to combine the effects of the levels and the dynamics of the signature genes' expression to determine the prognosis of osteosarcoma.

The methods of the disclosure depend on the detection of differentially expressed genes for expression profiling across heterogeneous tissues. Thus, the methods depend on profiling genes whose expression in certain tissues is activated to a higher or lower level in an individual afflicted with a condition, for example, cancer, such as osteosarcoma, relative to its expression in a non-cancerous tissues or in a control subject. Gene expression can be activated to a higher or lower level at different stages of the same conditions and a differentially expressed gene can be either activated or inhibited at the nucleic acid level or protein level.

Differential signature gene expression can be identified, or confirmed using methods known in the art such as qRT-PCR (quantitative reverse-transcription polymerase chain reaction) and microarray analysis. In particular embodiments, differential signature gene expression can be identified, or confirmed using microarray techniques or any similar technique (e.g., NanoString technique). Thus, the signature genes can be measured in either fresh or paraffin-embedded tumor tissue, using microarray technology. In this method, polynucleotide sequences of interest are plated, or arrayed, on a microchip substrate. The arrayed sequences are then hybridized with specific DNA probes from cells or tissues of interest.

The expression level of a signature gene in a tissue sample can be determined by contacting nucleic acid molecules derived from the tissue sample with a set of probes under conditions where perfectly complementary probes form a hybridization complex with the nucleic acid sequences corresponding to the signature genes, each of the probes including at least two universal priming sites and a signature gene target-specific sequence; amplifying the probes forming the hybridization complexes to produce amplicons; and detecting the amplicons, wherein the detection of the amplicons indicates the presence of the nucleic acid sequences corresponding to the signature gene in the tissue sample; and determining the expression level of the signature gene.

The expression level of nucleic acid sequences corresponding to a set of signature genes in a tissue sample can be determined by contacting nucleic acid molecules derived from the tissue sample with a set of probes under conditions where complementary probes form a hybridization complex with the signature gene-specific nucleic acid sequences, each of the probes including at least two universal priming sites and a signature gene-specific nucleic acid sequence; amplifying the probes forming the hybridization complexes to produce amplicons; detecting the amplicons, wherein the detection of the amplicons indicates the presence of the nucleic acid sequences corresponding to the set of signature genes in the tissue sample; and determining the expression level of the target sequences, wherein the expression of at least two, at least three, at least five signature gene-specific sequences is detected.

The present invention also pertains to a collection of isolated probes specific for osteosarcoma signature genes comprising at least two genes selected from the group consisting of AMZ2P1 , ASXL2, BIRC6, C11orf58, C1orf53, C5orf28, CCDC34, DCUN1 D2, DNMT3A, ECI1 , EIF3M, ESCO1 , FAM35DP, GAGE12D, GAGE12G, HADHA, HIST2H2AA3, HIST2H2AA4, MAPK1 , MEAF6, MYO1 B, MYOM3, NDP, PAIP1 , PCID2, PEG10, PEX26, POLR2C, PRR14L, RALGAPA2, SLC2A6, SLC7A4, ST7L, THAP9-AS1 , TIMP3, TTPAL, WNT4. In particular, the invention provides, for each of the collection of signature genes mentioned above, an appropriate collection of probes. Examples of such probes are provided in Table 3 below.

The invention includes compositions, kits, and methods for determining the prognosis of osteosarcoma for an individual from which a sample is obtained. A kit is any manufacture (e.g. a package or container) including at least one reagent, e.g. a probe or a collection of probes as above-described (e.g., nucleic acid probes of SEQ ID NO: 1-41 or a subset thereof), for specifically detecting the expression of a gene signature as described herein.

As shown in Example 2 below, the inventors demonstrated that a PPARy antagonist decreases the cell viability in 7 human osteosarcoma lines. This is a proof of principle that inhibitors of the PPARy pathway may constitute a breakthrough innovation in the armamentarium for treating osteosarcoma.

In another of its aspects, the present invention thus pertains to the use of a PPARy inhibitor as an antineoplastic treatment in a subject with osteosarcoma, especially for treating an adolescent or a young adult.

Examples of PPARy inhibitors which can be used according to the present invention include T0070907 (CAS N°313516-66-4), BADGE, FH 535, GW 9662, SR 16832, SR 202 and combinations thereof.

According to a particular embodiment, a PPARy inhibitor is used for treating a patient who has been identified as a poor responder to chemotherapy by a method of the invention, as above-described. In a particular embodiment, the PPARy inhibitor can be combined to another treatment, such as a treatment with antiangiogenesis agents (e.g., multityrosine kinase inhibitors and monoclonal antibodies), anti-osteoclasts (e.g., biphosphonates and anti-RANKL), etc.

Other characteristics of the invention will also become apparent in the course of the description which follows of the biological assays which have been performed in the framework of the invention and which provide it with the required experimental support, without limiting its scope.

EXAMPLES

The inventors reasoned that in osteosarcoma, gene expression at RNA level, in addition to giving access to 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 than phenotypic features.

In order to test this hypothesis, we performed RNA-sequencing (RNA- seq) of 79 primary osteosarcoma tumors sampled at diagnosis from patients accrued in and homogenously treated by the first line OS2006 trial (NCT00470223). These samples recapitulate the main clinical feature distribution of the whole cohort (Table 1). To study the transcriptome landscape of both osteosarcoma and TME cells, we first dissociated their respective transcriptome programs. To do so, we decomposed RNA gene expression matrix by independent component analysis (ICA, BIODICA implementation) ²³ with ICASSO stabilization procedure ²³ into 50 independent components (ICs), also called gene modules.

Table 1 : Patient primary tumor biopsy sample cohort explored by

RNA-seq analysis in OS2006 trial. Patients and tumor characteristics, treatment and outcome of the 79 patients with RNA-seq analysis in the first-line phase-ill therapeutic OS2006 trial.

Example 1: Immune infiltrate and tumor microenvironment transcriptional programs stratify pediatric osteosarcoma into prognostic groups at diagnosis

Materials and Methods

Population and Materials

Biological samples were prospectively collected for patients (up to 50 years) registered into the therapeutic approved French OS2006/sarcoma09 trial (NCT00470223) ⁶ . This study was carried out in accordance with the ethical principles of the Declaration of Helsinki and with Good Clinical Practice guidelines. A specific informed consent for blood and tumor samples was obtained from patients or their parents/guardians if patients were under 18 years of age upon enrolment. The information given to the patients, the written consent used, the collection of samples and the research project were approved by an independent ethic committee and institutional review boards. As part of the ancillary biological studies, RNA-seq, CGH array and WES were performed at 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 (Table 1) were comparable to the whole osteosarcoma population registered in the OS2006 trial. The sex ratio M/F was 0.55 and the median age of 15.43 years (4.71-36.83). Primary tumors were located in the limb (94.9%) and 20.2% had metastases at diagnosis. Chemotherapy was administered as per OS2006 trial; either MTX-etoposide/ifosfamide (M-EI) regimen in 86% of the patients or adriamycin/platinum/ifosfamide (API-AI) regimen for 12.6% of patients, and 34.1% received zoledronic acid according to the randomization arm. 77 patients had surgery of the primary tumor and poor histological response was observed in 24 % of the patients. 29 patients relapsed with a median delay of 1.55 years (range 0.09.-3.62). 22 patients died with a median delay of 2.25 years (range 0.07. -6.9). The median follow-up was of 4.8 years. Overall, the 3-year PFS and OS were of 65.82% and 82.27%, respectively.

Nucleic acids extraction

DNA and RNA were isolated using AHPrep DNA/RNA mini kit (Qiagen, Courtaboeuf, France) according to manufacturer's instructions. Quantification/qualification were performed using Nanodrop 2000 spectrophotometer (Thermo Fisher Scientific, lllkirch, France) and Bioanalyzer DNA 7500 (Agilent Technologies, Les Ulis, France).

RNA sequencing

RNA sequencing libraries were prepared with TrueSeq Stranded mRNA kit following recommendations: the key steps consist of PolyA mRNA capture with oligo dT beads using 1 pg total RNA, fragmentation to approximately 400pb, DNA double strand synthesis, and ligation of Illumina adaptors amplification of the library by PCR for sequencing. Libraries 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 were 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 in RNA-seq

Quality of stranded pair-ended RNA-seq libraries was evaluated with fastqc (a quality control tool for high throughput sequence data available on the Babraham bioinformatics pages of the website of the Babraham institute). Reads were mapped with RSEM using GRCh37 ENSEMBL mRNA dataset as reference sequences.

ICA was performed using BIODICA, one of the most performant implementations of fastICA including icasso stability analysis. Genes represented by less than 100 reads were filtered out from the ICA and gene expression matrix was log scaled and scaled by genes before decomposition in 50 components. Pearson correlations between ICs, to assess inter IC relationship, were calculated from the metagenes matrix. For each IC we defined contributive genes, as genes within 3 standard deviations of the mean of the corresponding IC metagene vector. For each component, gene set enrichment was performed using msigdb for positional enrichment and g:Profiler for functional enrichment. For each component, only the most significant enrichment detected among the GO, REACTOME or TF categories has been reported. For positional enrichment, only the cytobands with a minimum FDR of 1e ^-4 have been reported. When several cytobands from the same chromosomal arm were found enriched only the less significant p-values have been reported in Table 2. From the gene expression matrix, we selected only the contributive genes, within 2.5 standard deviations of the mean of the corresponding IC metagenes vector, to perform a network inference using Pearson distance to calculate edge lengths. Edges with distance above 0.25 were discarded. Network and subnetwork illustrations and analyses have been produced with Cytoscape v3.6.1. We used the REACTOME Fl viz plugin to detect clusters in the network with the spectral partition based network clustering algorithm and estimate for each cluster their functional enrichments in GO biological process terms.

Hierarchical classification of metagene matrix and corresponding heatmap was generated with pheatmap R package using Pearson distance and Ward’s construction method. Pls-da function from the mixOmics R package was used to select the most contributive independent components to the G1/G2 patient classification.

Kaplan-Meier OS and PFS curves were generated using the survminer and survival R packages. Log-rank test were used to compare survival distribution inbetween groups and to provide corresponding p-values.

Table 2: Functionnal annotation of the independent components by functional enrichment based on MSigDB and g: Profiler databases Differential mRNA expression was estimated with DESeq2 R package from 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 CNAs (inferred from CGH array in pair with RNA-seq samples) modeling the best each IC. Likewise, we detected dependencies between ICs corresponding to large chromosomal transcriptional modulation and copy number alterations from the same chromosomal region. Chromosomal region written in bold letters in Table 2 corresponds to the region, enriched in ICs, overlapping with the CNA contributing the most to the model, which suggest a dose-dependent relationship. Figure 5 illustrates those cases by the scatter plot of IC signal vs copy number log fold change of the copy number alteration with the highest coefficient in the regression model.

A logistic regression model using the glmnet R package was used to define a minimum gene signature able to discriminate G1 from G2 RNA-sequencing libraries (parameters: type.measure="mse", alpha=0.35, family- 'binomial"). Cross- validation strategy was used to select the best lambda. We then predicted with this signature G1 and G2 tumors from the 82 osteosarcoma RNA-seq dataset generated by the TARGET consortium as well as 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 OS and PFS curves were generated using the survminer and survival R packages. Log-rank test were used to compare survival distribution in-between groups and to provide corresponding p-values.

Oligonucleotide Array Comparative Genomic Hybridization (aCGH) assay.

In all experiments, sex-matched normal DNA from a pooled human female or male (Promega, Madison, Wl, USA) was used as a reference. Oligonucleotide 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 Alul and Rsal (Fermentas, Euromedex, France). The fragmented DNAs were labelled with cyanine Cy3-deoxyuridine triphosphate (dUTP) or Cy5-dUTP. Hybridization was carried out on SurePrint G3 Human CGH Microarray 4x180K (Agilent Technologies, Santa Clara, CA, USA) arrays for 24 hours at 65°C in a rotating oven (Robbins Scientific, Mountain View, CA) at 20rpm. The 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 pm resolution at 20°C in low ozone concentration environment. Data were extracted from scanned TIFF images using the Feature Extraction software (v11.5.1.1 , Agilent), along with protocol CGH_1105_Oct12. All further data treatments were performed under the R statistical environment in v3.4 (http://cran.r-project.org). Acquired raw intensities were transformed to Iog2(test/ref). Joint normalization of the whole cohort of raw CGH profiles was performed using the ‘cghseg’ package (v1 .0.2.1) using default parameters: a common “wave-effect” track was computed, then subtracted to all individual profiles through a lowess regression. As a second step, an individual normalization step was performed for each profile by subtracting pre-computed GC-content tracks through a lowess regression. Joint segmentation of the whole cohort of normalized CGH profiles was performed using a penalized least-square regression implemented in the ‘copynumber’ package (v1.16.0). To set a unique value to the gamma (penalty) parameter for the whole cohort, optimal gamma value was computed for each individual profile using 100 cross-validation folds, and the median of optimal gammas was chosen (median = 18, sd = 4.2, range = 11 :41). Joint segmentation resulted into 1604 segments. Individual centering of profiles was performed using an in-house method selecting the most-centered mode in the distribution density of probes’ Iog2(ratio) values. No calling of aberrations was performed.

Oligonucleotide aCGH analysis.

All genomic coordinates were established on the LICSC human genome build hg19. Hierarchical clustering of samples was performed under R on segmented data, using [[Euclidean I Pearson I Spearman]] distances and Ward’s construction method. Clustering of samples based on non-negative matrix factorization and spherical k-means were performed using the ‘NMF’ (v0.20.6) and ‘skmeans’ (v0.2.10) packages, respectively. Minimum common regions (MCR) analysis was performed using GISTIC2 (v2.0.22). Differential analyses were performed using non-parametric statistical tests (Wilcoxon for 2-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 association between the overall survival and the major clinical factors (sex, histological response, metastatic status , tumor size, treatment, chemotherapy pubertal status) and the G1/G2 signature, we considered stability selection with boosting using Cox model (Hofner et al., 2015: Controlling false discoveries in high-dimensional situations: boosting with stability selection; BMC Bioinformatics volume 16, Article number: 144).

The selection frequencies to 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 n _th to 0.9 (which can be relaxed, increasing the risk of false positive predictor selection). The per-family error rate (PFER), 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 n _th ) from the definition provided by Meinshausen and Buhlmann (Meinshausen and Buhlmann (2010); Stability selection; Royal Statistical Society 1369- 7412/10/72417J. R. Statist. Soc. B (2010)72, Part 4, pp. 417-473). These analyses were performed using the mboost R package.

Targeted RNA expression analysis

A direct gene expression analysis was conducted with a Nanostring standard custom approach and following the supplier recommendations. A specific panel of 41 probes of interest including five housekeeping genes (CNOT1 , EIF4G2, SF1 , SLC39A1 , and SLIRF4) was designed according to Nanostring. (Table 3 List of targets). RNA integrity was measured with a Fragment Analyzer system (RNA concentration, RNA Quality Number, Percentage of RNA fragments > 300nt), RNA concentrations and purity were controlled with a Nanodrop ND8000.

Table 3: List of targets

Due to the number of targets (<400), 100ng of total RNA were hybridized to Nanostring probes. Our cohort composed of 176 samples was splitted in two batches of hybridizations, including a No Template Control (NTC, water) sample and a Universal RNA sample (Agilent Technologies, P/N: 740000) per batch. Before hybridization, Nanostring positive and negative controls were added to samples as spike in controls. Probes and mRNA were hybridized 16h at 65°C prior processing on the NanoString nCounter preparation station to remove excess of 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. The Nanostring nSolver 4.0 software was used to control raw data, and to normalize data. Imaging QC criteria was higher than 93% (threshold: 75% of FOV). Binding density was monitored as Positive and Negative controls signals. Counts obtained for each target in NTC samples (1-9 counts) were deduced to corresponding target in samples of interest. Geometric average of Housekeeping gene signals was used to normalize RNA content. Comparison of Universal RNA were in the platform agreements (r2: 0.98792). Out of 176 samples, three samples with low performances were rejected from the primary analysis: one sample has a low binding density (lower threshold: 0.1 ; EX148_ARN010: 0.09) and two samples have a high mRNA content normalization factor (threshold: 20; EX148_ARN071 : 25.7 and EX148_ARN165: 870.6).

G1/G2 strata prediction from targeted RNA expression panel

To detect suspicious discrepancies, we first compared the Nanostring and RNA-seq estimation of the RNA expression of the 35 genes as well as the 5 housekeeping genes. Likewise, we filtered out 7 genes with a p-value of correlation (pearson method) higher than 10' ⁵ (Figure not shown, Fig 6A). From the 166 samples, 70 samples, already analyzed by RNA-seq and therefore classified G1 or G2, were chosen as training set to retrain our logistic model with elastic net regularization using glmnet on the 28 targeted RNA expression panel. Alpha was chosen to 0.8 and lambda by a cross validation strategy. Then we predicted the G1/G2 strata on the 96 remaining samples and Kaplan-Meier OS and PFS curves were generated using the survminer and survival R packages. Log-rank test were used to compare survival distribution inbetween groups and to provide corresponding p-values.

Results

1.1. Independent components recapitulate biological functions and their interations

We first tested whether specific ICs were separately correlated to clinical variables. Apart from IC2, associated with sex and obvious gene enrichment from Y chromosome (Figure not shown), no other ICs were significantly associated with clinical items after p-value adjustment.

We next characterized ICs by functional enrichment analysis using the most contributive genes of each IC. We observed 90% of the IC are associated significantly to either large chromosomal region over transcriptional modulation with msigdb database or specific cellular/molecular functions with g:Profiler (Table 2). Several ICs appeared specific to tumor clones with altered gene expression involving statistically significant large chromosomic regions. Some of these regions recapitulate known CNVs observed in osteosarcomas, as confirmed by a comparison to 70 CGH arrays produced in pair with our RNA-seq samples (Table 2, Fig 5). Contributive genes from the remaining 10 ICs recapitulate the transcriptional program of diverse TME cells or their interaction with tumor cells. We identified osteosarcoma functional gene modules belonging to bone microenvironment (osteoclasts/bone resorption IC25; osteoblast/ossification IC44), angiogenesis (IC13), immune response (IC41), but also neuron projection (IC48) and muscle (IC1). Pearson correlation study of the components revealed a proximity between specific ICs with two main groups of components (Figure not shown). Such IC clustering shows that osteosarcoma tumors share two specific IC patterns and also confirms strong relationship across TME, cancer cells and large chromosomal region modulations. To decipher these potential functional interactions, we inferred a co-expression network using only the genes contributing significantly to ICs (>3 standard deviation from the mean). Genes belonging to the same ICs demonstrate strong interconnection in the network, confirming the existence of a continuum between ICs latent variables and real measured gene expression (Figure not shown). In most cases, we identified matches between the ICs and functional subnetwork annotations performed with REACTOME Fl viz (Figure not shown), confirming the potency of ICA to decompose gene expression matrix from heterogeneous tumor into functionally relevant gene modules. In addition, annotation of the co-expression network helped us to classify genes from ICs into more precise subnetworks. IC41 , initially flagged as general immune response, has been subdivided into “adaptive immune response/neutrophil degranulation”, “Type 1 interferon response” and “adaptive immune response”. IC9 has been subdivided into “neutrophil degranulation/defensin” and “adipogenesis/peroxisome proliferator-activated receptor (PPAR) signaling pathway”. At a larger scale, the network is structured around several hubs of highly connected ICs/subnetworks, probably mirroring the active crosstalk between TME and tumor cells.

1.2. Stable independent components stratify patients by prognosis

We questioned whether ICs/subnetwork interconnections reflecting the TME composition were linked to clinical variables. We categorized tumors by hierarchical classification using stable ICs (Stability index>0.5) reflecting clinical annotations (Figure not shown). This unsupervised analysis defined two groups of tumors associated to different average “living status” and discriminating the cohort into 35 low risk tumors, coined G1 (8.5% Death Rate, DR), from a higher risk group G2 (43.1% DR; Figure not shown). We confirmed the significant different overall survival (OS), estimated by the Kaplan-Meier method, between G1/G2 groups using log-rank test (p-value=0.00042; Fig 1A). With a median follow-up of 4.8 years, the 3y-OS was of 100% and 67.8% for G1 and G2 respectively. Difference was also observed for progression-free survival (PFS) although less significant (p-value=0.042; Fig 1 B). We next challenged the prognostic value of our stratification with a multivariate proportional hazard model including known osteosarcoma prognostic factors. Model tested in this stability analysis identified the G1/G2 stratification at diagnosis as the most contributive factor to OS (with 68.5% inclusion frequency; Fig 1C). 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 1 D). As confirmed by a supervised partial least square discriminant analysis (Fig 1 E), IC39 best represents the G1/G2 stratification but several other ICs participate to the distinction between both groups (Fig 1 D, 1 E). First principal component (PCA1) confirmed high association of several ICs to the prognosis, when the second component (PCA2) has not been significantly associated to any clinical item. Therefore, although PCA2 might be associated to 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. These observations support our hypothesis that the complex traits related to our stratification and patient prognosis emerged from the interaction of several gene modules.

1.3. Stratifying 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 G1 versus G2 tumors Iog2 fold-change gene expression on the inferred networks (Data not shown).

We identified that the G1 group expresses genes from two IC41 subnetworks at higher level which are significantly enriched in genes involved in the “innate immune response/interferon 1 response” and “inflammatory response” (Fig 2). Interestingly, IC39, a complex component in terms of function and the most contributive to the G1 group, appears to be related to the epigenetic reprograming of the immune TME (e.g. HDAC11/HDAC6 negatively regulates IL10, KDM6A positively regulates IL6). These observations support conclusions from previous publications about favorable prognostic role of osteosarcoma anti-tumor immunity. In addition, some G1 tumors also express specific cancer testis antigens (CTA; IC26), known as a source of neoantigens. In contrast, the unfavorable prognosis G2 group exhibited connected sub-networks resuming other aspects of the TME previously described as related to poor prognosis and pro-metastatic phenotype such as osteoclast differentiation (IC25), tumor angiogenesis/VEGFR (IC13) and neutrophil degranulation (IC9). Two other isolated subnetworks were strongly upregulated in G2 group (Fig 2): i) The fibroblastic network (IC1) probably reflecting the myofibroblastic reprogramming of osteosarcoma stem cells, crucial for lung metastasis or the cancer-associated fibroblasts (CAF) contributing to mesenchymal-like phenotype and metastasis, and ii) the adipocyte/ PPAR signaling pathway (IC9) involved in inflammatory states and metabolism reprogramming favoring cancer progression. Combined, these G2 characteristics depict an unfavorable osteosarcoma TME prone to induce metastases and imply an early event locking the TME and tumor fate into a vicious cycle.

To complement this functional approach and test the validity of our findings, we returned to the gene expression matrix and performed a differential gene expression analysis between G1/G2 (Fig 2). Consistent with the above results, in G1 tumors we observed the re-expression of specific CTAs, including CTAs from the GAGE family clustered on Xp11.23, and others (including FMR1 NB, PASD1 , NLRP4, CT45A7, MAGEC3, PAGE2, MAEL, TDRD1 , IL13RA2) in addition to genes implicated in spermatogenesis (e.g. RHOXF2B, GTSF1 , RHOXF2, DAZ3). This re-expression evokes some major epigenomic changes at methylation level as required in the control of pre- meiosis or meiosis. In accordance, we observed a significant upregulation of some genes (adjusted p-value< 1.3x1 O' ⁴ ) involved in the RNA gene silencing pathway (e.g. MAEL, DDX4, TDRD1), and known to suppress transposable element expression during the meiosis through the piRNA pathway. Similarly, functional characteristics of G2 group tumors were confirmed by the differential gene expression analysis with upregulation of genes implicated in adipocyte differentiation, osteoclastogenesis and angiogenesis, including some link to other G2 functions and poor prognosis in osteosarcoma (e.g. PLA2G2A, PAQR3, HP).

Thus, all three complementary functional analysis methods were consistent, underlying the major contribution of TME to osteosarcoma progression, with innate immune response associated with good prognosis G1 tumors while angiogenesis, osteoclastogenesis and adipogenesis were associated with poor prognosis G2 tumors. This is congruous with the observed clinical efficacy of multityrosine kinase inhibitors with anti-angiogenic activity in relapsed osteosarcoma while the observed inefficacy of zoledronate in front line osteosarcoma treatment is thought to be partially linked to its action on the immune system.

1.4. Stratifying ICs are associated with distinct large chromosomic regions

In addition to these distinct TME compositions, we observed altered gene expression involving large chromosomic regions, probably related to CNAs dosage effects (Figure not shown, Table 2) and seemingly associated to each of these G1/G2 stratification groups. The most influential IC contributors to the G1 group (Figure not shown, Table 2) reflected known osteosarcoma cancer cell characteristics such as the cytoband 12q14.1 (CDK4, OS9) either associated with 4qter (IC34) or 6p21.1-22.1 (IC6: RLINX2, CDC5L, LIBR2), which contain genes involved in osteosarcoma tumorigenesis and response to chemotherapy. Other G1 group contributive ICs were rather 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 contributive gene of the IC19 was DAXX which, as part of the ATRX/DAXX/HistoneH3.3 complex, regulates telomere maintenance through alternative lengthening of telomeres. For the unfavorable G2 group, three main chromosomic regions were dysregulated on chromosome 6p (IC19), 8q (IC23) and 22q (IC33). The first two are well-known high copy gain or amplification events associated with osteosarcoma oncogenesis more frequent in recurrent/metastatic than in primary osteosarcomas, and previously linked to poor prognosis. Cytoband 6p is a recurrent amplified region in osteosarcoma dysregulating the expression level of the oncogene CDCL5. Cytoband 8q copy number gain is strongly suspected to participate in tumorigenesis through MYC-driven super-enhancer signaling and to be a prognosis factor in osteosarcoma. Surprisingly, chromosome 22q was not previously described as prognostic in osteosarcoma and may require further investigation. These different chromosome imbalances specific to each group might reflect different tumorigenic pathways.

Altogether, those ICs with dysregulated chromosomal regions emphasized a major contribution of CNVs to G1/G2 stratification, even if we cannot exclude that some global modulations might emerge from epigenomic changes. To estimate such potential contribution, we integrated our results with 70 CNV profiles (CGH array) paired with our RNA-sequencing samples. The differential analysis of the aCGH profile using GISTIC2.0 characterized four chromosomes with regions differentially and significantly altered between prognostic groups (adjusted p-value<0.1 , Fig 3A). Hence, gains of 2p21-23 and 22q11.22 cytobands and loss of 11 p12-p15 cytobands were detected in the G2 group while 13q12 region loss was more pronounced in G1 than G2 group. Among these regions, only the chr22:21859431-22318144 (8.27E-03) region was also identified at the expression level (IC33) as participating to G1/G2 stratification. This region contains PI4KAP2, RIMBP3B, RIMBP3C, UBE2L3, YDJC, CCDC116, SDF2L1 , MIR301 B, MIR130B, PPIL2, YPEL1 , MAPK1 , PPM1 F and TOP3B genes. The MAPK signaling pathway has been previously proposed as an important driver of the metastatic stage, whereas PI3K-Akt signaling pathway may participate to early and late stages during osteosarcoma evolution. We also noticed that chr8q cytoband, including MYC, was close to significance in aCGH profile analysis, supporting multiple reports about the MYC amplification involvement in osteosarcoma. Considering the number of ICs associated with chromosomal regions and copy number alterations, relatively few of them have been associated with G1/G2 by our analysis, supporting lower contribution of tumor cells to disease progression than TME composition.

1.5. Validation of the osteosarcoma RNA-seq prognostic signature Gene expression in oncology often shows high dispersion between samples and often 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.

In order to validate the robustness of our stratification in an independent cohort, we identified four gene signatures, based on respectively 37, 15, 8 and 5 genes, predicting G1/G2 tumors by machine learning from initial gene expression count matrix, with logistic regression regularized by elastic-net (Tables 4 to 7), as well as a 20 gene signature predicting G1/G2 tumors by machine learning from initial gene expression count matrix, with logistic regression regularized by LASSO (Table 8).

Table 4: 37 genes signature predicting G1/G2 tumors (obtained by machine learning from initial gene expression count matrix, with logistic regression regularized by elastic-net)

Table 5: 15 genes signature predicting G1/G2 tumors (obtained by machine learning from initial gene expression count matrix, with logistic regression regularized by elastic-net)

Table 6: 8 genes signature predicting G1/G2 tumors (obtained by machine learning from initial gene expression count matrix, with logistic regression regularized by elastic-net)

Table 7: 2 genes signature predicting G1/G2 tumors (obtained by machine learning from initial gene expression count matrix, with logistic regression regularized by elastic- net)

Table 8: 20 genes signature predicting G1/G2 tumors (obtained by machine learning from initial gene expression count matrix, with logistic regression regularized by LASSO) We tested the 15 genes signature shown in Table 5 to predict G1/G2 tumors from an independent cohort of 82 pediatric osteosarcoma tumors whose gene expression count table and paired clinical data are available via open access on the pages regarding the Osteosarcoma project, on the website of the Office of Cancer Genomics (National Cancer Institute). To check prediction validity, we compared OS between predicted G1/G2 and generated related log rank p-values (Fig 3B). Consistent with observations in the OS2006 cohort, results in this independent cohort demonstrated that predicted G2 tumors are significantly associated with worse prognosis than predicted G1 tumors (p-value: 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 in the OS2006 cohort. Another valid criticism to RNA-seq based signature, are about reproducibility in laboratory with different technical setup. We therefore designed a custom Nanostring panel based on the RNA-seq signature that we applied to 166 samples from the OS2006 cohort. 70 samples, already analyzed by RNA-seq, were used to retrained our elastic-net model. G1/G2 prediction of the 96 remaining tumors again confirmed that this stratification is associated to significantly different overall and progression free survival (Fig 3C and Fig 6B).

Finally, we used the RNA-seq based signature to estimate the validity of our stratification at relapse to interrogate the reversibility of this prognostic signature throughout disease progression. We observed similar proportions of G1/G2 tumors sampled at diagnosis from the patients in OS2006 RNA-seq cohort who experienced relapse and at relapse from patients in the MAPPYACTS trial (NCT02613962) (Fig 3D). This result supports the persistence of this stratification across disease progression.

1.6. Discussion

With unsupervised machine learning strategy, we defined the repertoire of gene components describing osteosarcoma tumor clones and TME. We observed that component interactions through co-expression stratify the cohort into good and poor prognosis tumors. Functional characterization of the components associates good prognosis tumors with specific innate immune expression and poor prognosis tumors with angiogenic, osteoclastic, and adipogenic activities; with distinct CNVs specific to each group.

These distinct functional characteristics enable treatment stratification in osteosarcomas, for instance immune modulation (e.g. mifamurtide) for G1 group tumors, and anti-osteoclastic and anti-angiogenic therapies for G2 group. We also identified underlying biological pathways of osteosarcoma, involving PPARy, piRNA or CTAs, which highlight new actionable targets. Our data suggest that early drastic genetic/transcriptomic perturbations in specific clones might influence TME, as well as tumor evolution, response to treatment and metastatic potential.

Finally, we confirmed the predictive power of a minimal prognostic signature of 15 genes within an independent cohort of 82 osteosarcoma primary tumors but also with reproducible Nanostring assay. This work paves the way to the development of prognostic tests, leading to personalized therapy in osteosarcoma, especially by helping clinicians to identify hard-to-treat patients at diagnosis.

Our data highlight that personalized therapies in osteosarcoma should not only be based on genetic abnormalities in the tumor itself (e.g. mutation/CNV), but also on RNA expression profiling to take into account the DNA abnormalities transcribed at RNA level and the TME landscape ²² .

Example 2: Modulation of the PPARy pathway - a new therapeutic target in osteosarcoma

Materials and Methods

MTS cell proliferation assay

HOS, 143B, MG-63, LI2OS and IOR/OS18 cell lines were seeded at 5,000 cells per well, Saos-2, Saos-2-B and IOR/OS14 cell lines were seeded at 10,000 cells per well in a 96-well plate containing a final volume of 100 pl/well and left to settle overnight in DMEM with 10% fetal calf serum.

The cells were treated with different drugs (PPARy agonists troglitazone (TGZ) or rosiglitazone (RGZ), or PPARy antagonist T0070907, dissolved in DMSO) at concentrations ranging from 0 pmol/l to 100 pmol/L. The control (untreated cells) get the same volume of DMSO when the cells are treated (10pl per 1 ml). Cell viability was determined 48 and 72 hours after exposure. Old medium was removed and a solution with MTS/new medium (20 pl of MTS solution - final concentration 0.33 mg/ml) (CellTiter 96 Aqueous One Solution Cell Proliferation Assay; Promega Corporation, Charbonnieres, France) was added. A set of wells were prepared only with MTS/new medium for background subtraction at the same time. An incubation from 1 to 7 hours at 37°C (cell line metabolism dependent) was performed follow by a colorimetric measurement at 490 nm in an automatic plate reader (Elx808; Fisher Bioblock Scientific SAS, lllkirch, France). Western Blot

HOS, 143B, MG-63, U20S, Saos-2, Saos-2-B IOR/OS18 and IOR/OS14 osteosarcoma cell lines were seeded into a 60mm plate dish and collected at approximately 80% confluence. The cells were collected in a lysis buffer (10ml of TNEN 5mM buffer, 1 protease inhibitor pill, 50pl of NaF and 50pl of Orthovanadate (phosphatase inhibitor)) and then obtained by performing an alternation of 5 cycles Frozen in nitrogen and thaw in a water bath at 37°C. Protein supernatants were collected after centrifugation at 4°C for 20min at 13,200rpm.

Protein quantification was performed with the kit BCA of ThermoFisher Scientific (Thermoscientific Pierce TM BCA protein Assay Kit), using a range of bovine serum albumin concentrations (BSA, Euromedex, 04-100-812-E, Souffelweyersheim, France). Absorbance was read at 570 nm using an automatic microplate reader (Elx808; Fisher Bioblock Scientific SAS, I llkirch, France).

Proteins (30pg/well) were separated by 4-20% polyacrylamide gel electrophoresis, tris-glycine extended (Mini-Protean TGX, Bio-Rad, CA, USA), and transferred to polyvinylidene fluoride (PVDF) membranes (Trans-Blot Bio-Rad, CA, USA) using a Trans-Blot Turbo transfer system (Bio-Rad Laboratories, CA, USA).

Membranes were saturated with a 5% BSA for 45 min at room temperature and an antibody against PPARy (5mg/mL; P5505-05B; US Biological, USA) or against AdipoQ (1 :1000; ab75989; Abeam, Cambridge, UK) was added and incubated at 4°C overnight. After, the membranes were washed several times with buffer (TBS 1 % and Tween 0.1 %) and incubated with a secondary antibody (in goat anti-rabbit IgG, 1 :5000; A9169; Sigma-Aldrich, St Louis, MO, USA), at least for 2 hours.

Membranes were revealed using the Clarity Western ECL Substrate kit (Bio-Rad Laboratories, CA, USA), and protein bands detected by chemiluminescence using the Bio-Rad ChemiDoc imaging system (Bio-Rad Laboratories, CA, USA).

A stripping was then performed (stripping buffer 62.5 mL Tris 0.5M, 50 mL SDS 20%, 3.5 mL B-mercaptoethanol, 500 mL H2O) in order to incubate the membranes with the p-actin antibody, directly conjugated to HRP (HRP conjugate; 1 :1000; #5125; Cell Signaling, MA, USA) and revealed as described before.

Results

We studied the therapeutic potential of PPARy pathway in vitro by evaluating cell proliferation under the effect of PPARy agonists and PPARy antagonists (MTS test). The PPARy antagonist T0070907 decreases the cell viability in 7 human osteosarcoma lines with an IC50 ranging between 9-18|JM (median IC50 20|JM). With the agonists, an effect is barely observed at the highest concentrations tested, which could be explained by a negative feed-back loop of PPARy through the expression of a PPARy negative dominant isoform.

Conclusion

Modulation of the PPARy pathway is a new therapeutic target in osteosarcoma.

Example 2: PPARy therapeutic potential in osteosarcoma: preclinical in vitro and in vivo experiments

Materials and Methods

Cells Culture

Human osteosarcoma cell lines HOS, HOS R/MXT, HOS R/DOXO, 143B, U2OS, Saos-2, Saos-2B, MG-63, IOR/OS14, and IOR/OS18, with different genetic background were cultured in Dulbecco’s modified Eagle medium (DMEM, Invitrogen, Saint Aubin, France) supplemented with 10% (v/v) fetal bovine serum (FBS, Invitrogen, Saint Aubin, France) at 37°C in a humidified atmosphere (5% CO2 and 95% air), under mycoplasma free conditions. shRNA and overexpression were performed on bioluminescent cell line (Luciferase/mKate2). qRT-PCR of PPARg expression

PPARg expression was evaluated by qRT-PCR. ARN extraction was performed with AHPrep DNA/RNA mini kit (Cat. No. I ID: 80204 Qiagen) according to manufacturer instructions. RNA (1 pg) was then reverse using M MLV reverse transcriptase (Invivogen ref 28025.013). Finally, 5pL of cDNA was mix with 6pL of nuclease free water, 1.5pM of primers at 10pM (PPARy-s TTGACTTCTCCAGCATTTCTAC (SEQ ID No: 42) + PPAR _Y -as

CTTTATCTCCACAGACACGAC(SEQ ID No: 43)) and 12.5pl of syber green master mix (ref K0223 Thermofisher scientific). Amplification of PPARy and GAPDH was performed with 1 cycle at 50°C for 2 minutes and 95°C for 10 minutes followed by 40 cycles of 95°C for 15 seconds, 60°C for 1min. Melting curve was performed at the end of the PCR (95°C for 15 s, 60°C for 1 min and 95°C for 15 s) to identify unique PCR products. Amplifications were monitored with ViiA 7 Real-Time PCR System. GAPDH was used as housekeeping gene. Calculation of the relative expression of each transcript was performed using the ₂ -AAct method.

Western blot

Osteosarcoma cell lines were seeded as mentioned before incubated at 37°C overnight and treated or not in the day after. Cell pellets were resuspended on 10OpI of lysis buffer (in 10ml of TNEN 5mM buffer add ¹ X protease inhibitor pill, 50pl of NaF and 50pl of Orthovanadate) 6, 24, 48 or 72h after. Proteins were extracted by frozen cell suspension in nitrogen and thaw in a water bath at 37°C (5X), follow by centrifugation at 13.200rpm at 4°C for 20min. Protein quantification was performed using the BCA protein Assay Kit (Thermoscientific PierceTM BCA protein Assay Kit) according to manufacturer instructions

Briefly, proteins (20pg- 30pg) were separated with 4-15% Mini-Proteam TGX stain free gel, (ref 4568086 Bio-Rad)and then transfered into nitrocellulose membrane (Ref 1704156 Bio-Rad) with the transblot Turbo transfer system (Bio-Rad). Membrane are then incubated overnight with PPARY primary antibody (Anti-PPARg polyclonal P5505-05B; US Biological). Membrane were washed ( 5x with wash buffer) and incubated for 2h with secondary antibody (Goat anti-rabbit a9169 sigma, 1/5000 All immunoblots images were performed using Bio-Rad ChemiDoc imaging system. Membrane were then incubated in striping solution for B- actin characterization (same procedure). Relative band intensities were performed using Image J software.

Compounds

Doxorubicin (DOXO), methotrexate (MTX), were purchased from Sigma Aldrich (Lyon, France), mafosfamide (MAF), T0070907, Rosiglitazone (RGZ) and Troglitazone (TGZ) from Clinisiences. (Nanterre, France). All compounds were solubilized in dimethyl sulfoxide (DMSO; Sigma Aldrich, Lyon, France), at 10 mM stock solutions, and stored at -20 °C. For in vivo all the compounds used were solubilized in PBS, except T0070907 solubilized in 2,5%DMSO+ 47,5%PEG 400+ 50% PBS.

In-Vitro Cell Viability Assays

Parental HOS, 143B, MG-63, IOR/OS18, resistant derived HOS R/MTX and derived PPARy shRNA or overexpressing PPARy cell lines were seeded at 5000 cells/well, parental , Saos2, Saos-2B, IOR/OS14 , resistant derived HOS R/DOXO and derived PPAR shRNA or overexpressing PPARy\ir\es were seeded at 10,000 cells/well in a 96-well plate in DM EM supplemented with 10% FBS for both assays.

The day after seeding, cells were incubated in the presence of a range of drug concentrations for 72h (from 0 to 100 pmol/L for DOXO, MAF; 0 to 500 pmol/L for MTX, T0070907, Rosiglitazone (RGZ) and Troglitazone (TGZ)). Cell viability was evaluated using the CellTiter 96 Aqueous One Solution Cell Proliferation Assay (MTS assay) (Promega, Charbonnieres, France), according to the manufacturer instructions. The half-maximal inhibitory concentration (IC50) was determined using the GraphPad Prism5 software (Graphpad Software Inc., California, USA).

Combination studies

Cells were treated with increasing concentrations of drugs either alone or in combination at their equipotent molar ratio concomitantly. Effects on cell number were determined by MTS assay according to the manufacturer instructions. The results were analyzed using the median effect analysis method (12) and by deriving the combination index (Cl), which was calculated at equipotent combined drug concentrations that inhibit growth at 50% (ED50). Exclusive Cl values were used to analyze combinations.

PDX models.

Experiments were validated by the CEEA26, Ethic committee (approval number: APAFIS#27183-2020091414229028 v3) and carried out under conditions established by the European Community (Directive 2010/63/UE). Animals were purchased at Gustave Roussy (Villejuif, France) and maintained in the respective animal facilities following standard animal regulation, health and care, and ethical controls. Osteosarcoma PDXs were established from relapsed patients by transplantation in immunocompromised NSG mice and considered as established when tumor growth was sustained after at least two in vivo passages. Further implantations will be performed by direct implantation of a tumor fragment from the previous passage, either fresh or conserved by soft congelation (frozen in FBS, +10% DMSO). Under anesthesia (3% isoflurane, 1.5L/min air), tumor samples were implanted in an orthotopic position, paratibially (~2mm3) between muscle and bone tibia after a 0.5 cm skin incision and a gentle activation of the periosteum (periosteum denudation). To avoid bone pain, an analgesic (buprenorphine at 0.3 mg/kg) will be applied in addition to the general anesthesia or when symptoms appeared. Clinical status, tumor uptake and tumor growth will be evaluated one to 3 times a week. Paratibial tumor was detected by palpation and tumor gross appearance (caliper measurements). The experiments lasted until tumors reached specific tumor volume -1500 mm3, significant weight loss, or difficulty to walk. Mice were then anesthetized and bone structure alterations were analyzed by CTscan imaging. Mice were euthanized at the endpoint and samples harvested and processed as describe below. In Vivo T0070907 Treatments of osteosarcoma PDX orthotopic models.

Groups of 8 animals bearing the same OTS orthotopic PDX model were treated from day 14 after tumor implantation with intraperitoneal (IP) injection (volume of 10 ml/kg) every day with either T0070907 (10mg/Kg/injection) or with vehicle (saline, control group). For the combination of T0070907 with methotrexate, 4 groups of X animals were treatment with T0070907 (10mg/injection) at J1-J4 or J1-J4 and J6-J7 and methotrexate (10mg/injection) at J5, for combo group and saline in the control group. Clinical status, tumor uptake and tumor growth were evaluated each twodays. Paratibial tumors were detected by palpation, tumor gross apparition (caliper measurements). Tumor CTscan imaging were performed once a week. Tumor leg, normal leg, lungs, spleen and liverwere harvested at time of sacrifice and conserved for future analysis (RNAseq, WES, Histology, single cell, spatial transcriptomic...).

In vivo CTscan imaging

IVIS SpectrumCT (Perkin Elmer, Courtaboeuf, France) was used for images acquirement. This system allows the primary tumor detection by X-ray tomography. CTscan imaging were performed under anesthesia with 3%(v/v) isoflurane.

Histology

Organs were fixed in a 4%(v/v) paraformaldehyde, and embedded in paraffin. Tissues were stained with hematoxilin-eosin-safranin (HES) for morphology. Paraffin sections were processed for heat-induced antigen retrieval (ER2 corresponding EDTA buffer pH9) for 20 min at 100°. Slides were incubated with a mouse monoclonal anti human Ki67 antibody (clone MIB1 ; 1 :20; Agilent Dako) or with Anti-PPARy polyclonal (P5505-05B ; US Biological) for 1h at room temperature. The nuclear signal was revealed with the Klear mouse kit (GBI labs). Slides were examined using light microscopy (Zeiss, Marly-Le-Roy, France) and a single representative whole tumor tissue section from each animal was digitized using a slide scanner NanoZoomer 2.0-HT (C9600-13, Hamamatsu Photonics). Histology was reviewed by a bone expert pathologist.

Background

The outcomes of adolescents/young adults with osteosarcoma have not improved in decades.

From the RNA sequencing of 79 osteosarcoma diagnostic biopsies, we identified stable independent components recapitulating the tumor microenvironment and clones (Marchais et al Cancer Research 2022, in press, see Example 1). Metagene unsupervised classification 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 favorable G1 tumors with innate immunity and unfavorable G2 tumors with angiogenic, osteoclastic and adipogenic activities as well as PPARy pathway upregulation.

PPARy (gamma isotype of Peroxisome Proliferator Activated Receptors) is a nuclear receptor implicated in different biological processes, adipogenesis, angiogenesis, and immunity (macrophage polarization) ⁷⁰ ’ ⁷¹ ; and have been variously implicated in cancer.

Both pro-and anti-cancer effect have been described depending on the tumor type. Consequently, both antagonists and agonists have been explored as potential anti-cancer therapies.

PPARy is usually considered as a tumor suppressor gene through its antiproliferative and pro-apoptotic and re-differentiation role in several cancers. Loss of expression ⁷² , mutation impairing ligand binding and thus its transcriptional activity ⁷³ , have been described in epithelial cancers. In these cases, synthetic agonists of PPARy might harbor anti-cancer activity. Oncogenic fusions PAX8- PPARy are described in thyroid cancers and PPARy ligand pioglitazone induce trans-differentiation toward adipocytic-like cells. Negative dominant isoforms might also be expressed in some cancers (e.g. ORF4 in colic cancer) ⁷⁴ .

At the opposite, PPARy activation in some cancer cells might induce an energetic switch toward fatty acid utilization and be implicated in cell proliferation and survival ³⁵ . Several target genes of PPARy are also implicated in tumor aggressiveness, such as angiogenesis. PPARy activation in a fatty acid rich microenvironment might favor metastatic development in some cancers (breast cancer, melanoma) ⁷⁵ . In that case, T0070907, a PPARy specific synthetic antagonist, might decrease metastatic occurrence ⁷⁵ .

In osteosarcoma, in vitro, PPAR-y agonist troglitazone (5pM) favor cell line survival through inhibition of the AKT-dependent spontaneous apoptosis ⁷⁶ . Higher concentrations had an opposite effect with troglitazone (100pM) having an antiproliferative effect in vitro and anti-tumor effect in vivo through pro-apoptotic et prodifferentiating effect, possibly due to negative retrocontrol loops with dominant negative isoforms ⁷⁸ ’ ⁷⁹ .

Known effects of PPARy activation are summarized in Figure 2 of Ahmadian et al. (2013 Nature Medecine; 19(5): PPARy signaling and metabolism: the good, the bad and the future). PPAR-y and targets in osteosarcoma biopsies at diagnosis: QS2006 cohort

Based on the RNAseq study we performed on OS2006 osteosarcoma cohort at diagnosis, PPARy signalling pathway is associated with the G2 poor prognostic group identified with our signature. RNA expression of several PPARy targets was correlated with other components linked to poor prognosis identified in this cohort pro- angiogenic, osteoclastic and adipocytic activity, in addition to PPARy activation. All this suggests pro-tumor and pro- metastatic activity of PPARy in osteosarcoma and a potential therapeutic role of PPARy antagonist in these patients.

PPAR-y in osteosarcoma: in vitro preclinical data

We used 8 osteosarcoma cell lines, all expressing PPARy (Western blot; Figure 4B) and their resistant counterparts to methotrexate or doxorubicin ⁷⁹ to evaluate in vitro anti-proliferative effect (MTS assay) of two PPARy agonists, troglitazone and rosiglitazone, and one synthetic specific antagonist T0070907, alone and in combination with chemotherapy used in osteosarcoma patients (methotrexate, doxorubicin, mafosfamide) (median effect analysis method) ⁸⁰ .

In vitro anti-proliferative effect of T0070907 in osteosarcoma cell lines.

Median IC50 of T0070907 for all 8 osteosarcoma parental cell lines was 20pM (range 9.3-37.4; MTS assay) while for the two agonists, IC50 was not reached with concentrations up to 100pM (Figure 4A). T0070907 IC50 in osteosarcoma cells were similar to IC50 of various histologies of adult cancers ⁸¹ . RNA sequencing of HOS cells exposed to T0070907 at 25pM (IC50) for 24h confirm induction of pro-apoptotic and osteoblastic differentiation programs (Figure 7). For the PPARy agonist (Rosiglitazone), only the adipocyte differentiation program was significantly up- regulated in treated osteosarcoma cells as expected, PPARy being a major regulator of adipogenesis.

In vitro synergy between of T0070907 with chemotherapy in HOS osteosarcoma cell lines, parental and resistant to methotrexate and doxorubicin.

We then tested the combination effect of T0070907 and chemotherapy routinely used in osteosarcoma (methotrexate, doxorubicin, mafosfamide) using MTS assays and the median effect analysis method ⁷⁸ in the osteosarcoma cell line HOS and their resistant counterpart to methotrexate and doxorubicin. Synergy was observed with all drugs in all cell lines and in the timing of administration of T0070907 (24h before, simultaneously, or 24h after chemotherapy), with a resistance index higher when T0070907 was administered first (Table"!). compound concentrations that inhibit osteosarcoma cell growth by 50% as determined by the median effect analysis method in HOS parental cell line and its resistant counterpart to methotrexate and doxorubicin, for PPARy antagonist T0070907 combined chemotherapy, (methotrexate, MTX; doxorubicin, DOXO; mafosfamide, MAF)

PPARy in osteosarcoma: in vivo preclinical data in PDX models

According to our initial hypothesis, PPARy antagonist might modify the composition of the tumor microenvironment of osteosarcoma with poor prognosis. To test this hypothesis, we treated with T0070907, alone and in combination with methotrexate, an osteosarcoma orthotopic (Paratibial) patient derived xenograft (PDX) model issued from a metastatic sample of patient at relapse MAP-217-PT. This model was chosen because of its expression of PPARy and its capacity to form lung metastasis. PPARy expression by IHC was observed in several PDX models (Fig. 8).

After an initial dose testing to find the concentration of T0070907 tolerable by the mice, the effect of T0070907 alone administered intraperitoneally (IP) at 10mg/kg/injection daily (Fig. 9A) or administered intraperitoneally (IP) at 10mg/kg/injection Day1-4 or Day1-4 and 6-7 combined with methotrexate at Day5 at 10mg/kg/injection (Fig. 9B), was evaluated in vivo by tumor growth of the primary tumor and analysis of the lung metastases at mice sacrifice.

T0070907 induced delayed tumor growth in MAP-217 PT PDX model (delay to reach 2.5 time the initial leg volume reflecting the tumor volume was 20% higher in T0070907 treated mice compared to control, Fig. 9A) and decreased lung metastases in number and size (quantification ongoing by the pathologists; Fig. 9B). In addition, T0070907 combined with methotrexate seems to have a greater effect than each drug alone in terms of primary tumor growth in vivo (Fig. 9C). Effect of the combo on metastasis is currently being quantified by the pathologist.

To further understand the effect of T0070907 alone and in combination with methotrexate on the tumor cells and tumor microenvironment, spatial transcriptomic is being performed.

Abreviations used in this text:

API-AI Adriamycin/platinum/ifosfamide

AUC Area under curve

BCA BiCinchoninic acid Assay

BSA Bovine Serum Albumin

CAF Cancer-Associated Fibroblasts

CGH Comparative Genomic Hybridization

CNA Copy Number Aberrations

CNV Copy Number Variation

CTA Cancer Testis Antigens

DASL cDNAmediated Annealing, Selection, and Ligation

DMEM Dulbecco/Vogt modified Eagle's minimal essential medium

DMSO Dimethylsulfoxide

DNA Deoxyribonucleic acid

DR Death Rate dllTP Deoxyuridine triphosphate

FOV Fields of view

HRP Horse Radish Peroxidase

IC Independent Components

ICA Independent Components Analysis

IL Interleukine

MAP Methotrexate, doxorubicin and cisplatinum

MCR Minimum common region

M-EI MTX-etoposide/ifosfamide

OLA Oligonucleotide Ligation Assay

OS Overall survival

OTS Osteosarcoma

PCA Principal Component Analysis

PFER Per-family error rate

PFS Progression-free survival

PPARy Peroxisome prol iterator-activated receptor

PVDF Polyvinylidene fluoride qRT-PCR Quantitative reverse-transcription polymerase chain reaction

RGZ Rosiglitazone

RNA Ribonucleic acid

RNA-seq RNA-sequencing

TGZ Troglitazone

TME Tumor microenvironment

WES Whole Exome Sequencing

WGS Whole Genome Sequencing REFERENCES

1. Perry, J. A. et al. Complementary genomic approaches highlight the PI3K/mTOR pathway as a common vulnerability in osteosarcoma. Proc. Natl. Acad. Sci. 111 , E5564-E5573 (2014).

2. Lorenz, S. et al. Unscrambling the genomic chaos of osteosarcoma reveals extensive transcript fusion, recurrent rearrangements and frequent novel TP53 aberrations. Oncotarget 7, (2016).

3. Trama, A. et al. Survival of European adolescents and young adults diagnosed with cancer in 2000-07: population-based data from EUROCARE-5. Lancet Oncol. (2016) doi : 10.1016/S1470-2045(16)00162-5.

4. Bielack, S. S. et al. Methotrexate, Doxorubicin, and Cisplatin (MAP) Plus Maintenance Pegylated Interferon Alfa-2b Versus MAP Alone in Patients With Resectable High-Grade Osteosarcoma and Good Histologic Response to Preoperative MAP: First Results of the EU RAMOS-1 Good Response Randomized Controlled Trial. J. Clin. Oncol. 33, 2279-2287 (2015).

5. Neyssa M Marina et al. Comparison of MAPIE versus MAP in patients with a poor response to preoperative chemotherapy for newly diagnosed high-grade osteosarcoma (EURAMOS-1): an open-label, international, randomised controlled trial. Lancet Oncol 1396-1408 (2016).

6. Piperno-Neumann, S. et al. Zoledronate in combination with chemotherapy and surgery to treat osteosarcoma (OS2006): a randomised, multicentre, open-label, phase 3 trial. Lancet Oncol. 17, 1070-1080 (2016).

7. Meyers, P. A. et al. Osteosarcoma: The Addition of Muramyl Tripeptide to Chemotherapy Improves Overall Survival-A Report From the Children’s Oncology Group. J. Clin. Oncol. 26, 633-638 (2008).

8. Ozaki, T. et al. Genetic imbalances revealed by comparative genomic hybridization in osteosarcomas. Int. J. Cancer 102, 355-365 (2002).

9. Man, T.-K. et al. Genome-wide array comparative genomic hybridization analysis reveals distinct amplifications in osteosarcoma. BMC Cancer 4, 1 (2004).

10. Atiye, J. et al. Gene amplifications in osteosarcoma-CGH microarray analysis: CGH Microarray Analysis of Osteosarcoma. Genes. Chromosomes Cancer 42, 158-163 (2005).

11. Selvarajah, S. et al. Genomic signatures of chromosomal instability and osteosarcoma progression detected by high resolution array CGH and interphase FISH. Cytogenet. Genome Res. 122, 5-15 (2008).

12. Fletcher. Identification of chromosomal aberrations associated with disease progression and a novel 3q13.31 deletion involving LSAMP gene in osteosarcoma. Int. J. Oncol. 35, (2009).

13. Kovac, M. et al. Exome sequencing of osteosarcoma reveals mutation signatures reminiscent of BRCA deficiency. Nat. Commun. 6, 8940 (2015).

14. Bousquet, M. et al. Whole-exome sequencing in osteosarcoma reveals important heterogeneity of genetic alterations. Ann. Oncol. 27, 738-744 (2016).

15. Chen, K. S. et al. A novel TP53-KPNA3 translocation defines a de novo treatmentresistant clone in osteosarcoma. Mol. Case Stud. 2, a000992 (2016). 16. Chiappetta, C. et al. Whole-exome analysis in osteosarcoma to identify a personalized therapy. Oncotarget 8, 80416 (2017).

17. Chen, X. et al. Recurrent Somatic Structural Variations Contribute to Tumorigenesis in Pediatric Osteosarcoma. Cell Rep. 7, 104-112 (2014).

18. Ribi, S., Baumhoer, D. & Lee, K. TP53 intron 1 hotspot rearrangements are specific to sporadic osteosarcoma and can cause Li-Fraumeni syndrome. Oncotarget 6, 7727 (2015).

19. Behjati, S. et al. Recurrent mutation of IGF signalling genes and distinct patterns of genomic rearrangement in osteosarcoma. Nat. Commun. 8, 15936 (2017).

20. Gambera, S. et al. Clonal dynamics in osteosarcoma defined by RGB marking. Nat. Commun. 9, (2018).

21. Petitprez, F. et al. Quantitative Analyses of the Tumor Microenvironment Composition and Orientation in the Era of Precision Medicine. Front. Oncol. 8, 390 (2018).

22. Quail, D. F. & Joyce, J. A. Microenvironmental regulation of tumor progression and metastasis. Nat. Med. 19, 1423-1437 (2013).

23. Kairov, II. et al. Determining the optimal number of independent components for reproducible transcriptom ic data analysis. BMC Genomics 18, 712 (2017).

24. Liberzon, A. etal. Molecular signatures database (MSigDB) 3.0. Bioinformatics 27, 1739-1740 (2011).

25. Raudvere, II. et al. g:Profiler: a web server for functional enrichment analysis and conversions of gene lists (2019 update). Nucleic Acids Res. 47, W191-W198 (2019).

26. Cheng, F. 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. 60, 44-53 (2014).

27. de Groot, A. E. & Pienta, K. J. Epigenetic control of macrophage polarization: implications for targeting tumor-associated macrophages. Oncotarget 9, (2018).

28. Buddingh, E. P. et al. Tumor-Infiltrating Macrophages Are Associated with Metastasis Suppression in High-Grade Osteosarcoma: A Rationale for Treatment with Macrophage Activating Agents. Clin. Cancer Res. 17, 2110-2119 (2011).

29. Gomez-Brouchet A 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 QS2006 phase 3 trial. Oncoimmunology. Aug 24, doi: 10.1080/2162402X.2017.1331193. (2017).

30. Salmaninejad, A. et al. Cancer/Testis Antigens: Expression, Regulation, Tumor Invasion, and Use in Immunotherapy of Cancers. Immunol. Invest. 45, 619-640 (2016).

31. Kelleher, F. C. & O’Sullivan, H. Monocytes, Macrophages, and Osteoclasts in Osteosarcoma. J. Adolesc. Young Adult Oncol. (2017) doi: 10.1089/jayao.2016.0078.

32. Yang, J. et al. Genetic amplification of the vascular endothelial growth factor (VEGF) pathway genes, including VEGFA , in human osteosarcoma: VEGFA Amplification in Osteosarcoma. Cancer 117, 4925-4938 (2011). 33. Mollinedo, F. Neutrophil Degranulation, Plasticity, and Cancer Metastasis. Trends Immunol. 40, 228-242 (2019).

34. Zhang, W. et al. Adaptive Fibrogenic Reprogramming of Osteosarcoma Stem Cells Promotes Metastatic Growth. Cell Rep. 24, 1266-1277. e5 (2018).

35. Liu, Y. et al. The Role of PPAR-5 in Metabolism, Inflammation, and Cancer: Many Characters of a Critical Transcription Factor. Int. J. Mol. Sci. 19, 3339 (2018).

36. Fratta, E. et al. The biology of cancer testis antigens: Putative function, regulation and therapeutic potential. Mol. Oncol. 5, 164-182 (2011).

37. McFarlane, R. J. & Wakeman, J. A. Meiosis-like Functions in Oncogenesis: A New View of Cancer. Cancer Res. 77, 5712-5716 (2017).

38. Ozata, D. M., Gainetdinov, I., Zoch, A., O’Carroll, D. & Zamore, P. D. PIWI- interacting RNAs: small RNAs with big functions. Nat. Rev. Genet. 20, 89-108 (2019).

39. Mintz, M. B. et al. An Expression Signature Classifies Chemotherapy-Resistant Pediatric Osteosarcoma. Cancer Res 8 (2005).

40. Ma, Z. et al. The tumor suppressor role of PAQR3 in osteosarcoma. Tumor Biol. 36, 3319-3324 (2015).

41. Arredouani, M. et al. Haptoglobin directly affects T cells and suppresses T helper cell type 2 cytokine release. Immunology 108, 144-151 (2003).

42. Kwon, J.-O. et al. Haptoglobin Acts as a TLR4 Ligand to Suppress Osteoclastogenesis via the TLR4-IFN-P Axis. J. Immunol. ji1800661 (2019) doi:10.4049/jimmunol.1800661.

43. Duffaud, F. etal. Efficacy and safety of regorafenib in adult patients with metastatic osteosarcoma: a non-com parative, randomised, double-blind, placebo-controlled, phase 2 study. Lancet Oncol. 20, 120-133 (2019).

44. Mannerstrdm, B. et al. Epigenetic alterations in mesenchymal stem cells by osteosarcoma-derived extracellular vesicles. Epigenetics 14, 352-364 (2019).

45. Lu, X.-Y. et al. Cell cycle regulator gene CDC5L, a potential target for 6p12-p21 amplicon in osteosarcoma. Mol. Cancer Res. 6, 937-946 (2008).

46. Martin, J. W. et al. Digital Expression Profiling Identifies RUNX2, CDC5L, MDM2, RECQL4, and CDK4 as Potential Predictive Biomarkers for Neo-Adjuvant Chemotherapy Response in Paediatric Osteosarcoma. PLoS ONE 9, e95843 (2014).

47. Hoelper, D., Huang, H., Jain, A. Y., Patel, D. J. & Lewis, P. W. Structural and mechanistic insights into ATRX-dependent and -independent functions of the histone chaperone DAXX. Nat. Commun. 8, 1193 (2017).

48. Lau, C. C. et al. Frequent amplification and rearrangement of chromosomal bands 6p12-p21 and 17p11.2 in osteosarcoma. Genes. Chromosomes Cancer 39, 11- 21 (2004).

49. Chen, D. et al. Super enhancer inhibitors suppress MYC driven transcriptional amplification and tumor progression in osteosarcoma. Bone Res. 6, (2018).

50. Tarkkanen, M. et al. DNA sequence copy number increase at 8q: A potential new prognostic marker in high-grade osteosarcoma. Int. J. Cancer 84, 114-121 (1999). 51. Wang, D. et al. Multiregion Sequencing Reveals the Genetic Heterogeneity and Evolutionary History of Osteosarcoma and Matched Pulmonary Metastases. Cancer Res. 79, 7-20 (2019).

52. Gaspar, N. et al. Results of methotrexate-etoposide-ifosfamide based regimen (M- El) in osteosarcoma patients included in the French OS2006/sarcome-09 study. Eur. J. Cancer 88, 57-66 (2018).

53. Piperno-Neumann, S. et al. Results of API-AI based regimen in osteosarcoma adult patients included in the French OS2006/Sarcome-09 study: Results of API- AI based regimen in osteosarcoma adult patients included in the French OS2006/Sarcome-09 study. Int. J. Cancer (2019) doi:10.1002/ijc.32526.

54. Shannon, P. Cytoscape: A Software Environment for Integrated Models of Biomolecular Interaction Networks. Genome Res. 13, 2498-2504 (2003).

55. Haw, R., Loney, F., Ong, E., He, Y. & Wu, G. Perform Pathway Enrichment Analysis Using ReactomeFIViz. in vol. 2074: 165-179. doi: 10.1007/978-1-4939- 9873-9_13. PubMed PMID: 31583638 (2020).

56. Le Cao, K.-A., Rohart, F., Gonzalez, I. & Dejean, S. mixOmics: Omics Data Integration Project. R package version 6.1.1. (2016).

57. Love, M. I., Huber, W. & Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15, (2014).

58. Love, M. I., Huber, W. & Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15, 550 (2014).

59. Friedman, J., Hastie, T. & Tibshirani, R. Regularization Paths for Generalized Linear Models via Coordinate Descent. J. Stat. Softw. 33, (2010).

60. Picard, F. et al. Joint segmentation, calling, and normalization of multiple CGH profiles. Biostatistics 12, 413-428 (2011).

61 . Nilsen, G. et al. Copynumber: Efficient algorithms for single- and multi-track copy number segmentation. BMC Genomics 13, 591 (2012).

62. Karolchik, D. The UCSC Genome Browser Database. Nucleic Acids Res. 31 , 51- 54 (2003).

63. Gaujoux, R. & Seoighe, C. A flexible R package for nonnegative matrix factorization. BMC Bioinformatics 11 , 367 (2010).

64. Hornik, K., Feinerer, I., Kober, M. & Buchta, C. Spherical k-Means Clustering. J. Stat. Softw. 22.

65. Mermel, C. H. et al. GISTIC2.0 facilitates sensitive and confident localization of the targets of focal somatic copy-number alteration in human cancers. Genome Biol. 12, R41 (2011).

66. Jiang P. et al. Genome-Scale Signatures of Gene Interaction from Compound Screens Predict Clinical Efficacy of Targeted Cancer Therapies. Cell.Syst. Mar 28;6(3):343-354.e5.(2018).

67. Liu X. et al. A prognostic gene expression signature for oropharyngeal squamous cell carcinoma. EBioMedicine. Nov;61 : 102805 (2020).

68. Vey, J. et al. A Toolbox for Functional Analysis and the Systematic Identification of Diagnostic and Prognostic Gene Expression Signatures Combining Meta-Analysis and Machine Learning. Cancer. Oct 21 ;11(10):1606 (2019).

69. Xia Y. et al. Genetic determinants of the molecular portraits of epithelial cancers. Nat.Commun. Dec 11 ;10(1):5666 (2019). 70. Christofides A, Konstantinidou E, Jani C, Boussiotis VA. The role of peroxisome proliferator-activated receptors (PPAR) in immune responses. Metabolism - Clinical and Experimental [Internet], 1 janv 2021 [cite 17 janv 2022];114.

71 Harmon GS, Lam MT, Glass CK. PPARs and lipid ligands in inflammation and metabolism. Chem Rev. 12 oct 2011 ;111(10):6321-40.

72. Niu S-F, Cui B-X, Huang J-Z, Guo Y. PPARy is correlated with prognosis of epithelial ovarian cancer patients and affects tumor cell progression in. :8.

73. Sarraf P, Mueller E, Smith WM, Wright HM, Kum JB, Aaltonen LA. Loss-of- Function Mutations in PPARy Associated with Human Colon Cancer. Molecular Cell. :6.

74. Sabatino L, Casamassimi A, Peluso G, Barone MV, Capaccio D, Migliore C, et al. A Novel Peroxisome Proliferator-activated Receptor y Isoform with Dominant Negative Activity Generated by Alternative Splicing. J Biol Chem. 15 juill 2005;280(28):26517-25.

75. Zou Y, Watters A, Cheng N, Perry CE, Xu K, Alicea GM, et al. Polyunsaturated Fatty Acids from Astrocytes Activate PPARy Signaling in Cancer Cells to Promote Brain Metastasis. Cancer Discov. dec 2019;9(12):1720-35.

76. Lucarelli E, Sangiorgi L, Maini V, Lattanzi G, Marmiroli S, Reggiani M, et al. Troglitazione affects survival of human osteosarcoma cells. Int J Cancer. 20 mars 2002;98(3):344-51.

77. Haydon RC, Zhou L, Feng T, Breyer B, Cheng H, Jiang W, et al. Nuclear Receptor Agonists As Potential Differentiation Therapy Agents for Human Osteosarcoma. Clinical Cancer Research, mai 2002;8:1288-94.

78. He BC, Chen L, Zuo GW, Zhang W, Bi Y, Huang J, et al. Synergistic Antitumor Effect of the Activated PPAR and Retinoid Receptors on Human Osteosarcoma. Clinical Cancer Research. 15 avr 2010;16(8):2235-45.

79. Marques da Costa ME, Marchais A, Gomez-Brouchet A, Job B, Assoun N, Daudigeos-Dubus E, et al. In-Vitro and In-Vivo Establishment and Characterization of Bioluminescent Orthotopic Chemotherapy-Resistant Human Osteosarcoma Models in NSG Mice. Cancers. 17 juill 2019;11 (7):997.

80. Chou T-C. Theoretical Basis, Experimental Design, and Computerized Simulation of Synergism and Antagonism in Drug Combination Studies. Pharmacological Reviews. 1 sept 2006;58(3):621-81.

81. Burton JD, Goldenberg DM, Blumenthal RD. Potential of Peroxisome Proliferator- Activated Receptor Gamma Antagonist Compounds as Therapeutic Agents for a Wide Range of Cancer Types. PPAR Research. 2008;2008:1-7.