Keywords: Image processing, statistical learning, brain tumours, data integration

Diffuse Intrinsic Pontine Glioma (DIPG) is a rare brain tumour located in the pons, mainly seen in children aged 5 to 7 years. It is considered one of the most aggressive paediatric tumours, with a survival rate of less than 10% beyond two years after diagnosis and a median overall survival of less than one year. DIPG is classified as a diffuse midline glioma (DMG), mainly characterized by a K27M mutation of the genes encoding the histone H3 protein and/or a loss of H3K27 trimethylation by overexpression of the EZHIP protein. The location of the tumour and its corresponding genomic alterations makes DIPG a completely different type of tumour than other high-grade gliomas.

This work proposes the integration of imaging data with genetic data in order to find biomarkers. First, we are interested in the extraction of the regions of interest of the images necessary for a radiomic study. Then, we propose a procedure for the integration of multisource data, which takes into account the complex graphs of interaction between genes. Finally, we apply our procedure to the available data in order to compare its performance with other existing models and to study the contribution of imagery and the graph to genetic data.

Radiomic analysis requires predefined regions of interest on available images. For our DIPG cohort, manual tumour segmentation was not feasible. Moreover, no database was created to train classical machine learning algorithms to automatically delineate tumor regions. This study focused on obtaining binary segmentations for DIPG using only FLAIR and T2w modalities, using models trained on glio-

List of Algorithms

RGCCA optimisation algorithm . . . . . . .

Titre : Machine learning pour l'intégration des données génomiques et d'imagerie appliquée à la neuro-oncologie Mots clés :traitement d'images, apprentissage statistique, tumeurs cérébrales, intégration de données Résumé :Le gliome infiltrant du tronc cérébral (Diffuse Intrinsic Pontine Glioma (DIPG)) est une tumeur cérébrale rare située dans le pons, principalement observé chez les enfants âgés de 5 à 7 ans. Elle est considérée comme l'une des tumeurs pédiatriques les plus agressives, avec un taux de survie inférieur à 10% au-delà des deux ans après le diagnostic et une médiane de survie globale inférieure à un an. Le DIPG est classé comme un gliome diffus de la ligne médiane (DMG), principalement caractérisé par une mutation K27M des gènes codant pour la protéine histone H3 et/ou une perte de la triméthylation H3K27 par surexpression de la protéine EZHIP. L'emplacement de la tumeur et ses altérations génomiques correspondantes fait du DIPG un type de tumeur complètement différent des autres tumeurs de haut grade.

Ce travail propose l'intégration des données d'imagerie avec les données génétiques afin de trouver des biomarqueurs. Dans un premier temps, nous nous intéressons à l'extraction des régions d'intérêt des images nécessaires pour une étude radiomique. Ensuite, nous proposons une procédure d'intégration des données multi-sources, qui prend en compte les graphes complexes d'interaction entre les gènes. Finalement, nous appliquons notre procédure sur les données disponibles afin de comparer ses performances avec d'autres modèles de la littérature et étudier l'apport de l'imagerie et du graphe aux données génétiques.

L'analyse radiomique nécessite des régions d'intérêt prédéfinies sur les images disponibles. Pour notre cohorte DIPG, la segmentation manuelle de la tumeur n'était pas disponible car elle ne fait pas partie de la routine clinique. De plus, aucune base de données spécifique n'a été créée pour entraîner des algorithmes d'apprentissage automatique classiques afin de délimiter automatiquement les régions tumorales. Cette étude s'est concentrée sur l'obtention de segmentations binaires pour le DIPG en utilisant uniquement les mo-dalités FLAIR et T2w, à partir de modèles entraînés sur des gliomes de haut grade. Nous proposons de combiner différents modèles simples de détection et de segmentation pour obtenir des résultats de segmentation satisfaisants.

En parallèle, un modèle d'intégration multiblocs prenant en compte des graphes complexes connus d'interactions entre les gènes a été développé et l'influence du graphe choisi sur la sélection des variables par le modèle a été étudiée. Nous proposons netSGCCA, un modèle combinant la Sparse Generalized Canonical Correlation Analysis (SGCCA) et la pénalité GraphNet. Nous avons appliqué notre modèle à l'ensemble de données multi-omiques TCGA-LGG (The Cancer Genome Atlas -Low Grade Glioma). Contrairement à Elastic-Net seul, la pénalité GraphNet est capable de sélectionner un ensemble raisonnable de gènes tout en offrant une interprétation biologique au niveau des voies biologiques et donc informative sur de potentielles cibles thérapeutiques. L'exemple sur l'ensemble de données TCGA-LGG montre la stabilité et la fiabilité de netSGCCA pour la sélection des variables d'intérêt.

Enfin, nous avons utilisé netSGCCA pour intégrer la radiomique et les données génétiques et l'avons appliqué à la tâche de prédiction de la survie. En raison de l'indisponibilité des données de survie sur notre cohorte DIPG, nous avons utilisé l'ensemble de données TCGA-LGG pour mener l'étude. Nous avons comparé netSGCCA avec d'autres approches de survie multi-blocs et des modèles construits en monoblocs. La netSGCCA s'est révélée être un modèle robuste capable de sélectionner des gènes connus dans le gliome de bas-grade et d'autres interagissants dans des voies biologiques pertinentes. Le bloc radiomique n'a pas fourni d'information supplémentaire au modèle. Cependant, les caractéristiques radiomiques extraites de la modalité T2, en mono-bloc sont des prédicteurs performants, ce qui représenterait un espoir pour les patients avec DIPG qui n'ont pas de biopsie. blastomas. Our proposition combines different simple detection and segmentation models to obtain satisfactory segmentation results.

In parallel, our work aims to build a multiblock integration model taking into account the intra-block correlation structure described in established complex graphs of gene-gene interactions (e.g. PathwayCommons). Moreover, our objective is to understand how the interaction graph influences the selection of variables. We propose netSGCCA, a model combining the Sparse Generalized Canonical Correlation Analysis (SGCCA) with the GraphNet penalty. We applied our model to the TCGA-LGG dataset. Unlike Elastic-Net alone, the GraphNet penalty was able to select a reasonable set of genes and gives an informative biological interpretation using biological pathway enrichment analysis. The example on the TCGA-LGG dataset shows the stability and reliability of netSGCCA for selecting variables of interest.

Finally, we used netSGCCA to integrate radiomics and genetic data and applied it to the survival prediction task. Due to the unavailability of survival data on our DIPG cohort, we used the TCGA-LGG dataset to conduct the study. We compared the results obtained with netSGCCA with other multi-block survival approaches and models built in mono-blocks. The netSGCCA has proven to be a robust model capable of selecting variables already linked to the pathology studied and interacting in relevant biological pathways. The addition of imageing did not improve the predictive ability of netSGCCA. However, the baseline results show that the radiomic features extracted from the T2 modality can be strong predictors.

I want to start my acknowledgements by thanking Doctor Jacques Grill and his team at the Gustave Roussy Institute. I deeply respect your commitment to paediatric glioma patients and your interest in various research angles. This thesis would not have been possible without your collaboration and the Imagine for Margo association, which funded this work.

I am also thankful to my thesis supervisors Vincent Frouin and Cathy Philippe. You gave me a unique opportunity to work on an interesting subject. The feedback I received always helped me improve and grow. Furthermore, I would like to thank all former and current members of the Brainomics team, who welcomed me among them and with whom conversation was always interesting.

I would also thank Neurospin and all its members for allowing me to work in such a diverse research laboratory.

I would also like to thank the people that evaluated my PhD work: Céline Lefebvre, Laurent Guyon, Arthur Tenenhaus, Irène Buvat and Anaïs Baudot. Thank you for taking some of your precious time to review my work. Both the questions raised and the discussion were very insightful. It is a deep honour to have you as my jury members. I would also acknowledge the profound impact of those around me during my whole life, especially my thesis preparation. My parents' support is critical for who I am today. I conclude by thanking my partner, Lina. Your curiosity towards my work and patience with me was a source of motivation. i To focus on the most interesting region, we only plotted the precision-recall scores for thresholds between 0.1 and 0.9. From the left to the right, using the FLAIR, using the 1 Recovering performances depending on configurations defined by the different cases defined by the vector u 2 and graphs. Corr is the correlation between the estimated components.

List of Figures

Precision, Recall and F1 correspond to the evaluation of u 2 against the computed weights.

Bold refers to highest values between netSGCCA and SGCCA. Low mean to variance ratio (c = 0.5). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2 Recovering performances depending on configurations defined by the different cases defined by the vector u 2 and graphs. Corr is the correlation between the estimated components.

Precision, Recall and F1 correspond to the evaluation of u 2 against the computed weights.

Bold refers to highest values between netSGCCA and SGCCA. High mean to variance ratio

6.4

The effect of pruning edges that connect genes selected when using the full PC graph.

Candidate and stable genes sets obtained when using the different pruned sub-graphs. . . we only kept the seven where the model did not fail. The table only shows radiomic features that were selected more than three times. The means and std were computed using only the samples where the feature was selected. . . . . . . . . . . . . . . . . . . . . . . . . . .

7.3

Selected radiomic features extracted from the FLAIR modality. The table only shows radiomic features that were selected more than 6 times. The means and std were computed using only the samples where the feature was selected. . . . . . . . . . . . . . . . . . . .

7.4

Selected radiomic features extracted from the CE T1 modality. The table only shows radiomic features that were selected more than 3 times. The means and standard deviation were computed using only the samples where the feature was selected. . . . . . . . . . .

7.5

Selected radiomic features extracted from all modalities. Out of 10 bootstrap samples, we only kept the eight with a convergent model. The table only shows radiomic features that were selected more than 3 times. The means and std were computed using only the samples where the feature was selected. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7.6

Survival prediction results obtained on the 419 patients from the TCGA-LGG dataset using mutation features. Validation results were obtained on three-fold cross-validation.

Test results were obtained using 10 bootstraps. Elastic-Net Cox was used . . . . . . . . . 7.7 Selected genes using the mutation matrix propagated in the PC graph. The table shows genes that were selected more than 5 times. The means and std were computed using only the samples where the feature was selected. . . . . . . . . . . . . . . . . . . . . . . . . . . 7.8 Selected genes using the mutation matrix propagated in the influence graph. The table shows genes that were selected more than 5 times. The means and std were computed using only the samples where the feature was selected. . . . . . . . . . . . . . . . . . . . . 7.9 Enrichr Enrichment results when using the five most selected genes using the influence applied on RNA data. The table shows features that were selected more than 5 times. The means and std were computed using only the samples where the feature was selected. . .

7.19

Enrichment results when using the most selected genes from the RNA block using netSGCCA.

A.1 Selected mutated genes using the raw mutation matrix. Out of 10 bootstrap samples, we only kept the 6 where the model did not fail.

Introduction

Gliomas in the WHO classification. Central Nervous System (CNS) tumours are neoplasms located in the brain or spinal cord tissue. Among them, gliomas have the distinctive feature of developing from a glial cell (astrocyte, oligodendrocyte). Glial cells support and protect nerve tissue by providing nutrients and oxygen to neurons; they produce the myelin sheath that makes the transmission of the electric signal or nerve impulse efficient. Gliomas constitute approximately 30% of all CNS tumours and 80% of all malignant brain tumours (Goodenberger and Jenkins, 2012). Until the 2007 edition of the World Health Organisation (WHO), CNS classification for operable tumours mainly relied on the histological comparison between the tumoural cells and their putative original cells [START_REF] Louis | The 2016 world health organization classification of tumors of the central nervous system: a summary[END_REF][START_REF] David N Louis | The 2021 WHO classification of tumors of the central nervous system: a summary[END_REF]. The WHO further categorised CNS tumours by their grade, reflecting the tumour aggressiveness, thus the patient's clinical evolution.

According to the 2007 WHO CNS classification of tumours, gliomas included Low-Grade Gliomas (LGG) and Glioblastomas. Low-Grade Gliomas comprised diffuse low-grade and intermediate-grade gliomas, graded II and III by the WHO. These tumours generally impact young, otherwise healthy patients and have an indolent course.

LGGs have a longer survival rate in comparison with highergrade gliomas. The lower grade gliomas are characterised by the IDH mutant and 1p/19q co-deletion [START_REF] Forst | Low-Grade Gliomas[END_REF]TCGA, 2015;[START_REF] Wong | Integrated proteomic analysis of low-grade gliomas reveals contributions of 1p-19q co-deletion to oligodendroglioma[END_REF]. On the other hand, Glioblastomas are classified as grade IV. Newly diagnosed patients with Glioblastoma have a median survival time of one year, with poor responses to treatments. This tumour is characterised by alterations in the pathways p53, Rb, receptor tyrosine kinases (RTK) and phosphoinositide 3-kinase (PI3K) signalling, among others [START_REF] Brennan | The somatic genomic landscape of glioblastoma[END_REF][START_REF] Mclendon | Comprehensive genomic characterization defines human glioblastoma genes and core pathways[END_REF].

Since the early 2010s, several research groups reported differences among the Grade IV gliomas, whether they affected adults (adult High-Grade Gliomas -aHGG) or children (pediatric High-Grade Glioma -pHGG). For example, the diagnostic from histopathology was poorly reproducible in pHGGs because they are heterogeneous [START_REF] Gilles | Pathologist interobserver variability of histologic features in childhood brain tumors: Results from the ccg-945 study[END_REF], and the therapeutic answer was very different for a supposed identical sub-type. The location of some specific pHGG in the brain is more often thalamic and infra-tentorial tumours which distinguished them from the aHGG [START_REF] Puget | Mesenchymal transition and pdgfra amplification/mutation are key distinct oncogenic events in pediatric diffuse intrinsic pontine gliomas[END_REF].

The 2016 edition of the WHO CNS classification introduced a paradigm shift in the diagnosis of CNS neoplasms. After the recommendations of the consortium cIMPACT which called for a better consideration of the molecular characteristics currently available, both histologic features and genetic alterations were incorporated into the diagnostic framework, classifying and grading brain tumours [START_REF] David N Louis | The 2021 WHO classification of tumors of the central nervous system: a summary[END_REF]. The WHO CNS classification has had a very recent update in 2021. Under the current version, adult and paediatric HGGs are divided based on an increasing emphasis on molecular markers [START_REF] Gaillard | Diffuse glioma[END_REF]. More specifically, four general groups of diffuse gliomas are recognised in the 2021 WHO classification: 1) adult-type diffuse gliomas, 2) paediatric-type diffuse low-grade gliomas, 3) pediatric-type diffuse high-grade gliomas, and 4) circumscribed astrocytic gliomas [START_REF] Osborn | The 2021 world health organization classification of tumors of the central nervous system: What neuroradiologists need to know[END_REF].

Pediatric Diffuse High-Grade Gliomas. Diffuse Intrinsic Pontine Glioma (DIPG) is a rare brain tumour located in the pons, mostly found in children between 5 and 7 years of age. It is considered one of the most aggressive paediatric tumours, with a survival rate of less than 10% beyond two years after diagnosis [START_REF] Paul | A clinicopathologic reappraisal of brain stem tumor classification[END_REF] and a median overall survival below one year [START_REF] Kenneth | Diffuse intrinsic pontine gliomas-current management and new biologic insights. is there a glimmer of hope?[END_REF]. The DIPG is categorised as a diffuse midline glioma (DMG), characterised mainly by a K27M mutation in genes coding for the histone H3 protein and/or a loss of H3K27 trimethylation through EZHIP protein over-expression [START_REF] Castel | Histone h3f3a and HIST1h3b k27m mutations define two subgroups of diffuse intrinsic pontine gliomas with different prognosis and phenotypes[END_REF]. The location of the tumour and its corresponding genomic alteration make the DIPG a particular type of tumour from other High-Grade Gliomas (HGGs) [START_REF] Louis | The 2016 world health organization classification of tumors of the central nervous system: a summary[END_REF].

DIPG grows fast. Its symptoms include eye problems, trouble with walking, muscle coordination or balance, and weakness in the arms or legs. Diagnosis of DIPG principally relies on Magnetic Resonance (MR) Images, but they are insufficient as the tumour can be confused with other similar tumours. Biopsies can be safely used for diagnosis [START_REF] Gupta | Prospective feasibility and safety assessment of surgical biopsy for patients with newly diagnosed diffuse intrinsic pontine glioma[END_REF]. However, it is rarely done

and has yet to be incorporated into the standard diagnosis protocol.

Currently, there is no curative treatment for patients with DIPG. Removing the tumour with surgery is not feasible due to its location. The tumour is in the brainstem, which helps control essential functions. The tumour can not be removed without damaging vital brain tissue. Most patients are treated with radiation therapy, which has been shown to increase survival and improve symptoms temporarily in most patients [START_REF] Gallitto | Role of radiation therapy in the management of diffuse intrinsic pontine glioma: A systematic review[END_REF]. Chemotherapy has been found challenging due to the blood-brain barrier. The blood-brain barrier is a highly selective semipermeable border that prevents solutes in the circulating blood from non-selectively crossing into the extracellular fluid of the central nervous system. However, some chemotherapy treatments are currently being examined in clinical trials [START_REF] Ho | Developing chemotherapy for diffuse pontine intrinsic gliomas (dipg)[END_REF].

Multiple clinical trials have been performed to investigate DIPG therapies (Rechberger et al., attributed to higher magnetic field strengths, using 3T magnets instead of 0.5T and 1.5T several years ago. [START_REF] Cantini | Benchmarking joint multi-omics dimensionality reduction approaches for the study of cancer[END_REF][START_REF] Herrmann | Large-scale benchmark study of survival prediction methods using multi-omics data[END_REF]. Unlike traditional machine-learning models that account for each block contribution separately, multi-block models allow taking into account the information shared among blocks and the interactions between variables across the blocks [START_REF] Philippe | Analyse intégrée de données de génomique et d'imagerie pour le diagnostic et le suivi du gliome malin chez l'enfant[END_REF].

The methodological work presented here is devoted to subtype 3 in WHO classification of CNS tumours, which we will name DIPG in this manuscript from now on. It aims at providing analysis tools of molecular and biomedical imaging profiles to inform treatment decisions to control the disease.

Introduction

Because of data sparse availability, our developments will leverage measurements in other subtypes when the category of measurements is similar across subtypes. The most salient point about these data is the small number of samples available, which, no matter how much effort is made, will remain limited given that this is a rare disease.

Open questions raised by multi-view datasets in rare tumours

Data about DIPG is scarce, and knowledge about the disease is limited. This situation can, however, be mitigated by the availability of a large amount of public data. These are of two kinds.

Firstly, there are data on other gliomas that are similar in some respects to DIPGs and contain the same types of measurements collected for BIOMEDE. These include data from the iconic The Cancer Genome Atlas (TCGA) 3 and The Cancer Imaging Archive (TCIA) 4 , projects which offer molecular and imaging data on many patients. These public datasets have the advantage of having a large number of available samples and extensive literature about them. Therefore, they can be used to study the behaviour and limitations of the methods used, and the results can be compared with existing literature. Additionally, while we know that DIPG is a unique tumour type, there exists, nevertheless, similarities with other gliomas that can be exploited to learn about the disease. The second kind of publicly available data is the curated and extensive knowledge about molecular data.

One can cite the Molecular Signature Database (MSigDB) 5 of the Broad Institute, which compiles collections of cancer signatures in the form of gene lists, and the PathwayCommons 6 , which compiles 1.1. Open questions raised by multi-view datasets in rare tumours 5 all the interactions described between different molecular entities in the form of interaction graphs.

These data can be used as prior knowledge of our methods and help the interpretation of the results.

It should also be noted that data scarcity is sometimes due to shortcomings in pre-processing methods adapted to data acquired under routine clinical conditions. The systematic and rational extraction of standard information from MR images requires numerous processing steps, such as the co-registration of images acquired with different sequences, their normalisation, harmonisation and quantification. Most of these steps have been implemented in clinical neuroscience work that operated on images without mechanical deformation, such as that caused by tumours. It is, therefore, necessary to revisit these approaches, and we have done this work on the data considered in our work. In contrast to these methods, whose adjustment does not constitute original work, the automatic delineation of DIPG/DMG type tumours constitutes an open research challenge. Indeed, the pHGG tumour cases are absent from the data compendiums used to train the models for automatic tumour delineation.

Overall, the presentation of the DIPG calls for methods to fully potentiate the images to generate the radiomics information we wish to study. This presentation of data also calls for models capable of integrating the different types of data and training paradigms that allow to take advantage of external information: these paradigms can consist in transfer learning -where a model trained on external data is only refined -or consist in penalising (or orienting) existing statistical models.

Non-black-box models are desired. In the data setting described above (collected from clinical trials combined with public data), statistical models are not only used to predict the outcome of diseases. Instead, they are expected to help to characterise the disease by finding a set of variables of interest explaining their outcome.

However, isolated variables associated with a pathology outcome usually do not, by themselves, give a biological meaning to this association. For example, at the gene scale, somatic mutation profiles from patients with the same pathology may differ from one to the other because different genes are mutated in the same pathway [START_REF] Lawrence | Discovery and saturation analysis of cancer genes across 21 tumour types[END_REF][START_REF] Le Morvan | Netnorm: Capturing cancer-relevant information in somatic exome mutation data with gene networks for cancer stratification and prognosis[END_REF][START_REF] Wood | The genomic landscapes of human breast and colorectal cancers[END_REF]. This leads to the need to identify networks and pathways that group variables that interact in complex patterns.

Several works have proposed the usage of complex graphs to identify sub-networks of variables of interest [START_REF] Zhang | Network-based machine learning and graph theory algorithms for precision oncology[END_REF]. Some works suggested using graphs as a post-analysis tool to identify the interactions between previously selected variables [START_REF] Kim | Identifying causal genes and dysregulated pathways in complex diseases[END_REF]Vandin et al., 2011a;[START_REF] Charles | Inference of patient-specific pathway activities from multi-dimensional cancer genomics data using PARADIGM[END_REF]. Others used the graphs as a pre-processing tool, such as smoothing variables over the interaction network [START_REF] Hofree | Network-based stratification of tumor mutations[END_REF][START_REF] Le Morvan | Netnorm: Capturing cancer-relevant information in somatic exome mutation data with gene networks for cancer stratification and prognosis[END_REF][START_REF] Rapaport | Classification of microarray data using gene networks[END_REF]. Graphs have also been integrated into statistical models as a penalty over the model parameters. This includes supervised models, such as survival models [START_REF] Zhang | Network-based survival analysis reveals subnetwork signatures for predicting outcomes of ovarian cancer treatment[END_REF] and regression models [START_REF] Li | Network-constrained regularization and variable selection for analysis of genomic data[END_REF], and unsupervised models, such as matrix factorisation models [START_REF] Zhu | Network-based cancer genomic data integration for pattern discovery[END_REF].

Overall, the expectations regarding statistical models used in clinical research projects concern the interpretability of these models and their ability to propose molecular hypotheses for treatment pathways or operating modes. This implies, for example, providing existing multi-block approaches with penalties so that the models selected by the learning procedure offer interpretable solutions: for example, we can try to constrain that the set of gene expression variables selected are "close" in the PathwayCommons graph and therefore belong to a same biological pathway. To contribute to the Introduction expected interpretability properties, graphical penalisation is scheduled on at least one block of the multi-block models. This penalty on models must scale to accommodate the block with a large number of variables, and it is necessary to have a framework to evaluate their contribution.

To provide interpretations even if few clinical annotations are available in a cancer sample, it is necessary to consider statistical models that are more than mere classifiers but also survival models.

Indeed, censored survival data is generally available. Yet it implies adapting the training framework of statistical models.

Finally, the question of the unitary and joint contribution of the different blocks arises in the heavy and complex experimental setting of collecting molecular and imaging data. It is necessary to evaluate these contributions under conditions of data scarcity. In addition, the question arises of the transfer learning approaches that could be constructed to move from one imaging + multi-omics dataset to another.

Manuscript outline

This work aims at studying the integration of radiomics with omics data in a multi-block analysis framework. Many issues are raised by the multi-block integration in DIPG. We addressed a few of them, highlighted in the previous paragraph. To achieve these tasks, we do not only consider the specific DIPG dataset (as it will be made available by the BIOMEDE clinical research). Instead, we consider several datasets of interest, namely LGG, Glioblastoma and DIPG. In some tasks, we used explicitly transfer learning, but in other, we implemented and evaluated methods to anticipate complete data annotation availability. Since radiomics are extracted from regions of interest in the images, in our case, the tumoural area, this work discusses brain tumour segmentation approaches.

The first part of this manuscript presents the basic definitions used and the related works. Then, it introduces the different datasets and the various preprocessing methods applied.

The second part is organised into three chapters. Chapter 5 presents our brain tumour segmentation approach, already published (Chegraoui et al., 2021a) that is demonstrated on DIPG data. In Chapter 6, we lay out our multi-block integration approach penalised with a graph a priori. The method is assessed on simulations and TCGA-LGG (TCGA, 2015) data. Finally, in Chapter 7, we study the potential of adding imaging to molecular data by applying our methods on the TCGA-LGG dataset, including imaging and molecular blocks. We show results obtained on the survival prediction with downstream analysis.

✽ ✽ ✽ ✽ ✽ ✽

Chapter 2

Notations and Definitions

In this chapter, the different notations and definitions of entities discussed in the thesis are introduced.

This chapter also includes the definition of classical methods discussed in this thesis.

Notations

Throughout this work, we use the following convention. Bold uppercase letters represent matrices, such as A. Bold lowercase represents vectors, such as a. Scalars are represented by lowercase letters, such as a. We denote with a star an optimal parameter (or vector) as a ⋆ . Data is represented in matrices, often denoted X, with n observations (rows) and p variables (columns). These data matrices are also called blocks. In this work, we are interested in the analysis of J blocks, which we denote X (1) , • • • , X (J) and their associated vectors are denoted w (1) , • • • , w (J) (β can also be used in some cases).

We define a graph as G = (E, V), where E is a finite set of nodes and V the set of edges. Unless specified otherwise, the graph is undirected. For sake of simplicity, the elements e 1 , e 2 , • • • of E can be referred to as 1, 2, • • • . The order of the elements is arbitrary but consistent. To the graph G, a binary adjacency matrix A is associated, where A i,j = 1 if the nodes i and j are connected, otherwise A i,j = 0. Each element of the graph has a degree counting its number of neighbours, which can be represented in a diagonal matrix D. The laplacian L associated with the graph G is defined as

L = A -D.
We denote the norm operator ∥ • ∥. The ℓ 1 norm of a vector β is defined as

∥β∥ 1 = |β i |. The norm ℓ 2 is defined as ∥β∥ 2 = β 2 i .
Let K positive-defined matrix, the norm ∥β∥ K = β ⊤ Kβ. For a real matrix M we define the Frobenius norm noted ∥M ∥ F = trace(M ⊤ M ). Recall that the trace function returns the sum of diagonal entries of a square matrix.

Data description

In this thesis, we are interested in two sources of data. First, imaging data was acquired through MR scans and described using radiomic features. Second, biological samples were obtained through biopsy, from which genetic data could be extracted. Here, we are mainly interested in gene expression, somatic mutation profiles and copy number variants (CNV).

The following sections will describe the different data used, starting with imaging and then genetic data.

Imaging Data

Images are acquired through an MRI scanner composed of a magnet, gradient coils, and radiofrequency coil. The magnet generates a magnetic field with strength in the order of a few teslas.

The gradient coils are responsible for spatially localising the signals. Finally, the radio-frequency coils generate and send radio-frequency signals in the order of a few microteslas. The scanner exploits the intrinsic spin property of protons to generate the images. MRIs mainly focus on protons in hydrogen in the water molecules. First, the magnetic field aligns the spin of the protons with its direction.

Radio frequencies are sent towards the protons, via the radio-frequency coil, which knocks them out of alignment. As the spin realigns with the magnetic field, they generate an electrical signal which fades through time due to the desynchronisation of spins. The decay of this electrical signal is characterised by its T2 relaxation time and is tissue specific. We can then measure the contrast between the decay-ing electrical signals at an operator-chosen Echo Time (TE). In parallel, as the protons realign with the magnetic field, their magnetism increases. This growth rate is tissue-specific and characterised by a T1 relaxation time. We can then measure the contrast between the increasing magnetism at an operator-chosen Repetition Time (TR). It should be noted that these two signals cannot be measured simultaneously. Acquiring the signal at TR gives the contrast between the different tissues and results in a T1-weighted image.

Choosing TE and TR results in different MR sequences (or modalities), which capture various characteristics of the tissue organisation (at a millimetric or submillimetric scale) or functional parameters. These sequences differ by how they contrast the different tissues or map functions in tissues.

Four modalities of MR images are commonly used for clinical diagnosis: The T1-weighted (T1w), the T2-weighted (T2w), Fluid Attenuated Inversion Recovery (FLAIR) and Contrast Enhanced T1weighted (CE T1w). The T1-weighted images are obtained with a short TE and a TR similar to the T1 relaxation time. The T1-weighted modality highlights fatty tissues and is best for observing anatomical structures such as the grey and white matter in the brain. CE T1w scans are T1w scans after the patient's infusion of a contrasting agent, usually, Gadolinium (Gad), a non-toxic paramagnetic agent. CE T1w images are especially useful in looking at vascular structures and breakdowns in the blood-brain barrier. The T2-weighted images are acquired with a TE similar to the T2 relaxation time and a long TR. T2-weighted images are bright on the fat tissues and fluids. The FLAIR acquisition parameters are similar to the T2-weighted. However, additional pulses are sent to the protons in the acquisition phase to nullify the signal from fluids. The contrast in the FLAIR images is similar to the T2-weighted, but the intensity of normal fluids is attenuated. Both T2-weighted and FLAIR sequences are utilised for tumour location. The GLCM counts the number of times two intensities are neighbouring each other. Formally, let M be the GLCM, M i,j is the number of times a voxel of intensity i is at maximum a distance δ from a voxel of intensity j. This matrix is normalised to reflect the joint probability distribution of image intensities. From the GLCM, 24 features are extracted, including Autocorrelation, Contrast and Joint Average.

Radiomics

The GLRLM counts the number of same intensity level sequences by length. Formally, let M be the GLRLM, M i,j is the number of sequences of intensity level i and length j occurring in the ROI.

From the GLRLM, 16 features are extracted, including Short Run Emphasis, Long Run Emphasis and Gray Level Non-Uniformity.

The GLSZM quantifies connex grey level zones with the same intensity level. Let M be the GLSZM, M i,j is the number of connex components of size j of intensity level i. Two voxels of the same intensity are neighbouring if they are adjacent in the image. From the GLSZM, 16 features are extracted, including Small Area Emphasis, Large Area Emphasis and Gray Level Variance.

The GLDM quantifies grey-level similarity in an image. For a parameter δ, two voxels are similar if their intensity levels g i and g j verify |g i -g j | ≤ δ. Formally, let M be the GLDM, M i,j is the number of times a voxel of intensity i has j similar voxels at a distance δ. From the GLDM, 14 features are extracted.

The NGTDM quantifies the difference between an intensity level value and the average intensity value of its neighbours. For each voxel of intensity level i, we compute the absolute value of the difference between i and the average intensity level of its adjacent voxels. Then, these differences are summed across all voxels of the same intensity level (s i ). Additionally, the normalised score p i = s i Σs j is computed for each intensity level. The two scores are used to calculate 5 radiomic features.

An illustration of the different texture matrices is given in figure 2.3. Multiple software can be used to extract the radiomic features such as LifeX [START_REF] Nioche | LIFEx: A Freeware for Radiomic Feature Calculation in Multimodality Imaging to Accelerate Advances in the Characterization of Tumor Heterogeneity[END_REF] and pyRadiomics [START_REF] Joost | Computational Radiomics System to Decode the Radiographic Phenotype[END_REF]. For the GLCM we only considered voxels with voxels in the right and left as neighbouring. Note that this matrix is symmetric, thus voxel couples are counted twice. For the GLRLM, we also used the same definition of the neighbourhood as for the GLCM. For all other matrices, the neighbourhood is given by the infinity norm. The radiomic features presented are only representative examples.

Genetic data

A biological sample from the patient can be extracted through a biopsy. This sample can then be analysed to obtain different generic profiles.

RNA sequencing

Gene expression measures the relative abundance of messenger RNA (mRNA) in a sample prepared from the tumour cells. This technique measures the expression level of a gene transcribed into an mRNA (in the form of one or several splicing variants for a given gene). Raw read counts correspond to the number of times an mRNA sequence was seen in the biological sample. To measure gene expression, RNA-seq is usually used. First, the mRNAs are isolated and then converted to a more stable double-stranded DNA and amplified. A sequencer reads the sequence from each fragment, which is then mapped to a gene and counted.

The read counts are said to be relative since they cannot be used to compare expression levels between samples due to the need to account for differences in transcript length and the total number of reads per sample. Multiple RNA-seq normalisations have been proposed. In this thesis, we used the Transcripts Per Million (TPM) normalisation. For a sample i and a gene j, the TPM matrix is defined as follows :

T P M i,j = q i,j /l j k q i,k /l k 10 6 (2.1)
where q i,j is the raw read count for the gene j of the sample i, and l j is the transcript length of gene j.

Somatic mutations

Somatic mutations are mutations observed in somatic cells only, as opposed to germline cells. These mutations can not be passed on to descendants and occur after conception. Mutation profiles do not allow for the identification of the subjects because they are obtained by subtracting a patient's genome from the genome observed in the tumour cells. In this thesis, somatic mutation profiles are represented as a matrix of n subjects and p genes where X i,j = 1 if a mutation occurred in the gene j for subject i.

Copy number variation (CNV)

Copy number variation (CNV) is the estimation of the variation in the copy number of genomic segments. The CNV is given by the ratio of the estimated concentration of the target gene to the estimated concentration of the reference gene. Sequencing can be used to detect genetic variants, which include CNVs.

Method Definitions

Machine learning (ML) is devoted to choosing and building mathematical models whose parameters are adjusted according to sets of available data. We define a few general terms used in ML that we will utilise in our work. Most statistical models define an objective function f (X; β), where X are the observed data, and β are the model parameters. The aim is to find the optimal set of parameters that minimises (or maximises) f . Machine learning approaches are traditionally divided into two categories.

The supervised models are used when a desired outcome to predict is available. The model learns a function g that maps the input X to the output Y . The unsupervised models have no labels given to the learning algorithm, leaving it on its own to find structure in its input.

This section defines several classical objective functions and models which are discussed in our work. We also present regularisation techniques that add constraints to the models. These constraints help the models discard unwanted noise and prevent it from closely fitting the observed samples without the ability to generalise on new observations. Finally, we discuss some metrics used to assess the model performances, allowing us to compare different models objectively.

Objective functions

Over the years, multiple objective functions (also called loss functions) have been defined. In this section, we describe the most commonly used objective functions.

Least-Square

The Least-Square objective function aims to find the best parameters for a mapping g(X; β), which will minimise the sum of squared errors between the observed data and the predicted data. The observed data come here as a couple (X; Y ). This loss function is formulated as follows:

f ((X; Y ), β) = ∥Y -g(X; β)∥ 2 2 (2.2)
Note that the least-square is only for supervised models. g can be a generative model that learns the best models that generate observed data X (In this case, the input X is also the approximated output). This loss function is generally used for regression or image denoising problems.

Log-likelihood

The likelihood f (X; β) is the probability of the realisation of the observed data X considered as a function of β. If the observed samples are assumed to be independent, then the likelihood function is the product of the likelihood of each sample. The log-likelihood is often used to transform the product into a sum. The statistical model aims at maximising the log-likelihood. The log-likelihood is often used for classification and survival models.

Regularisation

Three main possible benefits can arise from adding regularisation over the parameters of a statistical model. First, it can be added to ill-defined objective functions. For example, the loss function f

does not have a defined local or global minimum. Adding regularisation limits the search space and thus guarantees the existence of a solution. Additionally, regularisation can be used to prevent the statistical model from over-fitting the data, thus capturing unwanted noise and failing to generalise

Notations and Definitions

on new observations. Finally, regularisation helps retain a subset of the predictors, which helps interpretation.

In order to restrict the search space for the parameters or impose regularity over them, penalties can be added to statistical models. Three main penalties are widely used. First, the Ridge penalty [START_REF] Hoerl | Ridge regression: Biased estimation for nonorthogonal problems[END_REF], which introduces an ℓ 2 penalty over the model parameters, which helps shrink them. This penalty helps solve the multi-colinearity in the data and improves parameter estimation. The Least Absolute Shrinkage and Selection Operator (LASSO) [START_REF] Tibshirani | Regression Shrinkage and Selection Via the Lasso[END_REF] is another popular machine learning regularisation technique. This penalty introduces an ℓ 1 norm over the model parameters to promote sparsity, enhancing the model interpretability. The LASSO selects a few representatives among highly correlated variables in the observed data. Finally, ElasticNet [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF] is the linear combination of the Ridge and LASSO regularisations. The Elastic-Net performs variable selection while benefiting from the smoothness of the estimates introduced by the Ridge penalty.

Classical models

Supervised models

Having observed data or covariates X and target data Y , a supervised model finds the optimal parameters that map the covariates to the target. In this section, we present some commonly used supervised models.

Support Vector Machines Support vector machines (SVM) [START_REF] Cortes | Support-vector networks[END_REF] are classical ML methods. Given a set of labelled samples divided into two categories, an SVM training algorithm builds a model that assigns new samples to one category or the other by finding an optimal hyperplane that separates the two classes. The hyperplane defines a decision boundary between the different classes. Kernels can be used to learn high-dimension feature space. The SVM were initially introduced for two-classes classification tasks, but they can be extended to regression and survival tasks.

Partial Least Squares Partial Least Squares (PLS) [START_REF] Kramer | Chemometric Techniques for Quantitative Analysis[END_REF] is a linear regression model that maps the input X and target data Y by projecting both matrices into a new space. The underlying model is formulated as follows:

X = T P ⊤ Y = U Q ⊤ (2.3)
where P and Q are orthogonal matrices of loadings and T and U are the projections of X and Y respectively. The model aims at maximising the covariance between the scores T and U .

Neural networks Neural Networks [START_REF] Warren | A logical calculus of the ideas immanent in nervous activity[END_REF] are non-linear supervised ML models. Neural Networks are constructed as chains of functions composed together. For example, we might have three functions f (1) , f (2) , and

f (3) connected in a chain to form f (x) = f (3) (f (2) (f (1) (x)))
. The function f is a three-layers neural network, and each f (i) represents a network layer. Generally, the functions f (i) take the form f (i) (x) = g(W i x + b i ), where W i and b i are the layer parameters and g a non-linear function such as the Hyperbolic tangent.

Multiple Neural Network classes have been proposed, and we are especially interested in Convolutional Neural Networks (CNNs). CNNs are specialised neural networks for processing data with a known, grid-like topology of the input data, such as images. A CNN learns the parameters of a set of filters in each layer, which are used to convolve the input and pass its result to the next layer.

After passing through convolutional layers, the input data becomes abstracted to feature maps. Additionally, neural networks can be built as Encoder-Decoders, a Neural Network structure that learns to reconstruct the input data. Internally, it has a hidden layer representing the latent space from which the input data were generated. These structures are commonly used for image segmentation and dimensionality reduction.

Neural Networks constitute a flexible family of models that can be designed and adapted to a wide range of tasks, including regression, classification, segmentation and survival prediction. However, complex and deep Neural Networks are associated with a large number of parameters that must be estimated. Therefore, a populous dataset is required in order to meaningfully learn using Neural Networks.

Survival models

Survival analysis is the estimation of the expected duration of time until one event occurs. For example, survival analysis in the oncological context includes the estimation of the duration from the initiation of treatment to the occurrence of disease progression or death (progression-free survival).

Generally, instead of trying to estimate the event time directly, most survival models estimate a hazard function. A hazard function relates the passage of time and other covariates with the probability of an event occurring. Unlike typical regression problems, in which all training output variables are known, survival problems differ in that only some of the events in the training dataset are observed while the others are censored. In particular, several patients might not have a death date but only the last consultation date.

Cox-model

The Cox proportional hazard [START_REF] Cox | Regression models and life-tables[END_REF]) models the hazard function as the multiplication of a time-dependent baseline function (λ 0 (t)) and an exponential risk function which depends only on the individual's covariates. The hazard function (associated with the risk r), at a time τ and individual's covariates x, is given as:

h(x, τ ) = λ 0 (t)exp(x T β) r(x) = exp(x T β) (2.4)
The parameters β can be estimated without considering the time factor. Given n subjects with covariates x 1 , • • • , x n and survival times (t 1 , c 1 ), • • • , (t n , c n ) with t i being event times and c i being censoring times, we want to maximise the probability that individual j has an event at τ j = min(t j , c j ) compared to the other individuals still at risk at τ j . We note R j = i ∈ [1, n]|τ i ≥ τ j , and we formulate our target likelihood for individual j:

Notations and Definitions

L j = exp(x j β) i∈R j exp(x i β) (2.5)
When multiple individuals have the same event time, [START_REF] Breslow | Analysis of survival data under the proportional hazards model[END_REF] proposed regrouping them as having the same risk at τ j , and the reasoning shifts from the individual to the event time. Given D j individuals having the same survival time τ j , the likelihood function becomes for τ j :

L j = j∈D j exp(x j β) i∈R j exp(x i β) (2.6)
The parameters β are obtained by maximising the log-likelihood of observing all events:

L = log L j (2.7)
When covariates exceed the number of individuals, the log-likelihood is penalised using ElasticNet (Simon et al., 2011a). Given α and r two positive numbers, with r ∈ [0, 1], the loss function is defined as:

L penalised = -L + α r||β|| 1 + 1 -r 2 ||β|| 2 2 (2.8)
Random-forest Random forest [START_REF] Tin | Random decision forests[END_REF] is a machine-learning model based on the ensembling of multiple decision trees. Each tree is built by recursively splitting the dataset according to a split criterion. Individual trees are ensured to be uncorrelated by building each ones on a different bootstrap sample of the original data and by evaluating the split criterion only for a randomly selected subset of the variables.

In the case of survival prediction, multiple splitting criteria have been proposed. We chose to use the Log-rank test. This is a non-parametric hypothesis test to compare the survival distributions of two samples. Log-rang test estimates, under the null hypothesis H0, the probability that two samples of individuals have such an extreme difference in hazard functions. The chosen variables must minimise the Log-rank probability.

Unsupervised models

Having observed data X, unsupervised models focus on learning patterns in the data without a target at hand. This section focuses on some of the most commonly used dimension reduction models.

Non-negative Matrix Factorisation Non-negative Matrix Factorisation (NMF) [START_REF] Lee | Learning the parts of objects by non-negative matrix factorization[END_REF] Data is then projected into a smaller dimension where most variance is explained. This is done by diagonalising the estimated covariance matrix (X ⊤ X), which is a symmetric semi-positively defined matrix.

Canonical Correlation Analysis Canonical Correlation Analysis (CCA) [START_REF] Horst | Generalized canonical correlations and their applications to experimental data[END_REF] is a method for inferring information from cross-correlation matrices. Having two observation matrices X (1) and X (2) , the model projects the blocks into a latent space where their correlation is maximised. Formally, the model finds the optimal projectors w (1) , w (2) that maximise the function:

corr(X (1) w (1) , X (2) w (2) ).

Metrics

The performance of statistical models can be assessed through multiple methods. Here we present the main metrics used in this thesis.

Binary classification metrics

Having a binary prediction Ŷ , we aim to assess its agreement with the binary ground truth Y . The two can be viewed as sets of positive labels. Five metrics are commonly used.

The recall measures the fraction of positive labels retrieved by the prediction and is defined as :

Recall( Ŷ ) = | Ŷ ∩ Y | |Y | (2.9)
The precision measures the fraction of correctly predicted positive labels among all predicted labels.

P recision(

Ŷ ) = | Ŷ ∩ Y | | Ŷ | (2.10)
The dice or f-score is an aggregation of the precision and recall metrics. It is defined as

Dice( Ŷ ) = 2 × | Ŷ ∩ Y | | Ŷ | + |Y | (2.11)
For binary classification, the model returns a probability map for each class. Thus, these probability maps are first thresholded (typically at 0.5) to obtain the prediction. Recall, Precision and Dice all depend on the threshold. The Area Under Curve (AUC) is a metric independent of the threshold by integrating the precision-recall curve. A 0.5 value indicates random results.

Survival metrics

The concordance index (C-index) measures the survival model ability to correctly rank the survival times based on the individual risk scores. The C-index generalises the area under the curve, considering the censored data. A higher C-index is better, with 0.5 indicating random results and 1 reflecting a perfect ranking of the individuals. Let δ i be a binary variable indicating if the survival time of 20

Notations and Definitions

individual i is non-censored (δ i = 1 if survival time is not censored, 0 otherwise), η i his estimated risk score, and T i his observed event time, the C-index is given by:

C-index = i j 1 T j <T i • 1 η j >η i • δ j i j 1 T j <T i • δ j with: η i , the risk score of a patient i 1 T j <T i = 1 if T j < T i else 0 1 η j >η i = 1 if η j > η i else 0 (2.12) ✽ ✽ ✽ ✽ ✽ ✽
Chapter 3

Related Works

In this thesis, we are interested in integrating radiomic with multi-omic data. This integration is intended to reveal associations between non-invasive data and molecular data, which can improve our knowledge of the pathologies at hand. As a first step, it is necessary to position our project in the context of its related works.

We start by reviewing the various studies involving radiomic analysis. We are particularly interested in the tumours we want to study, namely Lower-Grade Gliomas, Glioblastomas and DIPG.

Since radiomic feature analysis relies on a precise Region of Interest identification, we will review brain tumour segmentation methods proposed in the literature. We are interested in the best-performing methods in international challenges. Later, we review different methods proposed to integrate multisource data, distinguishing between supervised and unsupervised approaches. Finally, we review several strategies to help understand and interpret the built models.

Related Works

Chapter Outline Contents

Radiomics

Several studies have shown the effectiveness of radiomic features in tumour analysis. These studies include but are not limited to, the classification of tumour types or subtypes and disease prognosis.

Often, the radiomic features are extracted from a region of interest in the image and then used in a machine learning or statistical model to make a prediction.

Classification of tumour types or subtypes using radiomics

Malik et al. ( 2021) studied the ability of radiomic features to distinguish between LGG and Glioblastomas. On 74 patients ( 32LGG and 42 Glioblastoma), the authors report an AUC score between 0.72 and 0.96, depending on the feature selection method and the classification algorithm. The most recurrent selected features are texture ones, especially from the CE T1w. However, their results show a lack of stability in variable selection: most features are chosen only by one method and on less than 70% of the runs.

Kong et al. ( 2019) studied the O 6 -methylguanine-DNA methyltransferase (MGMT) methylation status in patients with primary gliomas as classification target. On 107 patients (59 methylated and 48 unmethylated), the authors used radiomic features extracted from positron emission tomography (PET) images to classify the methylation status. They report an AUC of 0.86 on the test set.

Comparatively, the clinical data achieved an AUC of 0.69, and the combination of radiomic and clinical data gave an AUC of 0.85. Most variables of interest were first-order radiomics. to clinical and key molecular data for disease stratification. The study included 181 glioblastoma patients and reports that radiomics alone perform similarly to clinical data and better than molecular data. Combining the information from multiple sources reduced the prediction error compared to each single-source model. The highest accuracy model was built using data from all sources (i.e., clinical + molecular characteristics + radiomic features). Their approach identified 8 radiomic features that were the most important for predicting the outcome. These features included texture features (6/8 features derived from FLAIR and T2w) and volumetric features (2/8 features derived from CE T1w).

Survival Prediction using radiomics

Radiomics on the glioblastoma

Radiomics on the LGGs

Radiomics have also been used on LGGs. For instance, [START_REF] Liu | A radiomic signature as a non-invasive predictor of progression-free survival in patients with lower-grade gliomas[END_REF] studied the progression free survival of LGGs using radiomic features extracted from T2w MR images and compared the results obtained with clinical and histological data. The study found that first-order features are the most relevant for the problem at hand, and these radiomic features outperformed the other data sources. [START_REF] Wagner | Radiomics of pediatric low-grade gliomas: Toward a pretherapeutic differentiation of braf-mutated and braffused tumors[END_REF] studied the ability of radiomic features to discriminate between two mechanisms of implication (gene fusion or gene mutation) of the B-Raf proto-oncogene, serine/threonine kinase (BRAF) in pediatric LGGs. The study extracted the radiomic features from the FLAIR sequence on 115 patients and reports that the most important features were texture related, extracted from wavelet-transformed images. These features were successfully used to stratify the patients with a ROC-AUC of about 0.75.

Radiomics on the DIPG

Some studies, although few in number, have examined the DIPG using radiomics. For example, Tam et al. ( 2021) studied the survival prediction using radiomics from the T1w and T2w on 153 patients.

The study isolated five features as the most relevant. These features included three features from T1w and two features from T2w. All features were intensity and texture based on the wavelet-filtered images. The authors report that radiomics outperformed the clinical model that used sex and age at diagnosis as variables (a C-index of 0.55 for radiomics and 0.51 for clinical). When clinical features were combined with radiomics, the model performance increased but not significantly over radiomics alone. The performance of the radiomics features considered separately from sequences T1w and T2w

was lower compared to their combination. Additionally, their performances were not significantly better than the random predictions, with an average C-index of 0.51 and 0.55 when using the T1w only and the T2w only, respectively.

The reviewed articles show that radiomic are promising predictors for various tasks on different tumours. However, they also suggest stability issues, probably due to the important co-linearity among radiomic variables. The studies do not agree on the most critical modalities or predictive features.

No two studies have reported the same variables of interest. We note, however, that these studies are done in comparatively small studies, and further studies must be done to confirm the results.

Reproducibility of radiomics

Despite the successes of radiomic analyses in oncological studies, several factors influencing the stability and reliability of the feature extraction procedure have been identified in the literature [START_REF] Van Timmeren | Radiomics in medical imaging-"how-to" guide and critical reflection[END_REF]. For MR scans, these factors include image acquisition, resolution, reconstruction, preprocessing, and region of interest (ROI) delineation. [START_REF] Um | Impact of image preprocessing on the scanner dependence of multi-parametric mri radiomic features and covariate shift in multi-institutional glioblastoma datasets[END_REF] studied the effect of the MRI field strength and manufacturer on texture radiomics.

The authors studied 50 patients from The Cancer Genome Atlas (TCGA)-GBM. These patients met the following inclusion criteria: newly diagnosed GBM, pre-surgery status, and availability of multimodal MRI, including the T1w, FLAIR and CE T1w. The authors compared 32 scans acquired with a 1.5T with 13 with 3T scanners (they discarded other magnetic field strengths). Using the Wilcoxon test, the authors found that around 14% of GLSZM features and 7% of GLCM features were manufacturer dependent; and around 28% of GLSZM and 20% of GLCM features were field-strength dependent. The authors conclude that there is a bias introduced during the image acquisition. [START_REF] Molina | Lack of robustness of textural measures obtained from 3d brain tumor mris impose a need for standardization[END_REF] investigated the effects of the MR image resolution, slice thickness and dynamic grey-level range on 16 texture radiomic features extracted from Glioblastoma scans. The study was done on 20 patients who had T1w pre-treatment scans available. For each patient, the authors compared two matrices of the same slice: 432x432 (raw matrix) and 256x256 (standard MRI matrix size). They compared the features using the coefficient of variation, a standardised measure of dispersion defined as the ratio between the standard deviation and the mean of a series of data. The authors found that the investigated radiomic features are not robust to resolution and slice thickness.

Only the entropy extracted from the co-occurrence matrix resisted the robustness test regarding the dynamic range.

The dependence of texture indices with the image pre-processing (multiple sequence image alignment, intensity normalisation, level binning and ROI generation) was studied on about 30 patients with DIPG by [START_REF] Goya-Outi | Computation of reliable textural indices from multimodal brain mri: suggestions based on a study of patients with diffuse intrinsic pontine glioma[END_REF]. A task of classification of control tissue (white matter) versus tumour tissue was considered. Although the study reports AUCs above 0.80, it also shows that image preprocessing impacts the results, particularly binning, which is highly dependent on inter-site image harmonisation. [START_REF] Brynolfsson | Haralick texture features from apparent diffusion coefficient (adc) mri images depend on imaging and pre-processing parameters[END_REF] studied the effect of five imaging and pre-processing parameters on the GLCM texture features. These parameters are: the noise level, the resolution, the ADC map construction method, the quantisation method, and the number of grey levels in the quantized images.

The authors found that most features are affected by the investigated parameters using two datasets, glioblastoma and prostate cancer. Only the b-values used for constructing the ADC maps in the glioma data set had no significant effect on any feature. [START_REF] Moradmand | Impact of image preprocessing methods on reproducibility of radiomic features in multimodal magnetic resonance imaging in glioblastoma[END_REF] studied the effects of bias correction and noise reduction on the radiomic features. Overall, the total number of high-robustness features extracted from the necrosis region (30.6%) was higher than the number of features extracted from oedema regions (20.2%), enhancement regions (19.2%), and active tumour regions (17.3%). Additionally, the average number of highly reproducible features for baseline comparison with bias field correction (23.2%) was higher than for baseline comparison with noise correction (21.4%). Interestingly, performing the bias field correction before noise correction increased the stability (22.5%) compared with the noise correction before bias correction (20.4%). [START_REF] Poirot | Robustness of radiomics to variations in segmentation methods in multimodal brain mri[END_REF] studied the effect of segmentation on radiomic features. They used four classical segmentation software to compare radiomic features extracted from several brain compartments. The results show that despite a high segmentation agreement between the different methods (all above 0.75%), the radiomic features present significant differences. This is especially true for texture features extracted from T1w scans. [START_REF] Peerlings | Stability of radiomics features in apparent diffusion coefficient maps from a multi-centre test-retest trial[END_REF] studied the reproducibility of radiomic features. Scans were acquired twice within seven days under similar conditions, and their radiomic features were compared using the Concordance Correlation Coefficient (CCC). The results show that only around 25% of the features were reproducible. Nevertheless, radiomics features achieved comparable reproducibility after waveletfiltering, which can alleviate boundary inhomogeneity. The authors also found that features extracted from 3T images were more stable (32%) compared to 1.5T images (25%). No significant difference in stability was observed when comparing the scanner manufacturers; however, the stable features have only around a 58% overlap. Finally, 122 features were identified as stable across multiple cancer types (ovarian, lung and liver cancer).

The usage of radiomic features extracted from MR scans is still emerging compared to CT and PET scans. Thus, studies of the different stability issues, solutions, and protocols to alleviate these instabilities are still developing. [START_REF] Cattell | Robustness of radiomic features in magnetic resonance imaging: review and a phantom study[END_REF][START_REF] Van Timmeren | Radiomics in medical imaging-"how-to" guide and critical reflection[END_REF] present multiple works on the topic. The presented works suggest that radiomic features can be good predictors and can help in learning about various diseases. However, these features are not stable or reproducible to be used as biomarkers for tumour types.

Brain tumour segmentation

Image segmentation is a key step for radiomic analysis. Identifying the Region of Interest (ROI)

from the MR scan is important to extract meaningful descriptors of the tumour at hand. Given an input volume from one or multiple modalities, automatic brain tumour segmentation refers to the differentiation between the tumour (and its sub-compartments) from healthy brain tissue. This is done by statistical models that map each voxel into a set of sub-compartment labels.

In this section, we discuss the various automatic brain tumour segmentation methods. We focus on the best-performing models reported in the Brain Tumour Segmentation (BraTS) challenge (Menze et al., 2015). The challenge is focused on the segmentation of Glioblastoma and LGGs. The challenge has occurred since 2012 and reflects advances in the automatic brain tumour segmentation problem.

Using custom features

Random Forests have been one of the early successful models for brain tumour segmentation. For instance, [START_REF] Criminisi | Context-sensitive classification forests for segmentation of brain tumor tissues[END_REF] successfully used a Random Forest model to achieve a top rank in the BraTS challenge 2012. Their model segments the image into background, tumour and oedema by independently assigning a label to each voxel. First, given an intensity, they estimate its probability of belonging to each class. Then, these probability maps were combined with the original MR scans and used as input to their random forest. Their method achieved a Dice score of 0.7 on the whole tumour and 0.2 on the tumour core.

Reza and Iftekharuddin (2013) also used a Random Forest as an automatic brain tumour segmentation model. From the available MR scans, they computed the difference in voxel intensities between the different modalities, and they extracted the fractal Piece-wise Triangular Prism Surface Area (PTPSA) [START_REF] Islam | Multifractal texture estimation for detection and segmentation of brain tumors[END_REF] and textons [START_REF] Leung | Representing and recognizing the visual appearance of materials using three-dimensional textons[END_REF]. Then they combined these features with the MR scan intensities and trained their Random Forest. This method was among the highest-ranked models in BraTS 2013, with a reported Dice of 0.83 on the whole tumour and 0.72 on the tumour core. [START_REF] Dvorak | Structured prediction with convolutional neural networks for multimodal brain tumor segmentation[END_REF] proposed to redefine the segmentation output. The authors separated the segmentation problem into three binary segmentation sub-tasks: the whole tumour, the tumour core and the enhancing tumour segmentation. For each sub-task, they divided the label image into binary patches. This is done to extract a set of possible local appearances of each label, which they reduced into a few templates using k-means. Finally, the authors used Convolutional Neural Networks (CNN) to map the original image patch into its closest template, and all overlapping predictions of a neighbourhood are averaged. These predictions are thresholded at 0.5 to obtain the final prediction label. Using this method, the authors obtained a Dice score of 0.83, 0.75 and 0.77 on the whole tumour, the tumour core and the enhancing tumour, respectively.

Using Deep-learning

Since 2015, almost all best-performing models in the BraTS challenge have used Convolutional Neural Networks (CNN). Here, we discuss the multi-path networks and encoder-decoder framework. These two classes of CNNs have been prominent in recent research about automatic brain tumour segmentation.

Multi-path networks

A neural network path is defined as a chain of data processing. Most Multi-paths neural networks divide the input image into small patches and aim at finding the label of the central pixel or voxel at the centre of each patch. Using multiple paths allows the analysis of the patch at different scales.

Thus utilising local details and global information to make a prediction about the target voxel.

Related Works

Encoder-decoder networks

Using multi-path networks to propose a classification of the central pixel ( voxel ) of a patch is computationally challenging. There is a lot of computation redundancy, increasing training and inference times. Furthermore, the accuracy of the prediction can be heavily impacted by the patch and its source quality and size. To tackle these issues, models have been proposed to perform end-to-end segmentations. These methods encode the whole image into a latent representation, which is then decoded into the segmentation mask of the same size as the input.

The most popular encoder-decoder architecture is the UNet [START_REF] Ronneberger | U-Net: Convolutional networks for biomedical image segmentation[END_REF]. UNet is (2020) successfully used V-net for brain tumour segmentation. The authors report a dice score of 0.88 on the whole tumour and 0.8 on the tumour core when tested during the BraTS 2018 challenge.

a

Tackling issues with brain tumour segmentation

Brain tumour segmentation corresponds to a classification task with very unbalanced classes. This means that not all classes are represented equally: the tumour generally occupies 15% of the brain, of which the oedema represents generally 70% of its appearance. Additionally, using promising neural networks requires a large set of training data, which is often hard to satisfy. Multiple approaches have been proposed to tackle these issues. Here we focus on three main solutions: (a) training the models on multiple tasks (multi-task approaches), (b) customising the loss function, and (c) using priors.

Multi-task approaches

Multiple works introduced auxiliary tasks that are chained to enhance the model performance. The underlying idea is that tasks share information. Thus, the easier task helps solve the harder task. In the case of neural networks, this is done by making all the tasks share the earlier extracted features.

Using priors

Some works added prior knowledge to help their model learn helpful features for the brain tumour segmentation task. For instance, since the tumours appear differently on the FLAIR and CE T1w, [START_REF] Longwei | Three pathways u-net for brain tumor segmentation[END_REF] proposed using two models, each using a different modality. The first used the FLAIR to segment the whole tumour, and the second used only the CE T1w to segment the tumour core. When tested on the Brats 2018, the authors reported a Dice score of 0.85 and 0.72 on the whole tumour and core, respectively. where to look for. It must be noted that the bounding boxes have to be defined beforehand.

Multi-block approaches

Let X (1) , X (2) , • • • , X (J) be J blocks, representing measures of different sources on the same n samples.

Multi-block analysis frameworks are well established by numerous works that proposed to adapt classical statistical models for the multi-block setting. These methods often involve adapting the loss functions or the learning strategy. In this section, we discuss different methods that we think of interest.

Unsupervised models

Different approaches have been proposed to analyse multi-block structured data in an unsupervised manner. These models generally aim at a joint dimensionality reduction, thus getting a more comprehensive overview of the studied context. Here, we present the most relevant unsupervised methods from our perspective.

Models built on the NMF

The Non-negative Matrix Factorisation (NMF) has been extended to the multi-block framework. Joint NMF (jNMF) [START_REF] Zhang | Discovery of multi-dimensional modules by integrative analysis of cancer genomic data[END_REF] defines multiple NMF problems with a shared factor matrix across all blocks. This method is equivalent to the classic NMF with the concatenation of the blocks into a single matrix.

Integrative NMF (iNMF) [START_REF] Yang | A non-negative matrix factorization method for detecting modules in heterogeneous omics multi-modal data[END_REF] is another method based on NMF. Similarly to the jNMF, iNMF defines a set of NMF problems. However, instead of only using a single shared

Related Works

factor matrix, the latent factor matrix is decomposed into a shared and a block-specific factor matrix.

jNMF and iNMF can be seen as a simultaneous clustering of observations and variables.

The Multi-Omics Factor Analysis [START_REF] Argelaguet | Multi-omics factor analysis-a framework for unsupervised integration of multi-omics data sets[END_REF]) (MOFA) is another unsupervised method for decomposing the blocks into the products of a shared factor matrix and block-related weight matrices. The MOFA defines a probabilistic Bayesian problem, where prior probabilistic distributions are imposed on all the estimated variables. Instead of the least-square approximation used by the NMF, MOFA uses a variational inference to estimate the parameters. To make the model more interpretable, the authors propose adding sparsity constraints by using Automatic Relevance Determination (ARD) and spike-and-slab priors on the factor and weight matrices. This is done by the multiplication of the priors on the factor and weight matrices by a set of random variables having a Bernoulli distribution; each prior will select the variables of interest (or make them zero).

Models built on the PCA

The Multiple Factor Analysis (MFA) [START_REF] Escofier | Multiple factor analysis (afmult package)[END_REF]) is an extension of the classical factor analysis method to the multi-block framework. The MFA can be described as a Principal Component Analysis (PCA) of the concatenated blocks. Each block is weighted using the largest eigenvalue from its associated correlation matrix. The weighting ensures a single block variance or dimension does not dominate the global PCA first component.

The Joint and Individual Variation Explained (JIVE) [START_REF] Lock | Joint and individual variation explained (jive) for integrated analysis of multiple data types[END_REF] decomposes each block into the sum of two low-rank matrices. The first matrix represents the joint structure matrix, which contains common information between all the blocks. The second matrix represents the individual information of each block. The joint and individual matrices are required to be column-orthogonal.

JIVE can be seen as two successive PCA decompositions: the first decomposition extracts information from the global structure (using the concatenation of all the blocks), and the second PCA is done on the residualized individual blocks.

Models built on the CCA

The Canonical Correlation Analysis (CCA) is a statistical tool to find shared information from two blocks by maximising their correlation in a latent space. Multiple methods have been proposed

to extend the CCA to take into account more than two blocks. These methods are referred to as Generalized CCA (GCCA). Formulations of the GCCA include SUMCOR, MAXVAR, among others [START_REF] Horst | Generalized canonical correlations and their applications to experimental data[END_REF][START_REF] Kettenring | Canonical analysis of several sets of variables[END_REF]. The GCCA has been further expanded to include regularity or sparsity in the model Regularized GCCA (RGCCA) [START_REF] Tenenhaus | Regularized generalized canonical correlation analysis[END_REF] and Sparse GCCA (SGCCA) [START_REF] Tenenhaus | Variable selection for generalized canonical correlation analysis[END_REF]. The RGCCA and SGCCA will be detailed in later chapters.

Another method is the Multiple CO-inertia Analysis (MCOA) [START_REF] Bady | Multiple co-inertia analysis: a tool for assessing synchrony in the temporal variability of aquatic communities[END_REF], which extends the Co-inertia Analysis to a set of blocks. This is done by the optimisation of the covariance between several individual ordinations and a shared reference ordination.

Other methods have been proposed to analyse multi-blocks requiring special block configurations.

For example, the Tensorial independent component analysis (tICA) [START_REF] Teschendorff | Tensorial blind source separation for improved analysis of multi-omic data[END_REF]) is an extension of the PCA for tensorial data. However, this algorithm is designed for datasets with blocks that share the same samples and variables. The Data Fusion algorithm [START_REF] Zitnik | Data fusion by matrix factorization[END_REF] 3.3. Multi-block approaches 31 requires relationship matrices between pairs of blocks. A comparison between several discussed methods has been presented by [START_REF] Cantini | Benchmarking joint multi-omics dimensionality reduction approaches for the study of cancer[END_REF]. The study concludes that the different benchmarked algorithms have different results depending on the studied issue. For instance, for clustering tasks, iNMF performed the best as designed for this task, while JIVE and the RGCCA had the best performance among the set of methods not intrinsically designed for clustering. Furthermore, MCIA, JIVE, MOFA, and RGCCA were the most efficient for detecting latent features associated with survival.

Supervised models

Multiple approaches have been proposed to extend classic regression models to account for separate blocks, each having distinct variables. Two main strategies can be identified: combining multiple models, each learned from a single block, and learning a unique joint model from the different blocks.

Combination of multiple models strategies

Priority-Lasso [START_REF] Klau | Prioritylasso: a simple hierarchical approach to the prediction of clinical outcome using multi-omics data[END_REF]) is a multi-block regression model based on multiple regression models, each trained on a separate block. Priority-Lasso requires user-defined order of the blocks, which defines their priority to the model. The algorithm sequentially fits regression models on each block, using the obtained results from the previous block as an offset. The fitted models can be over-optimistic on the training set, making the offset capture variability not contained by the model.

To tackle this issue, the algorithm uses a cross-validation approach. On an acute myeloid leukaemia dataset, it showed similar results in terms of prediction accuracy compared to a Lasso model. The authors also showed a dependency of the model on the user-defined priority order.

Sequential and Orthogonalized-PLS (SO-PLS) [START_REF] Jørgensen | Combining designed experiments with several blocks of spectroscopic data[END_REF]) is a regression model based on the Partial Least Square (PLS) regression. Given two blocks, the algorithm starts by fitting the output to the first block using a PLS regression. Then, it uses the extracted PLS scores to orthogonalise the second block before using it in a second PLS regression. Finally, the model combines the two computed PLS scores and uses a classical least square regression to predict the output. The authors claim that the impact of the blocks order is minimal regarding the prediction. The method can be easily expanded to multiple blocks.

Parallel Orthogonalised-PLS (PO-PLS) [START_REF] Måge | Regression models with process variables and parallel blocks of raw material measurements[END_REF] is another multi-block regression model based on the PLS. Instead of sequentially regressing the different blocks, the authors propose to extract shared information first. This is done by computing PCA scores from the blocks, then using the Canonical Correlation Analysis (CCA) (or the Generalised CCA when more than three blocks are involved) on the PCA scores to obtain the canonical coefficients. Only the components having a canonical correlation exceeding 0.95 are kept. The original PCA scores are then orthogonalised against the canonical coefficients. Finally, a PLS regression model is built using the common and orthogonalised components.

Structured penalties strategies

Structured penalties have been proposed to integrate multiple blocks into machine learning models.

For example, Integrative LASSO with Penalty Factors (IPF-LASSO) [START_REF] Boulesteix | IPF-LASSO: Integrative L 1 -Penalized Regression with Penalty Factors for Prediction Based on Multi-Omics Data[END_REF] proposed

Related Works

using a weighted sum of the ℓ 1 norms of the coefficient vectors of each modality as a penalty. The proposed penalty adapts the sparsity for each block, which allows adjusting the proportion of variables of interest selected from each modality. However, this comes at the cost of the number of parameters needed to tune, which is proportional to the number of blocks at hand.

Random forests have been expanded to take into account the multi-block structure of the data in an approach called BlockForest [START_REF] Hornung | Block forests: random forests for blocks of clinical and omics covariate data[END_REF]. The classic random forest algorithm takes a small sample of the variables on which it computes a criterion and makes a splitting decision.

The imbalance of the number of variables from each block causes an over-representation of variables from larger blocks (in terms of the number of variables), regardless of their link to the outcome. To tackle this issue, the authors propose five sampling strategies. In VarProb, each variable has a blockspecific variable sampling probability. With SplitWeights, the split criterion is weighted according to the block. Using BlockVarSel, the algorithm samples a fixed number of variables from each block. In RandomBlock, each split is done on a single, randomly selected block. Finally, with BlockForest, a random set of blocks is chosen before BlockVarSel is applied.

Other methods based on classic machine learning strategies, such as gradient boosting, have been proposed to take into account the block structure of data. Neural networks have also been used for supervised multi-block learning, and a study of these methods can be found at [START_REF] Leng | A benchmark study of deep learning-based multi-omics data fusion methods for cancer[END_REF].

A benchmark of some of the discussed methods has already been done by [START_REF] Herrmann | Large-scale benchmark study of survival prediction methods using multi-omics data[END_REF] in the context of survival prediction with multi-omics. When tested on the survival prediction task, the study reports that BlockForest was the best-performing model among all tested models. However, the study also reports a high variability of the model performances when tested on different tumour types.

For instance, Priority-Lasso was the best-performing model when tested on liver, lung and ovarian cancers, but its performance also dropped below the random results on other cancers.

Variable selection

The usefulness of statistical models is not limited to their interference performance but also their ability to exhibit variables of interest. Variable selection is an essential regularisation tool to fit models when the number of studied variables far exceeds the number of samples available, which is the most frequent case when analysing biological data. As a side-effect, this allows better interpretability of a built model and the identification of the most critical variables, which opens further investigation into them. In this section, we will discuss the different proposed methods to perform variable selection.

Sparsity-based penalties

One of the most common approaches for variable selection is introducing a sparsity constraint, which refers to setting the coefficients of non-important variables to zero. Most of these methods are adaptations of the Least Absolute Shrinkage and Selection Operator (LASSO) [START_REF] Tibshirani | Regression Shrinkage and Selection Via the Lasso[END_REF]. Almost all of these methods focus on respecting or exploiting dependencies among the studied variables.

The Fused Lasso is another adaptation of the Lasso penalty proposed by [START_REF] Tibshirani | Sparsity and smoothness via the fused lasso[END_REF].

The penalty was designed for ordered data such as time series. The penalty aims at selecting a continuous sequence of variables. The penalty takes the form

λ 1 ||β|| 1 + λ 2 |β i -β i+1 |.
The authors proposed to estimate the order from the data when there is no a priori order known. For example, they organised similar variables near one another using hierarchical clustering and applied it to gene expression data. [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF] introduced the Group Lasso penalty, which allows a model to jointly select (or discard) all variables within the same group. The groups are predefined and reflect a priori knowledge.

The penalty takes the form λ 1 J j=1 ||β j || K j with K j a symmetric positive definite matrix and the estimates β have been separated into J disjoint groups.

The Octagonal Shrinkage and Clustering Algorithm for Regression (OSCAR) has been proposed by [START_REF] Bondell | Simultaneous regression shrinkage, variable selection, and supervised clustering of predictors with OSCAR[END_REF] and aims at simultaneously selecting variables while grouping them into predictive clusters. OSCAR encourages the grouping of highly correlated covariates. This is done by combining a Lasso penalty and a pair-wise ℓ ∞ . The penalty takes the form

λ 1 ||β|| 1 + λ 2 j<k max(|β j |, |β k |).
While the two previously discussed penalties rely on a priori knowledge, OSCAR infers the groups from the data.

Clustered Lasso, by She (2010), is another adaptation of the Lasso penalty. Similar to OSCAR, the penalty aims to group similar variables of interest into clusters without a priori knowledge. The penalty is built as a combination of the Lasso and a pairwise absolute difference between all the coefficients. The penalty takes the form λ 1 ||β|| 1 + λ 2 j<k (|β j -β j |). While Clustered Lasso and OSCAR seem similar, they differ in their treatment of anti-correlated variables. The authors found that while the predictive performance of Clustered Lasso was unsatisfying, they were able to use the estimated coefficients to order the covariates and run a successful Fused Lasso.

While some of these methods allow for the injection of a priori knowledge, they are not suited for complex interactions between the covariates, as it is the case between genes. Meanwhile, the grouping performed by the other methods can only be inferred from the data and do not necessarily reflect previously established interactions.

Graph based penalties

Biological interaction graphs are the result of many years of biological research. These graphs depict complex interactions between entities such as proteins or genes. Multiple propositions have been made to include these graphs in statistical models. Graph integration aims at smoothing the estimates over the network nodes and thus making the models more interpretable. Here we only focus on methods that add graph knowledge as a penalty to guide variable selection and some of their applications.

Graph penalty

Li and Li (2008) introduced a graphical penalty into their regression problem to smooth the estimated parameters over their graph. Instead of using an ElasticNet, the authors proposed to combine a Lasso penalty with a graph penalty that takes the form λ 1 ||β|| 1 + λ 2 ||β|| 2 L with L the normalised Laplacian associated to their graph. The authors demonstrated the grouping effect. In their study, the authors applied their model to glioblastoma data using genetic data (microarray). They showed that the performances of the graph constraint were similar to the ElasticNet while selecting more variables, which allowed them to exhibit sub-networks in the graph. [START_REF] Grosenick | Interpretable whole-brain prediction analysis with graphnet[END_REF] proposed another variant of the graphical constraint, which they called GraphNet. Instead of using the normalised Laplacian, they proposed using the Laplacian of the graph.

Their aim was to obtain sparse but structured solutions by combining structured graph constraints with a global sparsity-inducing prior that automatically selects important variables. The authors applied their method to functional MRI data and a brain connectivity graph. They found the same grouping properties of the penalty as previously discussed.

GraphNet assumes that connected nodes in the graph should have similar coefficients, which is biased against negatively correlated variables. [START_REF] Du | Structured sparse canonical correlation analysis for brain imaging genetics: an improved GraphNet method[END_REF] proposes to alleviate this issue by using the absolute value of the coefficients. The authors combined their proposed penalty with a CCA model and called it the Absolute value GraphNet Structured CCA (AGN-SCCA). Applied to imaging-genetic data from an Alzheimer dataset (Alzheimer's Disease Neuroimaging Initiative (ADNI) database), this penalty outperformed similar methods in extracting the canonical correlations. Additionally, it yielded sparser results but identified relevant genetic and imaging markers. [START_REF] Du | Detecting genetic associations with brain imaging phenotypes in Alzheimer's disease via a novel structured SCCA approach[END_REF] proposed two new graph-based penalties to capture better structures from the data. The first penalty, which they called Fused pairwise Group Lasso (FGL), takes the form

λ j β 2 j + β 2 j+1 .
The second penalty, Graph guided pairwise Group Lasso (GGL), takes the form λ (j,k)∈V β 2 j + β 2 k with V the set of edges of the graph. Both penalties are not sensitive to the correlation sign. The authors combined the penalties with the CCA model and compared them with other penalties, including the AGN, on the ADNI dataset. They obtained better canonical correlations with FGL.

Applications

GraphNet has been applied using different statistical models in various biostatistical contexts. [START_REF] Watanabe | Disease prediction based on functional connectomes using a scalable and spatially-informed support vector machine[END_REF] compared Fused Lasso, ElasticNet and GraphNet when applied to a Support Vector Machine (SVM) model. They tested the models on functional MRI data from a schizophrenia dataset (Center for Biomedical Research Excellence (COBRE)) and a spatial parcellation graph. Their results show that Fused Lasso and GraphNet both select spatially continuous regions, with the GraphNet selecting more variables. The increase in model complexity did not come at the cost of prediction accuracy. The classification accuracy with ElasticNet, GraphNet and Fused Lasso are comparable with similar sparsity levels. [START_REF] Kim | GraphNet-based imaging biomarker model to explain levodopa-induced dyskinesia in Parkinson's disease[END_REF] studied the GraphNet penalty combined with a regression model and applied it to Parkinson's disease. They compared GraphNet with Group Lasso and Lasso. Their results show that the graph penalised model resulted in a more sparse model than the Group Lasso and a less sparse model than the Lasso. The GraphNet penalty model outperformed the two other models in predicting the severity of the disease. [START_REF] Zhu | Network-based cancer genomic data integration for pattern discovery[END_REF] added the graph penalty to their Sparse Singular Value Decomposition (SSVD) model, which is called Sparse Network SVD (SNVD). The authors applied their method to integrate a prior gene interaction network from a protein-protein interaction network and gene expression data to identify underlying gene functional modules. Real cancer genomic data show that most co-expressed modules are significantly enriched on GO [START_REF] Ashburner | Gene ontology: tool for the unification of biology[END_REF]Consortium, 2020) [START_REF] Castel | Histone h3 wild-type dipg/dmg overexpressing ezhip extend the spectrum diffuse midline gliomas with prc2 inhibition beyond h3-k27m mutation[END_REF], leading to the same epigenetic dysregulation as H3K27M. Furthermore, the oncogenesis of these tumours seems to present general known mechanisms implicating the genes TP53, EGFR, PDGFRA, ACVR1, Src and mTOR pathway. The PHRC BIOMEDE and BIOMEDE 2.0 is an unprecedented effort in France and Europe that aims at collecting multi-omics and imaging data. The objectives are both clinical and fundamental: 1) to test new chemotherapies assigned based on patients' individual molecular profiles and to 2) to build a well characterized multi-modality measurement dataset to allow for a better understanding of the biological underpinnings of this tumour, to define new biomarkers, notably non-invasive ones based on imaging features, and derive some new hypnotical pharmaceutics in order to improve disease control.

During the duration of the PhD preparation, a significant fraction of the data could finally be accessed, but not the whole dataset. Without reconsidering the directions of the methodological work we chose to extend the datasets to tumours that have long been considered as a comparative context

for DIPGs even though it is now known that these tumours are of different origin and prognosis. We have chosen high grade gliomas (glioblastoma) and low grade gliomas Our work involves three different tumour types: Glioblastoma, LGG and DIPG. Glioblastoma and

LGG were also chosen for their publicly available data as well as the multiple previous works about them. This allows us to benchmark our methods and compare our results with earlier findings. Finally we used the DIPG to apply our techniques and report the results. These results are given in terms of prediction or tumour characterisation.

This chapter describes the different datasets used with their various preprocessing steps. We start by defining the data sources and various acquisition parameters. Then we lay out the preprocessing pipeline for the radiomic feature extracting and the tumour segmentations. Finally, we describe the preprocessing applied to omic data.

Glioblastoma

To study glioblastoma, we used the TCGA-GBM dataset [START_REF] Brennan | The somatic genomic landscape of glioblastoma[END_REF][START_REF] Mclendon | Comprehensive genomic characterization defines human glioblastoma genes and core pathways[END_REF]. The dataset comprises 602 patients from different sites, all diagnosed with glioblastoma. We also used the public image dataset BraTS'19 [START_REF] Bakas | Advancing the cancer genome atlas glioma MRI collections with expert segmentation labels and radiomic features[END_REF][START_REF] Bakas | Identifying the best machine learning algorithms for brain tumor segmentation, progression assessment, and overall survival prediction in the brats challenge[END_REF]Menze et al., 2015) (last access January 2020), which comprises 254 patients diagnosed for High Grade Gliomas. For each individual, 4 MRI volumes corresponding to T1w, CE T1w, T2w and FLAIR were available.

These volumes were acquired with different clinical protocols and with various scanners from eight institutions and originated from different studies.

From the HGG set, we isolated the 102 patients diagnosed with Glioblastoma and belonging to the TCGA-GBM sub-cohort (and only kept 97 patients having manual segmentation provided). From now on, the HGG dataset refers to the (254 -102 = 152) patients from BRaTs'19 HGG deprived of TCGA-GBM. The HGG dataset will be divided into HGG train and HGG val and used for the training and validation sets for our segmentation method.

Lower Grade Gliomas

To study Lower-Grade Gliomas (LGG), we worked with the TCGA-LGG dataset (TCGA, 2015) obtained from The Cancer Genome Atlas (TCGA) project (http://cancergenome.nih.gov). We used the data as they have been made available on https://www.openml.org/ by [START_REF] Herrmann | Large-scale benchmark study of survival prediction methods using multi-omics data[END_REF].

Data was collected between 1992 and 2013 and comprises five groups of variables (five blocks), including Clinical records, Copy Number Variants (CNV), micro RNA expression (miRNA), gene expression (RNAseq), and somatic mutations, obtained for 419 patients. Clinical information, including all TCGA-LGG patient demographics, is summarised in Table 4.2, and the estimated survival probability using Kaplan-Meier in Figure 4.2.

Additionally, we used the MR scans of 108 patients obtained from the BRaTs'19 dataset, 76 of which are publicly available. We used the 76 initially obtained images to assess our segmentation methods and all 108 images for radiomic extraction and data integration. All four MRI modalities corresponding to T1w, CE T1w, T2w and FLAIR were available for each individual. Similarly to TCGA-GBM, these volumes were acquired from five different institutions using various clinical protocols and scanners. All the cases have the T2w and FLAIR modalities available. Manual delineation of the whole tumour was achieved for all test cases and checked on the most relevant modality (for some patients, ASL and diffusion MR images were used) [START_REF] Raphael Calmon | Multimodal magnetic resonance imaging of treatment-induced changes to diffuse infiltrating pontine gliomas in children and correlation to patient progression-free survival[END_REF][START_REF] Calmon | Radiogenomics of diffuse intrinsic pontine gliomas (DIPGs): correlation of histological and biological characteristics with multimodal MRI features[END_REF]. These segmentations were obtained on images with a 4 mm 3 resolution. Informed consent of the use of clinical and radiological data was obtained and the protocol was approved by the local ethics committee of the reference institution.

The BIOMEDE project is an initiative led by the French Society of Paediatric Oncology and the European consortium "Innovative Therapies in Children with Cancer (ITCC)". The cohort includes 189 patients diagnosed with DIPG from multiple European institutions. Available data includes MR scans, gene expression and somatic mutation profiles. No delineation of the tumours is available.

The PREBIOMEDE cohort was only used to assess our segmentation method. As of 09/12/2022, no clinical data was accessible for the BIOMEDE patients.

Data availability

The complete data availability description can be found in Table 4.3. We used the TCGA-GBM, TCGA-LGG and PREBIOMEDE to assess the performance of our segmentation model (see chapter 5). Since the mutation and image require prior preprocessing, we only used the CNV, miRNA and RNAseq data to study variable selection in a multi-block framework from the TCGA-LGG dataset (see chapter 6). We used RNAseq, mutations and images for the radio-genomics integration (see chapter 7). Survival data were unavailable for BIOMEDE, and TCGA-GBM did not have sufficient data which precluded their use in the latter study.

Preprocessing

To use the different data at hand, each block underwent preprocessing steps. This section details the preprocessing pipeline applied to each block.

Image preprocessing

Radiomic analysis utilises intensity distributions to extract features. In order to obtain comparable intensity distributions within the studied cohorts, several preprocessing steps have been applied to each image. First, different MRI sequences of the same patient were skull-stripped and co-registered to each other. Then, image intensities were bias field corrected and normalised into a common distribution, allowing comparison across different patients. Finally, multiple radiomic features were automatically calculated using previously extracted regions of interest (ROI).

For the segmentation, we used the bias field corrected images. Then value harmonisation was obtained with a min-max intensity normalisation. Our choice of normalisation is justified by our intent to enhance tumour intensities, thus making our detection and segmentation tasks easier. The images were min-max normalised using the 5% and 95% percentiles to discard outliers. Out-of-range values were capped as in Eq 4.1 where v and v are respectively the original and normalised grey level

Skull-stripping

Skull-stripping aims at finding a 3D brain mask in the MR scan, which removes the skull. We used the FSL BET [START_REF] Smith | Fast robust automated brain extraction[END_REF] software. The algorithm first estimates the brain volume using the intensity distribution. This is done by removing outlier intensities from the image (2% and 98% percentile), then binarising the image using a threshold that lier at 10% of the intensity range. This allows the algorithm to estimate the brain centre of gravity (CoG) and its spherical radius. Starting from a spherical tessellated surface centred on the CoG with half the radius estimated, the surface iteratively expands until it reaches a stationary state.

Co-registeration

The co-registration problem aims at aligning a moving image on a reference volume. It is done by finding the best geometric transformation that maximizes the similarity between the transformed moving image and the reference volume. We used rigid co-registration to align each patient's different scans to their reference modality, which defines a transformation with six degrees of freedom (three translations and three rotations). In our case, the reference modality is the T1w, which is often the highest-quality image in terms of resolution and contrast.

We used the FSL FLIRT [START_REF] Jenkinson | A global optimisation method for robust affine registration of brain images[END_REF][START_REF] Jenkinson | Improved optimization for the robust and accurate linear registration and motion correction of brain images[END_REF] for the co-registration.

The software utilises a hybrid local-global optimisation algorithm on the same image at multiple resolutions. We used the normalised mutual information as our cost function, defined as

C(X, Y ) = H(X,Y ) H(X)+H(Y )
where H is the standard entropy and X and Y are two volumes with the same dimensions.

Bias field correction

Bias field correction aims at removing unwanted low-frequency intensity nonuniformity from the images. We used the N4ITK algorithm [START_REF] Tustison | N4itk: Improved n3 bias correction[END_REF]. As the N3 bias correction algorithm (Sled et al., 1998), the bias field is modelled with a multiplicative effect (i.e. I biased (x) = I unbiased (x)F (x)

with I the image and F the field at a voxel x). Their approach consists of finding a smooth, slowly varying gaussian field that maximises the frequency of the unbiased image. N4 improved the optimisation routine and the B-spline approximation method over its predecessor. The implementation used was provided by Advanced Normalisation Tools (ANTs) Software https://github.com/ANTsX/ANTs.

White-stripe normalisation

Intensity normalisation aims at reducing discrepancies between intensity distributions within tissue classes (white matter, grey matter ...) across the subjects. We used the white stripe normalisation [START_REF] Shinohara | Statistical normalization techniques for magnetic resonance imaging[END_REF]. The method focuses on matching the distributions between the different patients on the white matter. The intensity distribution in other tissues is adjusted accordingly:

I normalised = I -µ W M σ W M (4.2)
With I the 3D image. To estimate the mean (µ W M ) and standard deviation (σ W M ), the distribution of intensities within a 4 cm rectangle at the centre of the head is considered. This rectangle contains a large amount of white-matter voxels. The distribution is smoothed using B-splines of order four. The µ W M corresponds to the last peak in the smoothed distribution for the T1w modality and the highest peak for the T2w and FLAIR modalities. These intensities correspond to estimated whitematter voxels. To estimate σ W M , the set of voxels having intensities within the range f -1 (f (µ W M )-τ ) and f -1 (f (µ W M ) + τ ) are considered. With f the smoothed estimated intensity distribution and τ is a quantile tolerance, which takes the value 0.05 as recommended by the original paper.

We used the hybrid white-stripe to normalise the image, which takes the intersection of estimated white-matter voxels from multiple modalities of the same patient. We used in-house implementation of white-stripe, using python language (Van Rossum and Drake, 2009).

Radiomic extraction

We extracted features from the original images and eight wavelet-transformed images. This resulted in 846 ( (19 + 24 + 16 + 16 + 14 + 5) × 9 (original + 8 wavelet transformed images)) radiomic features extracted from each modality. Note that we did not use shape-related features. We used the pyradiomics software [START_REF] Joost | Computational Radiomics System to Decode the Radiographic Phenotype[END_REF] to extract the radiomic features.

Mutation data preprocessing

Statistical analysis of somatic mutations is challenging. Most tumour mutations are only present in a small sample of patients. The mutated gene may or may not be impactful on the pathology. Statistically, deciding the relevance of these genes from the small observed samples is infeasible. Additionally, mutations in genes may lead to different results. Some mutations lead to the dysfunction of the gene, which promotes the cancerous behaviour of the cells. Other mutations on some genes may impact other genes in the same biological pathway, even if the latter are not mutated. Finally, other mutations may not cause tumoural cells but rather be a consequence of the rapid cell division typical of cancer.

Starting from the matrix of somatic mutations (samples × genes), we used network smoothing proposed by [START_REF] Vanunu | Associating genes and protein complexes with disease via network propagation[END_REF] which diffuses the mutation information in the matrix according to the network. This procedure normalise the somatic mutation profiles under the hypothesis that mutation or dysfunction in a gene leads to the disturbance of the behaviour of other genes it interacts with. Using a gene-gene interaction graph, the mutation propagation is given by the equation:

X t+1 = αX t D -1 2 AD -1 2 + (1 -α)X 0 (4.3)
Where X t is the somatic mutation profile at iteration t, with X 0 the initial somatic profiles.

D is the diagonal matrix of the gene degrees (number of neighbours in the graph ). A is the binary adjacency matrix indicating if two genes are neighbours in the graph. α is a tuning parameter between 0 and 1, governing the distance that a mutation signal is allowed to diffuse through the network during propagation. We computed the propagation t times until convergence was found for the Frobenius norm.

Gene-gene interaction graphs usually describe complex interactions between genes. Suppose the case where a gene implicated in multiple biological pathways is mutated. We hypothesise that the influence on each of the biological pathways (and on each gene involved in these pathways) is lesser compared to a gene implicated in fewer pathways. The distance between two genes in the graph does not reflect their influence on each other. To solve this issue, we used the influence graph computed from the gene-gene interaction graph proposed by Vandin et al. (2011b). The method uses the heatdiffusion model on the gene-gene interaction graph to compute the influence. Let L γ = L + γI, where L is the Laplacian of the interaction graph, and γ represents the diffusion rate. The influence scores are given by L ′ = L -1 γ and L ′ i,j gives the influence of gene i on the gene j. The influence graph is defined as G = (E, V ), where E is the set of genes, and for each gene, its neighbours are the n genes with the highest influence score on it (by default n = 10). Note that the graph is undirected.

In our proposed thesis, we decided to preprocess the somatic mutation profiles using both the original Pathway Commons gene-gene interaction and the influence graph derived from it. We used γ = 10 -2 to compute the influence graph. For the propagation algorithm, we set α = 0.5, as recommended for Pathway Commons by [START_REF] Hofree | Network-based stratification of tumor mutations[END_REF].

✽ ✽ ✽ ✽ ✽ ✽ Part II

Contributions

Chapter 5

Brain Tumour Segmentation

Tumour lesion segmentation is a key step to study and characterise cancer from MR neuroradiological images. Nowadays, numerous deep learning segmentation architectures have been shown to perform well on the specific tumour type they are trained on (e.g. glioblastoma in brain hemispheres). But, a high performing network heavily trained on a given tumour type may perform poorly on a rare tumour type for which no labelled cases allows training or transfer learning. Yet, because some visual similarities exist nevertheless between common and rare tumours, in the lesion and around it, one may split the problem into two steps: object detection and segmentation. For each step, trained networks on common lesions could be used on rare ones following a domain adaptation scheme without extra fine-tuning.

This study evaluates the impact of adding an object detection framework into brain tumour segmentation models, especially when the models are applied to different domains. We identify object detection as a simpler problem that can be injected into a segmentation model as an a priori, and which can increase the performance of the segmentation models. All models were trained on glioblastoma, and they were evaluated on their performances on the glioblastoma, LGG and DIPG. This work has been the subject of a journal article (Chegraoui et al., 2021a).

Availability: Source code is freely available at https://github.com/neurospin-projects/ 2021_hchegraoui_DetectSegmBIOMEDE/tree/main

Introduction

Learning from images classically requires large cohorts for which tumours are finely delineated. The rarity of DIPG added to the fact that segmentation is not part of the clinical routine procedure makes it difficult to obtain robust statistical classifiers or predictors. Automatic tumour segmentation based on transfer learning could theoretically alleviate this problem, circumventing the small number of manual delineations which prevent directly and efficiently training these segmentation models.

In general, DIPG tumours have a central location involving more than 50% of the pons [START_REF] Warren | Diffuse intrinsic pontine glioma: poised for progress[END_REF]. On MRI scans and due to its infiltrating nature, the tumour appears as an intrinsic expansion of the brainstem and not as a distinct foreign mass compressing the pons. However, the tumour is not always restricted to the pons and it can infiltrate other compartments of the central nervous system such as the cerebral peduncles and supratentorial midline or the cerebellum [START_REF] Hankinson | Interpretation of magnetic resonance images in diffuse intrinsic pontine glioma: A survey of pediatric neurosurgeons -Clinical article[END_REF].

The deformation of the pons induced by the tumour and its infiltrating nature makes its detection and delineation non-trivial. On the T2w scans, the tumour presents a hyper-intense signal while it appears hypo-intense with indistinct tumour margins on T1w scans. Enhancement following gadolinium injection (T1Gd) is inconstant and often absent. Finally, the tumour is relatively homogeneous on FLAIR modality [START_REF] Warren | Diffuse intrinsic pontine glioma: poised for progress[END_REF]. Figure 5.1 exhibits a DIPG tumour on different modalities. Because rare tumours present some visual similarities with common tumours, in the lesion and around it, one may split the tumour delineation problem into two steps: object detection and segmentation. For each step, trained networks on common lesions could be used on rare ones following a domain adaptation scheme without extra fine-tuning. Our work suggests different architectures to solve the segmentation and detection tasks and two combination strategies of the two tasks. We assessed the performances of our strategies in three different configurations: i) using the same tumour as in training, ii) using a different tumour located in the supratentorial region as in training and finally iii) using a different rare tumour located in the brainstem, an unseen region during the training.

Detection-segmentation combination strategy

Combination strategies

In order to obtain a robust segmentation of brain tumours, we combined proven object-detection models and segmentation models. Considering we could not directly learn from a few labelled DIPG examples, we decided to train our models with HGG train examples. The HGG and DIPG tumours present both similarities and differences (see Figure 5.1). Tumour intensities have comparable characteristics, while the ages of the patients, tumour locations and image qualities differ. We hypothesised that in a restrained zone around the tumour, HGG and DIPG present enough visual similarities to allow the training of a segmentation model from the sole HGG data and which would be able to segment both types of tumours reliably. We used an object detection model to define these restrained zones around the tumour and bypass the dissimilarities between the two cohorts.

We chose You Only Look Once (YOLO) (Jocher et al., 2020;[START_REF] Redmon | You only look once: Unified, real-time object detection[END_REF] as our objectdetection framework. For the segmentation, we benchmarked UNet [START_REF] Ronneberger | U-Net: Convolutional networks for biomedical image segmentation[END_REF] and

Bounding-Box UNet (BB-UNet) [START_REF] Jurdi | BB-UNet: U-net with bounding box prior[END_REF]. It must be noted that both UNet and BB-UNet receive the whole 2D images as input for the training and inference; additionally, BB-UNet receives also an a priori bounding-box used internally to (non-exclusively) focus the learning segmentation process. Consequently, we examined two different procedures for combining the object-detection and segmentation. We called our first procedure Parallel YOLO UNet (pYU). In pYU, both YOLO and UNet are trained independently. In the inference phase, YOLO-generated bounding-boxes are merely used to mask UNet predictions, thus eliminating all segmented voxels outside the bounding-box. In our second approach, called Sequential YOLO BB-UNet (sYBBU), YOLO and BB-UNet are trained independently, but during inference, YOLO-generated bounding-boxes are provided as additional input to BB-UNet. Similarly to pYU, we used the bounding-boxes to mask the segmentation output. Figure 5.2 illustrates the two approaches, pYU and sYBBU, sketched in their inference stage.

Final masking

Let SEG, BBOX and GT be the sets of voxels belonging to the predicted segmentation, predicted bounding-box and ground truth respectively. Our approaches introduce a masking phase. Masking (versus no masking) will affect precision/recall scores in an anticipated direction, if we make an assumption that will be checked in our results: let v be a voxel and P a probability:

P (v ∈ GT |v ∈ (SEG ∩ BBOX)) ≥ P (v ∈ GT |v ∈ (SEG \ BBOX)) (5.1)

Detection-segmentation combination strategy 53

Under this assumption, it is more likely to find a true positive inside the bounding-box than outside.

Thus, it follows that:

P recision(SEG ∩ BBOX) ≥ P recision(SEG)

Recall(SEG ∩ BBOX) ≤ Recall(SEG)

(5.2)

Masking will always result in a decrease in the recall (amount of the tumour detected). However, we assume that the decrease in the recall will be outweighed by the increase in precision.

(

Input data

We decided to train all our models on 2D slices as we can extract a greater number of different training examples from each 3D volume. Furthermore, we trained all our models only on axial slices since these slices have the best resolution in all the studied cohorts.

In this work, we only used the T2w and the FLAIR modalities. Our choice is justified by several considerations. First, HGG and DIPG MRI scans present similar local tumour patterns mostly in these two sequences. Furthermore, we are only interested in a binary segmentation, and the modalities which best reveals all the different compartments of the tumour are the T2w and the FLAIR. Finally, the DIPG is representative of DIPG data obtained in clinical context, which contains many subjects with missing sequences, but most subjects have either the T2w or the FLAIR available.

Ensembling the inferences

Each model is trained twice : i) on the T2w and ii) on the FLAIR modalities. This makes our approach resilient to missing data. Depending on the data at hand, we retain the inference obtained from the single modality available or we combine the two inferences. In this latter case, we merged the predictions using a weighted average. As described by [START_REF] Shahhosseini | Optimizing ensemble weights and hyperparameters of machine learning models for regression problems[END_REF], we propose to find the optimal weights ω * according to the optimisation problem, where N is the number of validation

cases, ŷF LAIR v (resp. ŷT 2w v
) is the model's confidence scores for the FLAIR (resp. T2w) on a voxel, label(ŷ v ) is the thresholded confidence score and takes the values 0 (for non-tumoural voxels) or 1 (for tumoural voxels) :

with ŷv (ω) = ω × ŷF LAIR v + (1 -ω) × ŷT 2w v ω * = argmin ω - 1 N N k=1 v∈DISCR GT v × log( ŷv (ω)) + (1 -GT v ) × log(1 -ŷv (ω)) s.t. 0 ≤ ω ≤ 1 (5.3)
To find the optimal weight, we used the HGG val set. To better fit our needs, we adapted the optimisation problem by only considering the DISCR = {v|label(ŷ

F LAIR v ) ̸ = label(ŷ T 2w v
)} (i.e. the set of voxels where T2w output label and FLAIR output label are different). This is justified by the fact that the weighted average does not change the predicted label if both models predict the same label.

Ensembling was used both after the detection inferences and after the segmentation inferences.

More precisely, we used the confidence scores of the YOLO models to compute an ensembled boundingbox prediction, and the UNet and BB-UNet scores for an ensemble segmentation prediction. For the segmentation, the ensembling phase is prior to the bounding-box masking.

We used Brent's method [START_REF] Brent | Algorithms for minimization without derivatives[END_REF] to solve the optimisation problem (Eq. 5.3) and SciPy's optimisation package [START_REF] Virtanen | SciPy 1.0: fundamental algorithms for scientific computing in python[END_REF].

Reusing off-the-shelf networks

You Only Look Once (YOLO)

You Only Look Once (YOLO) framework is a multi-scale object detection neural network. YOLO is designed to detect multiple objects of different classes on natural images. The input image is divided into an S ×S grid of cells. The cell where the centre of the object falls into is responsible for predicting the bounding-box and class of the objects. Each cell predicts B bounding-boxes, a confidence related to the existence of an object in each bounding-box, and conditional class probabilities related to the object instances.

Reusing off-the-shelf networks 55

We used the YOLOv5, implemented by Ultralytics [START_REF] Jocher | ultralytics/yolov5: v4.0 -nn.silu() activations, weights & biases logging, pytorch hub integration[END_REF]. This network was pretrained on 416 × 416 images from the Common Objects In Context (COCO) dataset [START_REF] Lin | Microsoft coco: Common objects in context[END_REF].

We used transfer-learning to fine-tune the model parameters for the tumour detection task. We resized our 250 × 250 input images to the 416 × 416 dimension using zero padding. The model was trained to detect the smallest bounding-box around the tumour. Since YOLO is pre-trained on RGB images, we transformed our grey-scale images into RGB images by copying our input image into the three channels.

For the hyper-parameter B, we kept its default value 3. Because YOLO makes predictions with a multi-scale approach, S took successively 3 values in which the prediction is made. For 416 × 416 input images, the values taken were 13, 26, 52. Furthermore, we were interested in detecting only one object, namely the tumour (using the -config file requested by the software).

Finally, starting from the pretrained model, we fine-tuned it for 100 epochs, using an initial learning rate of 0.001. Other default training parameters were kept. Main parameters are listed in the Table 5.1 We assume that tumours are 3D-connected-component volumes. However, as the detection model took axial 2D images as input, there was thus no guarantee to obtain a connected component object in a plane perpendicular to the axial plane. The model might miss the tumour on some slices of the volume, or detect tumours on isolated slices. To overcome this issue, we used a morphological closing of the bounding-boxes, followed by an opening, along the perpendicular axis to the axial plane, with a kernel size of (1,1,6) voxels.

UNet and BB-UNet models

For the segmentation step of the tumours, we studied two models. Firstly, we used UNet, a fully convolutional neural network, which is classically used for biomedical image segmentation [START_REF] Ronneberger | U-Net: Convolutional networks for biomedical image segmentation[END_REF]. Similar to an auto-encoder, it has two paths: an encoding path consisting of the stacking of convolutions, non-linear activations and max-pooling; and a decoding path which consists of convolutions, non-linear activations and transposed convolutions. Skip connections are used between each encoding layer and its symmetric decoding layer. Our network is 5 levels deep on each path. We used the rectified linear unit (ReLU), defined as ReLu(u) = max(u, 0), as the non-linear activation function.

Secondly, we used the BB-UNet (Rosana et al., 2020) model, whose architecture is similar to that of the UNet, except that it takes a binary bounding-box mask as additional input. We added to the UNet a bounding-box path parallel to the encoder path. The binary mask follows similar transformations to the main image. At each skip connection, we carry out an element-wise multiplication between the encoded image and encoded bounding-box. The role of these bounding-boxes is to discourage the network to look beyond it. Our BB-UNet models were trained using the ground truth bounding-boxes obtained as the smallest bounding-boxes comprising all the tumour mask (in 2D). As stated before, we used YOLO-predicted bounding-box during the inference phase.

For both models, the last layer has a soft-max activation, and we used a binary cross-entropy as a loss function, defined as :

Loss(ŷ, y) = -y • log ŷ + (1 -y) • log (1 -ŷ) (5.4)
where y an ŷ are the ground-truth labels and the network confidence score matrices respectively.

We implemented both UNet and BB-UNet using Pytorch (Paszke et al., 2019). The neural networks were trained on 250 × 250 grey-scale images for 100 epochs, with an initial learning rate of 0.001 and the Adam optimiser [START_REF] Diederik | Adam: A method for stochastic optimization[END_REF]. Given the size of the used dataset (≥ 9000), no data augmentation was used.

Deepmedic

We compared all our results to a reference, patch-based brain lesion segmentation network, namely Deepmedic [START_REF] Kamnitsas | Efficient multi-scale 3d CNN with fully connected CRF for accurate brain lesion segmentation[END_REF]. Deepmedic is an 11-layers deep, double-pathway, multi-scale, 3D

CNN. Deepmedic achieved state-of-the-art results on brain tumour segmentation on BraTS'15, and it is continuously updated. We trained Deepmedic on mono-channel twice, with T2w and FLAIR, to make the results comparable with our approaches. Contrary to our models, input images for Deepmedic were normalised using z-scores to remain in line with the network procedure.

We used the implementation from https://github.com/deepmedic/deepmedic. We kept the default values for the hyperparameters as proposed by the original paper, including the number of layers, filters, learning rate, and optimiser.

Experimental designs 57

Experimental designs

To conduct our experiments, we divided the HGG dataset into 90% training set and 10% validation set. We tested all the models on the TCGA-GBM dataset, the LGG dataset and 30 patients (71 sessions) of DIPG. Table 5.2 sums up the dataset sizes. To assess the performance of our approaches on the different test datasets, we used the provided segmentation labels to compute precision and recall, alongside the Dice index. These metrics were measured after 3D reconstruction of the binary masks. We note that the object-detection outputs are also binary masks.

On the TCGA-GBM dataset, we performed a correlation analysis between the ensembled boundingbox performance and the ensembled segmentation performance, in order to establish the impact of the object-detection step on the final segmentation. Furthermore, since BB-UNet models were trained with the ground-truth bounding-boxes while inference was performed using YOLO predicted boundingboxes, we analysed the impact of the used bounding-boxes on the prediction performance of the networks.

On the DIPG dataset, we compared the object-detection performance with a generic bounding-box around the pons. This bounding-box was manually extracted from a template [START_REF] Fonov | Unbiased average age-appropriate atlases for pediatric studies[END_REF] with an enlargement of approximately 50% on each side. Figure 5.3 summarises the experimental design chosen to evaluate the methods.

Benchmark results

Object-detection results

Table 5.3 and 5.4 give the results of the detection phase on the TCGA-GBM and LGG datasets.

Overall, both the FLAIR and the T2w obtain a very high recall and a relatively low precision score.

The merging of both modalities helps further improve the recall and the stability of the predictions (lowering the standard deviations) while lowering the precision. Low precision scores were expected in this phase since the predictions are piece-wise squares while tumour shapes are complex meshes.

Therefore, the precision score depends heavily on the tumour shape and orientation. One must also note that a tumour generally occupies around 7% of the brain, in the studied dataset, which impacts the precision score. To choose the best model, it is important to remember that the main objective of this phase is to generate priors for a segmentation. It is therefore imperative to reliably detect the whole tumour (implying high recall), even if it comes with lower precision. (±.123) .940(±.094)

The detection framework achieved better performance in the TCGA-GBM dataset than in LGG.

This was expected since our model was solely trained to detect high-grade-gliomas. Even if the performance was degraded for the LGG dataset, this decrease is moderate, especially when comparing the results of the ensembled model. This shows that the object-detection model is able to detect different types of tumours that occur in the same tissues of the brain.

Morphological opening and closing showed a minimal effect. However, as these effects were always positive on both precision and recall, we kept them in our detection process.

Segmentation results

YOLO bounding-boxes, obtained in the previous phase, were used during the segmentation. Each segmentation model uses bounding-boxes obtained from the same input image and modality. The bounding-boxes used for the segmentation are all post-processed by the morphological transformations.

Segmentation Results on TCGA-GBM Table 5.5 describes the results obtained for the segmentation of TCGA-GBM. A voxel is considered tumoural if its confidence score ŷv is above 0.5. As expected, precision scores were considerably higher than during the detection phase, however, this came with a decrease in recall. The UNet models performed poorly on TCGA-GBM compared to the other models. Indeed the mean Dice index ranged from 0.70 to 0.78, which is below the other models' averages, with values greater than 0.84. This was mainly due to their low precision scores. A deeper look into the results

showed that UNet segments healthy bright spots of the brain the same way it segments the bright spots indicating the presence of the tumour, especially oedema. On the other side, UNet gave comparable results in the recall metric (i.e. percentage of the tumour detected).

Moreover, we can see a clear improvement after UNet segmentations were masked with the predicted bounding-boxes, i.e the pYU model. Mean precision scores were increased by nearly 15% in all configurations, while the standard deviations were reduced by nearly half. These results suggest that most of the false positives are outside of the bounding-boxes, which is in line with the assumption stated by the equation 5.1. These improvements came with a slight decrease in the recall, of around 1%. We consider this decrease is minor compared to the benefits of masking with bounding-boxes in precision.

Furthermore, bounding-boxes also have a positive effect when they are used as inputs in the BB-UNet models. There is an increase in mean precision and a decrease in standard deviations of precision and Dice index. After masking, the results of UNet and BB-UNet are very similar, with BB-UNet coming slightly ahead, with an improvement in precision between 1% and 3%, and a recall that remained similar among all the models. Figure 5.6 exhibits clearly that models using the boundingboxes perform better, especially sYBBU models. To compare the sYBBU model as the best approach using bounding-box with UNet that ignores them, we computed the AUC of the mean precision-recall graph. We obtained 0.91, 0.90 and 0.93 on the FLAIR, T2w and ens.(FLAIR, T2w) respectively, on the sYBBU model. Meanwhile, on the UNet model, we obtained 0.80, 0.83 and 0.90 using the FLAIR, T2w and ens.(FLAIR, T2w) respectively. The precisions and recalls in segmentation using the bounding-boxes are strongly correlated with the respective precisions and recalls of the bounding-boxes detection, by a correlation score between 0.6 and 0.8 (see Tables 5.6 and 5.7). Of note, UNet results are also positively correlated with the object detection results (correlation between the precisions is 0.40 and correlation between the recalls is 0.66), though not as strongly as in the other models. This suggests that part of the performance is related to the images themselves, and some tumours are especially hard or easy to detect or segment for any model due to image quality or tumour visual characteristics. However, the overall performance is strongly dependent on YOLO's ability to detect the whole tumour. This is shown on the Dice metric, which indicates a strong correlation between bounding-box recall and the Dice of models using the bounding-boxes, ranging from 0.61 to 0.72. This reinforces the strategy consisting in promoting recall over precision during the detection phase in order to obtain overall high performance.

Overall, FLAIR-based models perform better than T2w-based models. It appears that the FLAIR may reflect the diffuse characteristics of the tumour better, while in the T2w images, the intensity distribution of voxels inside the tumour is not as distinguishable from other bright regions of the brain.

However, the ensembled models always perform better, across all configurations, and have equal or lower standard deviations. When computing the optimal weights to merge the models, we found ω = 0.77 for pYU and ω = 0.50 for sYBBU. However, the gain of an optimized weighted average, as opposed to a basic average, was below 1% (results not shown). The weighted average improves the log-likelihood, but with little impact on the accuracy after binarisation.

Table 5.8 shows the differences in metrics when ground-truth bounding-boxes were used for the FLAIR in BB-UNet. BB-UNet with ground-truth bounding-boxes was unable to detect, on average, 10% of the tumour. When YOLO bounding-boxes were used, a 6% decrease in recall was observed.

This exhibits that two-thirds of the missed voxels are inherently related to BB-UNet and not to errors in YOLO bounding-boxes. Given these results, we can say that YOLO bounding-boxes are not the prevailing source of errors and they are sufficient to be integrated into our detection-segmentation approach. To focus on the most interesting region, we only plotted the precision-recall scores for thresholds between 0.1 and 0.9. From the left to the right, using the FLAIR, using the T2w, and ens.(FLAIR, T2w)

Segmentation Results on LGG Table 5.9 shows the segmentation results obtained on the LGG dataset. Overall, the proposed models exhibit comparable results with those obtained on the TCGA-GBM dataset but show an average drop in the Dice metric of 0.05. This reduction was expected since the networks were solely trained on High Grade Gliomas and were not readapted for the Low Grade

Benchmark results

63

Glioma cases. Unlike the proposed models using bounding-boxes, UNet showed poor performances on the LGG dataset. However, the pYU model shows an improvement in the overall results by increasing precision by 30% at the cost of a mean decrease of 10% of the recall. We obtained an AUC score of 0.80, 0.87 and 0.89 for the FLAIR, T2w and ens.(FLAIR, T2w) respectively when using the sYBBU.

Comparatively, we obtained an AUC score of 0.71, 0.73, 0.80 for the FLAIR, T2w and ens.(FLAIR, T2w) respectively when using UNet only. The precision-recall curves are given in Figure 5.7. The proposed procedures outperformed the Deepmedic network. On average, the Dice metric was between 6% and 9% lower for Deepmedic compared to our models. This exhibits the robustness of our strategy. Similar to the TCGA-GBM dataset, the Deepmedic model seems to prioritise high precision (the highest of all the models) over recall. From the 30 DIPG patients, 71 sessions were available obtained at different follow-up visits. Table 5.11 shows the inference detection results obtained on 62 out of the 71 DIPG sessions. The detection step failed to identify the tumour region (tiny bounding-box with recall < 0.015) in 9 sessions (13% of the sample). An example of a failed detection is presented in Figure 5.8. Table 5.10 shows detection results obtained on all 71 test sessions. Overall, the FLAIR exhibited significantly better results than the T2w, especially for the recall. Bounding-boxes obtained from the FLAIR show robust results. On average, 66% of the tumour is captured, with a mean precision equal to 61%. While the T2w alone failed to give significant results, ensembling T2w bounding-boxes with FLAIR boundingboxes improves the recall by 7% while lowering its precision by 17% on average. Both FLAIR and ens.(FLAIR, T2w) bounding-boxes performed better than the generic bounding-box around the pons, which gives a 63% recall with a 20% precision (results not shown). This shows that, even if the location of the tumour is known beforehand, the problem remains non-trivial because of the infiltrating nature of the tumour and its tendency to deform the surrounding tissue or structures (cerebellum, spinal .192 .596 ± .203 .599 ± .194 ens.(FLAIR,T2w) .688 ± .192 .622 ± .203 .618 ± .197 Table 5.13 shows pYU and sYBBU segmentation results using FLAIR detection masks. Table 5.12 shows segmentation results with the combined ens.(FLAIR,T2w) detection masks. Overall, across the configuration reported in Table 5.13, the mean Dice index for segmentation results is 61% (with 95% CI 0.56 to 0.66), which is satisfying considering the difficulty of the problem. An example of the segmentations obtained is presented in Figure 5.9. .733 ± .194 .598 ± .244 .622 ± .206 T2w .423 ± .193 .701 ± .231 .496 ± .190 .697 ± .195 .576 ± .233 .596 ± .200 ens.(FLAIR,T2w) .437 ± .201 .719 ± .235 .509 ± .192 .711 ± .195 .603 ± .243 .616 Since the detection and segmentation phases can be done independently, we computed the performance of segmentation on the T2w, using FLAIR and ens.(FLAIR, T2w) detection masks. Indeed T2w detection, and segmentation using T2w did not fail. However, its results were still below FLAIR ones. On the T2w, pYU performance exhibits a dependence on the detection performance. Indeed, the pYU segmentation model failed to discriminate between tumoural voxels and healthy tissue ones, thus the segmentation results follow the detection performance. This is not the case for sYBBU, using ens.(FLAIR, T2w) detection masks, which have lower precision scores, and did not impact the segmentation model as much as the pYU. Looking at the Dice measurements, FLAIR and ens.(FLAIR, T2w) have similar performances whichever the bounding-boxes and the model used. However, FLAIR tends to have a higher precision while ens.(FLAIR, T2w) has a better recall. The choice between the two approaches should be made in regard to the application.

Segmentation results on DIPG

Discussion

65

On the FLAIR, the Deepmedic network was outperformed by the detection model, and therefore obviously outperformed by the segmentation models that use the FLAIR mask. Deepmedic also failed to detect any tumour region in the same 9 cases excluded earlier.

Discussion

Our study proposes two detection-segmentation combination strategies that allowed us to obtain better results than the tested state-of-the-art networks (UNet and Deepmedic) on both BraTS'19, an openly available HGG and LGG dataset, and DIPG, a cohort of a rare paediatric tumour. Our strategies were able to segment the DIPG lesion while only training the models on the HGG cohort and without re-adapting the networks to the new tumour type. It was necessary to use this domain adaptation since we did not have access to enough annotated DIPG data nor a complete dataset to fine-tune each of the networks used.

Throughout this work, the FLAIR modality consistently appeared as the most important modality for any segmentation model, aiming at delineating globally the tumour lesion without distinguishing between its multiple compartments. It is therefore not surprising that our detection-segmentation algorithms prefer to rely on the FLAIR sequence. Moreover, the FLAIR modality has also been found as the most relevant for oncologists and features extracted from the FLAIR scans have shown the best results for survival analysis and tumour characterization for a range of tumours [START_REF] Kickingereder | Radiomic profiling of glioblastoma: Identifying an imaging predictor of patient survival with improved performance over established clinical and radiologic risk models[END_REF]. Specifically in DIPG, Castel et al. [START_REF] Castel | Histone h3f3a and HIST1h3b k27m mutations define two subgroups of diffuse intrinsic pontine gliomas with different prognosis and phenotypes[END_REF] identified differences in FLAIR index according to the type of histone mutated. Our segmentation, which is based on FLAIR imaging and produces a FLAIR-mostly derived delineation, produces regions of interest that appear to be relevant.

In addition to that it appears that, even if the T2w did not perform as well as expected for the DIPG dataset, its presence always helped the proposed segmentation models.

Our proposals consist of procedures implicating multiple different and distinct models. Having different models, trained separately, has several advantages. The models had very different architectures, and therefore, could have different weaknesses and strengths, which can be complementary. In the DIPG case, even when the T2w detection failed, we were able to use the trained T2w BB-UNet model efficiently using the alternate FLAIR bounding-boxes. This possibility allowed us to circumvent the differences between glioblastoma and DIPG.

Our proposed approaches consist in combining multiple models, each model is relatively small.

The inference time for each 3D example is around 5 seconds on an Nvidia Titan X, including 2.5 seconds for detection and 2.5 seconds for segmentation. Meanwhile, the training phase took roughly 4 hours each. Comparatively, training Deepmedic took 3 days and had an inference time of 3 minutes per 3D example, on similar machine and software configurations.

Most recent segmentation efforts have focused on developing deeper and more complex models.

While these solutions can be suitable for tumour lesions for which large curated and well-documented datasets, there is no indication that they can be easily adapted to small cohorts of rare tumours, such as DIPG, with missing data and heterogeneous quality. We found that Deepmedic, trained with four 5.7. Conclusion 67 modalities (FLAIR, T2w, T1w, T1wgd), performed exceptionally well for HGG, with an average Dice of 0.9, but fails on DIPG with an average Dice of 0.3. In our proposition, the segmentation model is not fully dependent on the object detection performance, given bounding-boxes can be obtained from other input images. This allows us to use the best bounding-boxes assessed during a quality check. Our results are in line with the work of [START_REF] Isensee | nnu-net: a self-configuring method for deep learning-based biomedical image segmentation[END_REF], which found that recent efficient very complex and deep networks cannot necessarily be easily fine-tuned for rare oncological lesions segmentation problems with few training examples and, promoted the UNet architecture.

Our study presents several limitations. The ground-truth segmentations obtained on the DIPG are done on thick slices of 4mm 3 , which negatively bias the results obtained even if it does not question the magnitude of these results. Furthermore, we have only one set of rare tumour data, further studies should investigate the robustness of the method using other rare tumours. Additionally, throughout this study, the used networks are considered black boxes. An ablation study could be made to investigate the limitation of the networks. We also did not investigate what the networks learnt and how they make the inference. Finally, we only focused on binary segmentation using either the FLAIR or T2w, further studies should investigate multi-compartment segmentation possibly using other modalities.

Conclusion

This paper addresses the problem of rare tumour types, for which no database can be built to train a deep neural segmentation network. Our work shows that state-of-the-art segmentation methods perform poorly when applied on test cohorts on which they were not trained. We propose to combine different simple models of detection and segmentation to allow us, not only to improve UNet performance but also to obtain satisfying results on a cohort that contained differences compared to the training dataset regarding, among others, patient age, image quality and tumour type.

Although all the sets presented in the paper concern cerebral tumours, the differences between an adult brain (in the case of the HGGs and LGGs) and children's brains (in the case of DIPG) give rise to challenges during inference. We think that using a set of a wider range of brain tumour types in children might help solve this issue. Additionally, the paper does not explore alternatives to the object detection framework YOLO. Work should be done to compare it to other algorithms, especially the ones dedicated to medical imaging and not only natural images. Lastly, other detection-segmentation strategies, such as the weak supervision paradigm, can be explored and compared to the proposed approaches.

Overall, the bounding boxes give useful a priori knowledge for segmentation purposes. We have already established that these bounding boxes, when correctly detecting the tumour, can be used for radiomic analysis with similar performances to radiomics extracted using fine delineation of the tumour Chegraoui et al. (2021b). In this work, we were able to obtain satisfying segmentation for the DIPG. These segmentations and performance will allow us to perform further clinical work to characterise this rare pathology using radiomics.

✽ ✽ ✽ ✽ ✽ ✽

Chapter 6

Multiblock Integration Method

In the medical field, statistical models are often used to discover and identify variables of interest that influence the development or status of the studied pathology. These statistical models not only predict the outcome but also links a set of variables of interest to it. However, genes interact in complex patterns. Some genes regulate others, and some genes are co-expressed with each other.

These interaction patterns are described in gene-gene interaction graphs and can be used in multiple ways in the biostatistical field. Including these graphs in the statistical models is one way to use them, highlighting not only the variables of interest but also the biological interaction of interest.

Furthermore, a new class of statistical models emerged as multiple measurements became available for the same patients from various sources. These statistical models differentiate between the different sources of data in separate blocks. The goal was expanded to understand the interactions between the different sources of data.

This work focuses on the addition of a graph penalty to a Regularized Generalized Canonical

Correlation Analysis (RGCCA) model. The chapter starts with a mathematical definition of the RGCCA. Then we describe the variable selection process of the RGCCA. Later we present the RGCCA with a graph penalty.

Throughout the chapter, we define our models on a set of J blocks X (1) , X (2) , • • • , X (J) . Each block consists in data from the same sample of n patients and p (j) variables. The blocks are assumed to be centred and scaled.

Availability

Regularized Generalized Canonical Correlation Analysis

The RGCCA was proposed by [START_REF] Tenenhaus | Regularized generalized canonical correlation analysis[END_REF] and extended the classical Canonical Correlation Analysis (CCA) to a set of J blocks. The model aims at extracting shared information between the blocks using the following optimisation problem :

argmin

w (1) ,w (2) ,••• ,w (J) J j,k=1 j̸ =k -c j,k g 1 n w (j)⊤ X (j)⊤ X (k) w (k) s.t.w (j) ∈ Ω (j) RGCCA j = 1, • • • , J (6.1)
The coefficients c j,k are equal to 1 or 0 and define if the blocks k and j are connected. They can be seen as the coefficients of the adjacency matrix of an a priori design graph of interactions among the blocks.

The function g can be any real continuous differentiable function. Authors of the RGCCA propose the identity g(x) = x, the absolute value g(x) = |x| or the square function g(x) = x 2 . Only the identity penalises the negative correlation between the blocks.

Ω (j)
RGCCA defines a set of constraints on each block weights. Here Ω (j)

RGCCA = {w (j) ∈ R p (j) ; w (j)⊤ M j w (j) = 1}, where M j = τ (j) I p (j) + 1-τ (j)
n-1 X (j)⊤ X (j) . τ (j) is a shrinkage parameter between 0 and 1. τ (j) = 1 means the projectors (w (j) ) must have a unit norm, while τ (j) = 0 means the projections (X (j) w (j) ) must have a unit variance and the optimisation problem 6.1 aims at maximising the correlation between the blocks.

The weights w (j) are estimated by maximising the sum of the correlations between pairs of latent components. Cycle block coordinate ascent, described in Algorithm 1, is used to solve the optimisation problem defined in 6.1. This algorithm updates each block weight in turn, keeping the others fixed (line 8 in Algorithm 1). Using the Lagrangian multipliers, the solution of the optimisation problem for block X (j) , if all other block parameters are considered fixed, is given by the Equation 6.2.

Algorithm 1 RGCCA optimisation algorithm

Require: X (1) , X (2)

, • • • , X (J) , τ (1) , τ (2) , • • • , τ (J) , design matrix C Ensure: Optimal w (1)⋆ , • • • , w (J)⋆ 1: for j ∈ 1, • • • , J do 2: w (j)
0 is a random vector or the first singular vector of X (j) 3:

Ensure w (j) 0 ∈ Ω (j) RGCCA 4: end for 5: repeat 6:

for j ∈ 1, • • • , J do 7: z t+1 j = j k=1 c j,k g ′ w (j)⊤ t X (j)⊤ X (k) w (k) t+1 X (k) w (k) t+1 + J k=j+1 c j,k g ′ w (j)⊤ t X (j)⊤ X (k) w (k) t X (k) w (k) t 8: w (j) t+1 = 1 n M -1 j X (j) z j ∥ 1 n M -1 j X (j) z t+1 j ∥ 9:
end for 10: until ∥w

(1) t+1 -w (1) t ∥ ≤ ϵ, • • • , ∥w (J) t+1 -w (J) t ∥ ≤ ϵ w (j) = 1 n M -1 j X (j) z j ∥ 1 n M -1 j X (j) z j ∥ s.t. z j = J k=1 c j,k g ′ w (j)⊤ X (j)⊤ X (k) w (k) X (k) w (k) (6.2)

Sparse Generalized Canonical Correlation Analysis

As defined in equation 6.1, the RGCCA does not allow for the identification of variables of interest among a large set of variables in each block. To tackle this issue, the Sparse Generalized Canonical Correlation Analysis (SGCCA) was introduced [START_REF] Tenenhaus | Variable selection for generalized canonical correlation analysis[END_REF]. The authors extended the RGCCA by adding a sparsity constraint on the model weights, using an ℓ 1 penalisation. The set of constraints for the SGCCA becomes Ω (j)

SGCCA = Ω (j)
RGCCA ∩ {w (j) ∈ R p (j) ; ∥w (j) ∥ 1 ≤ s (j) }, where s (j) is a user-chosen parameter. Additionally, the authors also fixed the shrinkage parameter τ (j) in the RGCCA constraints to 1 for all the J blocks. Thus, the RGCCA constraints becomes Ω (j)

RGCCA = {w (j) ∈ R p (j) ; ∥w (j) ∥ 2 = 1}.
Since ∥w (j) ∥ 2 ≤ ∥w (j) ∥ 1 ≤ p (j) ∥w (j) ∥ 2 holds for all vectors, Ω (j) SGCCA is not empty if and only if s (j) ≥ 1. Additionally, the ℓ 1 is active if and only if s (j) ≤ p (j) .

Using the Lagrangian multipliers, the authors demonstrated that the optimal solution has the form

w (j) = S( 1 n X (j)⊤ z j ,λ (j) ) ∥S( 1 n X (j)⊤ z j ,λ (j) )∥ 2
, where S denotes the soft-thresholding operator and λ (i) is chosen such that ∥w (j) ∥ 1 ≤ s (j) . The optimisation algorithm is similar to the one introduced in Algorithm 1, where the update of w (j) , in line 7, takes the new form.

netSGCCA

This section describes a modification to the SGCCA that injects prior graphical knowledge into the model. Given a graph G = (E, V), where E is the set of a block X (g) variables, we aim at imposing the graph structure on the block coefficients. The objective is to reflect the basic intuition that interacting variables (neighbouring variables in the graph) should contribute similarly to the model.

Following the work of [START_REF] Guigui | Network regularization in imaging genetics improves prediction performances and model interpretability on alzheimer's disease[END_REF], we propose to extend the SGCCA to netSGCCA by adding a GraphNet penalty on one block, denoted X (g) , relaxing the ℓ 2 equality constraint and using a gradient descent method to optimise the problem. The extension of our proposal to penalise more than one block is straightforward. Given the graph G = (E, V), we introduce its Laplacian L G (or its normalised Laplacian) into the SGCCA optimisation problem stated, which then becomes:

argmin

w (1) ,w (2) ,••• ,w (J) J j,k=1 j̸ =k -c j,k cov X (j) w (j) , X (k) w (k) + γ G λ max w (g)⊤ L G w (g)
s.t. w (j) ∈ Ω (j)

netSGCCA j = 1, • • • , J (6.3)
The binary symmetric design matrix C = (c j,k ) is the same as in the RGCCA optimisation problem.

The constraints sets Ω

(j) netSGCCA = {w (j) ∈ R p (j) ; ∥w (j)⊤ ∥ 2 ≤ 1} ∩ {w (j) ∈ R p (j) ; ∥w (j) ∥ 1 ≤ s (j) }.
As in the SGCCA optimisation problem, the s (j) are the sparsity parameters.

netSGCCA 73

The Karush-Kuhn-Tucker conditions state that if the inequality constraint of an optimisation problem is active, the optimal solution necessarily satisfies the corresponding equality constraint. As explained by [START_REF] Witten | A penalized matrix decomposition, with applications to sparse principal components and canonical correlation analysis[END_REF], if the chosen sparsity constraints s (j) lead to an active ℓ 2 inequality constraint, the constraint sets Ω netSGCCA and Ω SGCCA are equivalent. We kept the relaxation of the ℓ 2 equality (i.e. we use the inequality in our formulas) as it leads to a convex optimisation problem with a convex set of constraints.

GraphNet

The Laplacian matrix associated with a graph G = (E, V) is defined as L G = D -A, where D is the diagonal degrees matrix, and A is the graph adjacency matrix. Using these matrices and for any vector w of size p, the GraphNet penalty takes the form :

w ⊤ L G w = p j,k=1 (k,j)∈V (w j -w k ) 2 (6.4)
The Normalised Laplacian associated with the graph is defined as

L n G = D -1 2 L G D -1 2 .
Then, with d j the degree of the node j, the GraphNet penalty becomes :

w ⊤ L n G w = p j,k=1 (k,j)∈V ( w j d j - w k √ d k ) 2 (6.5)
Looking at these formulations, we can see that GraphNet explicitly attempts to give similar weights to variables connected in the graph. Using the normalised Laplacian, one can expect variables with more neighbours to have higher weight values, which can be more stable. The higher weights for higher connected variables can be explained because, supposing two variables j, k with d j ≥ d k , GraphNet will try to impose

w j = √ d j √ d k w k , with √ d j √ d k ≥ 1.
GraphNet uses a symmetric, semi-defined positive matrix. Therefore, the largest eigenvalue λ max of this matrix defines the upper-bound of the penalty (and its gradient), i.e., w ⊤ L G w ≤ λ max |w| 2 2 . The hyper-parameter γ G , introduced in the equation 6.3 defines the importance of GraphNet in the objective function. By dividing the GraphNet penalty by its largest eigenvalue, γ G is comparable between different Laplacians, which will make clearer the benchmarks of graphs we will present thereafter.

The GraphNet penalty has a grouping effect, which means that the weights of linked variables are brought closer. To prove the grouping effect, we can rewrite the equation 6.3, using the lagrangian multipliers λ (j) 1 and λ (j) 2 respectively associated with the ℓ 1 and ℓ 2 constraints for each block: argmin

w (1) ,w (2) ,••• ,w (J) h(w (1) , w (1) , • • • , w (J) ) h(w (1) , w (2) , • • • , w (J) ) = J j,k=1 j̸ =k -c j,k cov X (j) w (j) , X (k) w (k) + γ G λ max w (g)⊤ L G w (g) + J j=1 λ (j) 1 ∥w (j) ∥ 1 + J j=1 λ (j) 2 ∥w (j) ∥ 2 2 (6.6)
For two variables u and v, only connected to each other in the graph and positively correlated, we can show (see Appendix B) that:

|w (g) u -w (g) v | ≤ Q 2(1 -corr(X (g) u , X (g) v )) 2 γ G λmax + λ (g) 2 (6.7)
Q is a constant only dependent on the correlation matrix of the other blocks. The proof can be found in the annex, and it follows the demonstration established by [START_REF] Li | Network-constrained regularization and variable selection for analysis of genomic data[END_REF] for the grouping effect of the GraphNet in the context of the regression.

Optimisation

As proposed by [START_REF] Löfstedt | Structured variable selection for regularized generalized canonical correlation analysis[END_REF], we used the Fast Iterative Shrinkage-Thresholding Algorithm (FISTA) [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF] to find the optimal solution to our optimisation problem. The FISTA is built on the Iterative Shrinkage-Thresholding Algorithm (ISTA), which is designed to solve optimisation problems of the form:

argmin x {F (x) = f (x) + g(x)} (6.8)
where the functions f and g verify: f is a smooth real convex function continuously differentiable with Lipschitz continuous gradient (i.e. ∥∇f (x) -∇f (y)∥ ≤ L(f )∥x -y∥ for every x, y), and g is a real continuous convex function and possibly non-smooth.

Instead of the classic gradient descent, the ISTA uses the equivalent proximal formulation. Starting from an initial point x 0 , the ISTA generates a sequence x k verifying:

x k+1 = prox g (x k -α k ∇f (x k )) (6.9)
where the proximal operator is defined as prox g (x) = argmin y {g(y) + 1 2 ∥x -y∥ 2 2 } and α k is a positive real number.

To accelerate the convergence rate of the ISTA, the FISTA uses Nesterov's accelerated gradient descent. The authors have shown that the FISTA converges at a rate O( 1 k 2 ), and it is described as follows, with θ 0 = 0:

θ k+1 = 1 + 1 + 4θ 2 k 2 β k+1 = 1 -θ k θ k+1 x ′ k+1 = prox g (x k -α k ∇f (x k )) x k+1 = (1 -β k+1 )x ′ k+1 + β k+1 x ′ k (6.10)
We solve the optimisation problem defined in 6.6 sequentially, i.e. we find the optimal solution for each block while fixing the parameters related to the other blocks. For a block i, fixing the other blocks parameters results in the following optimisation problem:

6.3. netSGCCA 75
F (w (i) ) = f (w (i) ) + g(w (i) ) f (w (i) ) = J k=1 k̸ =i -c i,k cov X (i) w (i) , X (k) w (k) + γ G λ max w (i)⊤ L G w (i) g(w (i) ) = λ (i) 1 ∥w (i) ∥ 1 + λ (i) 2 ∥w (i) ∥ 2 2 (6.11)
The function f and g satisfy the FISTA conditions and the optimisation algorithm is described in 2. Since f has a Lipschitz continuous gradient, the gradient descent step α k is fixed at α k = 1 L with L the Lipschitz constant of the gradient.

Algorithm 2 netSGCCA optimisation algorithm

Require: X (1) , X (2) , • • • , X (J) , L G , γ G , design matrix C Ensure: Optimal w (1)⋆ , • • • , w (J)⋆ 1: for j ∈ 1, • • • , J do 2: w (j)
0 is a random vector or the first singular vector of X (j) 3:

Ensure w (j) 0 ∈ Ω (j) netSGCCA 4: end for 5: repeat 6:

for j ∈ 1, • • • , J do 7: repeat 8: θ k+1 = 1+ √ 1+4θ 2 k 2 9: β k+1 = 1-θ k θ k+1 10: w (j) ′ k+1 = prox g (w (j) k -α k ∇f (w (j) k )) 11: w (j) k+1 = (1 -β k+1 )w (j) ′ k+1 + β k+1 w (j) ′ k 12:
until ∥w

(j) k+1 -w (j) k ∥ ≤ ϵ 13: end for 14: until ∥f k (w (1) ) -f k+1 (w (1) )∥ ≤ ϵ, • • • , ∥f k (w (J) ) -f k+1 (w (J) )∥ ≤ ϵ

The proximal operator

In the netSGCCA, the proximal operator can be written as proj g (x) = argmin y∈Ω netSGCCA { 1 2 ∥x-y∥ 2 2 }. This defines a projection into a convex set problem, which can be solved using the Dykstra projection algorithm [START_REF] Bauschke | Dykstra's alternating projection algorithm for two sets[END_REF].

The Dykstra algorithm is designed for the projection on the intersection of two convex sets. Let S and R be two convex sets, with proj S and proj R their respective projection operators, the Dykstra algorithm generates a sequence x k converging into the projection of x onto S ∩ R. Starting from x 0 = x and with p 0 = q 0 = 0, the algorithm is described as follows:

y k = proj S (x k + p k ) p k+1 = x k + p k -y k x k+1 = proj R (y k + q k ) q k+1 = y k + q k -x k+1 (6.12)
In the netSGCCA, S = {w (j) ∈ R p j ; ∥w (j) ∥ 1 ≤ s (j) } and R = {w (j) ∈ R p j ; ∥w (j) ∥ 2 ≤ 1}.

The proj S (x) is the soft-thresholding operator defined as:

(proj S (x)) i = sign(x i )max(|x i | -λ, 0) (6.13)
Where λ is a positive scalar chosen to satisfy the constrain set S.

The proj R (x) is defined as:

proj R (x) =    1 ∥x∥ 2 x ∥x∥ 2 ≥ 1 x otherwise (6.14)
The Dykstra algorithm converges to the point in the intersection between the two sets nearest to the starting point [START_REF] Bauschke | Dykstra's alternating projection algorithm for two sets[END_REF]. The netSGCCA algorithm is implemented in python and accessible at https://github.com/neurospin/netSGCCA. It depends on the parsimony library from which it inherited the optimisation framework pylearn-parsimony package https://github.com/ neurospin/pylearn-parsimony (de [START_REF] De Pierrefeu | Prediction of activation patterns preceding hallucinations in patients with schizophrenia using machine learning with structured sparsity[END_REF][START_REF] Lofstedt | Structured {Variable} {Selection} for {Regularized} {Generalized} {Canonical} {Correlation} {Analysis}, {The} {Multiple} {Facets} of {Partial} {Least} {Squares} and {Related} {Methods}[END_REF]. .

Application

In the previous part (see 6.3), we introduced netSGCCA, an extension of the SGCCA with an a priori graph. Graphs have already been used into non-multiblock statistical models as a penalty over the model parameters. Graphs were used in supervised models, such as survival models [START_REF] Zhang | Network-based survival analysis reveals subnetwork signatures for predicting outcomes of ovarian cancer treatment[END_REF] and regression models [START_REF] Li | Network-constrained regularization and variable selection for analysis of genomic data[END_REF], and also in unsupervised models, such as matrix factorisation models [START_REF] Zhu | Network-based cancer genomic data integration for pattern discovery[END_REF]. However, it is still unclear how these graphs interact with the models and how they impact the feature selection. Additionally, with the availability of different graphs encoding different information, one can wonder how the choice of the graph meaningfully impacts the results obtained.

In this part, we present a study of the netSGCCA. First, using simulated data, we investigate the ability of netSGCCA to recover variables of interest in various conditions on the graph and expected solutions sides. The aim here is not to compare the different conditions but to study the behaviour of the graph penalty. Then, on the TCGA-LGG dataset, we compare the grouping effect and stability of netSGCCA while considering different graphs with different properties and from various bio-informatics sources. Finally, we investigate the relationship between selected features and the disease outcome in the same real dataset using survival prediction and pathway enrichment analysis.

Simulated Data

We tested netSGCCA using two simulated blocks X 1 and X 2 with a graph penalty based on a graph G on the second block. Our simulation procedure followed Du's proposal [START_REF] Du | Detecting genetic associations with brain imaging phenotypes in Alzheimer's disease via a novel structured SCCA approach[END_REF]. We started by defining the vectors u 1 and u 2 of dimensions p 1 = 150 and p 2 = 100. Then we generated n = 80 samples for each row of the two blocks

x 1 |z ∼ N czu ⊤ 1 , Σ 1 (respectively x 2 |z ∼ N czu ⊤ 2 , Σ 2
), with z ∼ N (0, 1) a latent variable. We defined (Σ 1 ) kl = 0.1, and (Σ 2 ) kl = -0.9 × |u k -u l | + 0.9 if the variables k and l are adjacent in the graph G , and 0.1 otherwise. The variance of each vector is 1. c is a scalar defining the signal-to-noise ratio, in our case c takes the values 0.5 and 2.

Application 77

The vector

u 1 = (0, • • • , 0 60 , 1, • • • , 1 30 , 0, • • • , 0 60
) was used in all configurations.

We tested 12 configurations consisting of a vector u 2 and a graph. Three different cases of u 2 were used to simulate different interaction types between variables of interest. The first case,

u 2 = (0, • • • , 0 40 , 1, • • • , 1 20 , 0, • • • , 0 40
), represents a set a variables having similar contribution to the observations and all connected to each other. The second case,

u 2 = (0, • • • , 0 40 , 1, • • • , 1 10 , 0, -1, • • • , -1 10 , 0, • • • , 0 39 ),
represents two sets of variables having opposite contribution to the observations. And the third case,

u 2 = (0, • • • , 0 40 , 1, -1, 1 • • • , -1, 1 20 , 0, • • • , 0 40
), represents a two sets of variables having opposite contribution to the model and connected to each other via a non-selected variables.

Four different graphs were investigated, the path (where the edges are between subsequent variables), the star graph (where the 50th variable is connected to all the others), the union of the path and the star graph and finally, the complete graph. An illustration of the different cases can be found in Figure 6

Practical implementation remarks

We assume that the graph edges represent interactions between the variables, which are modelled as high correlations here. The way we stated the simulation has to be further discussed in terms of numerical implementation. Suppose a graph G 1 = ({x 1 , x 2 , x 3 }, {(1, 2), (2, 3}), a graph of three variables where the second variable is connected to the first and third, as illustrated in figure 6.2.

This scheme clearly states that x 1 and x 3 must be correlated. Thus, the defined Σ 2 is not a plausible correlation matrix 'as is', in this case.

More generally, the defined Σ 2 matrix is not guaranteed to be semi-definite. Instead, we used their nearest semi-definite matrices, relative to a weighted W Frobenius Norm, approximated by Higham algorithm [START_REF] Nicholas | Computing the nearest correlation matrix-a problem from finance[END_REF]. Having a symmetric matrix A ∈ R p×p , the Higham algorithm approximates :

γ (A) = min{||A -X|| s.t X ∈ U ∩ S} (6.15) With S = {Y ∈ R p×p s.t Y = Y ⊤ and Y ≥ 0} the set of symmetric positive semi-definite matrices. And U = {Y ∈ R p×p s.t Y = Y ⊤ and y ii = 1, i = 1, • • • , p} the set of
symmetric matrices with only ones in its diagonal. Higham used the Dysktra projection algorithm presented earlier 6.12.

The projection onto the set U, relative to the W norm, is defined as :

proj U (A) = A -W -1 diag(θ i )W -1 s.t (W -1 • W -1 )θ = diag(A -I) (6.16)
And the projection onto the set S, relative to the W norm, is defined as :

proj S (A) = W -1 2 ((W 1 2 AW 1 2 ) + )W -1 2 (6.17) Such that A + = Q ⊤ diag(max(λ i , 0))Q
, where the λ i and Q are respectively the eigenvalues and the eigenvectors of A. In our case, we chose to work with W = I. And θ = diag(A -I). And we refer to Σ1 (respectively Σ2 ) as the corrected correlation matrix.

Results

If we rewrite our simulation protocol, the expected covariance is :

z ∼ N (0, 1) ϵ 1 ∼ N (0, I p 1 ×p 1 ) ϵ 2 ∼ N (0, I p 2 ×p 2 ) x 1 = czu ⊤ 1 + Σ-1 2 1 ϵ 1 x 2 = czu ⊤ 2 + Σ-1 2 2 ϵ 2 then cov(w ⊤ 1 x 1 , w ⊤ 2 x 2 ) = w ⊤ 1 u 1 u ⊤ 2 w 1 (6.18)
Using the Cauchy-Schwarz, one can see that the expected covariance is maximum when w 1 = u 1 and w 2 = u 2 . Thus, the model performance can be assessed on its ability to recover the variables Table 6.1: Recovering performances depending on configurations defined by the different cases defined by the vector u 2 and graphs. Corr is the correlation between the estimated components. Precision, Recall and F1 correspond to the evaluation of u 2 against the computed weights. Bold refers to highest values between netSGCCA and SGCCA. Low mean to variance ratio (c = 0.5).

netSGCCA SGCCA γ G Corr Precision Recall F1 Corr Precision Recall F1
Case Path 10 -4 0.56 ± 0.05 0.2 ± 0.11 0.17 ± 0.1 0.18 ± 0.1 0.46 ± 0.04 0.22 ± 0.32 0.03 ± 0.03 0.04 ± 0.05 Star 1 0.6 ± 0.04 0.13 ± 0.07 0.3 ± 0.15 0.18 ± 0.09 0.48 ± 0.04 0.02 ± 0.09 0.01 ± 0.02 0.01 ± 0.04 Union 10 -3 0.58 ± 0.05 0.13 ± 0.07 0.26 ± 0.15 0.17 ± 0.1 0.47 ± 0.04 0.09 ± 0.25 0.01 ± 0.03 0.02 ± 0.12 Complete 10 -2 0.39 ± 0.04 0.0 ± 0.0 0.0 ± 0.0 0.0 ± 0.0 0.4 ± 0.04 0.0 ± 0.0 0.0 ± 0.0 0.0 ± 0.0 0.39 ± 0.05 0.0 ± 0.0 0.0 ± 0.0 0.0 ± 0.0 0.39 ± 0.05 0.0 ± 0.0 0.0 ± 0.0 0.0 ± 0.0 Table 6.2: Recovering performances depending on configurations defined by the different cases defined by the vector u 2 and graphs. Corr is the correlation between the estimated components. Precision, Recall and F1 correspond to the evaluation of u 2 against the computed weights. Bold refers to highest values between netSGCCA and SGCCA. High mean to variance ratio (c = 2). 1.0 ± 0.0 0.19 ± 0.07 0.31 ± 0.1 Complete 10 -4 0.36 ± 0.09 0.02 ± 0.07 0.07 ± 0.23 0.03 ± 0.11 0.39 ± 0.09 0.05 ± 0.22 0.01 ± 0.07 0.02 ± 0.1 of interest and maximise the covariance between extracted latent variables. To do so, we computed the precision, recall and F1 metrics between the true u 2 vectors and the weights w j estimated by the model. For each configuration, we chose the hyper-parameter γ G by running the model 20 times, with γ G ranging from 10 -4 to 10 4 each time. The best γ G was selected using the best average F1 score because our objective is mainly to recover the variables of interest. We also ran the model for each configuration without using a graph and compared the results. The sparsity value was fixed to √ 25

Case

for the first block (resp. √ 20 for the second block) in all runs.

Table 6.1 and Table 6.2 show the results obtained on the simulated data. It shows that when we have a high mean/variance ratio, the F1 scores are higher, which means that the models using the graphs focus on retrieving the underlying projector u 2 . This was expected since a low mean/variance ratio means noisier data.

Overall, the tables also show that using netSGCCA outperformed the SGCCA without a graph.

When c = 2 (mean to variance ratio), SGCCA selected very few variables, about 3, leading to a high precision but very low recall. In contrast, the graph penalisation allowed the model to select more variables, retrieving all the variables of interest and a high F1 score. However, this increase of the F1 score came with a slight decrease in the correlation between the estimated components, by around 2%. Additionally, when c = 0.5, the F1 score continues to show improvement when the graph penalisation is used, but to a lower degree. However, the correlation between the estimated components also increased by 0.13 on average.

Comparing each graph type, we can see that, when c = 2, the path graph recovers all the variables of interest perfectly, with an average F1 score of 1 in all cases. The star graph was also able to obtain a perfect recall but with much lower precision. This is because the star graph selects many more variables, about 45 on average. Knowing that the sparsity level is the same for all configurations, the hub in the star graph seems to spread the weights more into its neighbours compared to the path graph. However, the correlations between estimated components are comparable. When c = 0.5, the models seem to select the variables randomly, which is shown by an F1 score close to 0.2. Additionally, the weights do not seem to resemble the original u 2 . However, even this result is better than without the graph a priori, which only selected a couple of features and resulted in an F1 score close to 0. In this situation, by choosing a greater number of variables, the star graph performed better in the correlation score compared to the path graph. Additionally, the union of the star and path graph exhibited behaviour similar to the ones of path and star graphs. Finally, the models with the complete graph always failed to outperform all the other models. This result is expected since the complete graph does not contribute to bring any information.

If we fix the graph and the mean-variance ratio, for all the cases u 2 considered, we observe no significant difference in the precision of the variable selection process nor in the extracted correlations.

This observation holds for all graph types and mean-to-variance ratios. The correlations between neighbours in the graph did not change the selected variables.

Overall, netSGCCA seemed to outperform the SGCCA in terms of retrieving the variables of interest, in nearly all configurations tested. It appeared that it is through its properties and structures that the graph have an influence on the behaviour of the model. We seek to investigate these results on real oncological data in the next sections. To assess the behaviour of the proposed method on real data, we used the TCGA-LGG dataset.

TCGA-LGG analysis design and results

The dataset includes 419 patients. All patients have CNV, which has 57964 variables, miRNA with 645 variables, and RNA with 22297 variables.

netSGCCA (and SGCCA) is an unsupervised model that jointly reduces the dimension of the different blocks. We want learned features that not only maximise the correlation between different blocks but are also relevant to survival. [START_REF] Mark | Microarray gene expression data with linked survival phenotypes: diffuse large-Bcell lymphoma revisited[END_REF] previously demonstrated that solving a sparse Cox model is equivalent to solving a linear regression when the null deviance residuals replace the survival. The null deviance residuals are derived from the highly skewed martingale residuals by a normalising transformation analogous to deviance residuals for Poisson regression :

D i = sign( Mi ) -Mi -δ i log((δ i -Mi )/δ i ) ≈ δ i -Êi Êi (6.19)
Where δ i indicates whether or not an event is observed for subject i, Êi is the expected number of events at time t i (event time for individual i), and Mi = δ i -Êi is the martingale residuals.

Using this equivalency, [START_REF] Bastien | Deviance residuals based pls regression for censored data in high dimensional setting[END_REF] proposed to predict survival using the PLS regression model by replacing the outcome with the null deviance residuals. Similarly, we added the null deviance residuals into the netSGCCA (and later to SGCCA). We aim to make the model learn to correlate latent features extracted from the other blocks to survival. Note that this extra block is added only during training and does not intervene during the inference stage.

We applied the GraphNet penalty on the RNA block. Gene identifiers were mapped to their HUman Genome Organisation (HUGO) names with non-matching genes removed, leading to 19864 genes remaining in the RNA block.

Different gene-gene interaction graphs with different properties were used. Pathway Commons v12 (PC) [START_REF] Cerami | Pathway commons, a web resource for biological pathway data[END_REF] is an aggregation of multiple subgraphs from various sources, containing 15710 nodes and 841690 edges. The Molecular Signature Data Base (MSIGDB) [START_REF] Liberzon | Molecular signatures database (MSigDB) 3.0[END_REF][START_REF] Subramanian | Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles[END_REF] C3 regulatory target gene set is one of the subgraphs of the PC graph containing 10463 nodes and 82962 edges. The Kyoto Encyclopedia of Genes and Genome (KEGG) [START_REF] Kanehisa | KEGG: Kyoto encyclopedia of genes and genomes[END_REF] is also another subgraph of the PC graph, containing 776 nodes and 11963 edges. The full description of these graphs can be found in the Table 6.3. Since graphs do not include all genes, missing genes were added as isolated nodes, and genes mapped to the same HUGO name were duplicated in the graph with their edges.

We aimed to compare the gene selection abilities of the different graphs. First, using the PC graph, we established the differences between the raw and normalised Laplacian, in terms of the number of selected variables, as the γ G varies, and the stability of variable selection. Then, we compared the different graphs. Later, to establish the impact of the graph in the selection process and disentangle the graph penalty effect from the RGCCA model effect, we permuted the nodes of the PC graph and compared the results with the non-modified graph. We also removed edges from the PC graph to investigate the effect of the graph density and the importance of the edges between selected nodes.

Finally, we evaluated the application of the features extracted by using them to predict survival and discussed the biological pathways potentially involved through the selected genes. Details of the experimental design are presented in Figure 6.3.

In order to achieve the experimental design, we first isolated 15% of the patients (63 patients) as a test set and performed the training on the remaining 356 patients. We stratified the patients using the event status. The training set was split into five folds of equal sizes. This allowed us to compute a standard deviation for the c-index metric, and identify two sets of selected genes: the candidate selected genes (the union of gene sets selected in each fold) and stable selected genes (the intersection of gene sets selected in each fold). The test set was only used for the survival prediction in the final part of our study.

Throughout this work, c j,k = 1 for all j and k, with j ̸ = k. This means that each block is connected to all the other blocks in the netSGCCA model. Additionally, the ℓ 1 constraints have been fixed as the best yielding parameters -in terms of the c-index -on the five folds without GraphNet using a grid search.

Comparison between normalised and raw Laplacian

To compare the normalised Laplacian with its raw version, using the PC graph, we varied the γ G between 10 -4 and 10 4 . Figure 6.4(a) shows the evolution of the number of selected genes, for both Laplacians, as γ G increases. We can see that the higher the γ G , the more genes are selected; this is in line with previous findings in the related works. It also shows that, without a graph penalty, very few genes are selected (about 3 in each fold). The Figure 6.4(a) also highlights a similar behaviour between the raw Laplacian and the normalised Laplacian, in terms of the number of genes selected.

However, a closer inspection shows that the normalised Laplacian is more stable as γ G increases. To

show this, we defined the stability as the # of stable selected genes # of candidate selected genes and computed it for each γ G . Results are shown in Figure 6.4(b). Figure 6.4(c) shows the distribution of the number of folds in which a gene was selected, when γ G was fixed to 10 2 . In the case of the normalised Laplacian, around 51% of genes selected by the model were chosen in all five folds. In contrast, only 30% were selected in all five folds with the raw Laplacian. In fact, a large number of genes were selected only once.

We examined the degree distributions by considering the subset of selected genes across five folds in the full PC graph on the one hand and, on the other hand, in the PC sub-graph containing the selected genes. Figure 6.5(a) shows the degree distribution of selected genes in the full PC graph.

Both the normalised and raw Laplacian follow the same behaviour. The two distributions are also similar to the degree distribution considering all genes in the PC graph. Considering the proportion of isolated nodes in PC, they are underrepresented among the selected genes. This indicates that GraphNet does not favour genes with a high degree in the graph but discriminates against isolated nodes. However, GraphNet did not select neighbours in the graph. This is exhibited in Figure 6.5(b), which gives the degree distribution in the PC sub-graph of the selected genes. It shows a shift toward lesser degrees suggesting that genes are mostly selected because of their covariance similarity across the patients, but not because of their neighbourhood in the graph penalty (also shown later). The GraphNet penalty, in our configuration, seems to mainly smooth weights over similarly correlated genes but does not give rise to gene communities because of their adjacency in penalty graph.

( We looked at netSGCCA gene weight distribution according to the degree of the gene in the PC graph. Figure 6.6 shows that the gene degree does not influence the weight. However, isolated nodes tend to have significantly higher weight variance, ranging from -0.1 to 0.1 for the normalised Laplacian. This is expected as the more a node is connected, the more its weight is constrained. The same pattern can be shown for the normalised and raw Laplacian, except that the normalised Laplacian produced smaller weights in absolute value. Additionally, Figure 6.7 shows the weight difference between genes according to the distance between the genes in the graph (we used the shortest path in the graph). It shows that the closer the genes, the closer their final weights, going from 0.016 on average between neighbouring nodes to 0.017 if the shortest path between them is five, for the normalised Laplacian.

The same pattern can be seen on the raw Laplacian with higher values. Thus, even if the graph does not select sub-networks, it has a grouping effect.

Finally, the model seems to select genes by correlation, which is exhibited in Figure 6.8 which shows correlation distribution among selected genes by the normalised Laplacian. The model selected two distinct, negatively correlated groups. The correlation distribution between selected genes does not follow the distribution of the correlation among gene expression profiles in the whole LGG dataset.

The GraphNet penalty seems to partially override the ℓ 1 penalty and allow the model to select groups of highly correlated variables. Figure 6.7: Average absolute weight difference between selected nodes according the Dijkstra distance between them in the PC graph Figure 6.8: Correlation distribution of selected genes. Between stable selected genes ( Intersection of selected genes in the 5 folds ), Candidate selected genes ( Union of selected genes in the 5 folds ) and All the genes in the dataset.

Comparisons between different graphs

Since the normalised Laplacian allowed the selection of more stable gene sets, the following experiments will use only normalised Laplacians. We ran our model using the MSIGDB and KEGG graphs on the same five-folds and the same model hyper-parameters. We compared the number of genes selected by each graph and the overlap of selected genes between each graph and netSGCCA using the PC normalised Laplacian. More precisely, we compared the overlap between the stable selected genes and the candidate selected genes. Figure 6.9 shows the results obtained using the MSIGDB and KEGG Laplacian. Both graphs selected fewer genes than the PC graph, with the MSIGDB selecting more than the KEGG. However, there is a large overlap between the selected genes given each graph penalty.

This suggests that the graph structure strongly influences the grouping effect of the method, but also suggests that it does not substantially impact the set of genes selected. To better exhibit the influence of the graph semantic (information embedded in the graph), we randomly permuted gene labels within the PC graph, and using this new graph, we extracted genes in terms of stable or candidate genes. This was done ten times to check the validity of the results.

For the stable and candidate sets, we computed the Intersection Over Union metric (IOU) from each permutation and the selected genes using the PC graph. The different permuted graphs selected a similar number of candidate selected genes as the original PC graph, but the IOU was always below 75% as shown in Figure 6.10. Same applies for the stable selected genes. This illustrates that, for a biological graph of a given density, its semantics is weakly passed by the GraphNet penalty and has a limited influence on the selected genes. Table 6.4: The effect of pruning edges that connect genes selected when using the full PC graph. Candidate and stable genes sets obtained when using the different pruned sub-graphs.

Ref Figure 6.11: The evolution of the number of selected genes when the number of edges decreases (x-axis is the percentage of edges removed from 0% to 90%.)

To show the impact of the graph density on the number of selected genes, we removed edges randomly from the PC graph. We made sure that we removed edges from candidate selected genes 6.4. Application 87 using the normalised PC graph in the same proportions as from the whole edge set. Figure 6.11 shows the evolution of the number of selected genes as edges are removed. It exhibits a strong correlation between graph density ( # of edges # of all possible edges ) and the number of selected genes. This is the case for both candidate and stable selected genes. Additionally, selected genes after the removal of edges are all included in the original set of selected genes, for both candidate and stable. This shows that GraphNet penalty achievement for feature selection depends on the density of the graph it used.

We investigated the importance of the direct and the indirect paths between the genes selected by a model with the GraphNet penalty based on the full PC graph. To do this, we removed inner edges (edges connecting two selected genes), outer edges (edges connecting a selected and non-selected gene), and all edges of the selected nodes, which would make them completely isolated. We considered the candidate and stable selected gene sets obtained with the PC normalised Laplacian, which produced 822 candidate selected genes from the five folds of which 423 were stable selected genes. As shown in Table 6.4, removing the inner edges (around 1450 edges) did not change either the number nor the set of genes selected. Removing the outer edges (around 65285 edges) diminished the number of selected genes while keeping a large overlap. Making all the originally selected genes isolated nodes reduce the number of selected genes by more than a half. Even when removing either the inner or the outer edges, some path may remain between the selected genes, which helps the model to retrieve them. However, making them isolated lowers their probability of being selected. In the previous part, we presented results regarding the variable selected in the various blocks of our multi-block setting, but netSGCCA also yields components for each blocks. Using the components extracted from the CNV, RNA, and the miRNA by the netSGCCA method, we performed a survival prediction using a simple Cox model. We compared the results using the different graphs and a model without a graph constraint, and we present the results in Table 6.5. It shows that models using the GraphNet penalty have similar C-index scores compared to the model without a graph. When using the MSIGDB and KEGG, results are slightly better compared to the other models on the validation and test sets. This indicates that the GraphNet penalty did not weaken the ability of the estimated components to predict patients' survival. However, the model without GraphNet could not reliably extract genes, which hindered its ability to identify pathways of interest. SGCCA, without a graph, selected about three variables per fold. It is in line with the results found so far, as the ℓ 1 penalty tends to select few representatives among highly correlated variables. GraphNet forced the model to gather these correlated variables. Using the different graphs for the GraphNet penalty did not change the C-index noticeably, which is also expected as the variable selection has shown a substantial overlap.

Survival Prediction

Adding the graph made the signatures from other blocks less stable. The more stable RNA variables are selected, and the fewer and less stable variables are observed from the CNV block. As a result, no stable variables are obtained from the CNV block when a graph was used, but it must be balanced by the fact that for each fold, the model selected about 150 variables when the normalised Laplacian was used (compared to 250 when no graph was used). Genes function together in biological pathways. In order to find the most significantly represented pathways by the selected genes, we used Enrichr [START_REF] Edward Y Chen | Enrichr: interactive and collaborative HTML5 gene list enrichment analysis tool[END_REF][START_REF] Kuleshov | Enrichr: a comprehensive gene set enrichment analysis web server 2016 update[END_REF][START_REF] Xie | Gene set knowledge discovery with enrichr[END_REF]. Enrichr uses the Fisher exact test to compute the significance of the enrichment of a pathway. This is a proportion test that assumes a binomial distribution and independence for the probability of any gene belonging to any set. We used the MSigDB C6 collection to investigate associations between signatures and the 423 selected genes.

Results are presented in Table 6.6. We found pathways that have already been associated with low-grade gliomas. Genes associated with astrocytes in the set CAHOY ASTROGLIAL, have also been selected by the model. These genes have previously been studied for their link with brain tumour development [START_REF] Irvin | Intrinsic astrocyte heterogeneity influences tumor growth in glioma mouse models[END_REF][START_REF] Katz | Astrocyte-specific expression patterns associated with the PDGF-induced glioma microenvironment[END_REF]. ATF2 is known to promote invasion in malignant glioma [START_REF] Zhang | miR-622 suppresses proliferation, invasion and migration by directly targeting activating transcription factor 2 in glioma cells[END_REF]. Additionally, TGF-β (Transforming growth factor-beta) has been targeted to limit brain tumour growth [START_REF] Han | TGF-β signaling and its targeting for glioma treatment[END_REF]. Other pathways have been labelled by enrichment analysis but the corrected p-value was not significant. For example, the RAF Fusion has been associated with pediatric low-grade tumours [START_REF] Lind | Novel RAF fusions in pediatric low-grade gliomas demonstrate MAPK pathway activation[END_REF]. Mutations in KRAS, HRAS and NRAS are known in gliomas and are often concomitant with BRAF mutations and fusions [START_REF] Knobbe | Mutation analysis of the ras pathway genes NRAS, HRAS, KRAS and BRAF in glioblastomas[END_REF].

Discussion

Our work establishes some characteristics of netSGCCA, a data integration method that implements the GraphNet penalty. It shows, in the context of multi-view analysis, that GraphNet helps to group variables instead of selecting a few candidates, which is in line with previous results involving the same penalty. We also exhibit better interpretation capabilities for netSGCCA compared to SGCCA, as it allows the selection of sound and stable candidate variables within a block. When netSGCCA implements the GraphNet penalty using the normalised Laplacian of the a priori graph and not the raw, it leads to an even greater stability.

From the application of netSGCCA to a real multi-modal oncological dataset, we have derived presumably general observations. For the block being submitted to the GraphNet penalty, the similarity of the variable profiles is the primary driver for co-selecting variables, and the proximity of the variables (as nodes) in the graph is secondary. The graph seems to present candidates that are grouped according to their correlations. The density of the graph a priori used in the GraphNet strongly influences the final number of selected variables. The denser the graph, the more variables were selected. Yet, our results also suggest that the capacity of the netSGCCA at extracting variables of interest capitalises mainly on the variables initially selected by the SGCCA (without a graph).

Finally, the large overlap between selected variables in the penalised block when using graphs with equivalent density but with different semantics, suggests the netSGCCA model relies on the data first.

Regarding the LGG pathology, as it may be studied from the TCGA multi-modal dataset, we demonstrate two remarkable achievements of netSGCCA that overrides the current performance of other multi-modal integration frameworks. First, netSGCCA predicts survival very well. Large c-index values were obtained without using the medical data (no eCRF information other than survival). The results found are similar to reference values reported in a recent work which did consider medical data [START_REF] Herrmann | Large-scale benchmark study of survival prediction methods using multi-omics data[END_REF]. Second, using netSGCCA with the PathwayCommon graph penalty for the gene expression block, we took advantage of the stability of variable selection to propose a list of candidate genes and candidate biological pathways that explain the pathological outcome. Overall, this suggests that graph penalty in multi-modal analysis model like netSGCCA is able to bring original molecular biology insights into the pathology.

A limitation of this study is the small number of applications considered, namely one simulation and the TCGA LGG dataset. While a general observation is that GraphNet selects more variables than the classical Elastic-Net, discussions remain about its stability when there are multiple variables of interest with no correlation between them. It is a data-related problem, and the model performance on the LGG dataset is insufficient to give indications about its behaviour in such a case.

Conclusion

The present work focuses on the analysis of the GraphNet penalty available in the mulit-block netS-GCCA model and applied to the TCGA-LGG dataset. Contrary to Elastic-Net alone, GraphNet penalty is able to select a reasonable set of genes and yields informative biological interpretation from the pathway enrichment analysis. The example on the TCGA-LGG dataset exhibits the stability and reliability of netSGCCA for selecting variables of interest. However, it is important to note that we show that the co-selection of variables is not primarily influenced by the structure of the graph, but rather by its overall density. Therefore, an interpretation in terms of the paths read in the graph is elusive. Nevertheless, the method did extract genes that have Our strategy using netSGCCA and SGCCA is similar to the one introduced in the previous chapter.

We used joint dimensionality reduction, followed by a Cox model. In this strategy, we consider the null Cox deviance residuals as an extra block, as proposed by [START_REF] Bastien | Deviance residuals based pls regression for censored data in high dimensional setting[END_REF]. This allows taking into account survival information in the unsupervised dimensionality reduction. We note that this survival block is used only during unsupervised training and not by the Cox model.

This section presents different state-of-the-art multi-block models, which will be compared to netSGCCA and SGCCA in our survival analysis. All used models aim at maximising the log-likelihood function presented in equation 2.7. We compared netSGCCA and SGCCA with multi-block supervised survival models, which directly infer the risk of each patient. We consider X (1) , • • • , X (J) of sizes ) the different blocks of data obtained on the same set of patients.

n × p (1) , • • • , n × p (J
We tested two supervised models, Priority-Lasso [START_REF] Klau | Prioritylasso: a simple hierarchical approach to the prediction of clinical outcome using multi-omics data[END_REF] (see also 3.3.2.1) and Block-Forest [START_REF] Hornung | Block forests: random forests for blocks of clinical and omics covariate data[END_REF] (see also 3.3.2.2), whose features will be detailed below. We chose these two algorithms because, in Herrmann's benchmark study [START_REF] Herrmann | Large-scale benchmark study of survival prediction methods using multi-omics data[END_REF], these two methods obtained the highest c-index among the compared methods. We note that, even if Block-Forest does not perform variable selection, we added it to our study as it presented the best results on the various cancers benchmarked.

Priority-Lasso [START_REF] Klau | Prioritylasso: a simple hierarchical approach to the prediction of clinical outcome using multi-omics data[END_REF] proposed a multi-block adaptation of a linear model, which includes the Cox model. The algorithm requires a priority order of the blocks. This order is userchosen and may reflect hypotheses about the blocks performances, or their importance to the issue at hand.

Let X (1) , • • • , X (J) be the different blocks of data. Here we suppose that the blocks are already ordered, with X (i) having a higher priority than X (i+1) . The algorithm sequentially fits a model on each block and uses the results of the previously fitted model as an offset. Since the outcome Y is a part of the data used to estimate the model parameters for each block, the model tends to be more optimistic when computing the offset. This results in underestimating the influence of lowerpriority blocks on the outcome. To tackle this issue, the authors propose to compute the offset using a cross-validation method: first, they split the training set into K equal-size parts S 1 , • • • , S K ; then, for individuals in part S k ; they estimated the model parameters using all parts of the training set except

S k and finally, they computed their offsets. The Priority Lasso is described in Algorithm 3.

For this study, we used the Priority-Lasso as implemented in the R package prioritylasso [START_REF] Klau | prioritylasso: Analyzing Multiple Omics Data with an Offset Approach[END_REF]. The chosen block order was RNA, Mutation, and Radiomics.

BlockForest Variable sampling in random forests is biased toward blocks with the largest number of variables. As a result, the small blocks are explored a small relative number of times, regardless of their relationship to the target function. BlockForest adapted the variable sampling to the multi-block structure of data. The various adaptation procedures have already been presented in Chapter 3. We chose to work with the BlockForest sampling method. The procedure starts by choosing a subset of the available J blocks by selecting each block with a probability of 0.5 (this leads to 0.5 J probability

Radio-genomics integration and survival prediction

Algorithm 3 Priority-Lasso algorithm Require: X (1) ,

X (2) , • • • , X (J) , Y , S 1 , • • • , S K 1: of f set(1, • • • , n) ← 0 2: for j ∈ 1, • • • , J do 3: model (j) ← train_model(X (j) , Y , of f set) 4: for k ∈ 1, • • • , K do 5: S c ← complement of S k 6: model ← train_model(X (j) [S c ], Y [S c ], of f set) 7: risk ← predict_risk(model, X (j) [S k ]) 8: of f set[S k ] ← of f set[S k ] + risk 9:
end for 10: end for that all blocks or no blocks are selected; in the latter case, the sampling is repeated). Then from each chosen block, a fixed number of variables is selected (the number of selected variables is p (j) , where p (j) is the number of variables of the block j). The criterion is then evaluated on the selected variables from each block as in the classical random forest algorithm. We used the BlockForest as implemented in the R package blockForest (Roman Hornung, 2022).

We only tested the default parameters except for the parameter splitrule (splitrule = Log-rank).

Experimental design

To conduct this study, we used the TCGA-LGG dataset. First, we assess the performance of each block when all data is available. We call this stage the baseline study. This is done to estimate the performance of each block when maximum data is available. Then, we compare the mono-block and multi-block models on only 83 patients having all blocks available. For this study, we used the RNA, Mutation profiles and Radiomics as the three blocks of our study. Figure 7.1 illustrates the described stages.

We used the C-index metric to assess the model performances. The datasets used were split into 80% training set and 20% test set. The training set was used to choose the hyperparameters for each model, including the α and r for the ElasticNet Cox model, the sparsity constraints s (j) and λ G for SGCCA and netSGCCA, and the α for Priority-Lasso. To do so, the training set was split into three equal-sized sets, and a three-fold cross-validation procedure was performed. The chosen parameters are those that yielded the best C-index while selecting less than 10% of variables. The latter condition was added because we are interested in sparse models.

To assess the model performance and evaluate each model variable selection (and its stability), we used ten bootstrapped samples (of equal size to the original training set) with repetition from the training set. Each bootstrapped sample is used to train a model, and its performance is assessed on the isolated test set. The most selected genes are examined from Enrichment Analysis using Enrichr [START_REF] Edward Y Chen | Enrichr: interactive and collaborative HTML5 gene list enrichment analysis tool[END_REF][START_REF] Kuleshov | Enrichr: a comprehensive gene set enrichment analysis web server 2016 update[END_REF][START_REF] Xie | Gene set knowledge discovery with enrichr[END_REF] in order to discover significantly enriched biological pathways. We used the Kyoto Encyclopedia of Genes and Genomes (KEGG), MSigDB hallmarks and MSigDB oncogenic signatures collections. We only considered significantly enriched pathways with a False Discovery Rate (FDR) p-value correction. We have ensured that the overall ratio of censored and non-censored data is kept across all samples and splits. Figure 7.2 illustrates the entire hyperparameter choice and model assessment procedure. We used the Cox model implemented in the glmnet R package [START_REF] Friedman | Regularization paths for generalized linear models via coordinate descent[END_REF]Simon et al., 2011b). Table 7.2: Selected radiomic features extracted from the T2 modality. Out of 10 bootstrap samples, we only kept the seven where the model did not fail. The table only shows radiomic features that were selected more than three times. The means and std were computed using only the samples where the feature was selected. 7.3 shows the most selected radiomic features. Most selected features are texture features computed from the GLCM, and only the Correlation computed from the GLCM of the LHL wavelet transformed image was selected in all configurations. The selected variables have lower coefficients compared to the models using the T2 radiomics, which is explained by the optimal low ℓ 1 ratio chosen for the FLAIR radiomics. As with the T2 modality, none of the most selected radiomic features has the sign of their coefficient changed throughout the different runs. but the instability of the models using the CE T1 radiomic features. Three of the ten samples led to a C-index above 0.7 on the test set, while four scored below 0.43 and thus worst than a random ranking. The models that used CE T1 radiomic selected few variables (about 12), most of them first-order features, and only four were chosen by four or more models, as displayed in table 7.4. The coefficients of the selected variables have small magnitudes, and the models are not stable in their feature selection. Comparing the three runs with the highest C-index on the test set, we observe that selected variables do not overlap except for the "First Order Kurtosis" extracted from the LLL wavelet transformed image. As with the T1, the CE T1 does not show oedema, which can explain why it had difficulties finding robust features. All modalities radiomics. Finally, we trained Cox models with all four modalities that gave similar results to those when using the T2 modality, with an average C-index of 0.80 for the validation set and 0.74 for the test set. The models trained on radiomic features extracted from all modalities failed to converge twice and returned the null model. Interestingly, only one sample resulted in the model using the T2 and all modalities not converging. This is even though both models were trained in the same sample, which suggests a dependence on the initialisation of the optimisation algorithm. The chosen parameters were α = 0.085 and r = 0.716, which resulted in 15 variables selected on average. As shown in Table 7.5, the T2 "First Order Kurtosis" extracted from the LLL wavelet transformed image was again the most selected variable, with similar coefficients. This suggests that this radiomic is a good predictor for survival in the case of LGGs. The most selected variables are first-order variables. Overall, the T2 and FLAIR modalities performed better than the T1 and CE T1. Using all modalities gave similar performances to those obtained with the T2, but gave different results in terms of variable selection. We hypothesize that the stability issues in the feature selection process are due, on one hand, to the ElasticNet penalty, which can lead to inconsistent results, and on the other hand, to the nature of the radiomic features, which by design contains a lot of redundant information.

For the coming comparisons with the multi-block analysis, we used the combination of all modalities because the performance was not degraded while containing more information than the T2 modality alone.

Somatic mutations

We investigated the baseline performance of the Cox model using the mutation matrix. Table 7.6

shows the results obtained when the Cox model was trained on the original or pre-processed mutation profile using the 419 patients from the TCGA-LGG. a standard deviation of 0.07. On average, the model selected 876 genes, and 40 genes were selected in all six runs. The list of the 40 selected genes can be found in Table A.1. Enrichment analysis using the Enrichr revealed the oxytocin signalling pathway as the only significantly enriched (with an adjusted p-value of 0.025894) from the Kyoto Encyclopedia of Genes and Genomes (KEGG). This pathway has already been linked to brain tumours and specially LGGs [START_REF] Liu | The identification of key genes and pathways in glioma by bioinformatics analysis[END_REF][START_REF] Yang | Identification of core genes and screening of potential targets in glioblastoma multiforme by integrated bioinformatic analysis[END_REF] Mutations Propagated according to the PC graph. The models using the mutation matrix after propagation using the Pathway Commons (PC) graph showed similar results to the raw mutation matrix on the validation set, with a C-index of 0.74, but generalised better on the test set with a Cindex of 0.70 on average. The chosen parameters were α = 0.07 and r = 0.85, which resulted in 43 variables selected on average. Using graph propagation helped identify similarities between patients that were not observable beforehand. However, only five genes were selected in more than half the runs. The most selected genes were only mutated in three or fewer subjects, as shown in Figure 7.5.

Without propagation, highlighting these genes would have been infeasible. Even though the model identified genes linked to gliomas, no gene was selected consistently, and Enrichr could not find any significantly enriched pathway.

Table 7.7: Selected genes using the mutation matrix propagated in the PC graph. The table shows genes that were selected more than 5 times. The means and std were computed using only the samples where the feature was selected.

Selected an average validation C-index of 0.71 and 0.70 for the test set. Additionally, the average number of selected genes is comparable, with 47 genes obtained with the parameters α = 0.26 and r = 0.32.

However, only the gene PALMD was frequently chosen for both approaches, as shown in Tables 7.8

and 7.7. Interestingly, the gene IDH1 was also selected when using the influence graph. This gene is strongly related to the disease outcome, as discussed earlier. Again, the most frequently chosen genes are only mutated in a small number of patients in the dataset, which shows the importance of graph propagation. Additionally, using the influence graph highlighted multiple pathways, as shown in table 7.9. These pathways involve the genes IDH1 and EGFR, which have already been studied for their link to brain tumours.

Table 7.8: Selected genes using the mutation matrix propagated in the influence graph. The table shows genes that were selected more than 5 times. The means and std were computed using only the samples where the feature was selected.

Comparison mono-block and multi-block approaches

This part is devoted to the comparison of the monovs. multi-blocks models in terms of performance and interpretation. A reduced sample set is used here consisting of the patients that have all the data available for the multi-block study: 83 subjects have the Radiomics, somatic Mutations and RNA data available in the TCGA-LGG dataset. The different sample sizes are summerized in Table 7.11.

Mono-block

We re-ran Cox models using each block (Radiomics, Mutations, RNA) independently to compare the different mono-block with the multi-block approaches. The obtained results from the three folds cross-validation and ten bootstrapping results on the test set are presented in Table 7.12 and Figure 7.7.

Table 7.12: Mono-block results obtained on the 83 subjects having the radiomics, somatic mutations and RNA data available. For the Radiomics, the concatenation of all the radiomics extracted from the four modalities has been used; The somatic mutations with graph propagation with the influence graph was used; Finally, the concatenation of all the blocks was also investigated Validation Test Radiomics block alone 0.69 ± 0.05 0.56 ± 0.13 Mutations block alone 0.81 ± 0.08 0.74 ± 0.13 RNA block alone 0.75 ± 0.10 0.59 ± 0.20 All blocks Concatenated 0.83 ± 0.14 0.53 ± 0.07 Using the same strategy to tune the models and find the optimal parameters described in the previous section, we compared the radiomics obtained from the four modalities collated together, the somatic mutations profile with the influence graph propagation, the gene expression levels, and finally the concatenation of the three previous blocks.

Radiomics block alone.

Using the radiomic features, the Cox model had above random performances on the validation sets with a C-index of 0.69 on average, with a standard deviation of 0.05. However, the model could not generalise on the test set. Seven of the ten bootstrapped runs on the test set failed to converge. On the three remaining runs, two models had a C-index above 0.8, while the remaining one had a C-index of 0.48. The models had an average C-index of 0.71 with a standard deviation of 0.15. The learnt sparsity and regularisation parameters were α = 0.02 and r = 0.25, respectively. This led to 105 variables selected on average (by taking into account only the three runs which did not return the null model). 230 unique features were selected in at least one of the three runs. Among them, 16 were chosen in the three runs and are described in Table 7.13. Most of the selected covariates are texture features from the GLZM. Most of the consistently selected variables were not identified in the baseline models except for four which have been chosen in the concatenated model. The distribution of the selected radiomic features among the four modalities is shown in Figure 7.8. It displays that the models selected more T2 features in each run. However, the features extracted from the T1 modalities are the most consistently selected, even if the baseline T1 features had poor performances. None of the FLAIR features was consistently selected.

Mutation block alone. Somatic mutation profiles pre-processed with the influence graph showed better results compared to all the other studied data sources. On average, the Cox models had a C-index of 0.81 on the validation and 0.74 on the test sets. The mutation profiles had high and stable results, similar to those obtained with baseline study that used the complete patient sample (n = 419 vs. n = 83). Among the ten bootstrapped runs, nine achieved a C-index above 0.7, while the remaining one only reached a C-index of 0.42. The chosen parameters were α = 0.72 and r = 0.05, which resulted in 1440 variables selected during each run. The stability of the performances of the Cox models using somatic mutation profiles was not reflected in the stability of chosen variables. 4680 genes were selected in at least one run. Among them, 1225 were selected in more than half of the runs which is a key part of this study. The netSGCCA model was also able to select genes and pathways already linked to the development of gliomas and the prognosis of LGGs.

The radiomic integration, using the netSGCCA, confirmed the results obtained in the previous chapter. The graph penalty allows the selection of more variables, and the correlation is the main driver of the grouping effect. Furthermore, the netSGCCA and SGCCA had similar performances when the Pathway Commons (PC) graph was used to construct the penalty. This was also a previously obtained result. Yet, the baseline results demonstrate that radiomic features can be good predictors of the survival of Lower Grade Glioma. The T2w modality had especially high performances. This has already been observed in the literature. Unfortunately, we could not recover obtained results in the baseline when using the radiomics on the 83 patients having all blocks available. The First Order Kurtosis from the T2w modality was the best predictor among all radiomic features. A similar result has already been established in the chapter of related works [START_REF] Liu | A radiomic signature as a non-invasive predictor of progression-free survival in patients with lower-grade gliomas[END_REF], though the feature was computed on a different wavelet-transformed image.

Overall, data availability restriction hindered the predictive performance of the blocks. Furthermore, the Cox model became less stable when data from only 83 patients were used. Additionally, we expected an overlap between the variables selected on the baseline and on the restricted dataset, however results obtained show otherwise. The netSGCCA showed promising results in term of survival prediction and variable selection.Therefore, further studies are required to validate the results. It is also necessary to perform this study on a larger dataset.

Conclusion

The netSGCCA showed to be a robust model able to select variables already linked to the studied pathology and interacting in relevant biological pathways. Adding the imaging features hindered the ability of the netSGCCA to extract features that can be used for survival prediction. However, the baseline results show that radiomic features extracted from the T2w modality can be strong predictors.

Thus, larger datasets are required to further investigate the potential of radiomic features. The usage of graphs was shown to be impactful both when used to propagate somatic mutations and as a penalty in the netSGCCA. Our results show that multi-block integration under the netSGCCA procedure can uncover candidate genes and biomarkers, which can be investigated in later studies. Our framework could be used to study the DIPG and discover new relevant biomarkers.

✽ ✽ ✽ ✽ ✽ ✽

Chapter 8

Conclusion

Our work was motivated by the study of the Pediatric Diffuse High Grade Gliomas (DIPG), a rare pediatric high grade glioma. We aimed to design statistical methods able to integrate imaging and genetic data. These methods must be interpretable, thus giving new insights into the genetic signature of the tumour and its link to imaging features.

Radiomic analysis requires predefined Regions of Interest (ROI) on the images at hand. For our DIPG cohort, manual segmentation of the tumour was not feasible. Additionally, no database has been built to train classical machine-learning algorithms to delineate the tumoural regions automatically.

This study focused on obtaining binary segmentations for the DIPG using only the FLAIR and T2w modalities, using models trained on Glioblastoma. We propose combining different simple models of detection and segmentation to obtain satisfying segmentation results on a cohort that contained differences compared to the training dataset regarding, among others, patient age, image quality and tumour type.

In parallel, our work aimed at building a multi-block integration model which takes into account known complex graphs of interactions between the genes. Additionally, our focus was on understanding how the graph of interactions influences variable selection. We propose netSGCCA, a model combining the Sparse Generalised Canonical Correlation Analysis (SGCCA) and the GraphNet penalty. We applied our model to the publicly available TCGA-LGG dataset. Contrary to Elastic-Net alone, GraphNet penalty can select a reasonable set of genes and yields informative biological interpretation from the pathway enrichment analysis. The example on the TCGA-LGG dataset exhibits the stability and reliability of netSGCCA for selecting variables of interest. However, it is essential to note that we show that the co-selection of variables is not primarily influenced by the structure of the graph but rather by its overall density. Therefore, an interpretation in terms of the pathways selected from the graph is elusive. Nevertheless, the method did extract genes that have been found co-(de)regulated in other studies of low-grade gliomas and other brain tumours.

Finally, we used netSGCCA to integrate radiomics and genetic data and applied it to the task of survival prediction. Due to the unavailability of survival data on our DIPG cohort, we used the TCGA-

LGG dataset to conduct the study. We compared the results obtained with netSGCCA with other multi-block survival approaches and the mono-block built models. Our study focused on the variables selected by the model and compared the results with the available literature. The netSGCCA showed to be a robust model able to select variables already linked to the studied pathology and interacting 114 Conclusion in relevant biological pathways. Adding the imaging features hindered the ability of the netSGCCA to extract features that can be used for survival prediction. However, the baseline results show that radiomic features extracted from the T2w modality can be strong predictors. Thus, larger datasets are required to further investigate the potential of radiomic features.

Future works

Studying rare tumours with statistical models is challenging. Data is scarce due to the nature of the tumour. Larger datasets are required to obtain conclusive and robust results. We also observe a high variability in quality and a high number of missing data, which is due to the fact that the data are collected in a clinical setting. The difficulty of adhering to a strict protocol designed for statistical analysis is obviously understood. It is important to carefully handle these datasets, as the available processing tools and models are not robust to noisy, high-variability data.

Our work proposed a segmentation approach that allowed us to obtain ROIs on the DIPG. However, further work should focus on finely delineating the tumours and their sub-components. Furthermore, the quality of segmentation and its impact on the radiomic features must be assessed in detail.

Missing data is inherent to multi-block studies, as collecting data from multiple sources in a clinical context is definitely hard to achieve. This thesis did not explore handling missing data in multi-block integration, which resulted in using a relatively small dataset and would deserve future work.

Finally, this work proposed using radiomic features in a multi-block framework. However, these features were biased towards the acquisition tools, reconstruction algorithm and preprocessing methods. No definitive conclusion can be drawn about image signatures of the tumours without strict control of all variability sources. Future works should exploit, for example, the already available follow-up MRI data when the embargo on the clinical annotations of the clinical study is raised.

✽ ✽ ✽ ✽ ✽ ✽

Table A.1: Selected mutated genes using the raw mutation matrix. Out of 10 bootstrap samples, we only kept the 6 where the model did not fail. Ensuite, les caractéristiques de statistiques de second ordre (texture) servent à décrire des configurations de niveau de gris fines et locales. L'hypothèse derrière l'analyse radiomique est que l'IRM peut capturer des caractéristiques de texture invisibles à l'oeil nu, qui peuvent être associées à l'évolution d'une maladie.

Nous faisons l'hypothèse que combiner les caractéristiques radiomiques avec les données de biologie moléculaire permettra de construire des modèles riches et très informatifs sur la biologie des tumeurs cérébrales. L'objectif est de trouver de nouveaux biomarqueurs pronostiques hybrides, voire de caractériser la tumeur à partir des seuls descripteurs de l'imagerie, ce qui fournirait à terme un moyen • La première méthode repose sur la suppression des segmentations en dehors de nos boundingboxes, justifié par le fait que la probabilité qu'un voxel tumoral se trouve en dehors de la bounding-boxes est faible.

• La deuxième méthode positionne les bounding-boxes en entrée des modèles de segmentation. Le résultat de la segmentation est ensuite masqué comme dans la première proposition.

Nous avons étudié les segmentations obtenues par nos deux approches sur 3 jeux de données;

• TCGA-GBM, une cohorte de patients présentant glioblastome. Cette cohorte comporte la même tumeur, localisée dans la même région du cerveau que le jeu de données utilisé lors de l'entraînement.

• TCGA-LGG, une cohorte de patients avec un gliome de bas grade. Il s'agit de segmenter une tumeur différente, localisée dans la même région du cerveau que le jeu de données utilisé lors de l'entraînement.

• PREBIOMEDE, une petite cohorte privée de jeunes patients avec un DIPG et pour lesquels une segmentation experte est disponible.

Les résultats de nos propositions sont compétitifs par rapport à ceux des méthodes de l'état de l'art Nous avons proposé d'étudier l'impact de la pénalité de graphe sur un modèle multi-blocs, dans le cadre de SGCCA (Sparse Generalised Canonical Correlation Analysis) [START_REF] Tenenhaus | Variable selection for generalized canonical correlation analysis[END_REF].

NetSGCCA consiste en une pénalité GraphNet ajoutée à un modèle SGCCA. Nous avons étudié l'effet de la pénalité sur le modèle à l'aide de l'ensemble de données TCGA-LGG.

Sur les données TCGA-LGG, en utilisant les données mRNA, miRNA et CNV, nous avons étudié l'effet de la normalisation du Laplacien du graphe dans la pénalité GraphNet. Ensuite, nous avons comparé différents graphes d'interaction gène-gène. Ces graphes sont différents dans leurs structures et informations contenues. Enfin, nous avons modifié un graphe en permutant ses noeuds et en supprimant quelques arêtes afin de vérifier l'utilité de l'information contenue dans le graphe.

Notre étude s'est focalisée, dans un premier temps, sur la sélection de variable avec la pénalité de graphe et sa stabilité. Dans un second temps, nous avons vérifié l'utilité des variables sélectionnées à prédire la survie et leur lien préalablement établi avec la tumeur étudiée dans la littérature.

Nous avons exhibé que cette pénalité permet de regrouper des variables, qui étaient séparées avec une pénalité l1 classique, et nous avons observé une meilleure stabilité quand le Laplacien associé au graphe est normalisé. En comparant différents graphes, nous avons montré une corrélation positive entre la densité du graphe (définie comme le nombre d'arêtes) et le nombre de variables sélectionnées.

Les résultats obtenus en utilisant les différents graphes ont été cohérents. Malgré un nombre de gènes sélectionnés différents avec les différents graphes, les ensembles de gènes sélectionnés étaient inclus les uns dans les autres. En permutant les noeuds du graphe, nous avons montré que la sémantique du graphe a un effet limité sur la sélection des variables, comparé à la densité du graphe.

Notre étude a aussi montré que les variables sélectionnées exhibent des voies biologiques précédemment associés aux gliomes de bas grade dans la littérature. Ceci démontre l'apport du graphe sur l'interprétabilité des résultats. De plus, la netSGCCA donne des résultats équivalents à la SGCCA en termes de prédiction de survie, voire dans certains cas meilleurs.

C.2.3 Chapitre 7 : Intégration Radio-génomique

Les caractéristiques radiomiques ont été utilisées avec succès pour différentes tâches et dans divers types de tumeurs. Cela comprend le DIPG, les glioblastomes et les gliomes de bas grade (LGG).

C.3 Conclusion

2. 1

 1 Acquired T2 signal. (a) As the spins realign with the magnetic field, their electrical signal decreases. The decrease is characterised by the T2 relaxation time. (b) The T2 relaxation time is tissue dependant. Acquiring the signal at TE gives the contrast between the different tissues and results in a T2-weighted image. . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.2 Acquired T1 signal. (a) As the spins realign with the magnetic field, their magnetism increases. The increase is characterised by the T1 relaxation time. (b) The T1 relaxation time is tissue dependant. Acquiring the signal at TR gives the contrast between the different tissues and results in a T1-weighted image. . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.3 A description of how the different texture features are computed. In an example 4×4 image ROI, three gray levels are represented by numerical values from 1 to 3. The red and blue annotations are used to highlight the different voxels used in order to compute some matrix values. To illustrate the different matrices: For the GLCM we only considered voxels with voxels in the right and left as neighbouring. Note that this matrix is symmetric, thus voxel couples are counted twice. For the GLRLM, we also used the same definition of the neighbourhood as for the GLCM. For all other matrices, the neighbourhood is given by the infinity norm. The radiomic features presented are only representative examples. . . . . . 13 4.1 Kaplan Meier plot on the 561 TCGA-GBM patients with survival available. Duration is given in days . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 4.2 Kaplan Meier plot on the 561 TCGA-LGG patients with survival available. Duration is given in days. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 4.3 The different pre-processing steps applied to MR images . . . . . . . . . . . . . . . . . . . 43 4.4 Intensity distribution on the MRI image for each patient on the LGG dataset. Columns correspond to the modalities T1w, T2w, FLAIR, CE T1w. First row represent distribution of the white matter voxels, and second row for the grey matter. (a) Before white-stripe and (b) after white-stripe normalisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 5.1 Top row: Example of T2w MRI scan of a patient with DIPG tumour extending beyond the pons. (a) Axial, (b) Coronal and (c) Sagittal slices. Bottom row: Example of glioblastoma and DIPG MRI scans, in axial slices. (d) glioblastoma T2w, (e) glioblastoma FLAIR, (f) DIPG T2w and (g) DIPG FLAIR. Tumours appear on the images as hyper-signal. . . . two approaches pYU and sYBBU . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 5.3 Experimental design diagram . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 5.4 Details of the 6 different models that are trained for each modality. These 6 models will be used for inference either independantly, or serialised parallelized. . . . . . . . . . . . . 58 5.5 The trained models are used for inference either independently, or serialised parallelized. Finally results may be ensembled accross the modality. . . . . . . . . . . . . . . . . . . . . 58 5.6 Mean precision-recall graphs of the different proposed segmentations on the TCGA-GBM dataset. To focus on the most interesting part of the plot, we only plotted the precisionrecall scores for thresholds between 0.1 and 0.9. From the left to the right, Using the FLAIR, using the T2w, and ens.(FLAIR,T2w). . . . . . . . . . . . . . . . . . . . . . . . . 61 5.7 Mean precision-recall graphs of the different proposed segmentations on the LGG dataset.

  T2w, and ens.(FLAIR, T2w) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 5.8 FLAIR image of a DIPG patient from the PREBIOMEDE cohort. Detection Failed on the image. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 5.9 Segmentation results obtained on one case of DIPG are superimposed on the FLAIR background. The top row displays the complete patient images. Ground truth mask in blue. Yolo detection contours in red. sYBBU Segmentation with FLAIR in orange. sYBBU Segmentation with T2w in green. pYU Segmentation with FLAIR in purple. pYU Segmentation with T2w in yellow. The tumours are presented as a hyper-signal. . . . . . . . 66 6.1 Star and Path graphs with different u 2 values. Grey nodes correspond to 0 values, red for 1, and blue for -1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 6.2 Illustration of a case were Σ 2 does not reflect a valid correlation matrix . . . . . . . . . . 77 6.3 Details of the proposed experimental design for the TCGA-LGG dataset. The analysis step refers to the comparison between the different graphs, analysis of the variable selection when nodes are permuted or edges are removed. Application refers to survival prediction and enrichment analysis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 6.4 Comparison between Raw and Normalised Laplacian. (a) Evolution of the number of selected genes as γ G varies, using the raw and normalised Laplacian. (b) Evolution of the stability metric as γ G varies, using the raw and normalised Laplacian. (c) Number of selection co-occurrence of selected genes, with γ G = 10 2 . . . . . . . . . . . . . . . . . . . . 83 6.5 Degree distribution of selected genes by fold for Raw and Normalised Laplacians. (a) For each selected gene, we counted the number of its neighbours in the PC graph. The black line represents the density of the degree distribution of all genes in the PC graph. (b) For each selected gene, we counted the number of its neighbours among selected genes. . . . . 83 6.6 Box plot of the weight distribution of selected genes for Raw and Normalised Laplacians. x-axis is log2 of the gene degrees in full PC graph. . . . . . . . . . . . . . . . . . . . . . . 84 6.7 Average absolute weight difference between selected nodes according the Dijkstra distance between them in the PC graph . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 distribution of selected genes. Between stable selected genes ( Intersection of selected genes in the 5 folds ), Candidate selected genes ( Union of selected genes in the 5 folds ) and All the genes in the dataset. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.9 Venn diagram showing the overlap between genes selected by the PC graph and the MSIGDB and the KEGG Graphs. Diagrams on the right shows the gene sets resulting from the candidate selected genes, while left diagrams show the gene sets from the stable selected genes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.10 Box plot for the IOU metric between the genes selected by the PC graph and the permuted PC graph . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.11 The evolution of the number of selected genes when the number of edges decreases (x-axis is the percentage of edges removed from 0% to 90%.) . . . . . . . . . . . . . . . . . . . . 7.1 Diagram of the different run models and the different comparisons discussed in this chapter. 7.2 Variable selection and model assessment methodology. The same CV folds and bootstrap samples are kept for all compared models. . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.3 Distribution of the C-index values according to the modality of the radiomics. Obtained from the ten bootstrap runs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.4 Distribution of the C-index values according to the pre-processing of somatic mutation profiles. Obtained from the ten bootstrap runs. . . . . . . . . . . . . . . . . . . . . . . . . 7.5 Number of times the selected genes appear to be mutated in the original somatic mutation profiles, using graph propagation (left) with the PC graph, and (right) with the influence graph. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.6 Distribution of the C-index values obtained from the ten bootstrap runs for RNA. . . . . 7.7 Distribution of the C-index values obtained from the ten bootstrap runs for RNA. . . . . 7.8 Number of radiomic features selected from each modality (left) at least in one run (right) in all three runs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.9 (a) Number of selected variables from each block, in each successful run using the Cox model with all blocks concatenated. (b) Total number of variables from each block using the Cox model with all blocks concatenated . . . . . . . . . . . . . . . . . . . . . . . . . . 7.10 Distribution of the C-index values obtained from the ten bootstrap runs, for the four multiblock models tested. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.11 Correlation matrices between the 17 most frequently selected genes from the RNA block by the netSGCCA model. (left panel): the minimum correlation between the genes in the ten bootstrap samples; (center panel): the maximum correlation between the genes in the ten bootstrap samples; (right panel): the difference between the two previous matrices. . .List of Tables4.1 TCGA-GBM Patient Demographics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 4.2 TCGA-LGG Patient Demographics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 4.3 Number of samples available of the different datasets used . . . . . . . . . . . . . . . . . . 42 5.1 YOLO parameters used for training . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 5.2 Dataset sizes of the different training and testing sets . . . . . . . . . . . . . . . . . . . . 57 5.3 Detection results on TCGA-GBM with 97 test patients. Results present the mean ± standard deviation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 5.4 Detection results on LGG with 76 test patients. Results present the mean ± standard deviation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 5.5 Segmentation results on TCGA-GBM with 97 test patients. Unlike (BB-)UNet approaches, Deepmedic network (Kamnitsas et al., 2017) is trained on 3D volumes from the HGG train dataset. Results present the mean ± standard deviation. . . . . . . . . . . . . . . . . . . . 60 5.6 Correlation study. Correlation values between detection precisions and final segmentation results obtained on the ensembled bounding-boxes. . . . . . . . . . . . . . . . . . . . . . . 61 5.7 Correlation study. Correlation values between detection recalls and final segmentation results obtained on the ensembled models. . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 5.8 Comparison of segmentation performances when using the real bounding-boxes and predicted bounding-boxes for the FLAIR without post-processing. Results present the mean ± standard deviation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 5.9 Segmentation results on LGG with 76 test patients. Results present the mean ± standard deviation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 5.10 Detection results on 71 test cases from the DIPG set. . . . . . . . . . . . . . . . . . . . . . 63 5.11 Detection results on 62 sessions from the DIPG dataset excluding the 9 sessions where detection step failed to detect anything. Results present the mean ± standard deviation . 64 5.12 Segmentation results on 62 test sessions from the DIPG set, using ens.(FLAIR, T2w) detection masks. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 5.13 Segmentation results on 62 sessions from DIPG, using FLAIR detection masks. Results present the mean ± standard deviation . . . . . . . . . . . . . . . . . . . . . . . . . . .

  (c = 2). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3 Different sources of prior knowledge graphs . . . . . . . . . . . . . . . . . . . . . . . . . .

6. 5

 5 Performances in survival prediction, number of selected variables and pathways depending on the type of graph to constrain the model. . . . . . . . . . . . . . . . . . . . . . . . . . 6.6 Top gene sets from MSigDB C6 collection. In bold, gene sets with an adjusted p-value lower than 0.05. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.1 Survival prediction results (C-index) obtained on the 95 patients from the TCGA-LGG dataset using radiomic features. Validation results were obtained on three-fold crossvalidation on 76 patients. Test results were obtained using 10 bootstraps on 19 patients. Elastic-Net Cox was used. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2 Selected radiomic features extracted from the T2 modality. Out of 10 bootstrap samples,

  Data presentation. DIPG or DMG are rare tumours with median overall survival of less than one year. This situation motivates clinical research projects aiming to acquire rich and multiple molecular measurement datasets obtained at different levels of biological organisation. These datasets could open new therapeutic opportunities and be used to derive different disease subtypes. In DIPG, tumours are inoperable because they are highly infiltrative and located in the pons. Consequently, the disease evolution and treatment impact on the tumour (radiotherapy and possibly targeted chemotherapy) is monitored using MR imaging assessment.Most of these clinical projects are pharmaceutical -they intend to study promising protocols -and include ancillary data for further study. Two clinical trials (PHRC BIOMEDE 1 and PHRC BIOMEDE 2.0 2 ) have been recently designed and promoted by Gustave Roussy Cancer Campus. They aim to inform treatment decisions according to molecular biomarkers and test targeted chemotherapies assigned based on patients' individual molecular profiles. This is to improve disease control along with conventional radiotherapy. The BIOMEDE pharmacological trials also include an ancillary study that brought additional data consisting of multi-omics from the biopsy, imaging acquired longitudinally, image biomarkers resulting from image processing, and clinical data. Even if the data of the ancillary study are not yet wholly available, they foreshadow the context of this type of study with multi-block data. In addition to their multi-block presentation, they include a significant amount of missing data.

  Figure 2.1: Acquired T2 signal. (a) As the spins realign with the magnetic field, their electrical signal decreases. The decrease is characterised by the T2 relaxation time. (b) The T2 relaxation time is tissue dependant. Acquiring the signal at TE gives the contrast between the different tissues and results in a T2-weighted image.

Figure 2

 2 Figure 2.3: A description of how the different texture features are computed. In an example 4×4 image ROI, three gray levels are represented by numerical values from 1 to 3. The red and blue annotations are used to highlight the different voxels used in order to compute some matrix values. To illustrate the different matrices:For the GLCM we only considered voxels with voxels in the right and left as neighbouring. Note that this matrix is symmetric, thus voxel couples are counted twice. For the GLRLM, we also used the same definition of the neighbourhood as for the GLCM. For all other matrices, the neighbourhood is given by the infinity norm. The radiomic features presented are only representative examples.

  [START_REF] Kickingereder | Radiomic profiling of glioblastoma: Identifying an imaging predictor of patient survival with improved performance over established clinical and radiologic risk models[END_REF] studied the performance of radiomics extracted from MR images when used to predict survival in 119 glioblastoma patients. The authors report a c-index of 0.65 when predicting overall survival with radiomics alone and 0.70 when clinical data were added to the radiomics.Comparatively, clinical data reached a c-index of 0.64. In this study, all variables of interest were texture features extracted from the FLAIR modality.[START_REF] Kickingereder | Radiomic subtyping improves disease stratification beyond key molecular, clinical, and standard imaging characteristics in patients with glioblastoma[END_REF] analysed the added value of radiomics (extracted from MR images)

  [START_REF] Havaei | Brain tumor segmentation with deep neural networks[END_REF] proposed a two-path 2D CNN-based model to classify the central pixel of each patch. Their approach only uses the axial view as it has the best resolution. The model combines the results of two simultaneously trained CNNs, each using a different patch size. The chosen sizes were 33x33 and 65x65. The output of the first CNN is concatenated with the second CNN input, thus treated as additional channels of the input patch. The local path uses 7x7 convolution kernels, while the global path uses 13x13 convolution kernels. Finally, the concatenation of the feature map is fed to the last output layer of the model to make a prediction. The authors achieved second rank in the 2015 BraTS challenge with an average reported dice of 0.88 on the whole and 0.79 on the core tumour.[START_REF] Kamnitsas | Deepmedic for brain tumor segmentation[END_REF] proposed DeepMedic an 11-layers deep, multi-scale 3D CNN. The model is built using two parallel paths, each processing the input at a different scale. The low-resolution and high-resolution paths take as input 3D patches of sizes 19x19x19 and 25x25x25, respectively. The network segments a small patch of 9x9x9 at the centre of the input data. The authors also propose to extend the original model with residual connections (adding the results of earlier layers as input to some later layers) to help train the network by facilitating signal preservation. The model was one of the top-performing networks during BraTS 2016, with a reported dice of 0.89 on the whole tumour and 0.76 on the tumour core. The addition of the residual connections had a small improvement in the results.Similarly toHavaei's model, Ben naceur et al. (2020) proposed CNN architectures designed to predict a central pixel class using the 2D patch from the image. They propose to use overlapping patches to help the model utilise local and global information. Thus, they used five predictions per patch, where the prediction of the adjacent patches influences the results obtained on the current patch. The authors report a dice score of 0.90 on the whole tumour and 0.83 on the core tumour when tested on the BraTS 2018 data.

  symmetric Encoder-Decoder CNN. The model has skip connections between the encoding part and the decoding part. These skip connections directly pass features from a latent encoding layer to its symmetric decoding layer, which helps the model recover details.[START_REF] Casamitjana | 3d convolutional neural networks for brain tumor segmentation: A comparison of multi-resolution architectures[END_REF] used the UNet architecture in the BraTS challenge 2015 and reported a dice score of 0.89 on the whole tumour and 0.76 on the core tumour.[START_REF] Isensee | Brain tumor segmentation and radiomics survival prediction: Contribution to the brats 2017 challenge[END_REF] proposed a modified UNet to perform brain tumour segmentation. To do so, the authors used feature maps from the decoding path to generate secondary segmentation maps.These are combined with the final segmentation, and the losses associated with these segmentation maps are added to the final loss function. The authors also used a dice loss function to tackle class imbalance and performed data augmentation. The model was one of the top-ranked during the BraTS 2017 challenge, with a dice score of 0.90 and 0.80 on the whole and core tumour, respectively.[START_REF] Milletari | V-net: Fully convolutional neural networks for volumetric medical image segmentation[END_REF] proposed V-net, a modified UNet. V-net defines the data flow as a chain of multiple stages operating on different resolutions, each stage comprising one to three convolution layers. V-net modifies UNet by adding the input of each stage in the encoding path to the output of the last convolutional layer of that stage in order to enable learning a residual function.Hua et al.

  [START_REF] Zhou | Learning contextual and attentive information for brain tumor segmentation[END_REF] used an architecture similar to UNet. The authors proposed to train their model on three tasks simultaneously: the binary segmentation of the whole tumour, the detailed 3.3. Multi-block approaches 29 segmentation of the tumour, and the binary segmentation of the enhancing tumour. When tested on the BraTS 2018 challenge, the authors report an average Dice of 0.88 on the whole tumour and 0.80 in the tumour core.[START_REF] Myronenko | 3d MRI brain tumor segmentation using autoencoder regularization[END_REF] proposed NV-net, an encoder-decoder architecture with an additional branch to the decoder part to reconstruct the original image. When tested on the BraTS 2018 challenge, the authors report an average Dice of 0.88 on the whole tumour and 0.82 in the tumour core.

  Rosana et al. (2020) proposed Bounding Box UNet (BB-UNet), a neural network model that expends the UNet model to consider prior bounding boxes. The proposed model incorporates priors through novel convolutional layers introduced at the level of skip connections to guide the model on

Figure 4 . 1 :

 41 Figure 4.1: Kaplan Meier plot on the 561 TCGA-GBM patients with survival available. Duration is given in days

Figure 4 . 2 :

 42 Figure 4.2: Kaplan Meier plot on the 561 TCGA-LGG patients with survival available. Duration is given in days.

  Figure 4.3 illustrates the preprocessing steps.

  Figure 4.4: Intensity distribution on the MRI image for each patient on the LGG dataset. Columns correspond to the modalities T1w, T2w, FLAIR, CE T1w. First row represent distribution of the white matter voxels, and second row for the grey matter. (a) Before white-stripe and (b) after white-stripe normalisation

Figure 5 . 1 :

 51 Figure 5.1: Top row: Example of T2w MRI scan of a patient with DIPG tumour extending beyond the pons. (a) Axial, (b) Coronal and (c) Sagittal slices. Bottom row: Example of glioblastoma and DIPG MRI scans, in axial slices. (d) glioblastoma T2w, (e) glioblastoma FLAIR, (f) DIPG T2w and (g) DIPG FLAIR. Tumours appear on the images as hyper-signal.

  Figure 5.2: The two approaches pYU and sYBBU

Figure 5 . 3 :

 53 Figure 5.3: Experimental design diagram

Figure 5 . 5 :

 55 Figure 5.5: The trained models are used for inference either independently, or serialised parallelized. Finally results may be ensembled accross the modality.

Figure 5 . 6 :

 56 Figure5.6: Mean precision-recall graphs of the different proposed segmentations on the TCGA-GBM dataset. To focus on the most interesting part of the plot, we only plotted the precision-recall scores for thresholds between 0.1 and 0.9. From the left to the right, Using the FLAIR, using the T2w, and ens.(FLAIR,T2w).

Concerning

  Deepmedic architecture, results on the FLAIR modality are slightly below results obtained with the proposed approaches. The Deepmedic model seems to prioritise high precision over recall. However, Deepmedic trained with the T2w failed to give any meaningful result, with an average Dice index of 0.08, which makes the T2w unusable in an ensembled model.

Figure 5 . 7 :

 57 Figure 5.7: Mean precision-recall graphs of the different proposed segmentations on the LGG dataset.To focus on the most interesting region, we only plotted the precision-recall scores for thresholds between 0.1 and 0.9. From the left to the right, using the FLAIR, using the T2w, and ens.(FLAIR, T2w)

  cord, thalamus). Careful inspection of the 9 cases with detection step failure indicates that the failure is mostly related to the tumours not being visible on the FLAIR and T2w MRI scans. Due to the low performance of object detection on T2w, we only used FLAIR and ens.(FLAIR, T2w) detection masks in the segmentation phase.

Figure 5

 5 Figure 5.8: FLAIR image of a DIPG patient from the PREBIOMEDE cohort. Detection Failed on the image.

Figure 5 . 9 :

 59 Figure 5.9: Segmentation results obtained on one case of DIPG are superimposed on the FLAIR background. The top row displays the complete patient images. Ground truth mask in blue. Yolo detection contours in red. sYBBU Segmentation with FLAIR in orange. sYBBU Segmentation with T2w in green. pYU Segmentation with FLAIR in purple. pYU Segmentation with T2w in yellow. The tumours are presented as a hyper-signal.

Figure 6

 6 Figure 6.1: Star and Path graphs with different u 2 values. Grey nodes correspond to 0 values, red for 1, and blue for -1

  ± 0.06 0.35 ± 0.08 Complete 10 -4 0.36 ± 0.09 0.02 ± 0.07 0.06 ± 0.23 0.03 ± 0.11 0.39 ± 0.09 0.05 ± 0.22 0.02 ± 0.09 0.03 ± 0

Figure 6

 6 Figure 6.3: Details of the proposed experimental design for the TCGA-LGG dataset. The analysis step refers to the comparison between the different graphs, analysis of the variable selection when nodes are permuted or edges are removed. Application refers to survival prediction and enrichment analysis.

  Figure 6.4: Comparison between Raw and Normalised Laplacian. (a) Evolution of the number of selected genes as γ G varies, using the raw and normalised Laplacian. (b) Evolution of the stability metric as γ G varies, using the raw and normalised Laplacian. (c) Number of selection co-occurrence of selected genes, with γ G = 10 2 .

  Figure 6.5: Degree distribution of selected genes by fold for Raw and Normalised Laplacians. (a) For each selected gene, we counted the number of its neighbours in the PC graph. The black line represents the density of the degree distribution of all genes in the PC graph. (b) For each selected gene, we counted the number of its neighbours among selected genes.

Figure 6 . 6 :

 66 Figure 6.6: Box plot of the weight distribution of selected genes for Raw and Normalised Laplacians. x-axis is log2 of the gene degrees in full PC graph.

Figure 6 . 9 :

 69 Figure 6.9: Venn diagram showing the overlap between genes selected by the PC graph and the MSIGDB and the KEGG Graphs. Diagrams on the right shows the gene sets resulting from the candidate selected genes, while left diagrams show the gene sets from the stable selected genes

Figure 6 .

 6 Figure 6.10: Box plot for the IOU metric between the genes selected by the PC graph and the permuted PC graph

Figure 7 . 1 :

 71 Figure 7.1: Diagram of the different run models and the different comparisons discussed in this chapter.

Figure 7 . 2 :

 72 Figure 7.2: Variable selection and model assessment methodology. The same CV folds and bootstrap samples are kept for all compared models.

Figure 7 . 3 :

 73 Figure 7.3: Distribution of the C-index values according to the modality of the radiomics. Obtained from the ten bootstrap runs

  Figure 7.4: Distribution of the C-index values according to the pre-processing of somatic mutation profiles. Obtained from the ten bootstrap runs.

Figure 7 . 7 :

 77 Figure 7.7: Distribution of the C-index values obtained from the ten bootstrap runs for RNA.

  totalement non-invasif de monitorer les DIPG ou les tumeurs cérébrales en général. Cette thèse est motivée par l'étude du DIPG inclue dans les études cliniques BIOMEDE dirigées par l'Institut Gustave Roussy. Ces gliomes pédiatriques sont rares (environ 40 nouveaux cas chaque année en France et 300 en Europe). Une fraction importante des données initialement prévue pour cette étude n'a pas été accessible. Nous avons décidé d'étendre les jeux de données étudiés aux gliomes de bas grade (LGG) et au glioblastome, pour lesquels des données multi-omiques sont publiquement disponibles grâce au projet The Cancer Genome Atlas (TCGA). Bien que ces tumeurs soient d'origine et de pronostic différents, nous utilisons ces données afin d'étudier les méthodes proposées et leurs limites. Ce travail propose l'intégration des données d'imagerie avec les données génétiques afin de trouver des biomarqueurs. Dans un premier temps, nous nous intéressons à l'extraction des régions d'intérêt des images nécessaires pour une étude radiomique. Ensuite, nous proposons une procédure d'intégration des données multi-sources, qui prend en compte les graphes complexes d'interaction entre les gènes.Finalement, nous appliquons notre procédure sur les données disponibles afin de comparer ses performances avec d'autres modèles de la littérature et d'étudier l'apport de l'imagerie et du graphe aux données génétiques.Le manuscrit de thèse est organisé en deux parties. La première partie est dédiée à la définition du contexte de l'étude. Elle débute par les définitions et notations utilisés dans toute la suite de la thèse (Chapitre 2). Ensuite, nous révisons les méthodes liées à notre étude (Chapitre 3). Enfin nous décrivons les données disponibles et utilisées durant la thèse avec les étapes de prétraitements utilisés (Chapitre 4). La deuxième partie de la thèse décrit notre contribution. Nous débutons par notre proposition pour la segmentation automatique de la tumeur (Chapitre 5). Ensuite, nous détaillons la méthode d'analyse multi-sources avec contrainte de graphe (Chapitre 6). Et enfin, nous appliquons nos méthodes sur les données disponibles (Chapitre 7). Dans la suite, nous résumons les trois chapitres de notre contribution.La segmentation d'images est une étape clé pour l'analyse radiomique. L'identification de la région d'intérêt (ROI) de l'IRM est critique pour extraire des descripteurs significatifs de la tumeur en question. Étant donné un volume (Image) d'entrée d'une ou plusieurs modalités, la segmentation automatique des tumeurs cérébrales fait référence à la différenciation entre les voxels appartenant à la tumeur (et sa structure) et ceux appartenant au tissu cérébral sain. Pour cela, on utilise des modèles statistiques de classification appliqués à chaque voxel de l'image.BraTS (Brain Tumour Segmentation challenge) est un défi international qui reflète l'avancement des méthodes de segmentation en oncologie. Depuis 2015, presque tous les modèles les plus performants du défi BraTS ont utilisé des réseaux de neurones convolutifs (Convolutionnal Neural Networks (CNN)). Notamment, l'architecture U-Net qui permet simultanément de reconstruire en sortie une segmentation de l'image d'origine. Nous proposons d'exploiter l'architecture U-Net[START_REF] Ronneberger | U-Net: Convolutional networks for biomedical image segmentation[END_REF] en nous appuyant sur la détection d'objets.Le défi consiste à obtenir des segmentations du DIPG sans entraîner des modèles sur cette tumeur car il n'y a pas assez de données disponibles, mais en les entraînant sur d'autres tumeurs pour lesquelles suffisamment de données sont publiquement accessibles. Nous avons décomposé le problème de délinéation des tumeurs en deux sous problèmes, i) de détection de tumeur dans des boîtes parallèlipidédiques (ou bounding-boxes) et ii) de segmentation des tumeurs dans ces boites. Nous proposons deux méthodes pour combiner un modèle de détection, YOU ONLY LOOK ONCE (YOLO)[START_REF] Redmon | You only look once: Unified, real-time object detection[END_REF], et le modèle UNET/BBUNET:

  sur TCGA-GBM et meilleurs sur BRaTS LGG et PREBIOMEDE. De plus, notre méthode permet de proposer des segmentations robustes sur le DIPG ce qui évite d'avoir recours à un expert.C.2.2 Chapitre 6 : Intégration multi-blocksDans le domaine de l'oncologie, des modèles statistiques sont utilisés pour la découverte de facteurs candidats qui influencent le développement de la pathologie ou sa réponse au traitement. Ceci est souvent réalisé en imposant des contraintes sur les modèles qui conduisent à une sélection des variables d'intérêt. Les contraintes basées sur des graphes a priori ont été utilisées dans la littérature comme un moyen d'améliorer la sélection des variables dans le modèle, ce qui fournit des modèles plus interprétables d'un point de vue biologique et/ou clinique. Cependant, les interactions entre les graphes choisis et le modèle sont mal caractérisées et on ignore comment ils impactent la sélection des variables. De plus, comme plusieurs graphes codant des informations différentes sont disponibles, on peut se demander comment le choix du graphe impacte les résultats obtenus.

  Cependant, il n'y a pas de consensus sur la modalité d'image qui est la plus informative ou le type de caractéristique le plus important. De plus, la recherche a également montré que l'ajout de caractéristiques radiomiques aux données moléculaires peut améliorer la prédiction de la survie. Mais ilC.3. Conclusion 125doit encore être intégré dans un cadre multibloc complet. Cela devrait nous aider à comprendre les interactions entre les caractéristiques radiomiques et génomiques.Ce chapitre est consacré à l'application de netSGCCA pour la prédiction de survie lorsque des données génétiques sont associées à des caractéristiques radiomiques. Nous portons une attention particulière à l'ensemble des variables sélectionnées et à leur lien avec la maladie étudiée. Étant donné que le cadre multi-blocs avec lequel nous travaillons nécessite la disponibilité des données pour toutes les sources pour tous les patients, il en résulte naturellement moins d'échantillons disponibles pour une telle analyse. Par conséquent, nous établissons les performances de base de chaque bloc lorsque toutes les données disponibles sont prises en compte avant l'intégration. Ensuite, nous comparons l'approche mono-bloc avec l'approche multi-blocs. Nous comparons également la netSGCCA avec les modèles les plus performants parmi tous les travaux de l'état de l'art. En raison de l'indisponibilité des données de survie sur l'ensemble de données DIPG, cette étude a été limitée à l'ensemble de données TCGA-LGG.La netSGCCA figurait parmi les modèles les plus performants de l'ensemble de validation et de test. Le C-index était stable pendant tous les essais, ce qui n'était pas le cas pour les autres modèles de Cox. Le modèle netSGCCA a également été capable de sélectionner des gènes et des voies déjà associés au développement des gliomes et au pronostic des LGG.Les performances de base démontrent que les caractéristiques radiomiques peuvent être de bons prédicteurs de la survie du gliome de bas grade (LGG). La modalité T2w montre des performances particulièrement élevées, comme déjà observé dans la littérature. Malheureusement, nous n'avons pas pu reproduire ces résultats lors de l'utilisation de la radiomique sur les 83 patients ayant tous les blocs disponibles.Dans l'ensemble, la restriction de la disponibilité des données a entravé les performances prédictives des blocs. Le modèle de Cox s'est montré également moins stable lorsque les données de seulement 83 patients ont été utilisées. De plus, la sélection de variables et les voies biologique identifiées dans l'étude de base et sur l'ensemble restreint ont donné des résultats significativement différents. Par conséquent, d'autres études sont nécessaires pour valider les résultats. Il est également nécessaire de réaliser cette étude sur un jeu de données plus important.

L

  'étude des tumeurs rares avec des modèles statistiques est difficile. Les données sont rares en raison de la nature de la tumeur. Cette rareté a deux conséquences. Des résultats peinent à être concluants et robustes et réclament des ensembles de données plus importants. La deuxième conséquence se traduit par des données collectées dans un contexte clinique, avec une grande variabilité de qualité et un nombre élevé de données manquantes. Tout en ayant conscience des difficultés que cela posent, un protocole plus strict conçu pour l'analyse statistique serait nécessaire. Notre travail a proposé une approche de segmentation qui nous a permis d'obtenir des ROIs sur le DIPG. La poursuite de ces travaux pourrait porter sur une délimitation fine des tumeurs et de leurs sous-compartiments. De plus, la qualité de la segmentation et son impact sur les caractéristiques radiomiques doivent être évalués en détail.

  graph. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.10 Survival prediction results obtained on the 419 patients from the TCGA-LGG dataset using RNA features. Validation results were obtained on three-fold cross-validation. Test results were obtained using 10 bootstraps. Elastic-Net Cox was used. . . . . . . . . . . . . . . . 7.11 Sample size comparison between the baseline setting and the current setting devoted to monomulti-block investigation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.12 Mono-block results obtained on the 83 subjects having the radiomics, somatic mutations and RNA data available. For the Radiomics, the concatenation of all the radiomics extrac-

ted from the four modalities has been used; The somatic mutations with graph propagation with the influence graph was used; Finally, the concatenation of all the blocks was also investigated . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.13 Selected radiomic features. The table shows features selected in convergent Cox models.

The means and standard deviation were computed using only the samples where the feature was selected. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.14 Selected genes using the somatic mutation profiles. The table shows genes that were selected in all 10 runs. The means and std were computed using only the samples where the feature was selected. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.15 Enrichment results when using the ten most selected genes using the influence graph. . . 7.16 Enrichment results when using the matrix with concatenated variables from all data sources.108 7.17 Multi-block results obtained on 83 subjects having the radiomics, somatic and RNA data available. For the radiomics, the concatenation of all the radiomics extracted from the four modalities has been used; The somatic mutations with graph propagation with the influence graph was used. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.18 Selected features using the netSGCCA, the normalised Pathway Commons Laplacian was

  Table only shows genes that were selected on all the runs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

A.2 Selected genes using the RNA data.

  Table shows genes that were selected in all runs. . . .

  The advances are not limited to hardware. Developments in software have led to increased imaging speeds, which reduced motion artefacts and screening costs, thus making MRI screening more available, accessible and affordable. Further details about MR images will be discussed in the Chapter In the case of WHO CNS subtype 3 diffuse high-grade gliomas, patients now undergo MRI exams for diagnosis but also for the follow-up of the answer to treatment. Recently[START_REF] Hugo | Decoding tumour phenotype by noninvasive imaging using a quantitative radiomics approach[END_REF] proposed to systematically extract textures of the computed tomography (CT) image in the different compartments of lung and head-and-neck tumours to study their potential prognostic power. This approach, which was generalised to other imaging types, including MRIs and was coined radiomics, allows the systematic exploration of the tumour macro/micro organisation

	2.
	DIPG and radiomics.

depending on the imaging sequence and resolution at hand. Practically, radiomic analysis includes all the extraction and analysis steps of a large number of features (200+) from a region of interest (ROI) on the available MRIs. These features include size and shape characteristics which describe apparent visual properties, then first-order statistics of the signal from the ROIs. Additionally, radiomic features include second-order statistics (texture) features that describe fine, local grey-level configurations. The hypothesis behind radiomic analysis is that MRI can capture texture features invisible to the naked eye, which can be linked to a disease outcome.

DIPG and high-throughput genomics. When a biopsy can be performed, the diagnosis of DIPG includes the use of advanced molecular immunohistochemistry tools. But for research purpose, highthroughput genomics is now often considered. Indeed in parallel to imaging progress, recent years have also known breakthroughs in molecular data acquisition technologies such as DNA/RNA Sequencing and Single Cell Sequencing, which increased the accessibility of these technologies. Coupled with the decrease in data storage costs and the rise in computational capabilities, these have made high-throughput multi-modal databases for hundreds of patients publicly available to study clinically relevant problems, especially in oncology.

DIPG and multi-view measurements. The recent evolution of data available in clinical research

for DIPG result in composite data obtained with various measurement tools from the same observations. It includes, but is not limited to, multi-omics data such as gene expression profiles, mutation profiles, copy number variants, clinical data from standardised electronic Case Report Forms, and imaging data. These data are regrouped into separate views or blocks, each block including multiple variables. Analysing these datasets required the development of new multi-view machine-learning and statistical models designed to integrate data

  Brain tumour segmentation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.2.1 Using custom features . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.2.2 Using Deep-learning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.2.3 Tackling issues with brain tumour segmentation . . . . . . . . . . . . . . . 28

3.1 Radiomics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.1.1 Classification of tumour types or subtypes using radiomics . . . . . . . . . . 23 3.1.2 Survival Prediction using radiomics . . . . . . . . . . . . . . . . . . . . . . . 23 3.1.3 Reproducibility of radiomics . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3.2 3.3 Multi-block approaches . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 3.3.1 Unsupervised models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 3.3.2 Supervised models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 3.4 Variable selection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 3.4.1 Sparsity-based penalties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 3.4.2 Graph based penalties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.1. Radiomics 23

4. Variable selection 35 hisa

  , 2000) pathways and correspond to dense sub-networks in the prior gene interaction network.

	✽ ✽ ✽
	✽ ✽ ✽	Chapter 4
	Datasets Description and Preprocessing

/KEGG (Kane-3.DIPG pediatric gliomas are rare (about 40 new cases each year in France and 300 in Europe) which make clinical research study very difficult. Unfortunately, they are resistant to classical therapy and only transiently sensitive to radiotherapy. Recent efforts from the DIPG community established that 90% DIPG/DMG share a common molecular anomaly, namely a point mutation on Lysine 27 of the histone H3 regulatory tail (H3K27M). In the remaining cases, it was recently demonstrated that tumour cells overexpress EZHIP

38 Datasets Description and Preprocessing Chapter Outline Contents 4

  .1 Datasets Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 4.1.1 Glioblastoma . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 4.1.2 Lower Grade Gliomas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 4.1.3 DIPG . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 4.1.4 Data availability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 4.2 Preprocessing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 4.2.1 Image preprocessing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 4.2.2 Radiomic extraction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 4.2.3 Mutation data preprocessing . . . . . . . . . . . . . . . . . . . . . . . . . . 46

Table 4

 4

	.1: TCGA-GBM Patient Demographics
	Sex	n (% of total)
	Male	366 (61)
	Female	230 (38)
	Total	602
	Age	Average (range), in years
		58.36 (11 -89)
	Survival	Average (std), in days
	Days to death	476.28 (529.01)
	Days to last follow-up	415.34 (482.37)
	Non-censored/Censored	447/149
	Table 4.2: TCGA-LGG Patient Demographics
	Sex	n (% of total)
	Male	285 (55)
	Female	230 (45)
	Total	516
	Age	Average (range), in years
		43.46 (14 -87)
	Survival	Average (std), in days
	Days to death	1205.80 (1121.42)
	Days to last follow-up	604.62 (846.96)
	Non-censored/Censored	92/422

Table 4

 4

	.3: Number of samples available of the different datasets used
		TCGA-GBM TCGA-LGG BIOMEDE
	# Images	102	108	89
	# Clinical	602	516
	# CNV	599	515
	# miRNA	5	514
	# RNAseq	160	512	69
	# Mutations	371	510	87
	# CNV & miRNA & RNAseq	5	419
	# RNAseq & Mutations & Images	20	90	28

Table 5

 5

	.1: YOLO parameters used for training
			Parameters
		lr0	0.00320
	Gradient	lrf	0.12000
		momentum	0.84300
		weight_decay	0.00036
	Training	warmup_epochs warmup_momentum	2.00000 0.50000
		warmup_bias_lr	0.05000
		box	0.02960
		cls	0.24300
		cls_pw	0.63100
	Loss function	obj	0.30100
		obj_pw	0.91100
		iou_t	0.20000
		anchor_t	2.91000
		degrees	0.37300
		translate	0.24500
		scale	0.89800
	Augmentation	shear	0.60200
		perspective	0.00000
		flipud	0.00856
		fliplr	0.50000

Table 5

 5

	.2: Dataset sizes of the different training and testing sets
	dataset	HGG train HGG val TCGA-GBM LGG	DIPG
	number of patients	142	15	97	76	30 (71 sessions)
	number of 2D images	9533	783

Table 5

 5

	.3: Detection results on TCGA-GBM with 97 test patients. Results present the mean ± standard
	deviation
				With morphological transformation
		Precision	Recall	Precision	Recall
	FLAIR	.577(±.104) .899(±.112) .599(±.105)	.927(±.075)
	T2w	.569(±.110) .880(±.152) .593(±.102)	.905(±.115)
	ens.(FLAIR,T2w) .511(±.103) .945(±.069) .527(±.105)	.956(±.059)

Table 5

 5

				With morphological transformation
		Precision	Recall	Precision	Recall
	FLAIR	.577(±.125) .849(±.203) .611(±.124)	.873(±.183)
	T2w	.581(±.122) .856(±.158) .610(±.135)	.883(±.136)
	ens.(FLAIR,T2w) .503(±.121) .926(±.121) .529

.4: Detection results on LGG with 76 test patients. Results present the mean ± standard deviation.

Table 5 .

 5 5: Segmentation results on TCGA-GBM with 97 test patients. Unlike (BB-)UNet approaches, Deepmedic network[START_REF] Kamnitsas | Efficient multi-scale 3d CNN with fully connected CRF for accurate brain lesion segmentation[END_REF] is trained on 3D volumes from the HGG train dataset. Results present the mean ± standard deviation.

	Architecture		Precision	Without masking Recall	Dice		Precision	With masking Recall	Dice
	UNet	FLAIR T2w ens.(FLAIR,T2w) .784 ± .260 .843 ± .103 .781 ± .192 .746 ± .295 .825 ± .100 .741 ± .215 .696 ± .275 .823 ± .163 .699 ± .230	pYU*	.902 ± .128 .813 ± .103 .845 ± .089 .858 ± .140 .812 ± .168 .813 ± .134 .914 ± .115 .830 ± .107 .861 ± .088
	BB-UNet*	FLAIR T2w ens.(FLAIR,T2w) .909 ± .114 .834 ± .094 .863 ± .087 .906 ± .120 .809 ± .091 .847 ± .087 .887 ± .090 .807 ± .127 .838 ± .096	sYBBU*	.921 ± .102 .807 ± .093 .854 ± .079 .901 ± .083 .806 ± .129 .843 ± .094 .925 ± .096 .835 ± .096 .869 ± .079
	DeepMedic	FLAIR	.913 ± .110 .774 ± .175 .820 ± .149

* Models using YOLO bounding-boxes.

Table 5 .

 5

		Precision Recall	Dice
	UNet	0.400	-0.083 0.334
	pYU	0.665	-0.073 0.422
	BB-UNet* 0.615	-0.037 0.363
	sYBBU	0.682	-0.033 0.357
	Table 5.7: Correlation study. Correlation values between detection recalls and final segmentation results
	obtained on the ensembled models.
		Precision Recall Dice
	UNet	0.167	0.665 0.319
	pYU	0.000	0.762 0.615
	BB-UNet* 0.189	0.791 0.684
	sYBBU	0.081	0.804 0.716

6: Correlation study. Correlation values between detection precisions and final segmentation results obtained on the ensembled bounding-boxes.

Table 5 .

 5 8: Comparison of segmentation performances when using the real bounding-boxes and predicted bounding-boxes for the FLAIR without post-processing. Results present the mean ± standard deviation.

	Precision	Recall	Dice
	YOLO Bounding-Boxes .906 ± .120 .809 ± .091 .847 ± .087
	Real Bounding-Boxes

.932 ± .072 .875 ± .073 .899 ± .043

Table 5

 5

	.9: Segmentation results on LGG with 76 test patients. Results present the mean ± standard
	deviation
	Architecture		Precision	Without masking Recall	Dice		Precision	With masking Recall	Dice
	UNet	FLAIR T2w ens.(FLAIR,T2w) .345 ± .275 .949 ± .062 .444 ± .303 .541 ± .360 .845 ± .139 .577 ± .296 .467 ± .305 .897 ± .172 .523 ± .295	pYU*	.863 ± .155 .768 ± .172 .792 ± .142 .785 ± .195 .831 ± .195 .766 ± .174 .871 ± .157 .800 ± .160 .814 ± .135
	BB-UNet*	FLAIR T2w ens.(FLAIR,T2w) .835 ± .161 .822 ± 0182 .815 ± .153 .844 ± .165 .773 ± .199 .790 ± .169 .828 ± .144 .797 ± .151 .801 ± .150	sYBBU*	.878 ± .144 .772 ± .201 .804 ± .165 .861 ± .132 .796 ± .181 .815 ± .148 .871 ± .141 .820 ± .183 .831 ± .152
	DeepMedic	FLAIR	.904 ± .146 .695 ± .254 .743 ± .215

* Models using YOLO bounding-boxes.

Table 5 .

 5 10: Detection results on 71 test cases from the DIPG set.

		Precision	Recall
	FLAIR	0.529 ± 0.265 0.580 ± 0.313
	T2w	0.436 ± 0.299 0.342 ± 0.341
	ens.(FLAIR, T2w) 0.395 ± 0.229 0.667 ± 0.305

Table 5 .

 5 11: Detection results on 62 sessions from the DIPG dataset excluding the 9 sessions where detection step failed to detect anything. Results present the mean ± standard deviation

			With morphological transformation
			Precision	Recall
		FLAIR	.606 ± .185	.664 ± .237
		T2w	.471 ± .285	.363 ± .334
		ens.(FLAIR,T2w) .447 ± .189	.734 ± .203
	Table 5.12: Segmentation results on 62 test sessions from the DIPG set, using ens.(FLAIR, T2w) detection
	masks.
			Precision	Recall	Dice
		FLAIR	0.73 ± 0.20 0.59 ± 0.24 0.62 ± 0.21
	pYU	T2w	0.66 ± 0.19 0.63 ± 0.24 0.61 ± 0.19
		ens.(FLAIR, T2w) 0.67 ± 0.20 0.63 ± 0.24 0.61 ± 0.20
		FLAIR	.725 ± .198 .609 ± .205 .627 ± .204
	sYBBU	T2w	.677 ±

Table 5 .

 5 13: Segmentation results on 62 sessions from DIPG, using FLAIR detection masks. Results present the mean ± standard deviation

	Architecture		Precision	Without masking Recall	Dice		Precision	With masking Recall	Dice
	UNet	FLAIR T2w ens.(FLAIR,T2w) .164 ± .172 .703 ± .234 .223 ± .173 .322 ± .235 .644 ± .225 .370 ± .210 .051 ± .048 .749 ± .232 .091 ± .077	pYU*	.753 ± .206 .571 ± .243 .611 ± .215 .680 ± .192 .616 ± .244 .606 ± .202 .695 ± .205 .622 ± .245 .609 ± .212
	BB-UNet*	FLAIR	.517 ± .220 .699 ± .236 .555 ± .195	sYBBU*

: Source code is freely available at https://github.com/neurospin/netSGCCA Chapter Outline Contents 6

  .1 Regularized Generalized Canonical Correlation Analysis . . . . . . . . . . . . . . . 71 6.2 Sparse Generalized Canonical Correlation Analysis . . . . . . . . . . . . . . . . . . 72 6.3 netSGCCA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 6.3.1 GraphNet . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 6.3.2 Optimisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 6.3.3 The proximal operator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 6.4 Application . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 6.4.1 Simulated Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 6.4.2 TCGA-LGG analysis design and results . . . . . . . . . . . . . . . . . . . . 80 6.4.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 6.4.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89

6.1. Regularized Generalized Canonical Correlation Analysis 71

Table 6

 6

		.3: Different sources of prior knowledge graphs
		# nodes # edges Diameter Radius % isolated nodes
	PC	15710	841690	6	4	0.20
	MSIGDB	10463	82962	6	4	0.46
	KEGG	776	11963	12	7	0.96

  . set is full PC graph Ref. set is full PC graph

			Candidate (822)	Stable (423)
		Edges removed
		in the pruned subgraphs	# of selected genes	# of selected genes
		Inner	822	423
	Candidate	Outer	586	262
		All	238	107
		Inner	822	423
	Stable	Outer	450	220
		All	239	108

Table 6 .

 6 5: Performances in survival prediction, number of selected variables and pathways depending on the type of graph to constrain the model.

		C-index		Stable variables selected in block No. of pathways
		Validation set Test set RNA MiRNA	CNV	from the RNA
	no graph	0.708 ± 0.122	0.588	0	1	208	0
	raw PC	0.692 ± 0.127	0.632	357	0	0	3
	normalised PC	0.709 ± 0.141	0.589	423	1	0	3
	normalised MsigDB 0.741 ± 0.092	0.625	232	1	0	2
	normalised KEGG	0.712 ± 0.110	0.690	106	1	0	2

Table 6

 6

	.6: Top gene sets from MSigDB C6 collection. In bold, gene sets with an adjusted p-value lower
	than 0.05.
	Term	Adjusted P-value
	CAHOY ASTROGLIAL	0.000384
	ATF2 UP.V1 DN	0.016688
	TGBF UP.V1 DN	0.041445
	KRAS.KIDNEY UP.V1 UP	0.104489
	RAF UP.V1 DN	0.260525

92 Radio-genomics integration and survival prediction Chapter Outline Contents 7

  been found co-(de)regulated in other studies of low-grade gliomas and other brain tumours. Future applications should focus on extending the results to other tumour types. Additionally, the multi-block model should enlarge the scope of multi-modality beyond molecular data. New data sources can be investigated, such as imaging data with their specific penalty resources, in order to increase prediction performance or unveil shared .1 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.1.1 Other Multi-block survival models from the state-of-the-art . . . . . . . . . 7.1.2 Experimental design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2 Baseline results using all per-block available data . . . . . . . . . . . . . . . . . . . 7.2.1 Radiomics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2.2 Somatic mutations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2.3 RNA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.3 Comparison mono-block and multi-block approaches . . . . . . . . . . . . . . . . . 7.3.1 Mono-block . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.3.2 Multi-block . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.4 Discusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	7.1. Methodology	93
	7.1 Methodology
	7.1.

1 Other Multi-block survival models from the state-of-the-art

Table 7 .

 7 3: Selected radiomic features extracted from the FLAIR modality. The table only shows radiomic features that were selected more than 6 times. The means and std were computed using only the samples where the feature was selected.

	Selected FLAIR radiomics	# times selected coeff	min	max
	wavelet.LHL_glcm_Correlation_flair	10/10	-0.0124 ± 0.0133 -0.0360 -0.0001
	wavelet.LHL_glrlm_RunEntropy_flair	9/10	-0.0047 ± 0.0035 -0.0124 -0.0008
	wavelet.LLL_glrlm_GrayLevelNonUniformity	8/10	0.0577 ± 0.0257	0.0253	0.1049
	wavelet.LLL_glcm_SumEntropy_flair	8/10	-0.0029 ± 0.0017 -0.0063 -0.0010
	wavelet.LLL_glcm_JointEnergy_flair	8/10	0.1065 ± 0.0604	0.0185	0.2133
	wavelet.LLL_firstorder_Uniformity_flair	8/10	0.0577 ± 0.0228	0.0282 0.0996
	wavelet.LLL_firstorder_Entropy_flair	8/10	-0.0028 ± 0.0017 -0.0064 -0.0010
	wavelet.HHL_glcm_Correlation_flair	8/10	-0.0504 ± 0.039	-0.1163 -0.0062
	wavelet.LLL_glcm_MaximumProbability_flair	7/10	0.0473 ± 0.0268	0.0158	0.0962
	wavelet.LLL_glcm_JointEntropy_flair	7/10	-0.0016 ± 0.001 -0.0034 -0.0004
	wavelet.HHL_glcm_InverseVariance_flair	7/10	0.1082 ± 0.099	0.0016	0.2910
	wavelet.HLH_glszm_SmallAreaEmphasis_flair	7/10	-0.0202 ± 0.0161 -0.0393 -0.0013
	wavelet.HLH_glszm_SizeZoneNonUniformity	7/10	-0.0111 ± 0.0085 -0.0203 -0.0003
	CE T1 radiomics. The Cox models trained on the radiomics from the Contrast Enhanced T1
	modality (CE T1) have similar results to the FLAIR, with a slightly lower average validation C-index
	of around 0.64 but a higher test C-index with an average of 0.56. However, as shown in figure 7.3,
	this test C-index does not reflect a higher ability of the model to generalised compared to the FLAIR,

Table 7 .

 7 4: Selected radiomic features extracted from the CE T1 modality. The table only shows radiomic features that were selected more than 3 times. The means and standard deviation were computed using only the samples where the feature was selected.

	Selected CE T1 radiomics	# times selected coeff	min	max
	wavelet.HHL_firstorder_Skewness_t1Gd	9/10	-0.0077 ± 0.0049 -0.0179 -0.0021
	wavelet.HHL_firstorder_Kurtosis_t1Gd	7/10	-0.0001 ± 0.0001 -0.0002 -0.0000
	wavelet.LLL_firstorder_Kurtosis_t1Gd	6/10	0.002 ± 0.0004	0.0013	0.0022
	wavelet.HHL_glszm_SizeZoneNonUniformity	5/10	-0.0877 ± 0.0698 -0.1787 -0.0036

Table 7 .

 7 5: Selected radiomic features extracted from all modalities. Out of 10 bootstrap samples, we only kept the eight with a convergent model. The table only shows radiomic features that were selected more than 3 times. The means and std were computed using only the samples where the feature was selected.

	Selected radiomics ALL MODALITIES	# times selected coeff	min	max
	wavelet.LLL_firstorder_Kurtosis_t2	7/8	0.2631 ± 0.1285	0.1080 0.4403
	wavelet.LLL_firstorder_Kurtosis_t1Gd	4/8	0.0461 ± 0.0442	0.0050 0.0973
	wavelet.HHL_firstorder_Skewness_t1Gd	4/8	-0.1052 ± 0.0996 -0.2531 -0.0379
	wavelet.LLL_firstorder_Uniformity_flair	4/8	1.2054 ± 1.1685	0.0655 2.7320

  Using the influence graph to propagate the somatic mutations gave similar performances to that obtained with the PC graph, with

	genes -Prop. Full # times selected coeff	min	max
	PALMD	8/10 144.1463 ± 116.9681	6.0754 296.3065
	CHEBI.7815	6/10 -32.797 ± 18.6147	-52.2561	-11.6946
	SRSF10	6/10 -0.9108 ± 0.713	-2.0581	-0.1349
	MED13	6/10 1.6335 ± 1.1587	0.5644	3.4945
	SLCO1A2	6/10 1.7941 ± 1.0129	0.6754	3.3610
	Mutations Propagated according to the PC influence graph.

Table 7 .

 7 9: Enrichr Enrichment results when using the five most selected genes using the influence graph.

	Selected genes -Prop. Influence # times selected coeff		min	max
	PALMD	10/10 26.5382 ± 18.8677 5.7802 51.9759
	IDH1	8/10 -0.2103 ± 0.1712 -0.4126 -0.0006
	CLTC	7/10 0.6579 ± 0.4416	0.0949	1.1603
	EGFR	6/10 0.3668 ± 0.3568	0.0871 0.8523
	CDR2	6/10 -0.1142 ± 0.1035 -0.2337 -0.0033
	Gene_set	Term	Adjusted P-value Genes
	MSigDB_Hallmark_2020 Protein Secretion		0.001891 CLTC;EGFR
	MSigDB_Hallmark_2020 PI3K/AKT/mTOR Signaling		0.001891 CLTC;EGFR
	MSigDB_Hallmark_2020 Glycolysis		0.004552 IDH1;EGFR
	KEGG_2021_Human	Central carbon metabolism in cancer		0.006356 IDH1;EGFR
	KEGG_2021_Human	Endocytosis		0.040868 CLTC;EGFR

Table 7 .

 7 10: Survival prediction results obtained on the 419 patients from the TCGA-LGG dataset using RNA features. Validation results were obtained on three-fold cross-validation. Test results were obtained using 10 bootstraps. Elastic-Net Cox was used.

	Validation	Test
	RNA 0.83 ± 0.02 0.79 ± 0.03

Table 7 .

 7 11: Sample size comparison between the baseline setting and the current setting devoted to monomulti-block investigation.

		Baseline study (sec. 7.2) This study (sec. 7.3.1 & 7.3.2)
	Radiomics	95	83
	Mutations	419	83
	RNA	419	83
	Blocks concatenated	-	83

Table A .

 A Tableonlyshows genes that were selected on all the runs. 2: Selected genes using the RNA data. Table shows genes that were selected in all runs. et une meilleure résolution spatiale. Ce qui a permis de caractériser les tumeurs de manière non invasive. L'analyse radiomique des IRM comprend toutes les étapes d'extraction et d'analyse d'un grand nombre de caractéristiques (plus de 200) d'une région d'intérêt (ROI) extraite des images cliniques. Ces caractéristiques décrivent la taille et la forme parmi d'autres propriétés visuelles apparentes des tumeurs. Puis des statistiques de premier ordre du signal d'intensité provenant des ROI sont calculées.

	122		Résumé en français
	Selected genes # times selected coeff CLTC 6/6 0.4685 ± 0.1167 Selected genes # times selected coeff EDA 6/6 SLC28A2.AS1 10/10 -0.0155 ± 0.0071 -0.0259 -0.0071 min max 0.2736 min max 0.5710 -0.0528 ± 0.0107 -0.0664 -0.0404 NLRP9 6/6 0.2261 ± 0.1934 0.0211 PLAA 10/10 -0.0553 ± 0.0203 -0.0880 -0.0258 0.5266 SLC45A4 6/6 0.1353 ± 0.1356 DZIP1L 10/10 0.0202 ± 0.0125 0.0041 0.0384 0.0348 0.3950 NAA11 6/6 0.1733 ± 0.1833 0.0221 BTG3 10/10 0.0102 ± 0.0056 0.0021 0.0208 0.4143 PHRF1 6/6 SFRP2 10/10 -0.0047 ± 0.0022 -0.0091 -0.0027 -0.1101 ± 0.1325 -0.3799 -0.0435 KCNS3 6/6 0.195 ± 0.1929 0.0229 COL8A1 10/10 0.0055 ± 0.0022 0.0016 0.0100 0.4788 PNMA3 6/6 0.167 ± 0.1359 ENTPD1 10/10 -0.0341 ± 0.0203 -0.0755 -0.0058 0.0020 0.3889 ZBED9 6/6 0.2607 ± 0.114 0.1019 SLC28A2 10/10 -0.0271 ± 0.0078 -0.0366 -0.0125 0.3695 NEBL 6/6 0.3167 ± 0.171 0.0526 FHDC1 10/10 -0.0089 ± 0.0047 -0.0157 -0.0022 0.5574 SLC11A1 6/6 0.4093 ± 0.2297 0.0117 WDR38 10/10 0.0035 ± 0.0017 0.0005 0.0057 0.5964 EDAR 6/6 0.5349 ± 0.0856 0.4519 RHOQ.AS1 10/10 0.0084 ± 0.0071 0.0020 0.0235 0.6780 SRSF10 6/6 LTV1 10/10 -0.0301 ± 0.0155 -0.0487 -0.0029 -0.3321 ± 0.2118 -0.6519 -0.1169 AFP 6/6 0.1447 ± 0.1228 0.0000 CASP9 10/10 -0.0261 ± 0.0058 -0.0341 -0.0182 0.3120 MED13 6/6 0.3441 ± 0.1756 0.1252 ELL3 10/10 0.0154 ± 0.0097 0.0014 0.0355 0.5231 IGHG4 6/6 0.2067 ± 0.1587 0.0524 H2BC11 10/10 0.0084 ± 0.0045 0.0018 0.0168 0.3725 OR1N2 6/6 0.0988 ± 0.1107 0.0202 CAAP1 10/10 -0.0355 ± 0.0173 -0.0623 -0.0157 0.3194 EGFR 6/6 0.2044 ± 0.1664 0.0248 NEK6 10/10 0.0093 ± 0.0067 0.0006 0.0187 0.4153 FLG 6/6 0.2363 ± 0.1303 MTAP 10/10 -0.0239 ± 0.0106 -0.0391 -0.0081 0.0959 0.4043 CLMN 6/6 NUAK2 10/10 0.0182 ± 0.0095 0.0029 0.0359 -0.2639 ± 0.1702 -0.5105 -0.1151 SLC4A3 6/6 0.2097 ± 0.1922 0.0091 SLC9B1 10/10 -0.0303 ± 0.0158 -0.0661 -0.0100 0.4483 ELK1 6/6 0.1748 ± 0.1876 TNFRSF11B 10/10 0.0062 ± 0.0036 0.0016 0.0132 0.0250 0.4886 SYBU 6/6 CCDC178 10/10 -0.0166 ± 0.0073 -0.0300 -0.0100 -0.0536 ± 0.0107 -0.0669 -0.0405 PDCD2L 6/6 LINC02198 10/10 0.0113 ± 0.005 0.0045 0.0188 -0.0529 ± 0.0107 -0.0665 -0.0403 IGHV3.43 6/6 LINC02044 10/10 -0.0138 ± 0.0101 -0.0356 -0.0013 -0.053 ± 0.0106 -0.0666 -0.0403 RP11.159L20.2 6/6 FLG.AS1 10/10 -0.0081 ± 0.0039 -0.0157 -0.0040 -0.0531 ± 0.0106 -0.0666 -0.0402 OR5B21 6/6 HLA.DQB2 10/10 0.006 ± 0.0027 0.0012 0.0092 -0.0534 ± 0.0105 -0.0667 -0.0401 OR52D1 6/6 EHBP1.AS1 10/10 0.0116 ± 0.0089 0.0010 0.0280 -0.0534 ± 0.0105 -0.0667 -0.0402 TBATA 6/6 CBWD4P 10/10 -0.027 ± 0.011 -0.0436 -0.0107 -0.0535 ± 0.0106 -0.0668 -0.0402 MASTL 6/6 MARCHF5 10/10 -0.0362 ± 0.0141 -0.0514 -0.0099 -0.0535 ± 0.0106 -0.0669 -0.0404 MIR133B 6/6 LINC00994 10/10 -0.0092 ± 0.0034 -0.0137 -0.0028 -0.0536 ± 0.0108 -0.0670 -0.0407 GPBP1L1 6/6 0.2042 ± 0.1514 0.0659 H1.9P 10/10 0.0064 ± 0.0018 0.0034 0.0098 0.4142 PRKAA1 6/6 NUTM2A 10/10 -0.0141 ± 0.0106 -0.0321 -0.0004 -0.0536 ± 0.0109 -0.0670 -0.0409 ORC4 6/6 ZNF648 10/10 0.01 ± 0.0049 0.0006 0.0168 -0.0535 ± 0.0111 -0.0671 -0.0411 ACBD3 6/6 GOLGA6L2 10/10 -0.0075 ± 0.0038 -0.0129 -0.0005 -0.0536 ± 0.0112 -0.0671 -0.0413 ALPK3 6/6 0.2049 ± 0.1899 TDRD12 10/10 -0.0149 ± 0.0084 -0.0316 -0.0046 0.0403 0.4929 PCDHGA5 6/6 0.2339 ± 0.1979 MLLT3 10/10 -0.0424 ± 0.0179 -0.0622 -0.0114 0.0396 0.5050 EPN3 6/6 0.2596 ± 0.1542 0.0450 PLEKHA2 10/10 0.0309 ± 0.0177 0.0090 0.0652 0.4117 PPP3R2 6/6 0.2117 ± 0.1243 0.0161 0.3088 KCNN2 10/10 -0.0106 ± 0.0076 -0.0232 -0.0002
	IDH1	6/6	-0.1428 ± 0.0543 -0.2377 -0.0916

Acknowledgements

Résumé en français

Les données manquantes sont inhérentes aux études multi-blocs car toutes les données doivent être disponibles pour tous les patients et toutes les sources. A ce fait, s'ajoute dans le cas du DIPG, la collecte de données dans un contexte clinique qui est difficile à réaliser. Cette thèse n'a pas exploré la gestion des données manquantes dans l'intégration multi-blocs, ce qui a entraîné l'utilisation d'un ensemble de données relativement petit. Cela pourrait aussi constituer une direction nouvelle de recherche.

Enfin, ce travail propose l'utilisation de fonctionnalités radiomiques dans un cadre multi-blocs.

Il faut cependant convenir, que les résultats obtenus ont pu être biaisées par les paramètres des séquences d'acquisition, les algorithmes de reconstruction et les méthodes de prétraitement. Des travaux supplémentaires considérant un vaste jeu de données permettant un contrôle strict de toutes les sources de variabilité permettrait de conclure plus avant sur les signatures d'image des tumeurs sans un contrôle strict de toutes les sources de variabilité.

.gov/ct2/show/NCT05476939 3 https://www.cancer.gov/about-nci/organization/ccg/research/structural-genomics/tcga 4 https://www.cancerimagingarchive.net/ 5 https://www.gsea-msigdb.org/gsea/msigdb 6 https://www.pathwaycommons.org/

//clinicaltrials.gov/ct2/show/NCT02233049 2 https://clinicaltrials

Multiblock Integration Method information between modalities.

✽ ✽ ✽ ✽ ✽ ✽

Chapter 7

Radio-genomics integration and survival prediction

Radiomic features have been successfully used for different tasks and in various tumour types. This includes the DIPG, Glioblastomas and Low-Grade Gliomas (LGG). However, it is unclear which image modality is the most informative or which feature type is the most important. Furthermore, research has also shown that adding radiomic features to molecular data can improve survival prediction. But

it has yet to be integrated into a full multi-block framework. This should help us understand the interactions between the radiomic and genomic features.

So far, we have proposed a procedure to extract the regions of interest required for the radiomic analysis. Additionally, we proposed the netSGCCA, a multi-block model that considers an a priori graph of interactions. We also established basic properties that help us understand the interaction between the multi-block model and the graph penalty.

The current chapter is devoted to applying netSGCCA for survival prediction when genetic data are associated with radiomic features. We pay close attention to the selected variables and their link to the studied disease. Since the multi-block framework we are working with requires data availability for all sources on all samples, it naturally results in fewer samples being analysed. Therefore, we establish the baseline performance of each block when all available data is considered before the integration.

Then, we compare the mono-block framework with the multi-block framework. We also compare the netSGCCA with the best-performing models in the related works. Due to the unavailability of survival data on the DIPG dataset, this study was restricted to the TCGA-LGG dataset.

Baseline results using all per-block available data

The multi-block analysis requires data availability for all subjects across all the blocks. This results in smaller sample sizes when using multi-block approaches compared to mono-block strategies. In order to compare the multi-block and mono-block approaches, we started by studying the predictive performances of each block when all available data is used for this block. This allows us to establish a baseline of the predictive power of each block, from which the multi-block study can be judged.

This section does not intend to compare the performances of the different blocks as they are studied across different sets of patients. However, it should show the major differences between the different experimental settings.

Radiomics

Table 7.1 shows the survival prediction results in terms of C-index using the radiomic features on 95 patients from the TCGA-LGG. Figure 7.3 displays the shape of the distributions of the C-index values revealed by the ten bootstrap runs. We discuss further the individual contribution of the radiomics brought by each modality. T1 radiomics. The Cox model could not give good predictions using the T1 modality. On both the validation and test sets, the C-index score remained around 0.5, equivalent to a random ranking of the patients. This result was expected as LGGs are manifested in small areas in the T1-weighted images. The regions of interest (ROI) used to extract the radiomics are dominated by oedema, which is not visible in the T1 weighted images, and provide non-informative radiomics.

T2 radiomics. The Cox model using the radiomics of the T2 weighted modality gave good results on the validation and test sets with an average C-index above 0.75. A closer inspection of the results obtained on the test set shows that 7 of the 10 bootstrap runs scored above 0.84, while the remaining 3 scored 0.5 due to convergence issues. This explains the standard deviation of 0.18 on the test set. The optimal parameters for the T2 radiomics were α = 0.03 and r = 0.85, which allowed the selection of, on average, 15 variables (except for the three non-convergent samples, where the null model has been returned). Looking only at the seven successful samplings, 64 of the 744 radiomics were selected at least once, and six variables were selected more than three times. Table 7.2 shows the number of times and the average coefficient of these six most selected variables. All features in the table kept the same PC graph propagation PC influence graph propagation Figure 7.5: Number of times the selected genes appear to be mutated in the original somatic mutation profiles, using graph propagation (left) with the PC graph, and (right) with the influence graph.

RNA

We investigated the survival prediction performances of the Cox model when using gene expression data obtained from RNA sequencing obtained on the 419 patients from the TCGA-LGG dataset.

Results obtained using RNA data have been the best yet, with a C-index of 0.83 on the validation set and 0.79 in the test set, on average, as shown in Table 7.10 and Figure 7.6. The best-performing parameters were α = 0.55 and r = 0.03 allowing the selection of 1104 genes on average from 19443 available genes. Among the 11078 uniquely selected genes, only 3079 were selected on over half the runs, and 38 were chosen in all runs (shown in Table A.2). Interestingly, none of these 38 genes coincides with any of the frequently selected genes obtained on the various somatic mutation models.

Using the frequently selected genes, we did not find significantly enriched pathways from the KEGG, MSigDB Hallmarks or the MSigDB oncogenic signatures collections. Radio-genomics integration and survival prediction and only 10 were chosen in all the runs. The ten genes selected in all the runs are shown in Table 7.14. Compared to the baseline, using the 83 individuals resulted in more genes selected. Additionally, all the consistently selected genes were not identified during the baseline phase. Interestingly, three consistently selected genes had coefficients which changed signs between the different runs. This adds inconsistency to the contribution of each gene to the disease outcome. Enrichment analysis using the ten most selected genes revealed four significantly enriched pathways from the KEGG collection, as shown in Table 7.15. These pathways were not identified in the baseline phase but have also been linked

to gliomas in previous studies [START_REF] Guo | Lipid metabolism emerges as a promising target for malignant glioma therapy[END_REF][START_REF] Huang | A prospective study of serum metabolites and glioma risk[END_REF][START_REF] Kim | Peroxisome metabolism in cancer[END_REF]. These pathways are linked to each other and to the previously identified genes in the baseline.

RNA block alone.

Using the gene expression profiles, the Cox model had good performances on the validation set, with a C-index of 0.75 on average. However, the results were lower on the test set across the ten bootstrapped runs. On average, the C-index on the test set was 0.59. As shown in Figure 7.7, three runs had performances below the random results with a C-index around 0.35.

Comparison mono-block and multi-block approaches 107

Table 7.14: Selected genes using the somatic mutation profiles. The table shows genes that were selected in all 10 runs. The means and std were computed using only the samples where the feature was selected.

Selected regularising parameters were α = 1.15 and r = 0.15, respectively. This resulted in 13 genes selected on average. However, the variable selection was unstable; most of the genes were selected only once, and only three were selected in more than three runs. Additionally, no gene was selected more than half the runs (five times), which explains the results instability. As a consequence, we do not present any genes from this block.

Radiomics, Mutations and RNA block concatenated. When using all data sources concatenated into a single matrix, the Cox model performed well on the validation set, with a C-index of 0.83.

However, as seen by the bootstrapped test runs, the models failed to generalise, which obtained an average C-index of 0.53. Among the ten runs, only three did not result in the null model, and their C-indices were 0.81, 0.47 and 0.86. The learnt sparsity and regularizing parameters were α = 2.53 and r = 0.01, respectively. This resulted in 1521 selected variables on average (only considering the three successful runs). As expected, most selected variables came from the Mutation block, which has shown the best result so far, as shown in Figure 7.9(a). We note that run_3 had the lowest C-index, which corresponds to the run where the Cox model using the mutations alone would not generalise.

Among all the variables selected, 262 were present on the three successful runs. Enrichment analysis on consistently selected variables from the mutation block highlighted the Phototransduction and Purine metabolism pathways.

Summary.

In general, removing subjects from the datasets considerably changes the performance of trained model especially as regards its generalization capabilities and selection stability. RNA data was the best-performing block in the baseline but performed poorly when only 83 subjects were kept from the whole dataset. Furthermore, mutation results gave similar performances between the baseline and the reduced mono-block model. However, the selected genes and their significant enriched pathways differed. For radiomic features, the results from the baseline and this mono-block study were similar, which was expected given the low number of removed subjects. Additionally, we remark that the concatenated results reflected the results obtained on the individual blocks.

The mono-block Cox model was generally unstable in terms of variable selection. Additionally, multiple runs failed to converge even with similar samples and parameters. The Cox model appears to be initialisation dependent, which weakens its reliability.

As an interim conclusion, these results question the reliability of the mono-block Cox model we obtained. It also raises questions about the interpretations drawn from this model that learns in a mere (concatenation) mono-block with the classical ElasticNet penalty.

Multi-block

In this section, we compared the different multi-block models using the 83 subjects having all data blocks available. The first two multi-block models we considered here are SGCCA, netSGCCA; They are extensions of the RGCCA framework that include respectively sparsity [START_REF] Tenenhaus | Variable selection for generalized canonical correlation analysis[END_REF] and graphical constraint developed in section 6.3. We compared these two models to two other approaches established independently: Priority-Lasso and BlockForest(see 7.1.1). The cross-validation and bootstrapped sets are the same as those used in the mono-block study. This is done to obtain comparable results. Similarly to the mono-block study, we used the RNA block, the radiomic features obtained from the four modalities and the somatic mutation profile processed with the influence graph 7.17. The C-index scores distribution on the test set is presented in Figure 7.10.

Table 7.17: Multi-block results obtained on 83 subjects having the radiomics, somatic and RNA data available. For the radiomics, the concatenation of all the radiomics extracted from the four modalities has been used; The somatic mutations with graph propagation with the influence graph was used.

Validation Test netSGCCA 0.88 ± 0.14 0.67 ± 0.10 SGCCA 0.88 ± 0.14 0.66 ± 0.12 Priority-Lasso 0.71 ± 0.11 0.54 ± 0.16 BlockForest 0.68 ± 0.06 0.72 ± 0.15 SGCCA and netSGCCA. Using the netSGCCA with the normalised Pathway Common graph penalty on the RNA, we obtained high performances on the validation set, with an average C-index of 0.88. However, the performance dropped to 0.67 for the test set, which is higher than most monoblock models but below the mono-block model based on the Mutations. On average, in each test run, the model selected 1800 variables. Among them, around 47 are selected from the Mutations block, 1752 from the RNA block and one variable from the radiomics.

Table 7.18 shows the 17 most frequently selected features from the RNA block. It exhibits the instability of the model in terms of variable selection. Two reasons can explain this. First, as shown earlier, the RNA block is not a good predictor for survival in this setting sample (n = 83, see also 7.3.1): Even if the GraphNet penalty allows more variables to be selected from this block, we can not expect the RNA block to yield stable selected variables. Second, the GraphNet penalty groups

Radio-genomics integration and survival prediction correlated features and using the bootstrap changes the correlation profiles between the genes, as shown in Figure 7.11. The correlation profile of the gene PGPEPE1 for example, is volatile across the folds. The instability of variable selection on the other blocks is somewhat inherited from the instability of the SGCCA, which has already been highlighted in the previous chapter.

Enrichment analysis results of the 17 most selected genes (RNA block) are reported in table 7.19.

The IL15 gene set has immune functions and has already been shown to improve anti-glioma activity [START_REF] Krenciute | Transgenic Expression of IL15 Improves Antiglioma Activity of IL13Rα2-CAR T Cells but Results Bibliography 137 in Antigen Loss Variants[END_REF]. The KRAS has been linked to increased risk of LGGs [START_REF] Guan | KRAS gene polymorphisms are associated with the risk of glioma: a two-center case-control study[END_REF][START_REF] Ryall | Pediatric low-grade glioma in the era of molecular diagnostics[END_REF]. PKCA has long been identified as a key regulator of the growth of malignant gliomas [START_REF] Arcos-Montoya | Progesterone receptor together with PKCα expression as prognostic factors for astrocytomas malignancy[END_REF]Baltuch et al., 1995;[START_REF] Cameron | Pkc alpha protein but not kinase activity is critical for glioma cell proliferation and survival[END_REF]. Table 7.18: Selected features using the netSGCCA, the normalised Pathway Commons Laplacian was applied on RNA data. The table shows features that were selected more than 5 times. The means and std were computed using only the samples where the feature was selected.

Selected radiomics # times selected coeff min max THRA 6/10 0.0001 ± 0.0015 -0.0023 0.0018 BID 5/10 -0.0009 ± 0.0023 -0.0037 0.0026 LINC02593 5/10 -0.0236 ± 0.0285 -0.0480 0.0180 DNM3 5/10 -0.0011 ± 0.0027 -0.0054 0.0017 MRTFA 5/10 -0.0012 ± 0.0019 -0.0041 0.0012 REPS2 5/10 -0.0017 ± 0.0032 -0.0074 0.0006 COL6A2 5/10 0.0004 ± 0.0023 -0.0029 0.0033 SPRY2 5/10 0.0018 ± 0.0053 -0.0050 0.0097 LINC00609 5/10 0.0027 ± 0.0428 -0.0402 0.0566 PGPEP1 5/10 0.0007 ± 0.0008 0.0001 0.0021 HS3ST3B1 5/10 0.0008 ± 0.0015 -0.0010 0.0030 CCNY 5/10 -0.0008 ± 0.0009 -0.0024 0.0001 MROH8 5/10 0.0005 ± 0.0019 -0.0012 0.0032 WASHC2A 5/10 -0.0009 ± 0.0019 -0.0028 0.0020 ABCA7 5/10 0.0011 ± 0.0023 -0.0010 0.0047 HOXC8 5/10 0.0024 ± 0.0041 -0.0027 0.0086 MAPT-AS1 5/10 -0.0249 ± 0.0462 -0.0733 0.0499

Discusion 111

Table 7.19: Enrichment results when using the most selected genes from the RNA block using netSGCCA. Further analysis revealed that removing the radiomic block from the multi-block netSGCCA model increased the C-index on the test runs, with an average score of 0.8. It should be noted that we did not run proper hyperparameters selection for these two models, thus the results are only indicative and not included in the result tables. On the other hand, using only the RNA and radiomic blocks yielded a lower average C-index of 0.62 on the test set.

SGCCA coupled with a Cox model had a similar performance as the one obtained with netSGCCA coupled with a Cox model on the validation and test sets, with a C-index of 0.88 and 0.66, respectively.

However, the model with SGCCA was more sparse on the RNA block by selecting, on average, only three variables. Variable selection on the other blocks did not meaningfully change between SGCCA and netSGCCA. These results confirm previous conclusions discussed in the last chapter.

Priority-Lasso and BlackForest. Using Priority-Lasso, we obtained an average C-index of 0.71 on the validation set. After a closer inspection of the obtained performances on the validation set, we observed an average C-index of 0.70 with a standard deviation of 0.14 using the RNA alone.

Adding the mutations did not significantly change the results, with an average C-index of 0.71 and a standard deviation of 0.15. Adding the radiomic features did not change the average C-index from 0.71; however, the standard deviation decreased to 0.11. Radiomic features helped the stability of the model. Priority-Lasso could not generalise the results obtained on the validation set into the test set with a mean C-index of 0.54 on the ten bootstrapped runs. Among the ten runs, three performed below 0.4, which is worst than average. The instability is reflected in the selected variables. In the ten runs, 50 unique variables were selected ( with an average of 7.6 variables per run ), and no variable was selected in more than six runs. Priority-Lasso only uses a Lasso penalty which does not add regularity to the model estimates, which explains the low number of selected variables and their instability.

Compared to the other multi-block models, BlockForest performed the best on the test set, with a C-index of 0.68 on the validation set and 0.72 on the test set. netSGCCa and SGCCA gave slightly lower but comparable results. However, Block Forest does not allow variable selection, which hinders the interpretability of the model and the extraction of genes of interest.

Discusion

The netSGCCA was among the best-performing models on the validation and test set. The C-index was stable during all test runs, which was not the case for the mono-block Cox models (see 7.3.1). The netSGCCA model was only topped by BlockForest among all multi-block models. However, it should be noted that the BlockForest is an unconstrained model and does not perform variable selection,

First Appendix

115

Second Appendix

We want to prove the grouping effect of GraphNet in the context of netSGCCA. Our proof is similar to [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF] and [START_REF] Li | Network-constrained regularization and variable selection for analysis of genomic data[END_REF]. We assume that two variables u and v are only connected to each other in the a priori graph. We also assume that u and v have similar contributions (i.e. their respective weights have the same sign). Give the lagrangian equation defined in 6.6:

argmin

Under the optimality conditions of w (g) (the weight of the block on which GrahNet is applied), we have:

For the variables u and v, let x u and x v be the columns v and u from the matrix X (g) . For these two variables, the equality becomes :

We subtract the two equations, with sign(w u ) = sign(w v ), which becomes:

Second Appendix

Using the Cauchy-Schwarz inequality, we can see that :

And x ⊤ v x u = ρ the estimated correlations between the variables u and v, we can show equation 6.7: