Anne-Claire Brunet

Rémy Burcelin

Jean-Michel Michel Loubes

Magali Champion

Radu M Neagoe

Maria Effernberger

Daniela T Sala

Florence Servant

Jeffrey E Christensen

Maria Arnoriaga- Rodriguez

Jacques Amar

Benjamin Lelouvier

Fabrice Gamboa

Herbert Tilg

Massimo Federici

Jose-Manuel Fernández-Real

Keywords: clustering

graphes . . . . . . . . . .

Chapitre 1 Introduction Générale 1 Contexte biologique 1.1 Les maladies métaboliques et cardiovasculaires

Ces dernières années, les progrès technologiques ont fortement marqué, à l'échelle mondiale, le mode de vie des populations notamment le rythme quotidien, la mobilité et la qualité d'alimentation. Ayant dû s'adapter à un environnement aussi changeant, les organismes, qui ne sont pas programmés biologiquement pour vivre dans ces nouvelles conditions sociétales, ont développé des maladies chroniques notamment métaboliques telles que l'obésité et une de ses conséquences graves, le diabète de type 2. Ce dernier est caractérisé par une augmentation de la concentration de glucose dans le sang le matin, à jeun ou suite à un repas, dû à un défaut de sécrétion et à un défaut d'action de l'insuline. Au-delà des aspects biochimiques et physiologiques, ces dérèglements métaboliques sont des facteurs de risque des maladies cardio-vasculaires, hépatiques et rénales et peuvent accentuer la sensibilité aux infections notamment virales. Ces maladies métaboliques augmentent le tonus inflammatoire de manière prolongée sur plusieurs années ce qui aboutit à des atteintes cellulaires pouvant conduire jusqu'à la mort cellulaire. Ainsi, l'élimination des cellules mortes et la réparation tissulaire dans le foie entraîne une fibrose progressive difficilement réversible. Cette fibrose peut aboutir à une insuffisance hépatique, une cirrhose et un cancer potentiel. Parmi les maladies hépatiques, nous retrouvons notamment la Non-Alcoholic Fatty Liver Disease (NAFLD), caractérisée par l'accumulation de graisses dans le foie qui, chez certains individus, peut altérer le fonctionnement des cellules hépatiques et provoquer une inflammation, la Stéatite Hépatique Non Alcoolique (NASH) (Figure 1.1).

Parmi les causes de ces évolutions pathologiques, la génétique exerce un impact majeur dans leur développement. Cependant, pour beaucoup de pathologies, de multiples facteurs liés à l'environnement entrent aussi en jeu et plus particulièrement la qualité et quantité alimentaire, la sédentarité et le stress. Les interrelations gènes/environnement à l'origine des maladies métaboliques sont d'autant plus difficiles à détecter que de nombreuses modifications moléculaires modulant l'expression des gènes peuvent interférer. Dans le cadre de l'obésité et du diabète de type 2, plusieurs études récentes ont révélé que le microbiote intestinal (l'ensemble des micro-organismes qui peuplent l'intestin) conciliait tous ces paramètres. Sous l'influence des conditions alimentaires, de la génétique et de l'environnement, la composition du microbiote intestinal évolue tout au long de la journée pour D'un point de vue statistique, afin de rendre compte de manière la plus générale possible du phénomène biologique sous-jacent, une première étape dans l'analyse des données consiste à réaliser une Analyse Exploratoire des Données (AED) cliniques, génétiques et microbiologiques. Elle peut être suivie d'analyses prédictives ou discriminantes dépendant de la nature des variables à étudier. Ceci permet à la fois d'avoir une vision globale d'un ensemble de données, d'en chercher les corrélations sans hypothèse a priori et d'étudier plus en profondeur le processus biologique en question. L'objectif est alors d'identifier des signatures biologiques, aussi appelées biomarqueurs, associés aux premiers stades de la maladie afin de poser un diagnostic précoce qui permettraient de prévenir, par une action thérapeutique, le développement de la maladie.

Dans la section suivante, les différents types de données exploités sont présentés.

Les données 1.2.1 Génomiques

Le génome représente l'ensemble du matériel génétique nécessaire au développement de l'espèce vivante. Il repose sur l'ADN, support de l'information génétique. Celui-ci est une molécule constituée de deux brins complémentaires en forme de double hélice composés d'une succession de nucléotides, comme représenté sur la figure 1.2. On trouve quatre nucléotides différents, notés A, G, C et T, du nom des bases correspondantes adénine, guanine, cytosine et thymine, dont l'ordre d'enchaînement (codes) détermine l'information génétique. Le génome joue un rôle crucial dans le fonctionnement de l'organisme en permettant la production de protéines qui assurent une multitude de fonctions au sein des cellules et des tissus. Les protéines sont des macromolécules composées d'acides aminés remplissant les fonctions vitales de la cellule qui sont aussi responsables de l'ensemble du métabolisme de l'ADN (synthèse, réplication, réparation). L'ADN détermine la synthèse des protéines par l'intermédiaire de l'acide ribonucléique (ARN). Comme pour l'ADN, l'ARN est un support moléculaire de l'information génétique mais il est constitué d'un seul brin. L'information est elle aussi encodée par quatre bases nucléiques, la thymine (T) étant remplacée par l'uracile (U). Il existe différents types d'ARN classés selon leur fonction. Certains ARN dit codants, aussi appelés ARN messagers, portent l'information génétique codant pour des protéines. Les autres ARN interviennent dans le fonctionnement de la cellule, la régulation de l'information et l'activité cellulaire. Par des techniques biochimiques, il est possible de déterminer le code génétique (partiel ou complet) d'un organisme donné : c'est le séquenc ¸age du génome.

Transcriptomiques

Le transcriptome est l'ensemble des ARN produits par transcription de l'ADN. Sa caractérisation et quantification dans une cellule donnée se fait par le biais de puces à ADN qui intéragissent avec les molécules d'ARN ou directement par séquenc ¸age de l'ensemble des molécules d'ARN. Ces techniques, très largement utilisées, permettent notamment d'identifier les gènes actifs, de déterminer les mécanismes de régulation d'expression des gènes et de définir des réseaux d'expression des gènes.

Métagénomiques

La métagénomique est la méthode d'étude des métagénomes, c'est à dire de l'information génétique (le génome) de l'ensemble des micro-organismes vivant au sein d'un environnement spécifique. Les micro-organismes regroupent l'ensemble des êtres vivants microscopiques, comportant les virus, parasites, champignons microscopiques, levures et bactéries. Cet ensemble d'organismes, le microbiote, est codé par le génome correspondant, appelé microbiome, c'est-à-dire l'ensemble des gènes du microbiote. La génomique est souvent, et par défaut, l'étude du génome des eucaryotes (organismes supérieurs complexes). La métagénomique englobe les génomes des microorganismes dans un milieu donné. Chez l'homme et la souris de laboratoire, il s'agit souvent des génomes bactériens.

Chaque bactérie a son propre génome, dont une partie très variable de l'information est partagée avec les autres bactéries d'un milieu donné. La mise en commun de l'ensemble des génomes, notamment bactérien, est appelé métagénome, même s'il s'agit de plusieurs génomes différents mais dont une partie de l'information peut être commune. Une analyse typique métagénomique bactérienne donne la composition de la taxonomie d'un microbiome bactérien c'est-à-dire les taxa notamment bactériens présents, leur abondance ainsi que leur diversité. Il existe deux grandes stratégies de séquenc ¸age en métagénomique : la stratégie globale (Figure 1. La métagénomique globale permet d'extraire, fragmenter et déterminer le code ADN de l'ensemble des génomes des micro-organismes présents dans le milieu étudié (la plupart du temps et par défaut bactérien). Les séquences obtenues, aussi appelées reads sont ré-assemblées bioinformatiquement afin de reconstruire les génomes bactériens d'origine. La métagénomique ciblée consiste à étudier un unique gène au lieu du génome complet. Ce gène doit être commun à plusieurs espèces tout en présentant des régions suffisamment variables afin de discriminer une espèce. En bactériologie, par exemple, le gène utilisé est celui de l'ARN ribosomique 16S. Ce gène possède des régions constantes (donc similaires) entre plus de 80 -90% des bactéries et d'autres régions variables très spécifiques d'une espèce bactérienne donnée. Les régions variables sur ce gène n'ont pas de rôle fonctionnel important et peuvent diverger au cours de l'évolution sous l'effet de mutations neutres. C'est ce qui permet de classer les microorganismes en taxons (famille, genre, espèce), sur la base de distances génétiques (Figure 1.4).

FIGURE 1.4 -Taxonomie des bactéries

Une fois le séquenc ¸age réalisé, à chaque read doit être assigné le nom d'une bactérie, on parle d'assignation taxonomique. En particulier, les séquences de nucléotides sont regroupées avec un seuil général de 97% de similarité dans des Unités Taxonomiques Opérationnelles (OTUs). Ces dernières sont ensuite comparées à des séquences de données de référence contenues dans des bases de données publiques pour identifier la taxonomie des bactéries en cause. Une fois la liste des taxons identifiée, leur abondance est estimée en fonction du nombre de séquences alignées sur leur référence. Le résultat final est représenté sous la forme d'un tableau de comptage à double entrée contenant le nombre de séquences par OTU et par échantillon.

2 Modélisation statistique des données

Caractéristiques des données

On considère un système composé de p variables quantitatives X = (X 1 , . . . , X p ) et n observations (X j 1 , . . . , X j n ) j∈{1,...,p} de ces p variables. Dans cette thèse, nous étudierons deux types de données : les données transcriptomiques, dont la littérature dans le domaine statistique est très riche et les données métagénomiques, encore assez peu connues, qui nécessitent des développements particuliers dus à leur nature très différente. Les données transcriptomiques se présentent sous la forme d'un tableau à double entrée représentant l'activité génomique, c'est-à-dire son expression en ARNm (ou quantité d'ARN messagers produite) de chaque gène par échantillon. Le nombre de gènes dans un tissu donné chez une centaine d'individus pouvant atteindre jusqu'à plusieurs dizaines de milliers, les données transcriptomiques sont caractérisées par la grande dimension (p >> n). Les données microbiologiques se présentent quant à elles sous la forme d'une matrice de comptage correspondant au nombre de reads de chaque espèce bactérienne observée par échantillon. Cette table de comptage présente plusieurs caractéristiques singulières. Il s'agit d'une matrice :

-creuse : elle contient un grand nombre de comptages nuls, pouvant varier de 70 à 90% ce qui peut être interprété comme l'absence ou la présence non détectée d'un OTU (échantillonnage trop faible), -de grande dimension dans le sens où p >> n, -dont les variables sont dépendantes en raison des interactions écologiques entre les espèces ou à l'appartenance à une même entité biologique, -avec une grande dispersion : variance d'un OTU généralement plus élevée que sa moyenne.

Transformation et normalisation des données métagénomiques

Les données de comptage sont des données très hétérogènes et de distribution asymétrique. Des modèles basés sur des lois discrètes comme la loi de Poisson ou la loi binomiale négative ont été proposés afin de les simuler le plus correctement possible [START_REF] Gary | Analysis of frequency count data using the negative binomial distribution[END_REF][START_REF] Linden | Using the binomial distribution to model overdispersion in ecological count[END_REF]. L'analyse statistique du jeu de données non transformé peut ne pas révéler l'information sous-jacente du fait de ses nombreuses caractéristiques singulières. Certains outils statistiques nécessitent une transformation des données afin de pouvoir être utilisés ou optimisés. C'est le cas notamment de la régression linéaire pour laquelle l'hypothèse de normalité est primordiale. D'autres outils comme l'Analyse en Composantes Principales présentent une certaine sensibilité aux trop grandes dispersions des données. Les méthodes de transformation des données se déclinent en deux grandes familles :

-la standardisation ou l'action de centrer-réduire les données pour diminuer la dispersion des données, -la normalisation des données qui consiste à rechercher la transformation adaptée au travers de laquelle elles suivront une loi normale. Parmi les normalisations les plus utilisées, on trouve la transformation log suivante : ∀i ∈ 1, n , Xi = log(X i + 1).

En revanche, celle-ci ne permet pas de stabiliser la variabilité du jeu de données. En effet, le nombre de lectures obtenues de chaque séquence, aussi appelé profondeur des séquences, peut varier d'un échantillon à l'autre, notamment à cause d'imprécisions techniques et technologiques. La comparaison de deux échantillons de profondeurs différentes risque alors de donner des écarts de proportions de gènes non nécessairement liés à une réelle distinction biologique. Pour les rendre comparables, la normalisation classique utilisée est la normalisation Total Sum Scaling (TSS) qui divise chaque comptage par le nombre total de variables observées par échantillon :

∀i ∈ 1, n , Xi = X i p j=1 X j i .
La matrice de comptage est ainsi convertie en matrice de proportion aussi appelée matrice d'abondance. La normalisation TSS engendre cependant un biais lors de l'estimation des abondances différentielles. Afin de réduire l'influence des OTUs sur-représentés dans une matrice de type creuse, [START_REF] Paulson | Differential abundance analysis for microbial marker-gene surveys[END_REF] ont développé une stratégie alternative, la normalisation par somme cumulée (CSS). Celle-ci consiste à remettre à l'échelle les OTUs en fonction d'un sous-ensemble (quartile) de taxons à faible abondance (relativement constants et indépendants), réduisant ainsi l'impact des taxons à forte abondance (dominants). Dans cette thèse, cette dernière méthode sera appliquée aux données de comptage log-normalisées avant toute analyse statistique.

Analyse de la diversité bactérienne

Afin de contrôler l'homogénéité des échantillons en terme de composition taxonomique, une étude de la diversité bactérienne est très souvent réalisée. [START_REF] Whittaker | Evolution and measurement of species diversity[END_REF] a introduit les deux types d'études suivantes :

-l'étude de la diversité α, dite locale, utilisée pour mesurer la diversité au sein d'un unique échantillon, -l'étude de la diversité β, utilisée pour estimer la diversité des espèces bactériennes entre les échantillons. La diversité α est mesurée à partir d'indices parmi lesquels : l'indice richesse qui mesure le nombre d'entités bactériennes (taxons) distinctes présentes dans l'échantillon mais donne parfois trop d'importance aux espèces rares. L'indice de [START_REF] Shannon | A mathematical theory of communication[END_REF], ou entropie relative permet par contre de prendre en compte à la fois le nombre d'espèces dans milieu et la répartition des individus au sein de ces espèces. Il est calculé comme suit :

∀i ∈ 1, n , H i = - p j=1 X j i p j=1 X j i log     X j i p j=1 X j i     .
Un autre indice très utilisé est l'indice de [START_REF] Simpson | Measurement of diversity[END_REF] qui mesure la probabilité que deux individus sélectionnés au hasard appartiennent à la même espèce :

∀i ∈ 1, n , S i = p j=1 X j i (X j i -1) p j=1 X j i p j=1 X j i -1
.

et ainsi donne plus de poids aux espèces abondantes qu'aux espèces rares. Finalement, un dernier indice utilisé est l'indice [START_REF] Chao | Nonparametric estimation of the number of classes in a population[END_REF] qui donne plus de poids aux faibles abondances (présentes seulement une ou deux fois) pour estimer le nombre d'espèces non observées et est calculé comme suit :

∀i ∈ 1, n , C i = p + F 1 (F 1 -1) 2(F 2 +1) ,
où F 1 respectivement F 2 correspondent au nombre d'espèces représentées une, respectivement deux fois chez chacun des individus.

Les différents indices présentés ci-dessus reflètent chacun des propriétés spécifiques de la diversité taxonomique. L'analyse de la diversité α se base donc sur la combinaison de tous ces indices.

La diversité β est quant à elle quantifiée à partir de mesures de distance permettant de comparer les échantillons deux à deux et de leur attribuer une valeur résumant leur ressemblance globale. Notons Ω = {1, . . . , n} l'ensemble des individus. On se propose alors de définir sur Ω×Ω différentes mesures d'éloignement entre deux individus. On introduit les notions de dissimilarité et de similarité suivantes :

Définition 1 Une dissimilarité d est une application de Ω 2 dans R + vérifiant pour tout (i, j) ∈ Ω 2 -d(i, j) = d(j, i),

d(i, j) = 0 ⇔ i = j, -i = j, d(i, j) ≥ 0.

Définition 2 Une similarité s est une application de Ω 2 dans R + vérifiant pour tout (i, j) ∈ Ω 2 -s(i, j) = s(j, i), -s(i, i) = S > 0, -s(i, j) ≤ S.

Une matrice de similarité peut se transformer en matrice de distance à l'aide de la transformation :

∀(i, j) ∈ Ω 2 , d(i, j) = (s(i, i) + s(j, j) -2s(i, j)) -1/2.

Les mesures de distance les plus utilisées pour estimer la diversité β sont les distances de [START_REF] Bray | An ordination of the upland forest communities of southern wisconsin[END_REF] et de Jaccard (1901). La distance de Bray-Curtis ou indice de dissimilarité de Bray-Curtis est utilisée pour évaluer la dissimilarité entre deux échantillons donnés, en terme d'abondance de taxons (niveau OTU, phylum, famille,...) présents dans chacun de ces échantillons. Elle est définie par :

∀i ∈ 1, n , ∀j ∈ 1, n , BC i,j = p k=1 |X k i -X k j | p k=1 (X k i +X k j )
.

La dissimilarité de Jaccard traite les espèces rares et abondantes de fac ¸on égale en comparant uniquement le nombre d'espèces en commun entre les échantillons. Elle est définie par :

∀i ∈ 1, n , ∀j ∈ 1, n , JC i,j = p k=1 1 { X k i >0, X k j =0 } +1 { X k j >0, X k i =0 } p k=1 1 { X k i +X k j >0 } .
L'étude couplée de la diversité bactérienne α et β permet d'obtenir des informations sur la composition et la diversité de la communauté microbienne étudiée. Elle ne permet cependant pas de décrire complètement l'effet de certains facteurs sur le microbiote et les interactions qui s'y opèrent.

Méthodes d'apprentissage statistique

Lorsqu'un système biologique est trop complexe pour être modélisé de manière déterministe et précise, l'apprentissage statistique propose un ensemble de méthodes et algorithmes pour décrire au mieux son comportement à partir d'une série d'observations. On distingue usuellement deux types de problèmes d'apprentissage : le non supervisé et le supervisé [START_REF] Vapnik | Statistical Learning Theory[END_REF][START_REF] Bishop | Pattern Recognition and Machine Learning (Information Science and Statistics)[END_REF][START_REF] Hastie | The Elements of Statistical Learning[END_REF], décrits dans les sections suivantes.

Méthodes d'apprentissage non supervisé

On suppose que l'on a un tableau de données représenté par une matrice X = (X j i ) i=1,...,n,j=1...,p où pour tout i ∈ {1, . . . , n}, j ∈ {1, . . . , p}, -la i-ième ligne X i := (X 1 i , . . . , X p i ) désigne les observations relatives à l'individu i, -la j-ième colonne X j := t (X j 1 , . . . , X j n ) désigne les observations relatives à la variable j. En apprentissage supervisé, les observations sont données sous la forme d'une entrée X sans aucune valeur cible Y . L'objectif est alors d'apprendre un modèle capable d'extraire les régularités présentes au sein des observations pour mieux visualiser ou appréhender la structure de l'ensemble des données. Un des principaux défis de la théorie de l'apprentissage statistique est posé par la grande dimension des données (la taille n de l'échantillon négligeable devant le nombre p de variables observées). Pour pallier à ce problème, la réduction de dimension se pose alors comme une étape importante avec des méthodes du type Analyse en Composantes Principales ou encore Analyse en Coordonnées Principales présentées dans les sous-sections suivantes.

Analyse en Composantes Principales

L'Analyse en Composantes Principales (ACP) est une méthode de réduction de dimension très classique pour l'étude exploratoire de tables de données en grande dimension [START_REF] Pearson | Liii. on lines and planes of closest fit to system of points in space[END_REF]. Elle est basée sur une transformation orthogonale qui convertit une famille de variables quantitatives en de nouvelles variables décorrélées, appelées composantes principales de dimension plus petite. Pour cela, un ensemble de k vecteurs (w 1 , . . . , w k ), appelés facteurs principaux, sont tout d'abord créés et ce, de manière successive de telle sorte à restituer le maximum de variance des données originales.

Munissons R p de la norme euclidienne et du produit scalaire euclidien. Alors, les facteurs principaux (w h ) h∈{1,...,k} sont définis comme solutions des problèmes d'optimisation (1.1) suivants : ∀h ∈ {1, . . . , k} , w h = arg max w∈R p Var(Xw),

(1.1) sous la contrainte w h = 1 et les contraintes d'orthogonalité ∀l ∈ 1, h -1 , t w t h XXw l = 0. Un ensemble de k composantes principales notées (T 1 , . . . , T k ) et définies comme des combinaisons linéaires des X j centrés sont ensuite créées :

∀h ∈ {1, . . . , k} , T h = Xw h .
Notons que ces composantes principales ont la caractéristique d'être orthogonales, gardant le maximum de variance des données originales et ordonnées de la plus informative T 1 à la moins informative T k . Ces nouvelles variables définissent des plans factoriels servant de base à une représentation graphique plane des variables initiales. L'interprétation de résultats se restreint généralement aux deux premiers plans factoriels, sous réserve que ceux-ci expliquent la majeure partie de la variance du nuage des variables initiales (Figure 1.5). i , . . . , X p i ) d'un jeu de données [START_REF] Kruskal | Multidimensional scaling by optimizing goodness of fit to nonmetric hypothesis[END_REF][START_REF] Cox | Multidimensional Scaling[END_REF]. Cette méthode est notamment très utilisée lors de l'analyse de la diversité β de la communauté microbienne étudiée. Elle construit une représentation euclidienne des observations dans un espace de dimension k << p à partir d'une matrice de dissimilarités/distances que l'on notera D, de terme général D ij = d(X i , X j ) où d est une mesure de dissimilarité définie en Section 2.3. Notons alors que la PCoA prend tout son sens lorsque les observations sont inconnues ou les distances entre les individus sont non-euclidiennes.

Rappelons la définition d'une matrice euclidienne : une matrice D est euclidienne s'il existe une configuration de vecteurs (X 1 , . . . , X n ) vérifiant :

(d(X i , X j )) 2 = X i -X j , X i -X j .
On introduit ensuite la matrice A de terme général :

∀i, j ∈ 1, n , A i,j = a(X i , X j ) = - d(X i , X j ) 2 2 .
Soit H la matrice de centrage, aussi appelée matrice de projection D-orthogonale, donnée par :

H = I n -1 T n 1 n D,
où 1 n est le vecteur dont tous les coefficients sont égaux à 1 et I n est la matrice identité de taille n × n. La proposition suivante permet alors de déterminer une configuration optimale euclidienne des distances originales entre les individus :

Proposition 1 Soit D une matrice de distance et B = HAH T une matrice centrée en lignes et colonnes.

-Si D est une matrice de distance euclidienne de points {X 1 , . . . , X n }, alors B se met sous la forme B = (HX)(HX) T et est appelée matrice des produits scalaires de la configuration centrée, -Réciproquement, si B est positive de rang p, une configuration de vecteurs admettant B pour matrice de produits scalaires est obtenue en écrivant sa décomposition spectrale B = U ∆U T .

Les différentes lignes de la matrice Z = U ∆ 1/2 permettent ainsi de définir p coordonnées principales (nouvelle configuration des distances), notées z 1 , . . . , z n ∈ R k qui minimisent une fonction de coût S(z 1 , . . . , z n ), appelée stress, définie par :

S(z 1 , . . . , z n ) = i =j (d(X i , X j ) -z i -z j ) 2 .
Remarquons que dans le cas euclidien, l'ACP et la PCoA sont directement connectées :

Proposition 2 Soit X la matrice des données sur laquelle est appliquée une ACP. Alors celle-ci fournit les mêmes représentations graphiques que la PCoA calculée à partir de la matrice de distances de terme général X i -X j . Si T désigne la matrice des composantes principales, alors les coordonnées principales sont les √ nT .

Introduisons à présent les méthodes d'apprentissage supervisé pour de la sélection de variable.

Méthodes d'apprentissage supervisé

En apprentissage supervisé, les observations d'entrée sont données sous la forme de couples entrée-sortie (X, Y ), où la sortie Y , appelée valeur cible, observée sur un même d'échantillon que l'entrée X, est à prédire. Il existe deux types de sous-problèmes en apprentissage supervisé :

-Régression : lorsque la variable cible à prédire est continue, -Classement ou classification : lorsque la variable cible est discrète. On suppose que Y est une observation bruitée de p variables explicatives X 1 , . . . , X p par le biais d'une fonction f inconnue :

Y = f (X) + ε,
où ε représente le bruit ou erreur de mesure. Il s'agit alors d'approximer f en commettant le moins d'erreurs possibles sur l'ensemble d'apprentissage tout en garantissant de bonnes prédictions pour des valeurs de Y non encore observées.

Un des principaux défis des méthodes d'apprentissage supervisé est aussi posé par la grande dimension, les méthodes de régression classiques présentant des limites. Pour contourner ce problème, une première solution consiste à utiliser des méthodes de régression pénalisée, une seconde possibilité est d'utiliser des méthodes réduction de dimension.

Régression Lasso

On se place dans le cadre particulier du modèle linéaire pour lequel la fonction f est linéaire. On dispose d'un tableau de données représenté par une matrice X = (X j i ) i=1,...,n,j=1...,p , et d'une variable aléatoire réelle à expliquer

Y = t (Y 1 , . . . , Y n ). Le modèle linéaire s'écrit : ∀i ∈ {1, . . . , n} , Y i = β 0 + β 1 X 1 i + • • • + β p X p i + ε i ,
où pour tout i ∈ 1, n , les X i := (X 1 i , . . . , X p i ) sont des vecteurs indépendants et de même loi, les paramètres β 0 , . . . , β p des paramètres inconnus supposés constants à estimer et (ε i ) 1≤i≤n est une variable aléatoire centrée, simulant la présence de bruit.

En notant les vecteurs ε = t (ε 1 , . . . , ε n ), β = t (β 0 , . . . , β p ) et X = (1 n , X 1 , . . . , X p ), le modèle s'écrit sous la forme matricielle :

Y = Xβ + ε. (1.2)
Dans le cadre de la grande dimension, la méthode Lasso [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] a pour objectif de sélectionner les variables les plus pertinentes. Il s'agit d'une version régularisée des moindres carrés avec une pénalité de type 1 .

Définition 3 Soit β 1 = p j=1
|β j |. L'estimateur Lasso de β dans le modèle (1.2) est défini par :

βLasso = arg min β∈R p ( n i=1 (Y i - p j=0 X j i β j ) 2 + λ p j=1 |β j |).
où λ est un paramètre de pénalité positif permettant de calibrer le nombre de coefficients mis à 0.

Ce modèle se réécrit sous la forme :

βLasso = arg min β∈R p ( Y -Xβ 2 2 + λ β 1 ).
La méthode Lasso présente cependant des limites lorsque les variables étudiées sont fortement corrélées. Présentons alors une seconde méthode permettant à la fois de pallier au problème de la grande dimension et à la multicolinéarité des variables explicatives.

Régression des Moindres Carrés Partiels

La PLS, introduite et développée par [START_REF] Wold | Multivariate analysis[END_REF] a pour objectif de rechercher un modèle de régression linéaire sur un ensemble de composantes décorrélées construites à partir de combinaisons linéaires des p variables explicatives centrées (X 1 , . . . , X p ). La construction des composantes est optimisée pour que celles-ci soient les plus liées à la variable Y à prédire au sens de la covariance empirique. Elles vérifient : ∀h ∈ {1, . . . , k} , T h = Xw h . où les w h , appelés facteurs principaux, vérifient les problème d'optimisation (1.3) suivants : ∀h ∈ {1, . . . , k} , w h = arg max

w∈R p Cov(Xw, Y ),
(1.3) sous la contrainte w h = 1 et les contraintes d'orthogonalité ∀l ∈ 1, h -1 , t w t h XXw l = 0. Un total de k composantes orthogonales T 1 , . . . , T k , appelées variables latentes sont ainsi créées. Il suffit ensuite d'effectuer la régression de Y sur les k variables latentes centrées. Se pose alors le problème du choix optimal du nombre de composantes k à conserver. Celui-ci est traditionnellement fixé par validation croisée. Basée sur une technique d'échantillonnage, elle consiste à diviser l'échantillon de taille n en sous-échantillons, le premier dit d'apprentissage et le second dit de validation ou de test. Le modèle est ensuite bâti sur l'échantillon d'apprentissage et validé sur l'échantillon de test avec un score de performance à définir (voir Section 3.2.3).

Lorsque le nombre de variables est important, on souhaite faciliter l'interprétation des données en sélectionnant les variables pertinentes et en mettant de côté les autres. Cependant, la régression PLS agit uniquement sur la compression (regroupement de variables qui se ressemblent dans un même axe) et ne répond pas à la question de sélection des prédicteurs. Dans ce contexte, Le [START_REF] Cao | A sparse pls for variable selection when integrating omics data[END_REF] a développé une méthode variante de la PLS, appelée sparse PLS, qui consiste à introduire une pénalisation lors du calcul des facteurs principaux w h , en mettant les moins pertinents à 0. Les variables latentes ensuite créées, dépendent seulement d'un sous-ensemble des prédicteurs originaux. Cette méthode s'inspire de celle du Lasso [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF]), qui adopte elle aussi une stratégie de pénalisation lors du calcul des coefficients de la régression linéaire.

Analyse Discriminante des Moindres Carrés Partiels

La régression PLS peut facilement s'adapter au cas de la classification supervisée (PLS-Discriminant Analysis) [START_REF] Barker | Partial least squares for discrimination[END_REF]. A la différence de la PLS, cette technique de classification nécessite des groupes d'appartenance de type qualitatif à m modalités (classes) des différents objets du jeu de données. Il suffit de générer un ensemble de m variables indicatrices Z 1 , . . . , Z m puis d'appliquer la régression PLS en considérant ces nouvelles variables comme des variables quantitatives.

Il existe aussi une version parcimonieuse de la méthode PLS-DA, appelée sPLS-DA [START_REF] Cao | Sparse pls discriminant analysis: biologically relevant feature selection and graphical displays for multiclass problems[END_REF], permettant la sélection des variables prédisant au mieux les variables réponses.

Techniques de validation croisée

Quand on a la possibilité d'avoir plusieurs choix de modèles, comme c'est le cas pour la régression Lasso, il se pose le problème du choix optimal du modèle le plus pertinent. La validation croisée est une méthode qui permet d'estimer la fiabilité d'un modèle. Il existe plusieurs techniques de validation :

-Le test-set validation : l'échantillon de taille n est divisé en un échantillon d'apprentissage, qui doit représenter plus de 60% de l'échantillon total et un échantillon de test, qui permet de valider le modèle, -La validation croisée à k blocs (k-fold cross-validation) : l'échantillon original est divisé en k échantillons de taille n k , puis on sélectionne un des k échantillons comme ensemble de validation tandis que les k -1 restants constituent l'ensemble d'apprentissage. On peut alors calculer une erreur. L'opération est répétée k fois pour que chaque groupe ne serve qu'une seule fois d'échantillon de validation. Enfin, on moyenne les k erreurs afin d'obtenir une estimation de l'erreur de prédiction. La pénalité qui minimise l'erreur est considérée comme la pénalité optimale.

-La validation croisée un contre tous (Leave-One-Out Cross-Validation) : il s'agit d'un cas particulier de la validation croisée à k blocs où k = n. Cela signifie qu'à chaque itération d'apprentissage-validation, l'apprentissage se fait sur n -1 observations et la validation sur l'unique observation restante. 1.7), -Un graphe non orienté est dit connexe si chaque sommet est accessible à partir de n'importe quel autre sommet. Autrement dit, si pour tout couple de sommets distincts (s i , s j ) ∈ V 2 , il existe un chemin entre s i et s j . X 1

Introduction à la théorie des graphes

. , s p ) ∈ V p (p ∈ N * ) tel que : • u = s 0 , v = s p , • (s i-1 , s i ) ∈ E pour tout i ∈ {1, . . . ,
G (V , E ) où V ⊂ V et E = {(u, v) ∈ E, u ∈ V et v ∈ V } (Figure 1.7), -Un sous-graphe partiel de G est un graphe G (V , E ) tel que V ⊂ V et E ⊂ {(u, v) ∈ E, u ∈ V et v ∈ V } (Figure
X 2 X 3 X 4 X 5 X 1 X 2 X 3 X 5 X 1 X 2 X 3 X 5 FIGURE 1.7 -Exemple de graphe G(V, E) (gauche), graphe induit par V = {X 1 , X 2 , X 3 , X 5 } (milieu) et graphe partiel défini par V = V \{(X 2 , X 3 ), (X 2 , X 5 ), (X 5 , X 1 )} (droite)
Définition 5 Une composante connexe d'un graphe non orienté G est un sous-graphe G de G qui est connexe.

Représentation d'un graphe

Un graphe G(V, E) peut être représenté par une matrice dite d'adjacence notée A à coefficients dans {0, 1} et de dimension p × p qui décrit la présence (ou absence) d'une arête dans le graphe :

∀i, j ∈ {1, p} , A ij = 1 si (i, j) ∈ E, 0 sinon .
Lorsque le graphe est pondéré, la matrice d'adjacence est remplacée par la matrice des poids W = (w i,j ) 1≤i,j≤p à coefficients dans R + tel que

w i,j = 0 si (i, j) / ∈ E w i,j sinon.
On parle alors de matrice d'adjacence pondérée.

Analyse de réseaux

L'analyse de réseaux est confrontée à deux problématiques de recherche très actives ces dernières années. La première concerne l'inférence de réseau : comment construire à partir des données un réseau dont les arêtes représentent des " liens directs" entre les variables ? La seconde problématique intervient lorsque le réseau est déjà construit ou donné : il s'agit de l'exploration de réseaux, au travers de méthodes de clustering de graphes, qui vise à en déterminer les principales caractéristiques. L'idée principale est de se focaliser sur des interactions entre groupes de noeuds densément connectés, dont les éléments constitutifs partagent des caractéristiques communes (similarités ou connections) dans un sens prédéfini. Ces groupes de variables, appelés clusters ou communautés, sont identifiés par les méthodes de clustering à partir d'un ensemble de données dont on ne connaît pas au préalable la structure. L'analyse des réseaux offre ainsi la possibilité de mettre en évidence des ensembles de communautés, qui ont de grandes chances d'avoir la même fonction (Figure 1.8). FIGURE 1.8 -Exemple de graphe non dirigé pondéré constitué de clusters

Définitions et propriétés d'un cluster

De manière très générale, les clusters forment des sous-groupes de noeuds densément connectés. Mais, cette définition n'est pas unique puisqu'il existe de nombreuses variantes utilisées dans la littérature. Trouver des communautés dans un réseau peut être une tâche difficile sur le plan computationnel, leur nombre étant inconnu et leurs tailles et densités hétérogènes. Pour faciliter leur détection, on peut souhaiter trouver des clusters satisfaisant des propriétés spécifiques.

La notion de coreset, outil permettant de restituer l'information contenue dans graphe en sélectionnant de petits sous-ensembles représentatifs, a connu un grand succès ces dernières années. Il se base sur la décomposition en cores d'un graphe qui attribue à chaque noeud un entier reflétant sa connectivité avec les noeuds qui l'avoisinent. On parle parfois de k-core décomposition lorsqu'il s'agit de décomposer le graphe en sous graphes dans lequel tous les sommets sont connectés avec au minimum k autres sommets du sous-graphe [START_REF] Seidman | Network structure and minimum degree[END_REF]. Cet outil permet ainsi de quantifier la pertinence, centralité des noeuds qui composent le graphe, notion cruciale dans le cadre de l'analyse de réseaux en grande dimension pour réduire la complexité du problème.

Dans cette thèse, nous utiliserons ces deux définitions comme base des algorithmes de clustering que nous avons développés. Rappelons dans un premier temps, les méthodes de clustering les plus citées dans la littérature.

Méthodes de clustering

On distingue deux grandes familles de méthodes de clustering : les approches hiérarchiques et celles par partitionnement. Les approches hiérarchiques produisent une séquence ordonnée de partitions emboitées dont le premier terme est la partition la plus fine qui ne contient que des singletons, et le dernier terme est la partition la plus grossière qui ne comporte qu'une seule partie. Elles conduisent à des résultats sous forme d'arbres hiérarchiques indicés, aussi appelés dendrogrammes (voir figure 1.9). -Les méthodes descendantes considèrent l'ensemble des observations du jeu de données et procèdent par division successive jusqu'à obtenir une partition formée de singletons, -Les méthodes ascendantes commencent avec la partition la plus fine i.e. chaque noeud est un cluster à lui tout seul. Une fois les deux clusters les plus proches trouvés, ils sont agglomérés en un seul cluster. Cette étape est répétée jusqu'à ce que tous les noeuds appartiennent à un seul cluster, constitué de l'agglomération de tous les cluters initiaux. On obtient ainsi un arbre binaire dont la racine correspond à la partition à une seule partie et dont les feuilles s'identifient aux différents singletons. Les différents noeuds intermédiaires correspondent à la fusion de deux parties. Ce mode de construction est appelé Classification Ascendante Hiérarchique (CAH) [START_REF] Ward | Hierarchical grouping to optimize an objective function[END_REF]. Contrairement aux approches hiérarchiques, les approches par partitionnement permettent de subdiviser l'ensemble des observations en un certain nombre de clusters fixé a priori. Pour cela, elles emploient une stratégie d'optimisation itérative qui consiste à générer une partition initiale, puis à chercher à l'améliorer en réattribuant les données d'une classe à l'autre. L'algorithme de partitionnement le plus connu est celui des k-moyennes (k-means) [START_REF] Macqueen | Some methods for classification and analysis of multivariate observations[END_REF][START_REF] Lloyd | Least square quantization in pcm[END_REF] qui consiste à répartir les observations

(X 1 , . . . , X n ) sur K ensembles C = {C 1 , . . . , C K } (k ≤ n) formant une partition C de E = {X 1 , . . . , X n } (telle que C k ∩ C l = ∅ , ∀k = l et K k=1 C k = E) de
fac ¸on à maximiser la similarité à l'intérieur des classes et à minimiser celle entre les classes :

n i=1 K k=1 (1 C k ) i ||X i -µ k || 2 , où ∀k ∈ 1, K , i ∈ 1, n , (1 C k ) i = 1 si x i ∈ C k , 0 sinon.
et µ k correspond à la moyenne du i ème cluster C k . D'autres méthode de clustering très connues ont été spécialement développées pour le clustering de graphe. Parmi les plus connues, nous retrouvons notamment le spectral clustering [START_REF] Shi | Normalized cuts and image segmentation[END_REF][START_REF] Ng | On spectral clustering: Analysis and an algorithm[END_REF] qui, en se basant sur les propriétés spectrales du réseau, permet d'identifier des structures de communautés au sein du réseau.

Contributions de cette thèse

Les travaux présentés dans cette thèse sont regroupés en quatre chapitres. Le chapitre II concerne l'étude d'un algorithme de core clustering pour la détection robuste de variables représentatives dans le cadre de données en grande dimension. Le chapitre III propose une nouvelle variante de l'algorithme du spectral clustering qui, permet de construire des clusters dans le cadre de graphes bruités. Le chapitre IV présente une étude statistique complète réalisée sur une cohorte de patients atteints de fibrose hépatique dont le but est de découvrir des signatures bactériennes impliquées dans le développement de cette maladie. Enfin, le chapitre V explore une problématique rencontrée dans le chapitre IV qui concerne l'application de méthodes d'apprentissage statistique à des jeux de données issus de deux cohortes différentes. Ce chapitre présente plusieurs approches d'apprentissage plus juste, basées sur des méthodes standards de réduction de dimension dont le but est d'expliquer la variabilité du jeu de données tout en limitant l'effet du biais engendré par les observations.

CORE-clustering et détection de variables représentatives en grande dimension

On se place dans le cadre d'un problème de clustering de graphe en grande dimension dont l'objectif est de détecter des groupes de variables très connectées formant des structures majeures du graphe. Une très grande variété de méthodes standards de clustering utilisées dans la littérature ont été développées en ce sens. Nombreuses d'entre elles requièrent de fixer en amont le nombre total de clusters désirés qui, dans le cadre de la grande dimension, entraîne souvent des instabilités. Ces dernières années, la notion de coreset est de plus en plus utilisée pour pallier à l'effet de la grande dimension en détectant de manière robuste de petits groupes d'observations représentatives dans des jeux de données en grande dimension. Dans le chapitre II de ce manuscrit est développé un nouvel algorithme, appelé CORE-clustering, permettant de détecter de manière robuste des variables centrales représentatives du jeu de données extraites de CORE-clusters, définis avec comme seuls paramètres d'entrée le nombre minimal de variables et le niveau minimal de similarité (voir Section 2.3) entre les variables de chaque cluster. Une version parcimonieuse alternative de cet algorithme, particulièrement adaptée à la grande dimension est aussi proposée.

Considérons un système complexe composé de p variables quantitatives et n observations de ces p variables (n << p) modélisé sous la forme d'un graphe G(V, E) non dirigé, pondéré et sans boucle où V = (V 1 , . . . , V p ) est l'ensemble des noeuds du graphe et E ⊆ V × V l'ensemble des arêtes reliant chaque paire de noeuds. Notons e i,j ∈ E l'arête reliant les noeuds V i et V j avec un poids w i,j calculé à partir de la valeur absolue du coefficient de corrélation de Pearson.

Commenc ¸ons tout d'abord par définir une mesure de similarité entre deux noeuds V i et V j du graphe. Rappelons pour cela qu'un chemin de longueur Λ est défini comme une liste d'indices {d 1 , . . . , d Λ } ⊂ {1, . . . , p} tel que

V i = V d 1 , V j = V d Λ ,
w d l ,d l+1 .
(1.4) Notons à présent P i,j l'ensemble de tous les chemins possibles connectant V i et V j . La cohérence entre deux noeuds c(V i , V j ) est définie en considérant le chemin P ayant la capacité maximale parmi tous les chemins possibles P i,j :

c(N i , N j ) = max P ∈P i,j cap(P ), = max P ∈P i,j min l=1,...,Λ-1 w d l ,d l+1 .
(1.5) Notons à présent S un sous-ensemble connecté de noeuds du graphe. La cohérence au sein de ce groupe notée c(S) est alors définie comme la cohérence minimale entre variables qu'elle contient :

c(S) = min (V i ,V j )∈S 2 c(V i , V j ).
(1.6)

Cette notion de cohérence permet de mesurer le niveau de similarité au sein d'un sous-groupe de noeuds, et donc de décomposer le graphe en sous-graphes dont les variables partagent de fortes similarités. Il s'agit de l'un des paramètres d'entrées de notre algorithme.

L'algorithme de CORE-clustering avec les paramètres d'entrée ξ (niveau de similarité minimal) et τ (nombre minimal de variables) a alors pour objectif de trouver U sous-ensembles de variables S = { S u } u∈{1,..., U } , où U n'est pas fixé, en optimisant :

S, U = arg max (S,U ) U u=1 c(S u ) , ( 1

Clustering de graphes par spectral clustering pénalisé

On se place toujours dans le cadre d'un problème de clustering de graphe dont le but est d'en détecter les structures sous-jacentes puis de les identifier. Un très grand nombre d'algorithmes ont été développés en ce sens et sont communément utilisés dans la littérature comme par exemple le k-means. Les algorithmes de spectral clustering, qui se basent de manière avantageuse sur les propriétés spectrales du graphe, ont également connu du succès dès le début des années 2000. Bien que très performants d'un point de vue théorique et algorithmique, ils restent cependant sensibles à la présence de bruit dans le graphe, de par l'utilisation du k-means, pour extraire des valeurs propres de la structure par bloc sous-jacente. Dans le chapitre III de ce manuscrit est présentée une version

∀k ∈ 1, K , ∀j ∈ 1, p , (1 C k ) j =
1 si le noeud j appartient au cluster C k , 0 sinon.

Alors, les K vecteurs propres v p-K+1 , . . . , v p associés aux plus grandes valeurs propres λ p-K+1 ≤ • • • ≤ λ p de la matrice d'adjacence A, correspondent, à constante près, à ces indicateurs. Considérons dans un second temps un graphe bruité dont les propriétés sont les mêmes que le graphe précédemment défini. Nous cherchons alors à détecter de manière robuste les indicateurs des K composantes connexes. Pour cela, nous montrons l'approche sparse suivante : l'indicateur de la première composante est solution du problème de minimisation suivant

1 C 1 = arg min (1,v)∈R p W v=-w v 1 ,
où w et W sont extraits des vecteurs propres correspondant aux p -K plus petites valeurs propres de A et définis par :

(v 1 , ..., v p-k ) =       b w T W T b       .
Les indicateurs des K -1 clusters restants sont calculés de manière similaire après avoir appliqué une procédure d'orthogonalisation de Gram-Schmidt pour garantir l'orthogonalité des indicateurs. Enfin, ces problèmes de minimisation sont résolus à l'aide de la méthode Lasso.

L'efficacité et le caractère robuste de cet algorithme sont prouvés au travers de nombreuses simulations, comparaisons et applications sur jeux de données réelles.

Cet article a fait l'objet d'une soumission dans le Journal Advances in Data Analysis and Classification et du développement d'un package disponible sur le CRAN (package l1spectral).

Découverte de signatures bactériennes chez des patients fibrotiques

Le chapitre IV de cette thèse est consacré à l'application d'un ensemble de méthodes statistiques sur un jeu spécifique de données réelles. Il s'agit d'une cohorte de patients atteints à divers degrés de Fibrose hépatique, maladie du foie qui, à un stade avancé, peut avoir des conséquences graves et irréversibles tel le cancer du foie. L'objectif de cette étude est d'identifier le profil clinique des patients puis des espèces microbiologiques impliquées dans le début du développement de la maladie afin de prédire et potentiellement de bloquer l'engagement du patient vers des états pathologiques graves en intervenant sur la causalité de la maladie.

L'ensemble des données étudiées provient de la fusion de deux cohortes contenant n 1 = 36 patients roumains et n 2 = 46 patients provenant de trois pays différents. Tous sont caractérisés par une variable d'intérêt, le stade de fibrose (classé sur une échelle à trois niveaux : F 0 pas de fibrose, F 1 fibrose légère et F 2 fibrose modérée) et ont un profil assez similaire (obésité sévère). Une première étude a été réalisée sur la table des données cliniques contenant q = 21 variables caractérisant chacun des n = 82 patients de la cohorte. Puis, nous nous sommes focalisés sur la table de comptage de p = 411 OTUs (Operational Taxonomic Units) contenant le nombre de séquences observées pour chaque OTU dans chaque échantillon. Cette table a été filtrée puis normalisée en appliquant la transformation CSS avant toute application de méthodes multivariées.

L'étude statistique réalisée sur cette cohorte se décompose en quatre étapes :

-Analyse exploratoire des données cliniques et bactériologiques (ACP, analyse α et β de la diversité bactérienne) pour dresser le profil global des patients à chaque stade de fibrose, -Analyse discriminante des données bactériologiques (LDA, PLS-DA et sPLS-DA) afin de déterminer les signatures bactériologiques impliquées dans l'évolution de la maladie, -Analyse discriminante plus juste et fair clustering (fair ACP, fair l 1 -spectral clustering, fairlet clustering) pour identifier des groupes d'OTUs discriminants tout en limitant le biais introduit par la fusion des deux cohortes, -Analyse discriminante des enzymes et pathways (PLS-DA et sPL-DA) afin de prédire les voies métaboliques fonctionnelles impliquées dans le développement de la fibrose. L'efficacité des diverses méthodes appliquées sur le jeu de données est prouvée au travers de nombreuses simulations et comparaisons mises en lumière sous la forme de figures et tableaux. La combinaison de l'ensemble des méthodes précédemment citées a permis d'identifier des signatures spécifiques du début du développement de la fibrose hépatique tant au niveau clinique que bactériologique.

Cet article a fait l'objet d'une soumission dans le Journal Microbiome.

Méthode fair de Régression des Moindres Carrés Partiels

Au cours de la précédente étude statistique, nous avons identifié des problématiques liées à l'origine des données. En effet, issues de l'agrégation de deux cohortes distinctes, elles sont caractérisées par une très grande variabilité phénotypique et une diversité bactériologique entre les deux groupes d'individus, qui vient interférer dans la découverte de signatures biologiques de la maladie. La précision des méthodes d'apprentissage standards, dépendante de l'information diffusée par les données, est très sensible à l'existence de biais. Pour remédier à cette problématique, nous avons proposé dans le chapitre précédent quelques adaptations de méthodes standards d'apprentissage statistique (exploratoires, prédictives et de clustering) qui, dans le chapitre V de cette thèse, font l'objet d'une recherche plus poussée.

Rappelons que l'objectif principal de cette étude est de comprendre l'ensemble des variabilités qui caractérisent le jeu données et à terme de détecter des biomarqueurs de la maladie. Notre étude s'est alors premièrement portée sur l'Analyse en Composantes Principales, méthode exploratoire souvent utilisée comme analyse préliminaire à toute autre analyse statistique, qui permet de créer une nouvelle représentation des données reflétant le plus justement possible les variabilités qui le composent. L'analyse de données biomédicales étant souvent centrée autour d'une variable clinique d'intérêt, nous nous sommes intéressés à une méthode standard d'apprentissage supervisé : la Régression des Moindres Carrés. Très proche de l'Analyse en Composantes Principales, elle créée une nouvelle représentation des données diffusant largement l'information des données originales tout en conservant un lien avec la variable d'intérêt. Cependant, ces deux méthodes principalement basées sur la notion de variance, sont beaucoup moins efficaces en la présence de biais. Pour cette raison, nous proposons dans ce chapitre deux approches d'apprentissage fair, basées sur les deux méthodes précédemment citées, afin de pouvoir expliquer l'ensemble des variabilités qui composent le jeu de données en limitant l'effet du biais. Notons S = t (S 1 , . . . , S n ) la variable, dite sensible, à l'origine du biais, ...,n,j=1,...,p chacune des observations qui composent le jeu de données. L'idée est alors de fixer un ou plusieurs critères permettant de dissocier la variable sensible des composantes créées par ces méthodes, afin d'éviter les interférences lors de la phase de prédiction. Cette notion d'indépendance est garantie par le rapport de corrélation et contrôlée par un seuil, déterminé à l'issue d'un processus de validation croisée. Dans ce chapitre sont présentés deux versions d'algorithmes fair, fair conditional PCA et fair conditional PLS dont l'objectif est de sélectionner un nombre optimal de composantes naturellement impliquées dans le développement de la maladie sans se soucier du biais généré par les données. L'efficacité de ces algorithmes est prouvée au travers de comparaisons et applications sur deux jeux de données réelles caractérisés par une variable clinique d'intérêt de type binaire et de type continue.

Y = t (Y 1 , . . . , Y n ) la variable d'intérêt et X = (X j i ) i=1,
Ce travail est encore en cours et ouvre à de nombreuses perspectives présentées en fin de chapitre.

Chapter 2 Detection of representative variables in complex systems with interpretable rules using CORE-clusters 1 Introduction

Discovering representative variables in high dimensional and complex systems with a limited number of observations is a recurrent problem of machine learning. Heterogeneity between the variables behavior and multiple similarities between variable subsets often make this process ambiguous. A convenient strategy to solve this task is to associate each representative variable of the complex system to a cluster of variables, and to model the relations between variables in a graph [START_REF] Liu | Graph-based data clustering via multiscale community detection[END_REF]. The complex systems are indeed typically modeled as undirected weighted graphs [START_REF] Boccaletti | Complex networks: Structure and dynamics[END_REF][START_REF] Newman | Networks: An Introduction[END_REF], where the nodes (vertices) represent the variables and the edge weights are a measure of the observed similarity between the variables of the dataset. The choice of a specific clustering algorithm over the wide variety of traditional methods often depends on the nature of the data (e.g. their structure or size). Determining the granularity level of the clusters is also a common issue in high dimensional data clustering. If the clustering granularity is high, some clusters have a high similarity rate between the nodes/variables they contain, but potentially many other clusters only contain noisy relations. A large amount of selected representative variables can then be meaningless. Conversely, a low granularity leads to few large clusters with high internal heterogeneous behaviors, which makes it hard to identify the representative variables of the system in each cluster. Importantly, this issue is particularly critical when the number of observations n is lower than the observations dimension p, because of the instability related to high dimensionality and high complexity.

As in k-means clustering algorithms (MacQueen, 1967), we will use in this paper the notion of cores to address with an interpretable strategy the choice of the granularity level. Based on a distance function, the k-means algorithms indeed form a controlled number of clusters. This notion of core is also used in [START_REF] Seidman | Network structure and minimum degree[END_REF]; [START_REF] Giatsidis | Corecluster: A degeneracy based graph clustering framework[END_REF], where the graph is partitioned into a maximal group of entities, which are connected to at least k other entities in the group. The method of [START_REF] Batagelj | Fast algorithms for determining (generalized) core groups in social networks[END_REF] is also related to the notion of coreness as it hierarchically calculates the core number for each node with a complexity in the order of O(p). Strong connections also exist between the issues we address and the notion of coresets [START_REF] Agarwal | Geometric approximation via coresets[END_REF]. This notion has recently gained a significant interest in the machine learning community as it deals with finding representative observations, and not variables, in large datasets. As described in [START_REF] Claici | Wasserstein measure coresets[END_REF], it can be used in supervised learning to reduce the size of large training sets. In a similar vein, it can also be used to robustify the generalisation of the trained decision rules [START_REF] Baharan | Coresets for robust training of deep neural networks against noisy labels[END_REF], for neural network compression [START_REF] Baykal | Data-dependent coresets for compressing neural networks with applications to generalization bounds[END_REF] or for unsupervised learning [START_REF] Bachem | One-shot coresets:the case of k-clustering[END_REF]. Note that [START_REF] Bachem | One-shot coresets:the case of k-clustering[END_REF] is also particularly close to core methods, as it specifically deals with k-clustering, i.e. finding at most k cluster centers. The proposed method however does not address explicit interpretability issues.

The challenge of high-dimensionality is clearly raised in [START_REF] Saeys | A review of feature selection techniques in bioinformatics[END_REF]; [START_REF] Zhao | Advancing feature selection research[END_REF]; [START_REF] Li | Feature selection: A data perspective[END_REF], where the authors proposed different feature selection techniques with an explicit regularization in order to speed up a data mining algorithm and to improve mining performance. In this spirit, [START_REF] Wu | Iterative spectral method for alternative clustering[END_REF] recently developed an approach based on an iterative spectral optimization technique that improves the quality, computation time and scalability to high dimension of an existing alternative clustering method (Kernel Dimension Alternative Clustering). In [START_REF] Chen | Crowdclustering with partitions labels[END_REF], the authors also used multinomial regression model to learn automatically the number of clusters, and then to limit strong assumptions required by the model in high dimension. Note that [START_REF] Yu | Feature selection for high-dimensional data: A fast correlation-based filter solution[END_REF] also defined a strategy for the detection of representative variables in high-dimensional data but did not explored a regularization strategy when the number of observations is much lower than the problem dimension. Those strategies require as well to make a decision about the number of clusters to determine. Motivated by the above issues, we propose a new formalism to robustly and intuitively estimate the representative variables of complex systems. This is first made through a graph clustering strategy for which the clusters do not necessarily cover all nodes/variables. This clustering strategy specifically estimates CORE-clusters, which are connected subsets of variables having (1) a minimum number of nodes/variables, and (2) a minimal similarity level between all their variables. The representative variables are then those having the lowest average distance to all other variables in each CORE-cluster. The detection of representative variables is therefore regularized using a control on the minimal COREcluster size and not the number of representative variables to be detected, or a LASSO derived penalty term. This point of view has been originally considered in [START_REF] Brunet | Gender shades: Intersectional accuracy disparities in commercial gender classification[END_REF][START_REF] Brunet | A new gene co-expression network analysis based on core structure detection (csd)[END_REF]. We present here a totally reformulated version of this initial idea, which makes fully interpretable the selection of the representative variables by introducing the notion of CORE-clusters. Fine algorithmic improvements, described in this paper, also make the original clustering algorithm more efficient. A greedy version of this original algorithm, which turns out to scale particularly well to high dimensional data, is additionally proposed. New results on synthetic data as well as comparisons with other methods now shed light on how the proposed strategy is robust and explainable. Finally, we now demonstrate the validity of our formalism on two high dimensional datasets representing the expression of genes and a road network.

Our methodology is described in Sections 2 and 3 and is then tested both on simulated and real data in Section 4.

Statistical Methodology

Graph-based representation of the observations

Let us consider a complex system of p quantitative variables X = (X 1 , • • • , X p ) and n observations (X j 1 , • • • , X j n ), j ∈ {1, • • • , p} of these p variables. The driving motivation of our work is to detect representative variables out of X when n p. As mentioned in Section 1, the detection of these representative variables is regularized using a graph-based approach. We then model the relations between the variables using an undirected weighted graph G(V, E), where

V = (V 1 , • • • , V p )
is the nodes set corresponding to the p variables, and E is the edges set. We also denote e i,j ∈ E the edge joining the nodes V i and V j with weight w i,j .

In order to handle the properties of application-specific similarity measures that can be encoded in the edge weights w i,j , we will consider in the remainder of the paper that all w i,j ≥ 0 and that the higher w i,j the closer the observed behavior of the variables X i and X j . The weights therefore represent a notion of similarity between the variables X i and X j . For instance, if the empirical correlations cor(X i , X j ) are measured between the pairs of variables X i and X j , with (i, j) ∈ {1, . . . , p} 2 , then w i,j = |cor(X i , X j )| can reasonably be used.

Coherence of a variable set

To define a notion of distance between two variables X i and X j , which are not directly connected in the graph, we use the notion of capacity [START_REF] Pollack | The maximum capacity through a network[END_REF][START_REF] Hu | The maximum capacity route problem[END_REF] of a path P between the corresponding nodes

V i and V j in G(V, E). A path P of a graph G from X i to X j of length Λ is a list of indices {d 1 , . . . , d Λ } ⊂ {1, . . . , p} such that X i = X d 1 , X j = X d Λ ,
and w d l ,d l+1 is known and is not equal to 0, for all l = 1, . . . , Λ -1. The capacity cap(P ) of path P is then the minimal weight of its edges, i.e. cap(P ) = min l=1,...,Λ-1

w d l ,d l+1 .
(2.1)

We also denote by P i,j the set of all possible paths connecting X i to X j . The coherence c(X i , X j ) between X i and X j is then defined by considering the path P having the maximum capacity among the paths of P i,j [START_REF] Hu | The maximum capacity route problem[END_REF]

, i.e. c(X i , X j ) = max P ∈P i,j cap(P ), = max P ∈P i,j min l=1,...,Λ-1 w d l ,d l+1 .
(2.2)

If the weight w i,j is known, its is interesting to remark that the coherence c(X i , X j ) is not necessarily equal to its value. For instance, both X i and X j may be very similar to a third variable X k but not similar to each other. From a computational point of view, the similarity in w i,j may also be unknown if the edge e i,j is not stored in a sparsified version of the complete graph. Since n p, pertinent relations may finally be lost in w i,j but recovered in c(X i , X j ) thanks to other relations that would be captured. We believe that these points are particularly important to define coherent variable sets in the complex data case.

We now denote by S a connected subset of the variable set X. The coherence c(S) of this variable subset is the minimal coherence between the variables it contains:

c(S) = min (X i ,X j )∈S 2 c(X i , X j ). (2.3)
If all the variables of S have a coherent observed behavior, then c(S) is high. The use of this notion on synthetic data is illustrated in Appendix 6.1. The coherence of a subset measures the strength of the variables it contains. The more coherent S, the more sense it makes to consider that its variables share common features measured by the chosen similarity. Decomposing the graph into maximal groups sharing a strong similarity, i.e. finding all the groups of variables with a large enough coherence is the core of the following subsections on CORE-clusters selection.

CORE-clusters

We recall that the goal of our formalism is to detect the representative variables of complex systems. In our formalism, each representative variable is extracted out of a CORE-cluster, defined as: Definition 1 A CORE-cluster S ξ,τ ⊂ X is a connected variable subset having a size higher than τ and a coherence c(S ξ,τ ) higher than a threshold ξ.

The parameters τ and ξ ensure that each representative variable has a non-negligible coherence with at least τ -1 other variables, which directly regularizes its selection: Large values of τ indeed lead to the detection of large sets of coherent variables. In that case, the representative variables are likely to be meaningful even if n < p. If τ is too large, each CORE-cluster may however contain several variables that would have been ideally representative. On the contrary, too small values of τ are likely to detect all meaningful representative variables but also a non-negligible number of false positive representative variables, especially if the observations are noisy or if n < p. A good trade-off, which depends on n, p and the level of noise in the observations has then to be found when tuning τ .

CORE-clustering

CORE-clustering consists in estimating an optimal set of CORE-clusters, so that the representative variables they contain explain as much information as possible in the observed complex system. We use the following definition:

Definition 2 CORE-clustering with parameters ξ and τ consists in finding U variable subsets S = { S u } u∈{1,..., U } , where U is not fixed, by optimizing:

S, U = arg max (S,U ) U u=1 c(S u ) , (2.4)
under the two constraints:

1. All S u are connected variable sets having a size higher than τ and a coherence c(S u ξ,τ ) > ξ. They therefore correspond to CORE-clusters and can be denoted S u ξ,τ . 2. There is no overlap between the clusters, i.e.

S u 1 ∩ S u 2 = ∅ for all (u 1 , u 2 ) ∈ {1, . . . , U } 2 .
It may first seem that U should be as high as possible, so that the union of all S u ξ,τ contains all the variables of X. Each CORE-cluster S u ξ,τ must however have a coherence higher than the threshold ξ. As illustrated in Appendix 6.1, the variables of X which are not coherent with at least τ other variables should ideally not be contained in any CORE-cluster as they would make drop their coherence. The CORE-clusters in S should then only contain pertinent variables so the optimal value for U is implicitly defined during the CORE-clustering procedure.

It is also important to remark that the potential number of subsets of X to find good CORE-cluster candidates is huge, even for moderate values of p. Moreover computing Eq. (2.4) is particularly demanding. The two optimization algorithms of Section 3 are then aggregative and divisive algorithms designed to optimize Eq. (2.4) without explicitly computing it.

Representative variables selection

We now present how each representative variable is extracted out of a CORE-cluster S ξ,τ . As mentioned Section 2.2 the pertinent relations between two variables X i and X j may be lost in w i,j since n p and recovered by their coherence c(X i , X j ) using other relations. In this example, the similarity s captures true positive and false negative relations and the notion of coherence makes robust the detection of relations. For the same reasons, it may however also capture false positive relations, so the CORE-clusters may contain undesirable variables. The variables captured by CORE-clusters using false positive relations should however be located at the cluster boundaries if the data are not too noisy. False positive connections are indeed less common and in average weaker than true positive connections in this case.

In order to limit the impact of the false positive relations, we then define the representative variables as the CORE-cluster centers. More specifically, each representative variable minimizes an average distance with the other variables of a CORE-cluster S ξ,τ . Instead of using distances based on the maximum capacity Eq. ( 2.2), we use a more standard notion of distance calculated as sums of weighted edges traversed by the optimal paths. This limits the phenomenon of having multiple variables with the same optimal distance due to the min-max strategy. The graph weights w i,j , which represent a similarity level, must however be converted into distances, which can be simply done by using d i,j = 1/(w i,j + ), where > 0. The representative variable of a CORE-cluster S ξ,τ is then the one that has the lowest average distance to the other variables of S ξ,τ . The impact of selecting the representative variables as CORE-cluster centers is discussed in Appendix 6.1.

General guidelines for the choice of ξ and τ

The selection of the representative variables directly depends on the parameters ξ and τ . As explained Section 2.3, a CORE-cluster S ξ,τ is indeed a connected variable subset having a size higher than τ and a coherence c(S ξ,τ ) higher than a threshold ξ. In order to estimate pertinent representative variables, we recommend to use the following guidelines: (1) First compute how the edge weights w i,j of the graph G(V, E) are distributed. The coherence c(S ξ,τ ) represents the weakest connection between the variables of S ξ,τ , so the threshold ξ should be relatively high regarding to the different values of w i,j . A value of ξ equal to the 80th percentile of the edge weights w i,j appears as reasonable.

(2) Choosing a suitable minimal amount of variables τ in S ξ,τ is more subtle as this choice both depends on the complexity of the relations expressed in G(V, E), and how the number of observations n is low compared with the observations dimension p. In all generality, tuning τ as equal to p/10 is reasonable as a first guess. (3) If no CORE-clusters are found with the initial parameters, the user may try to run again the CORE-clustering procedure with lower parameters ξ and τ .

From our experience, we recommend in all cases not using values of ξ lower than the 40th percentile of the edge weights or values of τ lower than 10. The CORE-clusters would be likely to contain variables with strongly heterogeneous behaviors or false positive connections in these cases. Note finally that when several connected CORE-clusters are identified with a given parametrisation, it is interesting to test whether stronger CORE-clusters would be locally found by using higher the values of ξ or τ . This is illustrated in Section 4.3.

Computational Methodology

Main interactions estimation

The very first step of our strategy is to quantify the similarity between the different observed variables. The similarity is first computed using the absolute value of Pearson's correlation and represented as a dense graph G(V, E), where V contains p nodes, each of them representing one of the observed variables, and E contains K E = p(p -1)/2 undirected weighted edges that model a similarity level between all pairs of variables. In this approach, the variables with no connection are associated with a correlation coefficient of zero. The algorithmic cost of this estimation can be O(n 2 p) but it can also be easily parallelized using divide and conquer algorithms for reasonably large datasets, as in [START_REF] Randall | Cilk: Efficient Multithreaded Computing[END_REF]. For large to very large datasets, correlations should be computed on sparse matrices, using e.g. [START_REF] Cysouw | R function cor.sparse[END_REF] to make this task computationally tractable.

CORE-clustering algorithms

Inspired by [START_REF] Hu | The maximum capacity route problem[END_REF] who solved the maximum capacity problem of [START_REF] Pollack | The maximum capacity through a network[END_REF] using optimal spanning tree, we estimate the CORE-clusters on optimal spanning trees. A spanning tree G(V, T ) is a subgraph of G(V, E) with no cycle and T ⊂ E. The maximum spanning tree of G is then the spanning tree of G having the maximal sum of edge weights. Using the maximum spanning tree to detect the CORE-clusters strongly limits the potential amount of node combinations to test while preserving the graph edges that are likely to be good candidates for the optimal paths of Eq. (2.2). Conversely, it is straightforward to show that the coherence of a variable subset in G(V, T ) is lower or equal to the coherence of the same variable subset in G(V, E). The edges of T are indeed a subset of E, so Eq. (2.2) between two variables X i and X j is lower or equal on T than on E. The CORE-clusters computed in G(V, T ) are therefore eligible CORE-clusters on G(V, E). By discussing the algorithmic cost of our algorithms, we will make clear that this reasonable domain reduction makes our problem scalable to large datasets. The impact of using maximal spanning trees on the measure of a cluster coherence is also discussed in Appendix 6.1 on synthetic examples.

Maximum Spanning Tree

Algorithm 1 Maximum spanning tree algorithm Require:

Graph G(V, E) with nodes V i , i ∈ 1, • • • , p and edges E k , k ∈ 1, • • • , K E . Require: Weight of edge E k is W (E k ).
1: Sort the edges by decreasing weights, so

W (E 1 ) ≥ W (E 2 ) ≥ • • • ≥ W (E K E ). 2: Assign label L(V i ) = i to each node V (i). 3: Initiate an edge list T as void. 4: for k = 1 : K E do 5:
We denote V i and V j the nodes linked by edge E k .

6:

if L(V i )! = L(V j ) then 7:
Add edge E k to the list T 8:

Propagate the label L(V i ) to the nodes that have label L(V j ).

9:

end if 10: end for 11: return Graph with a tree structure G(V, T ).

The maximum spanning tree of G is the simple and reliable modeling of the relationship between the graph nodes (only p -1 links). One of the most famous algorithm developed to find such trees is called Kruskal's algorithm [START_REF] Kruskal | On the shortest spanning subtree of a graph and the traveling salesman problem[END_REF]. The maximum spanning tree is built by adding step by step partial associations so that there will be no cycle in the partial graph.

We denote by G(V, T ) the resulting graph, where T has a tree-like structure. Details of the algorithm are given in Alg. 1. The algorithmic cost of the sort procedure (row 1) is O(K E log (K E )). Then, the for loop (rows 4 to 10) only scans one time the edges. In this for loop, the most demanding procedure is the propagation of label L(V i ), row 8. Fortunately, the nodes on which the labels are propagated are related to V j in G(V, E). We can then use a Depth-First Search algorithm [START_REF] Tarjan | Depth first search and linear graph algorithms[END_REF] for that task, making the average performance of the for loop O(K E log (p)).

CORE-clustering algorithm Algorithm 2 CORE-clustering algorithm

Require: Graph G(V, T ) with nodes V i , i ∈ 1, • • • , p; and edges T k , k ∈ 1, • • • , K T . Require: Weight of edge T k is W (T k ).
Require: Granularity coefficient τ and threshold ξ.

1: {Initiate the algorithm} 2: Sort the edges by increasing weights;

3: Assign label L(V i ) = i to each node V (i) and set CORElabel = -1. 4: {CORE-clusters detection} 5: for k = 1 : K T do 6:
We denote V i and V j the nodes linked by edge E k .

7:

if L(V i )! = L(V j ) then 8:
Propagate the label L(V i ) to the nodes that have label L(V j ).

9:

if number of nodes with label

L(V i ) ∈ {τ, • • • , 2τ -1} then 10:
Label CORElabel is given to the nodes with label L(V i )

11:

CORElabel = CORElabel -1 12: end if 13:
end if 14: end for 15: {Post-treatment of the labels} 16:

for v = 1 : V do 17: If L(V v ) > 0 then L(V v ) = 0 else L(V v ) = -L(V v ) 18: end for 19: {Filter the CORE-clusters S u ξ,τ s.t. c(S u ξ,τ ) < ξ} 20: for u = 1 : U do 21: If c(S u ξ,τ ) < ξ then Set L(V v ) = 0 to the nodes V v of S u ξ,τ . 22: end for 23: return Labels L.
In contrast with other clustering techniques, this CORE-clustering approach detects clusters having an explicitly controlled granularity level, and only gathers nodes/variables with a maximal path capacity. CORE-clusters are detected from the maximal spanning tree G(V, T ) by gathering iteratively its nodes V in an order that depends on the edge weights in T (increasing weight order). A detected CORE-cluster has a size higher than τ , where τ is the parameter that controls the granularity level. Thus, Alg. 2 first generates many small and meaningless clusters and parameter τ should not be too small to avoid considering these clusters as CORE-clusters. Then, the first pertinent clusters will be established using the edges of T with larger weights, leading to pertinent node groups. Finally, the largest weights of T are treated in the end of Alg. 2 in order to split into several CORE-clusters the nodes related to the most influential nodes/variables.

It worth mentioning that the Rows 20 to 25 of Alg. 2 are the only ones that require to compute the coherence c of the estimated CORE-clusters. Computing a coherence Eq. ( 2.3) is indeed demanding, so it is considered here as a post-treatment limited to pre-computed CORE-clusters. In practice it is also performed on the maximal spanning tree G(V, T ) and not on the whole graph G(V, E).

Remark too that the algorithmic structures of Alg. 1 and Alg. 2 are similar. However Alg. 2 runs on G(V, T ) and not on G(V, E). The number of edges K T in G(V, T ) is much lower than K E since T has a tree structure and not a complete graph structure. It should indeed be slightly higher than p [START_REF] Steele | Minimal spanning trees for graphs with random edge lengths[END_REF], which is much lower than K E = p(p -1)/2. Moreover the propagation algorithm (rows 9 and 11) will then never propagate labels on more than 2τ -1 nodes. The algorithmic cost of the sort procedure (row 1) is then O(K T log (K T )) and the average performance of the for loop O(K T log (τ )).

A greedy alternative for CORE-clusters detection

We propose an alternative strategy to the CORE-clusters detection algorithm: The edge treatment queue may be ordered by following decreasing edge weights instead of increasing edge weights. The nearest edges are then first gathered making coherent CORE-clusters as in Alg. 2, although one CORE-cluster may contain several representative variables. To avoid gathering noisy information, the for loop on the edges (row 6 of Alg. 2) should also stop before meaningless edges are treated. This strategy has a key interest: It can strongly reduce the computational time dedicated to Algs. 1 and 2. By doing so, Alg. 1 and modified Alg. 2 are purely equivalent to Alg. 3.

Again, the structure of this algorithm is very similar to the structure of the Maximum Spanning Tree strategy in Section 3.1. The algorithmic cost of the sort procedure (row 1) is O(K E log (K E )). Then, the loop rows 5 to 14 scans γ edges, where γ is the number of edges having a weight higher than ξ. In most cases, γ should be much lower than K E , which strongly limits the computational impact of this loop. Labels propagation in this loop (rows 8 and 10) is also limited to 2τ -1 nodes. The average performance of the for loop is therefore O(γ log (τ )).

Central variables selection in CORE-clusters

Once a CORE-cluster is identified, we use a straightforward strategy to select its central variable: the distance between all pairs of variables in each CORE-cluster is computed using a Dijkstra's algorithm [START_REF] Cormen | Introduction to Algorithms[END_REF][START_REF] Zan | Shortest path algorithms: An evaluation using real road networks[END_REF] in G(V, E). The central variable is then the one that has the highest average distance to all other connected variables in the CORE-cluster. As the computed CORE-clusters have less than 2τ nodes, the algorithmic cost of this procedure is O(τ 2 ) times the number of detected CORE-clusters, which should remain low even for large datasets.

Algorithm 3 Greedy CORE-clustering algorithm

Require: Graph G(V, E) with nodes V i , i ∈ 1, • • • , p and edges E k , k ∈ 1, • • • , K E . Require: Weight of edge E k is W (E k ).
Require: Granularity coefficient τ and threshold ξ.

1: Sort the edges by decreasing weights.

2: Define the number of edges γ having a weight higher than ξ.

3: Assign label L(V i ) = i to each node V (i). 4: Set CORElabel = -1. 5: for k = 1 : γ do 6:
We denote V i and V j the nodes linked by edge E k .

7:

if L(V i )! = L(V j ) then 8:
Propagate the label L(V i ) to the nodes that have label L(V j ).

9:

if number of nodes with label

L(V i ) ∈ {τ, • • • , 2τ -1} then 10:
Label CORElabel is given to the nodes with label L(V i )

11:

CORElabel = CORElabel -1 12: end if 13: end if 14: end for 15: for v = 1 : V do 16: If L(V v ) > 0 then L(V v ) = 0 else L(V v ) = -L(V v ) 17: end for 18: return Labels L.

Results

Core clustering of simulated networks

In this section, we compare our CORE-clustering algorithms with other standard methods on simulated scale-free networks. Such complex networks are indeed common in the data-science literature. They contain a little amount of highly connected nodes (hubs), that we will assimilate to representative variables, and many poorly connected nodes.

Experimental protocol

To generate the synthetic networks, we first simulated the profile (observations) of representative variables and then the profile of remaining variables around these hubs. We then considered

K different clusters of size p C 1 , p C 2 , . . ., p C K . A simulated expression data set X ∈ R n×p is then composed of p = p C1 + . . . + p C K variables.
In the cluster k, the observations are then simulated as follows: (1) Generate the observations x (1,k)

= (x (1,k) 1 , . . . , x (1,k) n
) of a representative variable using a normal distribution N (0, 1). (2) Choose a minimum correlation r min and a maximum correlation r max between the representative variable and the other variables of the predefined cluster. In this paper, we always used r max = 1 and r min = 0.5. (3) Generate the profiles x (j,k) , with j ∈ {2, . . . , p C k }, such that the correlation of the j-th profile with the profile of x (1,k) is close to

r j = r min + 1 -j p C (r max -r min ). For i ∈ {1, . . . , n}, we then use x (j) i = x (rep) i + (r -2 j -1) 1 2 (j)
i , where Note finally that an amount of 30 networks of each type was generated to assess the stability of our methodology. This will indeed make it possible to draw in Figure 2.1 the box-plots of the clustering quality for each type of network.

Measure of clustering quality

There exists various criteria to measure a clustering quality. External indices such as impurity and Gini indices measure the extent to which the clusters match externally supplied class labels. Internal indices like the modularity and the intra-cluster to inter-cluster distance ratio are also used to measure the quality of a clustering structure without any external information. Such criteria are however not suitable here, as we do not clusterize all variables but rather extract core structures that emphasize representative variables. Thus, we propose the following criterion for simulated data for which the block structure of the similarity matrix is known: Let X i , i ∈ {1, • • • , p} be the variables and C j (j ∈ [1, K]) be the ground-truth CORE-clusters of variables, typically on synthetic data. As in Section 2, we also denote S u ξ,τ (u ∈ [1, U ]) the CORE-clusters predicted by our algorithm. In order to evaluate the quality of the prediction, we compute a score R defined as:

∀u ∈ {1, • • • , U } , R u = max j∈[1,K] Card(X i ∈ C j ∩ S u ξ,τ ) , (2.5)
where

1 ≤ i ≤ p and R = 1 p U u=1 R u .
To compute this score, each R u is equal to 0 if there is no overlap between C j and any S u ξ,τ , and is equal to the number of variables in C j if a S u ξ,τ contains all the variables of C j . A score R equals to 1 then means that a perfectly accurate estimation of the C j was reached, and the closer to 0 its values the less accurate the CORE-clusters detection.

Results

We compared the standard and greedy CORE-clustering algorithms (Algs. 2 and 3) on these simulated datasets with two other graph-based clustering algorithms: the spectral clustering [START_REF] Shi | Normalized cuts and image segmentation[END_REF][START_REF] Ng | On spectral clustering: Analysis and an algorithm[END_REF], available in the R-package anocva, and Louvain method for community detection [START_REF] Blondel | Fast unfolding of communities in large networks[END_REF], available in the R-package igraph. Note that Louvain method requires as input parameter a graph modeling the dataset (the correlation matrix is transformed upstream into a graph) but not the final number of clusters. Boxplots of the computed scores R (see Eq. (2.5)) are shown in Figure 2.1. Note that in each boxplot, the dots represent the outlier scores which are either lower than q 0.25 -1.5(q 0.75 -q 0.25 ) or higher than q 0.75 + 1.5(q 0.75 -q 0.25 ), where q 0.25 and q 0.75 are the first and third quartiles of the scores, respectively.

The subplots of Figure 2.1(a-b) show that Alg. 2 is more robust than Alg. 3, and gives slightly better results than spectral clustering and Louvain method, when the level of noise is high. The same applies when the sample size decreases in the subplots of Figure 2.1(c).

Application to real biological data

We now present the results obtained on the classic Yeast dataset1 [START_REF] Spellman | Comprehensive identification of cell-cycle-regulated genes of yeast saccharomyces cerevisiae by microarray hybridization[END_REF]. With a total of about 1.3 × 10 6 weighted edges considered when representing the variables correlations in the graph G(V, E), the CORE-clustering procedure required about 160 and 3 seconds with Algorithms 2 and 3, respectively.

Yeast dataset

The well-known synchronized yeast cell cycle data set of [START_REF] Spellman | Comprehensive identification of cell-cycle-regulated genes of yeast saccharomyces cerevisiae by microarray hybridization[END_REF], includes 77 samples under various time during the cell cycle and a total of 6179 genes, of which 1660 genes are retained for this analysis after pre-processing and filtering. The goal of this analysis is then to detect CORE-clusters among the correlation patterns in the time series of yeast gene expressions measured along the cell cycle. Using this dataset, a measure of similarity between all gene pairs was measured with the absolute value of Pearson's correlation. A total of about 1.3 × 10 6 weighted edges are then considered when representing the variables correlations in the graph G(V, E).

Comparison of the two CORE-clustering algorithms

In order to compare the two proposed CORE-clustering algorithms described (Alg. 2 and Alg. 3) we tested them on the yeast dataset with τ = 30 and ξ = 0.75. We indeed empirically considered that CORE-clusters containing 30 to 59 variables are reasonable to regularize the problem, and that 0.75 is a threshold above which the absolute value of the correlation between two variables reasonably shows The computations were much faster using the greedy algorithm Alg. 3. It indeed required about 3 seconds. An amount of 11 CORE-clusters was found. To interpret this result, we computed the coherence c (see Eq. 2.2) between the 4 representative variables obtained using Alg. 2 and the 11 ones obtained using Alg. 3. Interestingly, Alg. 3 selected YLL026W=RV3 and YDL003W=RV5. Other variables very close to RV1, RV5 and RV6 (with c > 0.83 i.e. higher that the highest c within the CORE-clusters) were also selected: YHR219W, YLR103C, YLR276C and YLR196W. Results equal or close to those obtained with Alg. 2 were then obtained. The representative variables YDR418W, YML119W, YJL038C, YNL283C and YGR167W were additionally found. In this experiment, Alg. 3 therefore selected more representative variables and has then naturally a larger score (see Eq. (2.4)) than by using Alg. 2. However, it also obviously captured different representative variables that would be gathered in the same CORE-cluster using Alg. 2. The two algorithms have therefore slightly different properties but lead to coherent results.

Impact of the number of observations

In order to evaluate the stability of the results with respect to the number of observations, we tested again Alg. 2 with the same parameters, but by using only the 30 first observations of the yeast dataset out of the 77 observations. Interestingly, YDL003W=RV5 was selected and YML093W which is very close to RV6 (c > 0.86) was also selected. The two other representative variables found, YER190W and YLL026W, were however not similar to RV1 or RV3. The information lost in the 47 observation we removed therefore did not allowed to recover the influence of RV1 and RV3 on the complex system but the 30 remaining observations contained a sufficient amount of information to detect RV5 and RV6 as influent variables. This suggests that the strategy detects stable representative variables, even when the number of observations is very low compared with the dimension of the observations.

Comparison with spectral-clustering

In order to compare our CORE-clustering strategy with a standard clustering approach, we finally estimated representative variables in the Yeast dataset as the center of clusters estimated using spectral clustering. The standard version of the spectral clustering available in R was used. Its main parameter is the number of seeds η used in the k-means part of the spectral clustering. When using η seeds, an amount of 1 to η clusters (with more than one variable) are estimated using k-means and all variables are contained in a cluster. In order to fairly compare the spectral clustering and the CORE-clustering approaches, we then clusterized the variables of the Yeast dataset using η = {5, 30, 50, 70, 110}. Note that we only tested an η higher than 77, due to the fact that n = 77 observations are known. We tested η = 110 to evaluate the spectral clustering behaviour with η slightly higher than n.

An amount of {3, 3, 6} large clusters were obtained for η = {50, 70, 110}, respectively. For η = {5, 30} only a single cluster gathering almost all variables was also obtained. The average coherence of the estimated clusters was 0.44 for a standard deviation of 0.06. The highest coherence was 0.63 which is clearly lower than the considered threshold of ξ = 0.75 we used with the COREclustering. This makes ambiguous the interpretation of the role played by the representative variables.

Finally, it is worth mentioning that all representative variables (genes) obtained using Alg. 2 have a known physiological function and only two variables out of eleven have an unknown function by using Alg. 3. For the spectral clustering with η = 110 four selected variables out of six have known functions. For η = 70, three variables with unknown functions were selected. For η = 50, two variables out of three with known functions were selected. For η = 5 and η = 30 the single selected variable was the center of all variables with respect to the center definition in Section 3.3, and turns out to have a known function. It therefore appears in this experiment that the representative variables obtained using CORE-clustering are more interpretable than those estimated using spectral clustering.

Application to the U.S. road network

We now assess the CORE-clustering algorithms on the U.S. road network out of the 9th DIMACS Implementation challenge 2 . Our goal here is to discuss the pertinence of the detected CORE-clusters, and not specifically their representative variables, on a large scale and straightforwardly interpretable dataset. The graph contains here 2.4 × 10 7 nodes, each of them representing a crossing of the U.S. road/streets network, and 5.8×10 7 arcs representing the part of the roads/streets between two crossings. Interestingly the distance between the crossings is also associated to each edge of the graph.

To assess the CORE-clustering algorithms on the U.S. road network, we first transformed the distances between the crossings into weights that are higher and higher for increasingly closer crossings. A weight max ((4 × 10 3 -l)/(4 × 10 3 ), 10 -4 ) was specifically given to each edge, where l is its original length in feet.

We show Figure 2.3, the CORE-clusters 3 obtained using the greedy algorithm with τ = 5 × 10 4 and ξ = 10 -3 . This means that each CORE-cluster contains a road network of 5 × 10 4 to 10 5 crossings for which one can travel from any crossing to another one by only using streets of less than 3996 feet (about 1.2 km) between two crossing. As expected, the clusters represented in Figure 2.3 correspond to the 18 largest urban areas in the U.S. Note that two of them are made of three CORE-clusters (New-York City and Los-Angeles), and two other ones (Miami and Chicago) are made of two COREclusters. This is due to an algorithmic choice, discussed in Sections 3.2.2 and 3.2.3, which leads to the detection of CORE-clusters containing τ to 2τ -1 nodes by using the proposed CORE-clustering algorithms. As discussed in Section 2.6, an interesting strategy if connected CORE-clusters are found consists in running again the CORE-clustering algorithms with higher values of τ or ξ, to potentially detect more robust CORE-clusters : By reproducing this test with τ = 10 5 instead of τ = 5 × 10 4 , we now try to find CORE-clusters containing 10 5 to 2 × 10 5 crossings, only 5 urban areas are found: New-York City (2 CORE-clusters), Los-Angeles, Chicago, Miami and San-Fransisco. Now, by using τ = 5 × 10 4 and now ξ = 0.5, i.e. with distances between the crossings lower than about 2000 feet (about 0.6 km) instead of 3996 feet, only 3 urban areas are found: New-York City (2 CORE-clusters), Los-Angeles and Chicago.

What is interesting here in terms of interpretability is that we have been able to select specific subparts of the U.S. road network by explicitly controling the size of the clusters or a specific level of density in the network. This notion of interpretability can be further discussed by observing the CORE-clusters obtained in the New-York city urban area, as shown Figure 2. 3-(b-c-d). By using τ = 5 × 10 4 and ξ = 10 -3 , the three CORE-clusters split the densest parts of the New York City urban area into the New-Jersey state, Long-Island and the rest of New-York state. When having larger CORE-clusters, i.e. with τ = 10 5 , the New-Jersey cluster expands and the densest parts of the two New-York state clusters are merged. Now, when enforcing instead a stronger coherence inside of the clusters, i.e. with τ = 5 × 10 4 and ξ = 0.5, the New-Jersey and Long-Island CORE-clusters remain stable, but the Manhattan/mainland New York state cluster is not large and dense enough, so it is not captured as a CORE-cluster. Note that each CORE-cluster estimation required about 50 seconds and 2.4GB here without any parallelization.

Conclusion

Although complex systems in high dimensional spaces with a limited number of observations are quite common across many fields, efficient methods to treat the associated problem of graph clustering 3. Clusters represented using MI Map Tools: GeoPlotter: https://mobisoftinfotech.com/tools/ plot-multiple-points-on-map/ is an ambiguous task. Some of these techniques, based on assumptions in view of controlling the variables contribution to the global clustering, often do not allow to select the best graph partition.

In reply to this issue, we developed a formalism based on an original graph clustering strategy with specific properties. This formalism makes it possible robustly identify groups of representative variables of the studied system by tuning two intuitive parameters: (1) the minimum number of variables in each CORE-cluster, and (2) a minimum level of similarity between all the variables of a CORE-cluster. Its effectiveness was further satisfactorily assessed on simulated data and on real datasets. From a methodological perspective, an interesting research direction would be to mix Algorithms 2 and 3 into a single hybrid top-down and bottom-up optimization scheme. Our goal would be to scale well to very high dimensional datasets, as when using Algorithm 3, while being as robust as Algorithm 2 when potential CORE-clusters are coarsely identified. When the CORE-clusters found by either Algorithm 2 or 3 are very large, a stochastic strategy could also make faster the detection of their representative variables. Although this secondary part of our methodology has been addressed by using a standard Dijkstra's algorithm (see Section 3.3), it could be addressed by using an extension of [START_REF] Gadat | How to calculate the barycenter of a weighted graph[END_REF] where the optimal paths between all variables would not be pre-computed prior to the central variable detection. Application of our formalism in various fields such as gene regulatory networks, social networks or recommender systems would also be of interest. Our formalism is indeed sufficiently flexible to incorporate different types of similarity measures between the observed variables.

Note finally that our formalism was implemented in C++ and wrapped in a R package, which is freely available on sourceforge4 , so these results can be easily reproduced.

Appendices

Pertinence of the representative variable detection model

We illustrate in this section the notion of coherence defined in Section 2.2. We simulated an adjacency matrix that mimics the absolute values of the Pearson correlations between 15 variables. This symmetric matrix is shown Fig. 2.4(top-left) and represents graph G 1 following the method of Section 2.1. Each of its values is then equal to a similarity level encoded in w i,j between two variables X i and X j , where i and j are in {0, • • • , 14}. We can clearly remark that it contains two blocks of related variables

S 1 Ref = {X 0 , • • • , X 6 } and S 2 Ref = {X 8 , • • • , X 14
} and an independent variable X 7 . In addition, variables X 3 and X 11 are slightly more related to other variables in S 1

Ref and S 2

Ref ,respectively. They can then be considered as the most pertinent representative variables in these blocks. We also added to G 1 undesirable relations having an intermediate level between the variables {X 1 , • • • , X 5 } and {X 9 , • • • , X 13 } and saved the result in graph G 2 , as shown Fig. 2.4(topright). More quantitatively, the reference blocks S 1 Ref and S 2 Ref have inner similarities sampled following the normal law N (0.75, 0.1), the background ones are sampled following N (0., 0.1), and in Fig. 2.4, the undesirable relations are sampled following N (0.37, 0.1). In each reference block, the relation between the simulated reference variables and other variables are finally sampled following N (0.9, 0.1). The norm of these similarities are then considered and the similarities higher than one are set to one.

We will measure hereafter the coherence of different variable subsamples to illustrate this notion and make clear its interest in the CORE-clustering context. As the CORE-clustering algorithms proposed in Section 3 use maximal spanning trees, we will additionally discuss the corresponding coherences obtained on the maximal spanning trees of G 1 and G 2 , which are shown in Fig. 2.4 (bottom). The estimated representative variables will finally be given in all tests to make sure their estimation is robust.

S 1 T 1 S 2 T 1 S 1 T 2 S 2 T 2 S 1 T 3 S 2 T 3 S 1 T 4 S 2 T 4 G 1 0.
X1 T 1 X2 T 1 X1 T 2 X2 T 2 X1 T 3 X2 T 3 X1 T 4 X2 T 4 G 1 4

Definitions

Figure 2.4 -Adjacency matrices of the simulated graphs G 1 and G 2 of Appendix 6.1 and their maximum spanning trees (ST).

Illustration of the coherence

We give

= {X 0 , • • • , X 6 } and S 2 T 1 = {X 8 , • • • , X 14 }. (Test 2) Subsamples of S 1 Ref and S 2 Ref of size three, i.e. S 1 T 2 = {X 2 , X 3 , X 4 } and S 2 T 2 = {X 10 , X 11 , X 12 }. (Test 3) Samples S 1 Ref and S 2
Ref in which a variable was replaced with the independent variable X 7 , i.e. S 1

T 3 = {X 1 , • • • , X 7 } and S 2 T 3 = {X 7 , • • • , X 13 }. (Test 4) Samples S 1
Ref and S 2 Ref in which the variables X 4 and X 10 were swapped, i.e. S 1

T 4 = {X 0 , • • • , X 3 , X 10 , X 5 , X 6 } and S 2 T 4 = {X 8 , X 9 , X 4 , X 11 , • • • , X 14 }.
Let us first interpret the results in Table 2.1. The coherences first row on G 1 in tests T 1 and T 2 first show that similar coherences were obtained with CORE-clusters of size 7 and 3, although the coherences are slightly higher for smaller CORE-clusters. These coherences are also clearly higher to those obtained in tests T 3 and T 4 where all variables are not contained in the same simulated block. Interestingly, the results obtained on graph G 2 are similar to those obtained on G 1 . The only difference here is that the coherences of T 4 are slightly higher than in G 1 but still relatively low. Note that the tested CORE-clusters may lead to disconnected subgraphs when tested on the maximum spanning trees. Eq. (2.2) does not make sense in this case. When the simulated subgraphs were connected (12 cases out of 16), we obtained the same coherences on the maximum spanning trees and the whole graph. As discussed in Subsection 2.4, the coherence of a given CORE-cluster on a maximum spanning tree is indeed lower or equal to the coherence of the same CORE-cluster on the whole graph. We however computed the representative variables in all tested cases as the centrality measure makes sense even on disconnected subgraphs (see Subsection 3.3). The results in Table 2.2 show that the estimated representative variables are always in the reference sets S 1

Ref and S 2 Ref in the tested configurations. They also correspond to the simulated representative variables, X 3 and X 11 , in most cases or are very close to these variables. Note that slightly inaccurate reference variables were detected in S 1

Ref , and we can clearly see Fig. 2.4(top) that the influence of its simulated representative variable is less obvious than in set S 2 Ref .

Influence of the undesirable relations

We further study the influence of the undesirable relations between {X 1 , • • • , X 5 } and {X 9 , • • • , X 13 } in graph G 2 by simulating these relations with different strengths. In the previous subsection undesirable relations were sampled following N (µ, 0.1), where µ = 0.37. Here, we sampled 100 graphs G 2 for each strength µ ∈ {0.1, 0.40, 0.60, 0.80, 0.90}. For each graph, we then measured the portion of representative variable estimates in the true reference block of variables and the portion of representative variable estimates that are the ground truth representative variables.

Results are given in Table 2.3 and show that the representative variables detection is particularly stable in these tests, even for large values of µ. All estimated representative variables are indeed in the true block of variables except in Test 4 (where two variables of the reference sets are swapped) with µ = 0.9, i.e. with the same level of similarity as between the blocks representative variables and the other variables they contain. False estimations are however uncommon even in this case. Exact estimates of the representative variables are naturally less frequent as this test is more strict. They are however always clearly higher than random estimations which would have portions equal to 0.14. The estimates of pertinent representative variable therefore appears as robust, even with strong undesirable relations in the variable similarities and tested CORE-clusters which contain undesirable variables. ets). Each portion was computed on 100 simulated graphs (G 2 ) and their corresponding maximum spanning trees (ST (G 2 )). For each group of 100 graphs, a different strength µ of the undesirable relations is tested.

Test 1 Test 2 Test 3 Test 4 G 2 ST (G 2 ) G 2 ST (G 2 ) G 2 ST (G 2 ) G 2 ST (G 2 ) µ = 0.1 1. (0.

Introduction

Graphs play a central role in complex systems as they can model interactions between variables of the system. They are commonly used in a wide range of applications, from social sciences (e.g. social networks [START_REF] Handcock | Modeling social networks form sampled data[END_REF]) to technologies (e.g. telecommunications [START_REF] Smith | The integration of communications networks in the intelligent building[END_REF], wireless sensor networks [START_REF] Akyildiz | Wireless sensor networks: a survey[END_REF]) or biology (gene regulatory networks [START_REF] Davidson | Gene regulatory networks[END_REF], metabolic networks [START_REF] Jeong | The large-scale organization of metabolic networks[END_REF]). One of the most relevant features when analyzing graphs is the identification of their underlying structures, such as cluster structures, generally defined as connected subsets of nodes that are more densely connected to each other than to the rest of the graph. These clusters can provide an invaluable help in understanding and visualizing the functional components of the whole graph [START_REF] Girvan | Community structure in social and biology networks[END_REF][START_REF] Newman | Finding and evaluating community ctructure in networks[END_REF][START_REF] Abbe | Community detection and stochastic block models: recent developments[END_REF]. For instance, in genetics, groups of genes with high interactions are likely to be involved in a same function that drives a specific biological process.

Since the pioneering exploratory works in the early 50s, a large number of clustering methods have launched. Among them, partitioning algorithms, which include the well-known k-means [START_REF] Macqueen | Some methods for classification and analysis of multivariate observations[END_REF], classify nodes into a predefined number of groups based on a similarity measure and hierarchical clustering algorithms [START_REF] Hastie | The Elements of Statistical Learning[END_REF] build a hierarchy of clusters through dendrogram representations. More recently, spectral clustering algorithms, popularized over years by [START_REF] Shi | Normalized cuts and image segmentation[END_REF]; [START_REF] Ng | On spectral clustering: Analysis and an algorithm[END_REF], particularly draw the attention of the community research due to their speed, simplicity and numerical performances. As its name suggest, spectral clustering algorithms mainly use the spectral properties of the graph by (i) computing the eigenvectors of the associated Laplacian matrix (or one of its derivatives), which gives information about the structure of the graph, and (ii) performing k-means on it to recover the induced cluster structure. A large number of extensions of the original spectral clustering algorithm, as presented in [START_REF] Luxburg | A tutorial on spectral clustering[END_REF], have been proposed, with applications to different fields [START_REF] Zelnik-Manor | Self-tuning spectral clustering[END_REF][START_REF] Wang | Flexible constrained spectral clustering[END_REF][START_REF] Li | Spectral clustering in heterogeneous information networks[END_REF].

While spectral clustering is widely used in practice, handling noise sensitivity remains a tricky point [START_REF] Bojchevski | Robust spectral clustering for noisy data: Modeling sparse corruptions improves latent embeddings[END_REF], mainly due to the k-means algorithm, which is highly sensitive to noise. This issue has been considerably studied with extensions of the k-means to noisy settings so that it recovers the cluster structure in spite of the unstructured part of the input data [START_REF] Tang | Noise identification with the k-means algorithm[END_REF][START_REF] Pelleg | K-means with large and noisy constraint sets[END_REF]. More generally, the robustness of spectral clustering algorithms has recently been investigated for perturbed graphs derived from stochastic block models (SBM) [START_REF] Stephan | Robustness of spectral methods for community detection[END_REF][START_REF] Peche | Robustness of community detection to random geometric perturbations[END_REF]. In this paper, we develop an alternative method of the spectral clustering, called 1 -spectral clustering algorithm and based on Lasso regularization [START_REF] Tibshirani | Estimating the number of clusters in a data set via the gap statistic[END_REF]. Note that research papers have explored regularized spectral clustering to robustly identify clusters in large networks. Although [START_REF] Zhang | Understanding regularized spectral clustering via graph conductance[END_REF] and [START_REF] Joseph | Impact of regularization on spectral clustering[END_REF] show the effect of regularization on spectral clustering through graph conductance and respectively through stochastic block models. Equally, [START_REF] Lara | Spectral embedding of regularized block models[END_REF], shows on a simple block model that the spectral regularization separates the underlying blocks of the graph.

However, in our model, as in the spectral clustering algorithm, we carefully explore the underlying structure of the graph through the Laplacian matrix spectrum to cluster nodes. However, by directly promoting a sparse eigenvectors basis solution to an 1 -norm optimization problem, it does not require the k-means step to extract clustering structures, making it more robust in highly perturbed graph situations.

The paper is organized as follows: in Section 2, we introduce some preliminary concepts about graph clustering and more specifically spectral clustering. In Section 3 and 4, we present the 1spectral clustering we developed, from a theoretical and an algorithmic point of view. In Section 5, we finally show its efficiency and accuracy through experiments on simulated and biological real data set and compare it with state-of-the-art clustering methods.

Reminders about graph and spectral clustering 2.1 Graphs modeling and notations

This work considers the framework of an unknown undirected graph G(V, E), with no retroactive loop, consisting of p vertices V = {1, . . . , p} and a set of edges E ⊆ V × V connecting each pair of vertices. As usual, the graph G is represented by its associated adjacency matrix A = (A ij ) (i,j)∈E of size p × p, whose non-zero elements correspond to the edges of G:

∀(i, j) ∈ 1, p 2 , A ij = 1 if (i, j) ∈ E, 0 otherwise.
As G is undirected with no retroactive loop, the adjacency matrix A is symmetric with zero on its diagonal. Before turning to the next section, we recall some useful graph definitions.

Definition 3

The degree d i of a node i ∈ V of G is defined as the number of edges that are incident to i: d i = p j=1 A ij . The induced degree matrix D is then the p × p matrix containing (d 1 , . . . , d p ) on its diagonal and zero elsewhere:

D = diag (d 1 , . . . , d p ).
Definition 4 A connected component C of G is a subset of nodes from V such that each pair of nodes of C is connected by a path and there is no connection between vertices in C and outside C.

Connected components C 1 , . . . , C k are a k-partition of the set V of vertices if the three following conditions hold:

1. they are non-empty: ∀i ∈ 1, k , C i = ∅, 2. they are pairwise disjoints:

∀(i, j) ∈ 1, k 2 , C i ∩ C j = ∅,
3. their union form the set of all vertices:

k ∪ i=1 C i = V .
Definition 5 Let C 1 , ..., C k be a k-partition of the set of vertices V of G. Then, the indicators (1 C i ) i∈{1,...,k} of this partition are defined as the vectors of size p, whose coefficients satisfy:

∀i ∈ 1, k , ∀j ∈ 1, p , (1 C i ) j = 1 if vertex j belongs to C i , 0 otherwise.
In the present paper, we assume that the graph G is the union of k complete graphs, whose

set of vertices C 1 , . . . , C k form a k-partition of G. We denote by c 1 , • • • , c k their respective size ( k i=1 c i = p).
To simplify, we assume that the nodes, labeled from 1 to p, are ordered with respect to their block membership and the size of the blocks. From a matrix point of view, the associated adjacency matrix A is a k-block diagonal matrix of size p × p of the form:

A =                0 1 • • • 1 1 . . . . . . . . . . . . . . . . . . 1 1 • • • 1 0 c 1 0 . . . 0 0 1 • • • 1 1 . . . . . . . . . . . . . . . . . . 1 1 • • • 1 0 c k                . (3.1)

Graph clustering through spectral clustering

Graph clustering consists in grouping the vertices of the graph G into clusters according to its edge structure. Whereas some of the most traditional clustering algorithms are based on partitions (e.g. k-means) and hierarchies (e.g. hierarchical clustering algorithm), spectral clustering takes advantage of the spectral properties of the graph. A large number of spectral clustering algorithms exists in the literature. The most common version, presented in [START_REF] Luxburg | A tutorial on spectral clustering[END_REF] and recapped in Algorithm 4 below, uses the properties of the Laplacian matrix (Definition 6) to detect clusters in the graph.

Definition 6 Given a graph G, the Laplacian matrix L is defined as:

L = D -A,
where A is the adjacency matrix and D the degree matrix associated to G.

Algorithm 4 Spectral clustering algorithm

Require: G a graph, A its associated adjacency matrix, k number of clusters to build.

1: Compute the Laplacian matrix L = D -A.

2: Perform the spectral decomposition of L and store the k first eigenvectors

U := (u 1 , • • • , u k). 3: Cluster U with the k-means algorithm into clusters C 1 , • • • , C k. 4: return Clusters C 1 , • • • , C k.
By definition, the diagonal of L is equal to the degrees of the nodes. Moreover, in the ideal case where G has an underlying partition form with k connected components and a block diagonal adjacency matrix A, as given in Equation (3.1), the eigenvalue 0 of L is of multiplicity k and the associated eigenvectors correspond to the indicator vectors of the k components. These k components can then be recovered only by performing spectral decomposition of L. However, in the perturbed case, any perturbation caused by introducing and/or removing edges between and/or inside the components makes k -1 of the k eigenvalues 0 slightly larger than 0 and changes the corresponding eigenvectors. The final cluster structure is thus no longer explicitly represented. The spectral clustering algorithm then uses the k-means algorithm on these eigenvectors to discover the hidden underlying structure, which is hampered by perturbations.

Since the first development of the spectral clustering algorithm, it has been studied a lot and extended many times in different communities [START_REF] Hagen | New spectral methods for ratio cut partitioning and clustering[END_REF][START_REF] Hendrickson | An improved spectral graph partitioning algorithm for mapping parallel computations[END_REF][START_REF] Pothen | Graph partitioning algorithms with applications to scientific computing[END_REF][START_REF] Shi | Normalized cuts and image segmentation[END_REF][START_REF] Ng | On spectral clustering: Analysis and an algorithm[END_REF][START_REF] Zelnik-Manor | Self-tuning spectral clustering[END_REF] with powerful results. Refinements include the use of normalized versions of the Laplacian matrix, such as the symmetric and the random walk normalized ones [START_REF] Luxburg | A tutorial on spectral clustering[END_REF]. Nevertheless, the performances of the spectral clustering have shown to be very sensitive to perturbations, which often occurs when dealing with real data [START_REF] Bojchevski | Robust spectral clustering for noisy data: Modeling sparse corruptions improves latent embeddings[END_REF]. To provide more robustness with respect to perturbations, we thus developed the 1 -spectral clustering algorithm, described in Section 3.

An 1 -version of the spectral clustering algorithm

In this section, we describe the 1 -spectral clustering algorithm we developed as an alternative to the standard spectral clustering for clustering perturbed graphs. In this context, to ensure a good recovery of the connected components, the eigenvectors basis should be carefully defined. The key point is to replace the k-means procedure, which fails while the perturbation grows, by selecting relevant eigenvectors that provide useful information about the graph structure. As the spectral clustering algorithm, the 1 -spectral clustering focuses on the spectral properties of the graph.

Let G = (V, E) be a graph formed of k connected components, as defined in Section 2, and A its associated adjacency matrix. We denote by (λ i ) 1≤i≤p the p-eigenvalues of A, sorted in increasing order:

λ 1 ≤ ... ≤ λ p ,
and v 1 , ..., v p their associated eigenvectors. We define by V k the eigenspace generated by the k largest eigenvectors:

V k := Span(v p-k+1 , ..., v p ).

In the ideal case, where the graph is not perturbed, the indicators (1 C i ) 1≤i≤k of the connected components C 1 , . . . , C k correspond exactly to the eigenvectors of the Laplacian matrix L associated to the eigenvalue 0 of multiplicity k (see Section 2.2). As regards the adjacency matrix A, these indicators correspond this time to the k eigenvectors v p-k+1 , . . . , v p , associated to the k largest eigenvalues λ p-k+1 , ..., λ p . In the perturbed case, unlike the traditional spectral clustering, the 1spectral clustering algorithm does not directly use the subspace V k to recover the k connected components but computes another eigenbasis that promotes sparse solutions, as detailed in the next sections.

General 0 -minimization problem

v∈V k \{0} v 0 (P 0 )
has a unique solution (up to a constant) given by 1 C 1 .

In other words, 1 C 1 is the sparsest non-zero eigenvector in the space spanned by the eigenvectors associated to the k largest eigenvalues.

Proof We recall that, for all v ∈ R p , v 0 = | {j ∈ 1, p , v j = 0} |. Let v ∈ V k \ {0}. As (1 C j ) 1≤j≤p ∈ V k , v can be decomposed as v = k j=1 α j 1 C j where α = (α 1 , . . . , α k ) ∈ R k
and there exists j ∈ {1, ..., k} such that α j = 0. By definition of the 0 -norm, we then have:

v 0 = 1 α 1 =0 c 1 + • • • + 1 α k =0 c k , (3.2)
with c 1 ≤ ... ≤ c k the sizes of the k connected components. The solution of (P 0 ), which minimizes Equation (3.2), is thus given by setting α = (α 1 , 0, . . . , 0) with α 1 = 0.

Proposition 3 can then be generalized to iteratively find the indicators associated to the largest connected components introducing sparsity and orthogonality constraints. For i ∈ 2, k , let V i k refers to:

V i k := {v ∈ V k , ∀l = 1, . . . , i -1, v ⊥ 1 C l } . Proposition 4 Let i ∈ 2, k . The minimization problem arg min v∈V i k \{0} v 0 (P i 0 )
has a unique solution (up to a constant) given by 1 C i .

Solving (P 0 ) and (P i 0 ) 2≤i≤k is a NP-hard problem, which is not computationally feasible. To tackle this issue, the classical idea consists in replacing the 0 -norm by its convex relaxation, the 1 -norm, defined for all v ∈ R p as v 1 = 1≤j≤p |v j |.

In the next section, we show that the solution of the 0 optimization problems remains the same by replacing the 0 -norm by the 1 -norm, at the price of slight constraints on the connected components.

Relaxed 1 -minimization problem

From now on, we assume that we know one representative element for each component, that is a node belonging to each component, denoted by

(i 1 , ..., i k ) thereafter. Let Ṽk = {v ∈ V k , v i 1 = 1}.
Then, it is straightforward to see that the indicator vector of the smallest component is solution to the following optimization problem:

Proposition 5 The minimization problem arg min v∈ Ṽk v 1 (P 1 )
has a unique solution given by 1 C 1 .

Proof We recall that, for all

v ∈ R p , v 1 = p j=1 |v j |. Let v ∈ Ṽk . As (1 C j ) 1≤j≤p ∈ V k , v can be decomposed as v = k j=1 α j 1 C j where α = (α 1 , . . . , α k ) ∈ R k
and there exists j ∈ {1, . . . , k} such that α j = 0. By definition of the 1 -norm, we then have:

v 1 = |α 1 |c 1 + • • • + |α k |c k , (3.3)
with c 1 ≤ ... ≤ c k the sizes of the k connected components. The solution of (P 1 ), which minimizes Equation (3.3), is thus given by setting α = (α 1 , 0, . . . , 0) with α 1 = 1.

To simplify and without loss of generality, we assume that i 1 corresponds to the first node. We can then rewrite (P 1 ) as:

arg min

v∈R p-1 (1,v) T ∈V k v 1 .
Constraints in (P 1 ) can be converted into the following equality contraints: , v p-k ) the matrix formed by the eigenvectors associated with the p -k-smallest eigenvalues. We denote by w T its first row and W T the matrix obtained after removing w T from U k :

Proposition 6 Let U k := (v 1 , ...
U k := (v 1 , ..., v p-k ) =       b w T W T b       (3.4)
The minimization problem arg min

v∈R p-1 W v=-w v 1 ( P1 ) has a unique solution v * such that (1, v * ) T = 1 C 1 .
Proof Since A is symmetric, its eigenvectors form an orthogonal basis and, for all v ∈ V k , we have

U T k v = 0. Let (1, v) T ∈ V k .
Using Equation (3.4), we deduce that:

U T k 1 v = w + W v = 0.
The constraint in ( P1 ) is thus equivalent to the constraint in (P 1 ), which ends the proof.

Generalization of the relaxed 1 -minimization problem

Obviously, the indicator vector 1 C 1 alone is not sufficient to know the complete graph structure. However, Proposition 6 can be extended to find the remaining indicator vectors. To do so, as in Proposition 4, we add the constraint that the target vector is orthogonal to the previously computed vectors, which is done in practice by applying a Gram-Schmidt orthonormalization procedure (see Section 4 below for more details about the procedure).

4 The 1 -spectral algorithm

Global overview of the algorithm

In this section, we present a global overview of the 1 -spectral clustering algorithm we implemented to recover the components of any perturbed graph (see Algorithm 5 below). It is available as a R-package on GitHub at https://github.com/championcamille/l1-SpectralClustering. In the next paragraphs, some details about the algorithm and parameters setting are given.

Algorithm 5 1 -spectral clustering algorithm 1: Input: G a graph, A its associated adjacency matrix, k number of clusters to recover and (i j ) j∈{1,..., k} family of representative elements of each cluster. 2: Perform the spectral decomposition of A, sort the eigenvalues by increasing order and store the associated eigenvectors: V := (v 1 , ..., v p ). 3: for j = 1 to k do 4:

Define U k,j as the matrix that contains the p -k -j + 1 first columns of V : U k,j := (v 1 , ..., v p-k-j+1 ).

5:

Split U k,j into two parts:

w T := U i j k,j the i j -th row of U k,j and W T := U -i j k,j
the other rows of U k,j .

6:

Solve the 1 -minimization problem ( P1 ):

v * := arg min v∈R p-1 W v=-w v 1 . 7:
Recover the indicator of the j-th component:

1C j = (v * 1 , . . . , v * i j -1 , 1, v * i j , . . . , v * p ). 8: Update v j in V : v j ← 1C j . 9:
Perform Gram-Schmidt orthogonalization on V to ensure orthogonality between v j and the rest of the columns of V : V ← Gram-Schmidt(V ).

10: end for 11: Output: ( 1C j ) 1≤j≤ k the indicators of the k connected components.

Solving the 1 -minimization problem

This section is devoted to the resolution of the constrained 1 -optimization problem ( P1 ) (line 6 of Algorithm 5). To be simplified, it can be equivalently written as the following penalized problem:

arg min v∈R p-1 W v + w 2 2 + λ v 1 , (P Lasso )
where, for all v ∈ R p-1 , v 2 2 = p-1 j=1 v 2 j and λ > 0 is the regularization parameter that controls the balance between the constraint and the sparsity. Two methods are proposed thereafter to solve (P Lasso ).

Lasso solution enrg

The most traditional method to deal with such an 1 -minimization problem is the Lasso procedure, developed by [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF]. As for all regularizing methods, the choice of λ is of great importance. Here, especially, taking λ too large will lead to an over-constrained problem and a large number of nodes of G may not be clustered into components. In practice, K-fold cross validation, as implemented in the glmnet R-package, can be used to optimally set λ.

Thresholded least-squares solution enrg Another method consists in solving the least-squares problem:

v * := arg min v∈R p-1 W v + w 2 2 ,
and then thresholding v * given some predefined threshold t:

∀j ∈ 1, p -1 , v * j = 1 if v * j > t, 0 otherwise.
Of course, this thresholding step imposes sparsity on the solution. However, we can wonder if nodes with large coefficients should really be clustered together. In our model, the ideal parameters to recover (indicators of the components) do not take continuous values. Enforcing the coefficients of all representative elements to be equal to 1, under small perturbations, the coefficients of all other nodes belonging to the same components should then be close to 1. This specific behavior is underlined in Figure 3.1. In this example, we generated a graph G with 50 nodes, split into 5 connected components. We perturbed the structure of the graph by adding and removing edges with a probability P of 1%, 10%, 25% and 50%. We then solved (P Lasso ) to recover the first component only. As can be seen in Figure 3.1, the Lasso and thresholded least-squares solutions give almost the same results: for small perturbations (P ≤ 10%, at the top), the whole component is perfectly retrieved. For P = 25% (at the bottom left), all coefficients are tighter but both methods still work, forgetting only one node. As the perturbation becomes too large (P = 50%, at the bottom right), the selection of nodes belonging to the first component fails. Figure 3.1 -Evolution of the coefficients of v, solution of (P Lasso ), with respect to v 1 for different perturbations of the ideal graph (from top left to bottom right P = 1%, 10%, 25% and 50%). Red lines correspond to the coefficients belonging to the component we aim at recovering, in contrast with black ones. Dotted lines are related to the 1 -norm-threshold (vertical), associated with the Lasso solution, and the threshold on the value of the coefficients (horizontal), associated with the thresholded least-squares solution.

Optimally tuning the number of clusters

Traditional clustering algorithms, such as k-means, require the user to specify the number of connected components of the graph G to recover, which is, in practice, unavailable. Determining the optimal number of components k thus becomes a fundamental issue. A large number of methods have been developed in this sense: the hierarchical clustering for example looks for a hierarchy of components using dendrograms. The Elbow, average silhouette and gap statistic methods [START_REF] Tibshirani | Estimating the number of clusters in a data set via the gap statistic[END_REF] are also frequently used in addition to clustering techniques.

In our work, as proposed by [START_REF] Luxburg | A tutorial on spectral clustering[END_REF], we focus on the heuristic eigengap method, which consists in choosing k such that it maximizes the eigengap, defined as the difference between consecutive eigenvalues of the Laplacian matrix L. This procedure is particularly well-suited in a spectral context. Indeed, in the ideal case, perturbation theory ensures that there exists a gap between the eigenvalue 0 of multiplicity k and the next k + 1-th one. In the perturbed case, while being less strong, an eigengap still exists.

Finding the representative elements

In addition to the number of connected components, to run the 1 -spectral clustering algorithm, we need to know at least one representative element of each component. This assumption may be restrictive when working with real data. However, it makes sense in a large number of situations where clusters are chosen to classify nodes around specific elements of the graph.

To avoid an arbitrary choice of such elements, one solution consists in estimating them using a rough partitioning algorithm. Another solution is to explore the structure of the graph to find hubs of densely connected parts. In this work, this is done by computing the betweeness centrality score of all nodes. In graph theory, the betweeness score S b measures the centrality of a node based on the number of shortest paths passing through it:

∀ ∈ 1, p , S b ( ) = 1≤i,j≤p
# shortest paths from i to j # shortest paths from i to j passing through .

In practice, the representative elements of the k components are chosen to maximize this score. Note that the nodes with the highest betweeness scores should be those that connect the densest parts of the graph. The risk of clustering two nodes from different connected components may thus be high. To avoid this, we add a stabilization step to our algorithm. As soon as one of the nodes with the k highest scores is added to a component during the minimization step, it is removed from the list of potential representative elements. We then re-run the algorithm using the k nodes taken among the k + 1 ones with the highest scores, and so on until stabilization of the list of representative elements.

Numerical experiments

This section is dedicated to experimental studies to assess numerical performances of the 1spectral clustering algorithm through two datasets. First, we show that it behaves well on simulated data with a variety of different settings and in comparison with state-of-the-art spectral clustering methods. Then, using a gene expression data set from kidney cancer patients, we demonstrate the ability of our algorithm to discover relevant groups of genes that act together to characterize the disease.

Application to toy datasets

Numerical settings

To explore the capabilities and the limits of the 1 -spectral clustering algorithm with respect to state-of-the-art methods, we first considered simulated data, whose settings are detailed in the next paragraphs.

Simulated data set zoej

We generated random ideal graphs for a given number of nodes p (p = 20, 50 and 100) and a given number of connected components k (k = 2, 5, 10). The component sizes (c j ) 1≤j≤k were chosen in a balanced way: ∀j ∈ 1, k -1 , c j = p/k , with k j=1 c j = p. With a probability p in and p out of removing an edge from a component and of introducing an edge between two components varying from 0.01 to 0.5, we created multiple perturbed versions of the graph. All experiments were replicated 100 times each for better robustness.

Algorithm parameters oeruh

As some of the methods we compare with require the number of components to form, we focus on two versions of the 1 -spectral clustering: the one presented in Algorithm 5, for which the number of clusters and a list of representative elements are assumed to be known, and the self-tuned one, for which both of them are extracted from the graph, as explained in Section 4. The results being very similar, we choose to focus on the thresholded least-squares solution to solve the 1 -minimization problem ( P1 ) in Algorithm 5. The corresponding threshold parameter t is fixed using 5-fold cross-validation.

Comparison with state-of-the-art zrizi We compare the 1 -spectral clustering with three other graph-based clustering algorithms: first, the spectral clustering (Algorithm 4), which is available in the R-package anocva, then, SpectACl, which was developped by [START_REF] Hess | The SpectACl of non convex clustering: a spectral approach to density-based clustering[END_REF] with the aim of exhibiting both minimum cut and maximum density of the clusters. This algorithm can be viewed as a combination of DBSCAN, a density-based clustering algorithm [START_REF] Ester | A density-based algorithm for discovering clusters in large spatial databases with noise[END_REF] which is mainly used to identify clusters of any shape in a data set containing noise and outliers, and spectral clustering. SpectACl is publicly available on the Bitbucket platform as a Python code at https://bitbucket.org/Sibylse/spectacl/ src/master/. Both methods requiring the number of components to cluster as an input, we finally run the Self-Tuning Spectral Clustering from Zelnik-Manor and Perona (2005), available on GitHub as a Python code at https://github.com/wOOL/STSC. The latter is an improved version of the spectral clustering, in which the final postprocessing step (k-means) is removed and the structure of the eigenvectors is carefully analyzed to automatically infer the number of clusters. It is thus used to evaluate the performances of the self-tuned version of the 1 -spectral algorithm.

Performance metrics ieojr Performances are measured by comparing the learnt components with the true ones, which are obviously known in the context of simulated data. Among the large number of existing scores, we p=20 p=50 p=100 used the Normalized Mutual Information (NMI) score, for its ability to compare clusters that could be of different sizes. The closer to 1 the NMI score, the better the classification.

p in p out k = 2 k = 5 k = 2 k = 5 k = 10 k = 5 k = 10 0.01 0.01 1 1 1 1 1 1 1 0.1 1 1

Effect of the dimension and cluster sizes on perturbed graphs

First, we aimed at exploring the effect of the dimension and cluster sizes on the performances of the 1 -spectral clustering algorithm. For p ranging from 20 to 100 and k from 2 to 10, results, in terms of NMI scores, are summarized in Table 3.1. On the one hand, for fixed p and k, when the perturbations are small (p in , p out < 0.25), one may note that the 1 -clustering algorithm works well (it is clearly a favorable situation). On the other hand, a crude decrease in performance results can be observed as the perturbation grows (p in or p out ≥ 0.25). In that case, the perturbed graph is far from the original one, which makes hard the recovery of the components. As expected, this becomes even more significant while the dimension increases (from top left to bottom right for each value of p). For perturbations of 0.5 (last line), the NMI scores do not exceed 0.5, which means that the 1 -spectral clustering algorithm almost fails to recover the components. However, we must keep in mind that imposing a perturbation of 0.5 on a graph strongly affects its structure, with a probability of removing an edge inside a component and introducing an edge between components of 50%.

Performance results with respect to state-of-the-art

To give more credit to the 1 -spectral clustering algorithm, we also evaluated its robustness in comparison with the spectral and SpectACl algorithms (see Section 5.1.1) for clustering different perturbed versions of a graph with p = 100 and k = 10. For each perturbation, we generated 100 graphs and computed the clustering performances using NMI scores.

Results can be visualized in Figure 3.2, which also indicates the 50% confidence interval. As can be seen, the 1 -spectral and spectral clustering algorithms are very similar, especially for small perturbations (p out < 0.25). However, as the perturbations grow, the 1 -spectral clustering algorithm shows a smaller impact to noise sensitivity than the spectral one, being, almost ever, the best method. The results of SpectACl are oddly bad but this may be due to the fact that it was developed for clustering nonconvex shapes, which is beyond the scope of the present work. An interesting question is how the self-tuning version of the 1 -spectral clustering, for which the number of clusters and the representative elements are self-evaluated, compare with the Self-Tuning Spectral Clustering (see Section 5.1.1). For p = 20, k = 2 and perturbations ranging from 0.01 to 0.5, we generated 100 versions of the same graph. Results are given in Figure 3.3. At the top, the NMI scores indicate that the performances of the self-tuning 1 -spectral clustering (in red) decrease while the perturbations grow. On the contrary, the Self-Tuning Spectral Clustering (in blue) seems to be less sensitive to the increase of p in but provide bad results for p in < 0.25 and p out ≥ 0.25, even though it is a more favorable situation. Some of the lowest performances of the self-tuning 1 -spectral clustering can be explained by observing that this algorithm was developed in a sparse form, that is all nodes from a perturbed graph are not automatically clustered into components. In practice, this has huge consequences on the NMI score, which processes the non-classified nodes as wrongly-classified ones. When considering only the classified nodes, the NMI scores of the self-tuning 1 -spectral clustering can be found in green, with, of course, better results.

More generally, the NMI scores are particularly sensitive to the number of estimated clusters. At the bottom of Figure 3.3, for each perturbation, the estimated number of clusters of both methods, which should be close to k = 2, the true number of clusters, can be visualized. It is clear that the further to 2 the estimation is, the smaller the associated NMI scores. The self-tuning version of the 1 -spectral clustering was not optimized in this sense but it should be the key for an improvement and a stabilization of the performances.

Application to cancer data

This section is dedicated to the application of the 1 -spectral clustering algorithm on a real kidney cancer data set from The Cancer Genome Atlas (TCGA) project. After describing the data (Section 5.2.1), results are presented in Section 5.2.2 and followed by a discussion (Section 5.2.3).

The kidney cancer data set

The Cancer Genome Atlas (TCGA) is an american project from the National Cancer Institute (NCI) and the National Human Genome Research Institute (NHGRI), which was launched fifteen years ago with the aim of characterizing genetic mutations responsible for cancer using genome sequencing and bioinformatics methods. Since then, millions of data have been produced and made publically available. In this work, we focused on KIdney Renal clear cell Carcinoma, abbreviated to KIRC thereafter. KIRC is one of the most common types of cancer, usually affecting people (mainly men) around 60 years old. Even if the chances of surgical cures are good, KIRC is hard to detect with no early symptoms, which makes it even dangerous.

In this work, we extracted gene expression data for KIRC from the TCGA data portal http: //gdac.broadinstitute.org/. These data were produced using RNA-sequencing for a total number of 16, 123 genes and 532 cancer patients. After preprocessing by log-transformating, quantile normalizing the arrays and filtering genes based on variance, we only kept 75% of them, i.e. 12, 092 genes.

1 -spectral clustering algorithm on kidney cancer data

Applying the 1 -spectral clustering algorithm to cluster genes into components require the knowledge of an initial network that models the relationships between genes. The latter are usually represented through Gene Regulatory Networks (GRNs), which are directed graphs that connect genes based on regulations (activations/inhibitions). Here, we focused on co-regulated networks, a simplified version of GRNs, where no causality exists. Edges are thus undirected, representing co-regulations, or correlations in terms of expression between genes. To create such a network, we computed the correlation matrix, based on Pearson's correlation, between all pairs of genes and then thresholded the matrix by removing edges with correlation smaller in absolute value than 0.7. This network is made of 4,982 genes (see Figure 3.4 (a) for an overview of the network).

We then applied 1 -spectral clustering algorithm on the adjacency matrix associated with the co-regulatory network described above. The 4,982 genes were clustered into 186 components, from size 2 to 986, with an averaged number of genes of 27. These components are represented in Figure 3.4 (b), with different colors.

Clusters as hallmarks of kidney cancer

In this section, we investigate the biological hypotheses that can be deduced from the network. First, to assign biological meaning to each component of the network, we performed gene set enrichment analysis. To this aim, we used the databases GeneSetDB [START_REF] Culhane | GeneSigDB, a curated database of gene expression signatures[END_REF] and MSigDB [START_REF] Liberzon | The molecular signatures database (MSigDB) hallmark gene set collection[END_REF], restricted to hallmark (H), curated (C2), GO (C5), oncogenic (C6) and immunologic signatures (C7) gene sets, which include the gene sets most relevant to cancer gene expression profiles. Enrichments were evaluated by performing hypergeometric tests, corrected for multiple testing using the FDR [START_REF] Benjamini | Controlling the false discovery rate: a practical and powerful approach to multiple testing[END_REF]. Among the 186 identified components, four particularly drew our attention (see Figure 3.4 (c)). These components are described in details in the next paragraphs.

Transmembrane activity cluster

The first cluster we identified (Figure 3.4 (c), red cluster) is made of 43 genes. This cluster gathers genes involved in the same "transmembrane activity" pathway. Indeed, among the 43 genes of the component, 11 are members of the same group SLC of solute carriers transporters, which aims at facilitating the transport of substrates across membranes. This is confirmed by the gene set enrichment analysis we performed, with p-values ranging from 1.97 × 10 -6 from 1.89 × 10 -9 .

Epithelial-mesenchymal cluster The second cluster (Figure 3.4 (c), purple cluster) of seven genes is highly enriched in Epithelial-Mesenchymal Transition (EMT) pathways, a natural process that converts epithelial cells into mesenchymal phenotypes and is often altered in cancers. This cluster includes the gene SERPINH1, a known EMT-related gene, which has been identified as a potential biomarker of kidney cancers [START_REF] Qi | SERPINH1 overexpression in clear cell renal cell carcinoma: association with poor clinical outcome and its potential as a novel prognostic marker[END_REF].

T-cells associated cluster The third cluster (Figure 3.4 (c), blue cluster) includes 38 genes, mostly enriched in T-cells and inflammatory response associated pathways. To confirm this, we tested the correlation of the cluster expression (defined as the averaged expression across all genes from the cluster) with CD4+ and CD8+ T-cells, encoded by the genes CD4, CD8A/B, which play a major role in cancer immunotherapy [START_REF] Tay | Revisiting the role of CD4+ T cells in cancer immunotherapy -new insights into old paradigms[END_REF]. With correlations ranging from 0.46 (CD8B) to 0.86 (CD4) and associated p-values smaller than 10e -16 , we validate this relationship.

In addition, [START_REF] Kawashima | Tumour grade significantly correlates with total dysfunction of tumour tissue-infiltrating lymphocytes in renal cell carcinoma[END_REF] recently reported that both CD4+ and CD8+ T-cells were up-regulated in the patients with kidney cancer of high grade. We obtained the same results when comparing the cluster expression with low (grades 1 and 2) and high (grades 3 and 4) kidney cancer grades: the higher the expression, the worst the grade (p-values for t-test of 0.009, see Figure 3.5). .5 -Boxplots representing the association between the cluster expression (averaged expression across all genes from the cluster) and kidney cancer grades, which range from 1 (low grade) to 4 (high grade), grade x indicating that the grade could not be evaluated.

Liver-signature cluster In the last cluster (Figure 3.4 (c), yellow cluster), we found 96 genes, most of which belong to liver gene signatures, characterizing liver cancers. To a little extent, these genes are also associated with glutathione, an antioxidant that protects cells from important damages. [START_REF] Xiao | Glutathione metabolism in renal cell carcinoma progression and implications for therapies[END_REF] have recently shown that an increased level of glutathione is a hallmark of kidney cancer. To go further, we investigated whether the cluster expression could be used to predict survival in kidney cancer. To this aim, we used Cox proportional hazards modelling. Hazard ratios were used to report the direction of the survival effect and the Wald test was used to determine its significance. As shown in Figure 3.6, high gene expression is significantly associated with good survivals (Hazard ratio of 0.66, p-value of 0.009). This indicates that this cluster may be used as a prognosis of kidney cancers.

High expression Low expression

Conclusion

In this paper, we proposed a new spectral clustering algorithm, called 1 -spectral clustering, for detecting cluster structures in perturbed graphs. To tackle the noise robustness issue of the traditional spectral clustering, the k-means is removed and replaced by writting the indicators of the components as solutions of explicit 1 -constrained minimization problems. The performances of this algorithm are highlighted through numerical experiments, with competitive results when compared to state-ofthe-art. Nevertheless, many opportunities for further improvements can be considered. Firstly, from an algorithmic point of view, it would be interesting to better explore solutions for calibrating the optimal number of clusters and its representative elements. Secondly, future works include theoretical study of the eigenvectors stability, in order to validate the performances of the algorithm. A particular attention may be paid to the more general stochastic block model (SBM), where the edge probabilities depend on the community membership.

Introduction

Non-alcoholic fatty liver disease (NAFLD) is a common consequence of obesity and type 2 diabetes [START_REF] Dimaira | Liver fibrosis in the context of non-alcoholic steatohepatitis: the role of adipokines[END_REF][START_REF] Me | Non-alcoholic fatty liver disease: a systematic review[END_REF]. In NAFLD, the origin of inflammation and hepatocyte injury is related to dietary lipids, bile acids, adipokines and cytokines, to cite a few. Furthermore, gut microbiota seems to be one of the key players of NAFLD development [START_REF] Brandl | Intestinal microbiota and non-alcoholic steatohepatitis[END_REF][START_REF] Hoyles | Molecular phenomics and metagenomics of hepatic steatosis in non-diabetic obese women[END_REF]. Markers and receptors of microbiota-related injury features have been described in this disorder such as TLRs, NODS, and NLRP3 [START_REF] Denou | Defective nod2 peptidoglycan sensing promotes diet-induced inflammation, dysbiosis, and insulin resistance[END_REF][START_REF] Miura | Toll-like receptor 2 and palmitic acid cooperatively contribute to the development of non-alcoholic steatohepatitis through inflammasome activation in mice[END_REF][START_REF] Pierantonelli | Lack of nlrp3-inflammasome leads to gut-liver axis derangement, gut dysbiosis and a worsened phenotype in a mouse model of nafld[END_REF][START_REF] Roh | Toll-like receptors in alcoholic liver disease, non-alcoholic steatohepatitis and carcinogenesis[END_REF] as well as the activation of the innate and adaptive immune systems [START_REF] Bieghs | Innate immune signaling and gut-liver interactions in nonalcoholic fatty liver disease[END_REF]. In early sets of experiments, we initially showed that hepatic steatosis in the obese diabetic mouse was due to an increased circulating concentration of lipopolysaccharides (LPS) i.e. metabolic endotoxemia [START_REF] Cani | Metabolic endotoxemia initiates obesity and insulin resistance[END_REF]. Lipoproteins transport LPS [START_REF] Verges | Changes in lipoprotein kinetics associated with type 2 diabetes affect the distribution of lipopolysaccharides among lipoproteins[END_REF] to tissues, triggering the CD14/TRL4 pathway that increases liver inflammation and fat deposition [START_REF] Cani | Metabolic endotoxemia initiates obesity and insulin resistance[END_REF]. Gut bacteria were also reported to translocate through the intestinal tract to tissues [START_REF] Berg | Immunosuppression and intestinal bacterial overgrowth synergistically promote bacterial translocation[END_REF] such as the adipose depots and the liver, establishing a tissue microbiota as observed in rodents (Amar et al., 2011a;[START_REF] Garidou | The gut microbiota regulates intestinal CD4 T cells expressing RORgammat and controls metabolic disease[END_REF][START_REF] Pomie | Triggering the adaptive immune system with commensal gut bacteria protects against insulin resistance and dysglycemia[END_REF] and humans [START_REF] Sookoian | Intrahepatic bacterial metataxonomic signature in non-alcoholic fatty liver disease[END_REF]Amar et al., 2011b;[START_REF] Burcelin | Metagenome and metabolism: the tissue microbiota hypothesis[END_REF] which could trigger liver inflam-mation and the onset of fibrosis (Amar et al., 2011a). This mechanism activates immune cells, including Kupffer cells, to release various pro-inflammatory cytokines and chemokines [START_REF] Balmer | The liver may act as a firewall mediating mutualism between the host and its gut commensal microbiota[END_REF] damaging the surrounding tissues initi-ating fibrosis. This hypothesis is now largely supported by recent major advances in NAFLD research, which show gut and blood microbiota dysbiosis of patients with advanced stages of NAFLD [START_REF] Lelouvier | Changes in blood microbiota profiles associated with liver fibrosis in obese patients: A pilot analysis[END_REF][START_REF] Loomba | Gut microbiome-based metagenomic signature for non-invasive detection of advanced fibrosis in human non-alcoholic fatty liver disease[END_REF][START_REF] Schierwagen | Circulating microbiome in blood of different circulatory compartments[END_REF]. Hence, the identification of specific groups of translocated bacteria from dysbiotic gut microbiota could aid in the design of novel therapeutic strategies. To address this issue, we have sequenced and identified the bacterial 16S rDNA from liver tissue of a cohort of 36 Romanian, 17 Spanish, 19 Italians and 10 Austrian patients with different stages of liver fibrosis, notably at their early stages. We could design hypoth-eses regarding the putative causal role of liver microbiota in the development of liver fibrosis. We used this database to evaluate the efficacy of Principal Coordinate Analysis (PCoA) to visualize the different liver fibrosis group scores using Wilcoxon-Mann-Whitney statistical tests [START_REF] Cox | Multidimensional Scaling[END_REF][START_REF] Kruskal | Multidimensional scaling by optimizing goodness of fit to nonmetric hypothesis[END_REF]. Eventually, since the overall database of patients issued from different separated cohorts we anticipated some degree of heterogeneity of the overall cohort therefore, we adapted and developed a specific statistical approach i.e. 1 -spectral clustering with fairness. This approach establishes inter-relations between liver microbiota and low scores of liver fibrosis that allowed the identification of the translocated bacteria putatively causal to the disease and independent from the group size, the patient origins and sets of sequencing. Overall, we drew a European microbial profile of patients at early stages of liver fibrosis. /4.21). All participants provided informed consent prior to participation. The patients who gave their consent to perform a liver biopsy during the procedure were eligible. Exclusion criteria were serious liver diseases (eg hemochromatosis, alcoholic fatty liver disease, Hepatits B and Hepatitis C infection, chronic diseases, inflammatory systemic diseases, acute or chronic infections in the previous month, use of antibiotic, antifungal, antiviral drugs, proton-pump inhibitors, anti-obesity drugs, laxatives, excessive use of vitamin D supplementation, fiber supplements or probiotics or participation in a weight loss program or weight change of 3 kg during the previous 6 weeks, pregnancy or breastfeeding, or major psychiatric antecedents; neurological diseases, history of trauma or injured brain, language disorders, and excessive alcohol intake (≥ 40 g/day in women or 80 g OH/day in men) or intravenous drug abuse, and previous bariatric surgery.

Material and Methods

Subjects and Ethics

The cohort consists of 82 Caucasian patients where 34 were diagnosed with fibrosis stage 0 (F 0 ); 37 stage 1 (F 1 ) and 11 stage 2 (F 2 ), as diagnosed from histological analyses of liver biopsies (Table 4.1). The patients suffered from morbid obesity with a mean BMI 42.6 (±7.3). The mean waist circumference was 121.49 (±18.73) in male and 123.23 (±18.26) in female participants.

Liver biopsies and liver fibrosis diagnosis

Liver biopsies were performed during laparoscopic surgical bariatric procedures or via ultrasound guided liver biopsy. No energy devices were used for collecting the samples since hemostasis was done afterwards when the samples were extracted from the abdomen. Ultrasound (US) guided percutaneous liver biopsy (UPLB) was performed in 10 patients. In all patients, antiplatelet drugs and oral anticoagulation therapy was paused 1 week before UPLB was performed. One experienced physician (> 3000 US-exams and > 100 UPLB) performed the US-examinations with the Philips EPIQ 5 R (Philips Corporation, Amsterdam, The Netherlands). UPLB was performed using an 18 G Temno II semi-automatic tru-cut biopsy needle (Cardinal Health, Dublin, Ohio, USA). After UPLB, all patients were monitored for any signs of pain or clinically suspected bleeding by nursing staff over a 6-h period. If no serious complica-tions were evident, all patients would be discharged after the mandatory 6-h observation, a stable blood count and a normal ultrasound examination. All patients were follow-up in 2 weeks to review the results of the histology. All the samples were stored in a sterile container and kept at -80 • C until assayed. Furthermore, NAFLD was confirmed histologically by an independent pathologist.

Clinical assessments

Anthropometric measurement of each subject was performed by trained nurses in the morning after fast-ing for at least 8 h. Body height was recorded to the nearest 0.5 cm and body weight to the nearest 0.1 kg. BMI was defined as body weight (kilograms) divided by the square of body height (meters). Waist cir-cumference was measured in the horizontal plane midway between lowest rib and the iliac crest to the nearest 0.1 cm at the end of a normal expiration repeatedly in men and women by 3 trained nurses on 3 consecutive days. Blood pressure was recorded to the nearest 2 mmHg by a mercury sphygmomanometer with the arm supported at heart level after sitting quietly for 10 min. Fasting plasma glucose was meas-ured after fasting for at least 8 h. A standard oral 75-g glucose tolerance test was performed to measure 2-h postprandial plasma glucose. Hypertension was defined in accordance to the Guidelines of the Europe-an Heart Association or if the subject was taking medication for hypertension. Diabetes was diagnosed when fasting plasma glucose was ≥ 126 mg/dL (7 mmol/L), 2-h postprandial plasma glucose ≥ 200 mg/dL (11.1 mmol/L), and HbA1c ≥ 6.5% or if the subject was taking medication for diabetes.

Biochemical and molecular analyses

Plasma parameters: Biochemical analyses including total fasted plasma glucose, cholesterol, high-density lipoprotein (HDL) cholesterol, plasma liver enzymes i.e. aspartate aminotransferase (AST/GOT), alanine aminotransferase (ALT/GPT), gamma-glutamyl transferase (GGT), hematocrit and leukocytes were determined by Cobas 8000, (Roche, Basel, Switzerland) according to the manufacturer's specification. Elevated liver enzymes were defined as aspartate aminotransferase and alanine aminotransferase. HbA1c was measured by high-performance liquid chromatography (Bio-Rad, Muenchen, Germany) and a Jokoh HS-10 autoanalyzer.

16S rDNA sequencing and bioinformatic analysis

The V 3-V 4 hypervariable regions of the 16S rDNA were amplified by two steps PCR using v1 primers (Vaiomer) and sequenced using MiSeq Reagent Kit v3 (2 × 300 bp Paired-End Reads, Illumina, San Diego, CA, USA) as previously described [START_REF] Luche | Metabolic endotoxemia directly increases the proliferation of adipocyte precursors at the onset of metabolic diseases through a cd14-dependent mechanism[END_REF]. The MiSeq sequences were then analyzed using the bioin-formatics pipeline established by Vaiomer using FROGS v1.4.0 [START_REF] Escudie | Frogs: Find, rapidly, otus with galaxy solution[END_REF]. Briefly, after demultiplexing of the bar-coded Illumina paired reads; single read sequences are cleaned and paired for each sample independently into longer fragments. Operational taxonomic units (OTU) are produced with via single-linkage clustering and taxonomic assignment is performed in order to determine community profiles (generated by Blast+ v2.2.30+ against the Silva v128 Parc databank restricted to the bacterial kingdom).

Linear Discriminant Analysis (LDA) Effective Size (LEfSe)

The bacterial profiles were further compared between the three groups using LEfSe pairwise analysis with an alpha cut-off of 0.05 and an effect size cut-off of 2.0. The bacterial diversity analyses (alpha and beta diversity, MDS ordinations and taxonomic composition barplots) were generated using the 'Phyloseq (v1.14.0), 'vegan' (v2.4.0) and 'ape' (v3.5) packages under R environment v3.3.1. LEfSe analysis was performed on the OTU table using the online Galaxy interface to identify bacterial taxa that were differ-entially abundant in the three liver fibrosis groups [START_REF] Segata | Metagenomic biomarker discovery and explanation[END_REF]. Respective cladograms were generated with genus at the lowest level. Quantitative plots of differential features were generated from genus level percent relative abundance data showing means with standard deviation using GraphPad Prism 6 software. Using the LEfSe algorithm, bacterial taxa that were differentially abundant in analysis of liver fibrosis groups were first identified and tested using the Kruskal Wallis test.

Beta diversity analysis

The bacterial diversity (alpha and beta diversity) was analyzed and represented using the 'phyloseq' (v1.14.0), 'vegan' (v2.4.0), 'ape' (v3.5), and 'ggplot' (3.3.0) packages under R environment v3.5.1 with Chao, Inverse Simpson, Simpson and Shannon as indexes. The alpha diversity statistical significance was de-termined by Wilcoxon rank-test. The beta diversity was calculated for every pair of variables to generate a matrix of distance using Bray-Curtis, Jaccard, Unifrac, and weighted Unifrac indexes. From distance matrices, Multiple Dimension Scale (MDS) and hierarchical clustering were conducted for graphical representation of beta diversity. PERMDISP2 procedure was used for the analysis of multivariate homo-geneity of group dispersions. The Kruskall-Wallis test was performed to compare abundance across the three groups.

Multivariate analysis

To visualize the distribution of patients according to their clinical parameters, we performed a Principal Component Analysis (PCA) using 'FactoMineR' and 'factoextra' R packages. For the study of 16S rDNA diversity, we first filtered the less abundant OTUs to reduce the noise within the matrix before running the PCA. We eliminated those with abundance < 0.01. We then normalized the OTU table by using the Cumulative Sum Scaling normalization followed by a log transformation, using 'mixOmics' package [START_REF] Rohart | mixomics: An r package for 'omics feature selection and multiple data integration[END_REF]. To explore the metagenomic data and identify the largest sources of variation, another Principal Component Analysis was conducted. Also based on the projection of the dataset into a space of lower dimension and originally designed for regression we performed a Partial Least Square Discriminant Analysis (PLS-DA) and its sparse version (sPLS-DA) on the normalized OTU table count to predict and select the most discriminative features in the data that help to classify the samples according to the fibrosis variable (package 'mixOmics'). Since we observed the influence of the metagenomic data on the outcome, we used alternative method of classification such as random forest (package 'randomForest'). The random forest is built from a multitude of different decision trees and classifiers at training time thereby predicting and storing the predicted target outcome.

Cluster graphical analysis

The abundance matrix of OTUs can be modeled by a graph using 'PLNmodels' package under R where nodes represent OTUs and edges interactions between each pair of nodes. We developed an analysis in clusters i.e. the 1 -spectral clustering, implemented in R, a robust variant of the well-known spectral clustering that aims to detect the natural structures of a graph by taking advantage of its spectral properties. The adjacency matrix modeling the variable associations of the graph is used as an input of the 1 -spectralclustering algorithm. In front of the influence of the origin of the cohort on the graphical classification through clusters we applied "fair" technics with k-median clustering objectives. We identified k centers and assign each input point to one of the centers so that the average distance of points to their cluster center is minimized. In the fair-variant, the points are colored while the goal is to minimize the same average distance objective ensuring all clusters to have an approximately equal number of points of each color. This technique called 'fairtree' and developed in python takes as input the desired number of clusters, the desired cluster balance and the normalized table count.

Functional metagenomic prediction

Metagenome inference and predicted functional analysis were initiated by analysis of the OTU clustered 16S sequence count table data and the OTU representative sequences using the PICRUSt2 tool [START_REF] Douglas | Picrust2 for prediction of metagenome functions[END_REF] version 2.3.0b for each sample. The metagenome predic-tion process included four main steps: 1) The input OTU representative sequences were aligned against the PICRUSt2 reference alignment, 2) From this alignment, the input OTU were placed into the PICRUSt2 reference phylogenetic tree, 3) The metagenome functions were inferred by the hidden state prediction method using this phylogenetic tree. During this inference process, the abundance values of each OTU were normalized to their respective predicted 16S rDNA copy numbers and then multiplied by the respective gene counts of the target bacteria, 4) The predicted functions were mapped to the MetaCyc database to determine the minimum set of pathways present in the samples. The resulting core output was a list of enzyme functions (Enzyme Commission numbers) with predicted count data for each sample from step 3 as well as a list of MetaCyc pathways with predicted count data for each sample from step 4.

Data Availaibility Section

MiSeq 16S rDNA sequences were deposited under the primary accession number PRJEB41831 and a secondary number ERP125667 on December 9th 2020 with a release date on the 31st of December 2021.

Results

Graphical classification of the clinical variables by PCA

We aimed at identifying liver 16S rDNA profiles associated with the early onset of fibrosis. We aggregated together a library of liver biopsies from patients from four cohorts of different European countries. We first visualized the distribution of the patients according to the cohorts by performing a principal component analysis using the anthropomorphic and clinical data where the projection of the different clinical variables is represented (Figures 4.1 A,B). The ellipses calculated for each cohort show some degree of differential distribution suggesting that specific environmental factors have influenced the clinical outcomes. Interestingly, the Romania cohort was unifying all cohorts and could be used as a reference. In addition, we could detect numerous outlier patients from each cohort.

It is noteworthy that we voluntarily included all anthropomorphic and biochemical data, even if some were redundant and confounding, to remain within the frame of a non-a priory statistical approach. The age, diabetes and hypertension variables were the main drivers of the F 2 classification while HDL cholesterol and liver enzymes were drivers for the F 1 histological phenotype. These observations are supported by significant ANOVA tests (Table 4.1).

Analyses of the liver bacterial 16S rDNA ecology

To identify whether the graphical differences between the three liver fibrosis scores are associated with a differential liver bacterial DNA signature, we then performed PCA on the OTUs as entries in the database. The analysis using countries as groups shows that the different cohorts poorly overlapped suggesting the existence of specific environmental factors specific of each country cohort (Figure 4.2 A). Using the liver fibrosis scores as groups we could not clearly graphically discriminate the fibrosis scores since the distribution of the patients according to their OTU profiles were too scattered and seemed to be depending upon the largest Romanian cohort (Figure 4.2 B). To analyze differently the putative signatures according to the cohorts or the liver fibrosis scores, we studied the frequencies of the phylum and family taxonomic levels. The barplot analysis shows first a large degree of heterogeneity between all individuals at the phylum level (Figure 4.3) but still, we identified that the liver microbiota of the overall cohort was composed mostly of Proteobacteria, (> 75%) (Figure 4.4). Group comparisons showed that statistical differences were observed between the F 0 and F 1 groups for the Proteobacteria, Actinobacteria and Firmicutes phyla (Figures 4.5 A,B,C). At the family taxonomic level, the most prominent taxa were the Enterobacteriaceae and the Pseudomonadaceae which accounted for more than 50% of the overall taxa (Figure 4.6). Group comparisons showed that the Caulobacteraceae, Flavobacteriaceae and Propionibacteriaceae families were statistically different when comparing F 0 and F 1 (Figures 4.7 A,B,C).

To further identify whether liver fibrosis scores could be characterized by specific signatures we explored indexes of alpha and beta diversity of 16S rDNA in liver tissue. The data show that differences in abundances at the phylum, and family taxonomic levels were also associated with differences of the alpha diversity (Figure 4.8 A,B,C). Notably, the Observed, Shannon and Simpson indexes were significantly different between the F 0 and F 1 groups at the phylum and family levels. In addition to alpha diversity, we analyzed beta diversity and performed a principal coordinate analysis (PCoA) considering distances between variables (using Bray-curtis distance). The PCoA analyses showed that the F 0 group was distant from the two others which suggests a specific 16S rDNA signature (Supplementary Figure 4.9 A,B). It is noteworthy that outlier patients were also detected. Although, when analyzed together the three groups could not be classified clearly. The F 0 group differed graphically from the F 1 , F 2 groups suggesting a specific signature discriminating between F 0 and F 1 , F 2 . To determine if the ellipse centers of the F 0 group differs from the ellipse center of the other groups, a Permutational Multivariate Analysis of variance (PERMANOVA) followed by a Kruskall-Wallis test were performed and found a difference between F 0 and F 1 groups (p < 0.03). Along the same line of investigation, we performed different graphical representations such as heatmaps and Venn diagrams.

Identification of specific bacterial signatures

To identify the variables that are specific to Fibrosis scores we performed a first Venn diagram on the overall set of variables (Figure 4.10 A). Eighty-nine variables were common to all groups and considered as the core of the cohort while 21, 77, and 108 OTUs were specific of the F 2 , F 1 , F 0 groups, respectively. To isolate extremely rare variables and unbalanced distribution between groups we next considered only OTUs with more than 25% of non-zero counts and an average number of counts per group higher than 150 and similarly drew a second Venn diagram. We identified 12, 5, and 5 OTUs specific to F 2 , F 1 , and F 0 scores, respectively (Figure 4.10 B) and (Tables 4.2, 4.3, 4.4). To identify if these specific OTUs could be picked up using another approach we generated a heatmap where each OTUs was positioned while the fibrosis scores was fixed (Figure 4.11 A). We noted that the frequencies of the majority of OTUs equal 0 or are extremely low (< 0.01%) thereby, most of these variables do not bring information. Similarly, a minority of the variables of high frequencies were common to all liver fibrosis groups and did not provide discriminant information neither. Such OTUs could be considered as the core variable of liver microbiota. Conversely, a subset of OTUs could be considered as discriminant that was identified on a different heatmap following the removing of the non-informative OTUs (Figure 4.11 B).

To refine the identifications of the discriminant bacteria we performed a Linear discriminant analysis (LDA) coupled with effect size measurements (Figure 4.12,Figure 4.13,4.14). The data show that most of the discriminant information was identified when comparing between F 0 and F 1 . The Firmicutes, Flavobacteriaceae, Caulobacteraceae and Actinobacteria were specific to F 0 group and the Proteobacteria was specific to F 1 (Figures 4.4 A,B,C, Figures 4.6 A,B,C). On the boxplot the taxa enriched in patients with no fibrosis are indicated with a negative score and mild fibrosis enriched taxa are indicated with a positive score. We performed LEFSe between each score pairs and identified much less differences between F 1 F 2 suggesting that they could have a similar liver microbiota, as suggested in Figure 4.2 B despite the discriminant clinical variables identified in Figure 4.1 B. On these first sets of analyses, the number of fibrosis scores of each patient was too heterogeneous to perform a discriminant analysis (overfitting). As shown on Figures 4.15 A,B, 4.16 A,B they were almost no difference between F 1 and F 2 , therefore we merged F 1 and F 2 scores as F 1/2 group, increasing hence the number of patients of that group.

To validate the pertinence of such strategy we performed a partial least square discriminant analysis i.e. PLS-DA. To select the most discriminative features in the model we used its sparse version sPLS-DA based on a Lasso penalization. The number of variables to be selected per component involved in the visualization is optimized using leave-one-out cross-validation. On the sample plot (Figure 4.9 A), we observe a slight separation of the two fibrosis scores ellipses compared to the unsupervised PCA. From the most discriminant OTUs selected on each sPLS-DA component, a dissociation between the two groups can be visualized using a Clustering Image Map (CIM) (Figures 4.9 B, 4.10 A). The graphs show a clear classification of the patients based on the identified discriminant variables. Eventually, we calculated the ROC curve with all discriminant variables that shows an increased specificity and Altogether, some degree of graphical classification of the liver fibrosis score could be observed using the clinical database and the 16S rDNA database. However, in both instances the individuals appear to be still distributed according to the countries. Therefore, to overcome this issue we developed an ad hoc fairness statistical strategy allowing the classification of variables i.e. OTUs independently from the cohort.

Identification of clusters of cohort-independent 16S rDNA associated with different mild scores of fibrosis

In front of these numerous signatures and the influence of confounding factors such as the impact of the cohort set there is a need to identify clusters of variables specific to each liver fibrosis score but independent of the cohort impact. To this aim we considered three different fair approaches on the overall cohorts and then defined clusters of OTU variables independent from the cohort. The first fair approach consists in identifying principal components from the metagenomic dataset as signatures of the cohorts and removing them to generate a new dataset where no components would be cohort sensitive. To this aim we compared the largest cohort i.e. from Romania to the others. Principal components conditional distributions with respect to the cohorts were visualized (Figure 4.17 A). Then, we removed the principal components the most correlated with the cohorts when the absolute value of Pearson correlation was above a threshold. The remaining non-overlapping components are cohortinsensitive and used to identify the variables associated with the mild fibrosis score. Remarkably, more than 78% of the variation from the original data was still included into the selected principal components suggesting that the discriminant information was only marginally affecting our previous results. On this "fair" dataset we applied the standard random forest classification to predict fibrosis scores. From the variable importance plot, indicating the contribution of the variables to classify the data, we selected the 10 more predictive principal components and identified 3 significantly associated with fibrosis scores (Figures 4.17 The second fair clustering approach directly selects OTUs which are the most influenced by the cohorts and removes them from the analysis. The associated matrix is then modeled by a graph and subjected to a spectral clustering algorithm to which we applied an 1 penalty. The nodes represent OTUs and the edges show interactions between each pair of variables (Figure 4.18 A). Using this novel l 1 -spectral clustering algorithm we identified 5 clusters of OTUs among which 3 were significantly associated with the liver fibrosis scores (Figure 4.18 B).

Eventually, we performed the fair clustering method called "fair-tree". We used the 16S rDNA normalized table count to identify clusters with approximately equal number of patients from each cohort. Two of the three clusters found containing respectively 36 and 97 OTUs, were statistically significant when comparing F 0 versus F 1 scores (Figures 4.18 C,D).

To summarize all the identified OTUs significantly associated with the different low scores of fibrosis, we compiled them in (Tables 4.5,4.6,4.7,4.8) and identified their respective taxa. From the fair principal components identified, we only considered the five OTUs that contribute most to create each of these components. Then, from the Venn diagram we identified common OTUs signatures of low fibrosis scores from standard (sPLS-DA) and fair approaches (fair-tree, random forest, 1 -spectral 4.9) suggesting that most of the discriminant information could be due to these taxa. However, there is still most likely some information that this predominant family could be hiding. We therefore set a new mathematical strategy to exemplify the low frequency and meaningful bacteria by using the TF-IDF (Term frequency-inverse document frequency) approach.

Low frequency bacterial 16SrDNA gene contains classifying information

From the table count of all significant OTUs detected we generated a "word-cloud" (Figure 4.20 A,B) to visualize the most abundant TF-IDF transformed OTU counts, regardless of fibrosis scores when com-pared to those non-normalized. Cluster 2 emerged as the most important discriminant OTU (taxonomic identifaction=Bacteria|Proteobacteria|Gammaproteobacteria|Enterobacteriales |Enterobacteriaceae|Escherichia-Shigella) further confirming the important amount of information contained in the Proteobacteria.

Based on the identified specific signatures the next step was to generate hypotheses regarding their putative mode of action to the induction of the early events of liver fibrosis. We therefore performed predicted functional metagenomics.

Predicted functional metagenome pathways

To identify the pathways and enzymes involved in the early development of fibrosis, we run predicted functional metagenomics algorithms based on the fairness-selected bacterial taxa. The heatmap shows clusters of enzymes that are associated with the F 0 vs F 1/2 fibrosis scores (Figure 4.21 A). Eventually, sPLS-DA showed also a clear discrimination between the F 0 vs F 1/2 fibrosis scores. To evaluate the accuracy and sensitivity of our analyses as potential diagnostic tool, we drew a ROC and quantified the urea under curve with a score of 81.4% of accuracy (Figures 4.21 B,C). We performed a similar analysis on pathways and showed specific clusters also discriminately We then represented an listed all selected enzymes and pathways highly expressed in the two major discriminant components and (Tables 4.10,4.11). Three pathways were highly and negatively associated with the liver fibrosis score of F 1/2 when compared to the F 0 . We identified from the MetaCyc database that the preQ0 biosynthesis (PWY-6703), specific to Enterobacteriaceae such as E. coli, is involved notably in the synthesis of tetrahydrofolate and a class of nucleoside analogues that often possesses antibiotic, antineoplastic, or antiviral activities [START_REF] Dunphy | Separation and identification of menaquinones from microorganisms[END_REF][START_REF] Farrand | Physiological effects of menaquinone deficiency in bacillus subtilis[END_REF]) (Figure 4.24 A). In addition, two other pathways related to glucoryranose (PWY 6737) and glycogen (GLYCOCAT-PWY) degradation were identified probably providing energy to the main preQ0 biosynthesis pathway.

On the other hand, six major metabolic pathways were positively associated with the F 0 score from both components. One involves the glycolysis and pentose phosphate pathway (PWY-6629), while the 5 others are all involved in the menaquinones and demethylmenaquinones pathway (Figure 4.24 B). The low-molecular weight lipophilic components of the cytoplasmic membrane are considered vitamin K2 components that is found in most aerobic Gram-positive bacteria and are the main quinones that function as a reversible redox component of the electron transfer chain, mediating electron transfer between hydrogenases and cytochromes. Altogether the functional metagenomics prediction suggests that gram negative bacteria from the Proteobacteria family composed of preQ0 biosynthesis and glycolytic pathway are signature of F 1/2 fibrosis scores while the vitamin K2 biosynthesis pathway from gram negative bacteria such as Actinobateriaceae [START_REF] Iwata-Reuyl | Biosynthesis of the 7-deazaguanosine hypermodified nucleosides of transfer rna[END_REF][START_REF] Mccarty | The deazapurine biosynthetic pathway revealed: in vitro enzymatic synthesis of preq(0) from guanosine 5'-triphosphate in four steps[END_REF] would be signing F 0 liver fibrosis score.

Discussion

We here report a mathematical approach to identify a bacterial 16S rDNA signature in liver tissue and corresponding putative biochemical pathways in patients with low scores of fibrosis. Our main finding is that even low scores of fibrosis (F 0 vs F 1/2 ) can be classified by biomarkers from the Proteobacteriaceae family within the liver. The second observation is related to the importance of cohort heterogeneity in term of size and data variability which could be major confounding factors that must be taken into account in multi-centric clinical trials or database. We here present a mathematic approach that could help solving this major and common issue.

A gut metagenomics signature of liver fibrosis in humans has been recently described, suggestive of its causal role in the disease [START_REF] Loomba | Gut microbiome-based metagenomic signature for non-invasive detection of advanced fibrosis in human non-alcoholic fatty liver disease[END_REF]. However, such patients where mostly characterized by a high score of liver fibrosis questioning the putative causal role of the liver microbiota in the disease. We here focused our attention on low scores of liver fibrosis to putatively identify causal factors. We identified mostly gram negative bacteria and notably the Proteobacteria as signature of the F 1/2 liver fibrosis scores. Among the families the Proteobacteriaceae, Flavobacteriaceae, and Propionibacteriaceae were discriminating the low fibrosis scores from each other's. They synthesize LPS, a dramatically inflammatory suggesting a pathophysiological role in development of liver fibrosis, probably via the maintenance of a certain degree of immune vigilance. We further refined our analyses and mostly selected Enterobac-teriaceae family from the Proteobacteria phylum suggesting that the liver proinflammation observed dur-ing fibrosis would be due or associated with genera from the Enterobacteriaceae family [START_REF] Cani | Metabolic endotoxemia initiates obesity and insulin resistance[END_REF]. The Enterobacteriaceae encompass the genera Arthrobacter and Acinetobacter. The mechanisms through which such bacteria could induce inflammation might be linked to the unique structures of their LPS or peptidoglycans [START_REF] Farrand | Physiological effects of menaquinone deficiency in bacillus subtilis[END_REF]. Furthermore, since such bacteria are motile with flagella, one could also contemplate that the flagellar proteins are involved in the liver fibrosis process. However, data report that the TLR5 receptor of flagellin is rather associated with protection against metabolic syndrome, putatively ruling out this hypothesis [START_REF] Latino | Fine structure analysis of lipopolysaccharides in bacteriophage-resistant pseudomonas aeruginosa pao1 mutants[END_REF]. Through functional metagenomics production we identified the preQ0 biosynthesis pathway as a signature of F 1/2 fibrosis scores. Such pathway is notably identified from gram negative bacterial such as Proteobacteriaceae [START_REF] Kanehisa | Kegg for integration and interpretation of large-scale molecular data sets[END_REF][START_REF] Iwata-Reuyl | Biosynthesis of the 7-deazaguanosine hypermodified nucleosides of transfer rna[END_REF]. Conversely, the menaquinones and demethylmenaquinones pathways involved in K12 vitamin synthesis were the signature of the F 0 score. They are notably produced by the gram positive Actinobacteriaceae such as Bacillus subtilis [START_REF] Dunphy | Separation and identification of menaquinones from microorganisms[END_REF], therefore coherence with our metagenomics findings.

A major hurdle of aggregation of different cohort altogether is related to the heterogeneity of the size of the groups and of the diversity of the variables considered. Regarding invasive analyses such as liver biopsies the group size at completion of the inclusions could be different from what predicted during the calculation of power of the trial. Eventually, the distribution of the variables to be studied could be highly heterogeneous for a given disease. Altogether, we here faced several statistical challenges which are linked to liver fibrosis. The first major step preceding the microbial analysis was a prefiltering and then an adapted pathway to normalize the data to deal with their sparse nature. The package Mixomics (Le [START_REF] Cao | A sparse pls for variable selection when integrating omics data[END_REF] used for this study recommends CSS normalization on sparse OTU table counts that could prevent the bias included in the TSS normalization. In addition, it includes multivariate methods for microbiome studies and addresses its limits. In addition, we observed a strong impact of the cohort of origin since the largest cohort from Romania could discriminate the patients from the others based on the 16S rDNA OTU variables. The patients could even be classified by cohort when we used the clinical data as entries showing that this issue also has to be taken into account when analyzing the data. Mathematical approaches to overcome this issue are currently being developed however, little has been done regarding the handling of the 16S rDNA data now widely used by the scientific community that addresses the role of microbiota on diseases and notably liver diseases. Therefore, we here developed several approaches of fairness to overcome the classical cohort impact. Eventually, we noticed that two patients from the F 1 groups were distributed with the F 2 group. This ectopic distribution could be due to the extreme BMI (> 55) featuring a specific clinical phenotype. Conversely, a patient from the F 2 group was associated with the F 0 /F 1 distribution. This patient was characterized by his young age (< 40 years old) while the mean age of the F 2 group was of 54 years old.

The statistical approach required to properly analyze microbial data sets needed to be better fitted to the nature of the data. As a preliminary analysis we performed PCoA since better adapted than PCA to dissimilar and sparse data then followed by a sPLS-DA to identify subsets of 16S rDNA that are discriminatory for the liver fibrosis scores. PLS-DA aims to classify a data set according to the values of a qualitative variable by maximizing the covariance between linear combinations of the observed variables and the qualitative outcome. The sparse version, on the other hand, delivers variables per each component, only selected in the OTU dataset, that are the most discriminatory for the liver fibrosis scores. We focused our attention on the identification of the OTU frequencies within and across each group of patients and on the understanding of the importance that OTUs carry within and across the cohort. We found that the data set is mostly populated by a few high frequency OTUs. However, beside the level of information gained form this approach where overrepresented OTUs we identified cannot rule out that some more information could be obtained from OTUs rarely represented. Therefore, some information could be hidden in the low frequency OTUs. To test this hypothesis we introduced a new normalization approach called TF-IDF [START_REF] Wang | Prescription function prediction using topic model and multilabel classifiers[END_REF] originally developed for text mining, to attenuate the effects of the high frequencies OTUs in the data set. Furthermore, aside from the fibrosis scores, it reveals some new predominant taxa at the different taxonomic levels.

Conclusion

In conclusion, the first evidence of the existence of a liver microbiota opens alternate routes for novel therapeutic strategies since specific bacteria could be involved in the process of liver fibrosis. However, to generate information which could serve as a substratum to reach this aim, we here adapted predicted metagenomics and mathematical approaches to the original and novel nature of the tissue metagenomics data set. We here found that these data are constituted of high heterogeneity variables which are dominated by a few high frequency taxa such as Proteobacteria, signature of F 1/2 liver fibrosis scores, and Actinobacteria/Firmucutes, signature of F 0 liver fibrosis scores. These major taxa are masking information residing in the lower frequency taxa. Predicting metabolic pathways from selected 16S rDNA-based taxa revealed a role of folate metabolism in F 1/2 liver fibrosis scores while a role of vitamin K12 biosynthesis was characterizing F 0 liver fibrosis score. Altogether, the combined use of metagenomics, sPLS-DA, TF-IDF and fairness strategies appeared useful since we identified signatures specific to the lower scores of liver fibrosis i.e. at the onset of the disease.

Chapter 5 Fair conditional Partial Least Square 1 Introduction

Machine learning algorithms aim to find statistical patterns in a training dataset that can be generalized to the whole population. Here, the term generalization refers to the model's faculty to adapt and make accurate prediction of the outcome from previously unseen new data, which have been drawn from the same distribution as those used to build the model. Over the last past years, machine learning algorithms became efficient tools for the development of various technological applications, first mimicking human performances and then performing better. They are now widely used in many fields, among the most well-known: finance [START_REF] Trippi | Neural Networks in Finance and Investing: Using Artificial Intelligence to Improve Real World Performance[END_REF], education [START_REF] Mcarthur | The roles of artificial intelligence in education: Current progress and future prospects[END_REF], social media [START_REF] Zeng | Social media analytics and intelligence[END_REF], business [START_REF] Dirican | The impacts of robotics, artificial intelligence on business and economics[END_REF], robotics [START_REF] Dirican | The impacts of robotics, artificial intelligence on business and economics[END_REF][START_REF] Brady | Artificial intelligence and robotics[END_REF] automotive or industry [START_REF] Gusikhin | Intelligent systems in the automotive industry: applications and trends[END_REF][START_REF] Luckow | Deep learning in the automotive industry: Applications and tools[END_REF]. This list is not exhaustive and can be expanded to healthcare [START_REF] Morik | Medicine: Applications of Machine Learning[END_REF][START_REF] Hamet | Artificial intelligence in medicine[END_REF][START_REF] Yu | Artificial intelligence in healthcare[END_REF][START_REF] Sidey-Gibbons | Machine learning in medicine: a practical introduction[END_REF], which raised in recent years a growing interest with a wide potential of applications: radiology [START_REF] Hosny | Artificial intelligence in radiology[END_REF][START_REF] Pisarchik | From novel technology to novel applications: Comment on "an integrated brain-machine interface platform with thousands of channels" by elon musk and neuralink[END_REF], screening [START_REF] Bellemo | Artificial intelligence screening for diabetic retinopathy: the real-world emerging application[END_REF][START_REF] Mckinney | International evaluation of an ai system for breast cancer screening[END_REF], primary care [START_REF] Blease | Artificial intelligence and the future of primary care: Exploratory qualitative study of uk general practitioners' views[END_REF] or disease diagnosis [START_REF] Liang | Evaluation and accurate diagnoses of pediatric diseases using artificial intelligence[END_REF].

As Machine learning algorithms' accuracy highly depends on how the dataset conveys reliable information, handling algorithmic bias is a major concern [START_REF] Parikh | Addressing bias in artificial intelligence in health care[END_REF][START_REF] Decamp | Latent bias and the implementation of artificial intelligence in medicine[END_REF][START_REF] Tat | Addressing bias: artificial intelligence in cardiovascular medicine[END_REF][START_REF] Nelson | Bias in artificial intelligence[END_REF]. In addition, given the high inherent stakes in medical trials, researchers have to rethink how their algorithms are certified before being used at larger scales. To tackle this issue, it is particularly needed to understand where the bias comes from. First it may come directly from models that have been learnt from unbalanced data and will then fit the majority but behave badly for others. In a medical context, bias is particularly stringent since data depend on how it is collected. This is especially true when it comes from distinct cohorts, aggregated into a single one, where the nature of the phenotypes differs. Age, sex, ethnic origin are well known examples of variables that may introduce harmful bias to machine learning models. Yet, as quoted in [START_REF] Dias | Artificial intelligence in clinical and genomic diagnostics[END_REF], bias due to underrepresented population leads to the propagation of non-causal bias and thus bad prediction. It is the case in the diagnosis of Down syndrome which is influenced by the origin of the different populations, as pointed by [START_REF] Lumaka | Facial dysmorphism is influenced by ethnic background of the patient and of the evaluator[END_REF]. It may also come from variables with a higher variance in some populations or the presence of a confounder effect in a subgroup. An example of biased dataset which influences the prediction by detecting spurious correlations or confounders is presented in [START_REF] Esteva | Dermatologistlevel classification of skin cancer with deep neural networks[END_REF] as part of a skin cancer study. Finally, the evolution of the characteristics of the population, leading to shifts in the distribution of the patients, may jeopardize the mid or long term use of the machine learning algorithms, as described in [START_REF] Challen | Artificial intelligence, bias and clinical safety[END_REF]. To tackle the bias issue, new methods were developed in recent years, leading to the emergence of a new field called fairness [START_REF] Pérez-Suay | Fair kernel learning[END_REF][START_REF] Dwork | Fairness through awareness[END_REF][START_REF] Gillen | Online learning with an unknown fairness metric[END_REF][START_REF] Kearns | Average individual fairness: Algorithms, generaliza-tion and experiments[END_REF][START_REF] Kleindessner | Guarantees for spectral clustering with fairness constraints[END_REF], which is attracting increasing attention and is progressing in a wide range of applications (Joseph et al., 2016[START_REF] Joseph | Meritocratic fairness for infinite and contextual bandits[END_REF][START_REF] Berk | Fairness in criminal justice risk assessments: The state of the art[END_REF][START_REF] Brunet | Gender shades: Intersectional accuracy disparities in commercial gender classification[END_REF][START_REF] Siegel | Race-conscious but race-neutral: The constitutionality of disparate impact in the roberts court[END_REF][START_REF] Doherty | Information effect of entry into credit ratings market: The case of insurers' ratings[END_REF][START_REF] Johndrow | An algorithm for removing sensitive information: application to race-independent recidivism prediction[END_REF]. Studying the effect of some variables, which are more susceptible to change across time, and building fair models that mitigate their effect may improve the stability of machine learning algorithms and guarantee their performance. Clinical evaluation of such methods should include a careful study of the performances of the algorithms for each subgroup of the population, as pointed out by [START_REF] Kelly | Key challenges for delivering clinical impact with artificial intelligence[END_REF], and should be understandable by medical experts of the corresponding field. Hence, explainable models should be preferred over black box models.

In this chapter, we faced this problem when statistically analyzing two biomedical datasets: the Florinash genomic and the liver fibrosis microbial datasets, already presented in chapter IV. The common guideline of their statistical studies was to provide a new representation of the datasets that can be used to understand their total variability and, ultimately, to find, within the variables, biological signatures of the patients' liver pathology. Both were preliminary analyzed through an exploratory technique: Principal Component Analysis (PCA) [START_REF] Pearson | Liii. on lines and planes of closest fit to system of points in space[END_REF]. However, the generated graphic representations show some degree of differential distributions, suggesting the existence of specific environmental factors. The florinash dataset includes patients with singular phenotypic profile due to their nationality, whereas the liver fibrosis one stems from two distinct cohorts. As regard the induced bias, we must remove this intrinsic bias before further analysis using an unbiased algorithm. In this biological framework, an algorithm is unbiased when a variable S = (S 1 , . . . , S n ) T , also called protected or sensitive attribute, chosen prior to the analysis, does not interfere in the discovery of the biomarkers associated with the target Y = (Y 1 , . . . , Y n ) T from the observations X = (X j i ) i=1,...,n,j=1,...,p . We thus developed unbiased algorithms based on PCA exploratory analysis and on Partial Least Square (PLS) regression technique [START_REF] Wold | The collinearity porblem in linear regression. the partieal least squares (pls) approach to generalize inverse[END_REF][START_REF] Wold | Soft modelling by latent variables: The non-linear iterative partial least squares (nipals) approach[END_REF] to provide fair representations of the dataset related to the clinical variable of interest.

The most naive approach is to simply not include the biased features in the model but it is not enough to mitigate the bias. In fact, model's accuracy may be impacted by the lack of informative variables as some variables strongly correlated to the omitted features will propagate the bias they contain. In this work, fairness is guaranteed by controlling the independence between the orthogonal features built by standard PCA or PLS and the sensitive variables using a threshold term. The final insensitive transformed variables are selected with an optimal threshold through a cross-validation procedure.

The chapter is organized as follows: in Section 2, we introduce some preliminary concepts about standard PCA and PLS methods. Then, in Section 3 we present in general terms the fair strategy we developed, the theoretical and algorithmic development of the fair methods. As our work was driven by medical projects, we illustrated the efficiency and accuracy of the algorithms through experiments on the datasets.

2 Standard PCA and PLS

Principal Component Analysis

Principal Component Analysis (PCA) is a traditional multivariate statistical method that aims to reduce the dimension of a dataset by extracting few linear combinations of the variables without losing too much information [START_REF] Pearson | Liii. on lines and planes of closest fit to system of points in space[END_REF][START_REF] Jolliffe | Principal component analysis: a review and recent developments[END_REF]. This popular technique for dimensionality reduction is also commonly used as a preliminary step before applying any regression or classification technique to better visualize the variation that exists in the dataset.

Consider a system of p quantitative variables X 1 , . . . , X p and n observations X j := (X j 1 , . . . , X j n ) T , j ∈ {1, . . . , p} of these p variables. PCA creates a set of decorrelated variables, called principal components, by orthogonally transforming the p variables. Therefore, a set of p vectors (w h ) h∈{1,...,p} , called principal factors, is firstly found such that the maximal variance in the original dataset is preserved. Principal factors, are successively built as solutions of the following optimization problems: ∀h ∈ {1, . . . , p} , w h = arg max w =1

Var(Xw),

(5.1) under the constraint that for all l ∈ 1, h -1 , w T h X T Xw l = 0. Then, a set of p principal components T 1 , . . . , T p , defined as linear combinations of the centred variables X j , are built: ∀h ∈ {1, . . . , p} , T h = Xw h . Note that the principal components are ordered from the most informative T 1 to the least informative T p . One significant challenge in using PCA technique is to select the optimal number of principal components. To adress this issue, a large number of graphical methods have been developed in recent years. In what follows, we will note f PCA (w) = V ar(Xw) the objective function to optimize in PCA (equation 5.1).

Partial Least Square

Partial Least Square (PLS) is a technique closely related to PCA. While PCA reduces the dimension of a single dataset by iteratively projecting the data onto components of maximum variance, PLS reduces the dimension of a dataset by projecting the data onto components of maximum covariance with a second data set. It is a supervised technique developed especially for regression, discrimination and classification [START_REF] Wold | The collinearity porblem in linear regression. the partieal least squares (pls) approach to generalize inverse[END_REF][START_REF] Wold | Soft modelling by latent variables: The non-linear iterative partial least squares (nipals) approach[END_REF]. For an application to large omics datasets, we refer to the mixOmics package (Le [START_REF] Cao | A sparse pls for variable selection when integrating omics data[END_REF], dedicated to the development of multivariate methods. Note that a mathematical framework has been provided by [START_REF] Blazère | Pls: a new statistical insight through the prism of orthogonal polynomials[END_REF].

Fair conditional independence strategies

In this section, we present two direct fair strategies based on PCA/PLS techniques that aim to find fair components, containing enough variance to predict a clinical outcome. In the following subsections, we consider a response Y ∈ R n , p unprotected variables X 1 , . . . , X p and one sensitive variable S ∈ {0, 1} n inducing bias.

Overall strategy

Fairness is guaranteed in a post-processing step by imposing independence between the constructed orthogonal features and the sensitive variable using a correlation measure (see Figure 5.1 for a global overview). From all components created by PCA/PLS, we aim to select a set of variables both decorrelated to the variable causing bias and related to the clinical response. To this end, each correlation coefficient is thresholded to keep the less dependent features. The optimal threshold is selected using a cross validation procedure on the basis of the two following criteria:

-the MSE when predicting the outcome from a fitted linear model to guarantee a good prediction, -the total explained variance of the selected components to restitute enough variability of the original dataset. 2018) who showed how to find a single d-dimensional representation for all groups present in the input data and inducing bias, for which the maximum variance is optimized in a balanced manner. It is guaranteed by maximizing the minimum variance over the groups of the projection to a d dimensional subspace. [START_REF] Olfat | Convex formulation for fair principal component analysis[END_REF] published another version of fair PCA to reduce the projected data's dependence on a sensitive attribute through a Semi-Definite Programming. Here, we developed a fair algorithm, called fair conditional PCA algorithm, that directly uses the set of principal components to impose fairness.

A global overview of this algorithm we implemented to select principal components independent from the sensitive variable is presented in Algorithm 6. Fair conditional PCA algorithm relies on a standard PCA technique where a set of p principal components T 1 , . . . , T p is created. We choose to quantify fairness using a threshold µ applied on the correlation ratio η 2 T h /S [START_REF] Fisher | Statistical methods for research workers[END_REF] between each component T h and the sensitive variable S, which enables to characterize the association between a quantitative and a qualitative variable. In this study, the sensitive variable can only take K = 2 values. Hence, the empirical correlation ratio is defined as:

η 2 T h /S = K k=1 n k (T h k ) 2 n i=1 (T h i -T h ) 2 ,
where n k is the number of observations in category k and T h k is T h mean in category k. Then, components with a correlation coefficient C h (h ∈ {1, . . . , p}) lower than a threshold µ are preselected. Once 80% of the cumulative explained variance of the pre-selected components is reached, the remaining pre-selected components are removed. To select the best sample set of fair principal components that are representative of the original variabilities and good predictors of the outcome, we calibrated the threshold coefficient using a cross-validation procedure. Obviously, PCA does not deal with prediction, but fitting a linear model between the fair components and the clinical outcome enables to quantify the predictive power of the fair components. Fair components are then identified by determining an optimal threshold as the one that minimizes the MSE (or 10% of the minimum when the MSE always decreases with the number of components) after predicting the clinical outcome. In addition to the predictibility of the final fair components, we checked their representativity of the variance by calculating the cumulative explained variance. In summary, the final insensitive principal components are selected by computing standard PCA algorithm, thresholding the correlation ratio and predicting the outcome.

Algorithm 6 Fair conditional PCA algorithm Require: Centered matrix X of predictors, range of penalty coefficients µ, sensitive variable S.

1: Solve arg max

w =1 f P CA (w), 2: Extract principal components T h = Xw h , h ∈ {1, . . . , p}, 3: Compute the correlation ratio C h = η 2 T h /S , h ∈ {1, .
. . , p}, 4: Run a cross-validation procedure for calibrating the threshold µ, 5: Return μ the optimal penalty coefficient, 6: Select the fair components P h , h ∈ {1, . . . , k}, 7: Output: (P h ) 1≤h≤p the fair principal components.

Fair PLS based on correlation measure

In recent years, fair regression has been extensively studied. It deals with the prediction of a target Y while guaranteeing the notion of fairness with respect to a protected variable S. Several research studies aim at finding predictors that exhibit some form of independence from the protected variable. [START_REF] Kamishima | Fairness-aware classifier with prejudice remover regularizer[END_REF] seeks to fit a probabilistic model that satisfied statistical independence while [START_REF] Pérez-Suay | Fair kernel learning[END_REF] uses the notion of correlation in a RKHS to capture statistical independence. Finally, gouic et al. (2020) seeks a fair estimator using optimal transport. Yet, there have been no or really few research on the notion of fair Partial Least Square.

As Partial Least Square technique is somewhat similar to Principal Component Analysis, we developed an extension of the Fair conditional PCA algorithm (see Algorithm 7). This new algorithm, called fair conditional PLS algorithm relies on standard PLS technique that creates latent variables T 1 , . . . , T p which variance and correlation with the clinical outcome is maximal. Fairness is then imposed by thresholding the correlation ratio between each latent variable and the sensitive variable. Then, variables with a correlation coefficient C h (h ∈ {1, . . . , p}) lower than a predefined threshold µ and with the higher explained variance are chosen. To select the optimal number of latent variables independent from the sensitive variable i.e. to calibrate the optimal threshold coefficient, we used a cross-validation procedure.

Algorithm 7 Fair conditional PLS algorithm Require: Centered matrix X of predictors, outcome Y , range of penalty coefficients µ, sensitive variable S. 1: Solve arg max w =1 f P LS (w), 2: Extract latent variables T h = Xw h , h ∈ {1, . . . , p}, 3: Compute the correlation ratio C h = η 2 T h /S , h ∈ {1, . . . , p}, 4: Run a cross-validation procedure for calibrating the threshold µ, 5: Return μ the optimal penalty coefficient, 6: Select the fair components P h , h ∈ {1, . . . , k}, 7: Output: (P h ) 1≤h≤p the fair components.

Application to real data

Florinash dataset

The first medical dataset we used belongs to the project FLORINASH, which proposes an innovative research concept to address the role of intestinal microfloral activity in Non-Alcoholic Fatty Liver Disease (NAFLD). Hepatic steatosis, often observed in obese and diabetic patients, is a preliminary stage to NAFLD. The studied cohort [START_REF] Hoyles | Molecular phenomics and metagenomics of hepatic steatosis in non-diabetic obese women[END_REF] is made of obese patients featured with hepatic steatosis, a disease that results from an accumulation of genetic disorders. In this context, our objective is to find genes related to one indicator of the insuline-resistance, called clamp, that will help us to better understand the development of the disease.

The florinash dataset includes a 57×17827 gene expression matrix from the liver of 57 patients who come from Spain and Italy, of which 13370 genes are retained for this analysis after preprocessing and filtering. As this dataset gathers patients with two different profiles (genetically and phenotypically), it introduces a bias in the analysis.

Let X be the gene expression matrix, Y the clamp variable and S the following variable:

∀i ∈ {1, . . . , n} , S i = 0 if patient i comes from Italy 1 if patient i comes from Spain.

To overcome the bias issue, we applied algorithms 6 and 7 to identify, from the genomic dataset, country-insensitive components linked with the insuline-resistance variable.

Results

This section is dedicated to experimental studies conducted on florinash dataset to assess numerical performances of fair conditional PCA/PLS algorithms. We show how to select the final fair components, by finding a reasonable trade-off between fairness and high explained variance through a cross-validation procedure to select the optimal fair threshold µ among a large range of values. The choice of µ is of great importance because taking µ too large will lead to the selection of almost the total number of components. The results of fair PCA experiments can be visualized in Figure 5.2 which represents the Mean Squared Error, the number of components removed after thresholding and the explained variance according to the threshold values. These figures reflect the following tendency: a high threshold leads to a better prediction of the outcome and a high explained variance. With the trade-off introduced in the previous subsections, we get an optimal threshold of µ = 0.05, marked with a vertical line, which enables to select twenty principal components explaining 50% of the variance of the original data and predicting the outcome with a MSE equals to 0.75.

To assess the independence between the selected components and the sensitive variable, we can compare the plots of Figure 5.3, representing the marginal distributions conditioned to the country variable of the first four unfair and fair principal components selected. Note a separation between the marginals of the unfair components whereas this effect is attenuated or nonexistant for the fair components. We can observe the same tendency on the individuals factor maps (see Figure 5.4) where the ellipses, representative of the patient nationality, overlap on each other when the individuals are projected on the fair principal components. Obviously, these results are directly dependent on the value of the optimal threshold, which can hardly be decreased as it would lead to bad predictions of the outcome. The results of the fair PLS experiments can be visualized in Figure 5.5 representing the Mean Squared Error, the number of components removed after thresholding and the associated explained variance according to the threshold values. We get almost the same graphical trends than in fair PCA. The optimal threshold, marked in a red line, minimizes the MSE value with about 30 fair latent variables, explaining 50% of the variance of the original data and predicting the outcome with a MSE around 0.80. Note that standard PLS prediction accuracy on this dataset is already relatively low. We can deduce that it does not seem to be suitable for prediction.

Liver fibrosis dataset

Liver fibrosis is closely related to diabetes, obesity and cardiovascular diseases. It occurs when healthy tissues of the liver become too scarred to work well. Recent studies proved that gut microbiota can help understanding the development of such metabolic diseases. In this context, through a statistical analysis, our objective is to understand the pathological impact of gut microbiota on the development of liver fibrosis.

Fibrosis dataset includes a 82 × 411 microbial data count matrix from the liver of 82 patients of which 411 bacteria are retained for this analysis after preprocessing and filtering. This dataset stems from the fusion of two sets of analysis. The first set is derived from a Romanian cohort whereas the second one collected clinical and microbial information on Italian, Spanish and Austrian patients. From a primary analysis, we realized that the fusion of both cohort thereby introduced a bias in the study.

Let X be the microbial data count matrix, Y and S the following variables:

∀i ∈ {1, . . . , n} , Y i = 0 if patient i is in the early onset of liver fibrosis 1 if patient i is in a more advanced stage of liver fibrosis ∀i ∈ {1, . . . , n} , S i = 0 if patient i comes from Romania 1 otherwise.

To overcome this bias, we apply algorithms 6 and 7 to identify, from the microbial dataset, country-insensitive principal components that are linearly linked with the binary fibrosis variable.

Results

This section is dedicated to experimental studies conducted on liver fibrosis microbial dataset to assess numerical performances of fair conditional PCA/PLS algorithms where the clinical outcome to predict is a binary variable. The results of the fair PCA experiments can be visualized in Figure 5.6 representing the Mean Squared Error, the number of components removed after thresholding and the associated explained variance according to the threshold values. These figures reflect that the MSE decreases until a certain threshold is reached. Hence, the optimal threshold that guarantees good predictions of the outcome and good representativity of the variabilities is the one that minimizes the MSE, represented with a red line. About 15 fair principal components explaining 30% of the variance of the original data and predicting fibrosis with a MSE around 0.3 were selected. Then, to assess the independence between the selected components and the sensitive variable, we compared the plots of Figure 5.7 representing the marginal distributions conditioned to the country variable of the first two unfair and fair selected principal components. Note a separation between the marginals of the unfair components whereas this effect is attenuated or nonexistant for the fair components. We can observe the same tendency on the individuals factor maps (see Figure 5.8) where the ellipses, representative of the patient nationality, overlap on each other when the individuals are projected on the fair principal components. Consequently, the selected principal components are effectively independent from the countries. The results of the fair PLS experiments can be visualized in Figure 5.9 representing the Mean Squared Error, the number of components removed after thresholding and the associated explained variance according to the threshold. We get almost the same graphical trends than in fair PCA. The minimal MSE value is reached with a threshold in the range of 0.025 where about 30 fair latent variables are selected, explaining around 40% of the variance of the original data and predicting the outcome with a MSE lower than 0.30. Note that this fair PLS strategy seems to be better suited to this dataset than the PCA one as the selected country insensitive components represent almost half of the dataset variability and guarantee a low MSE. However, when comparing unfair and fair PCA/PLS procedures, we can observe a loss of accuracy, in terms of MSE, in the order of 0.1.

Results of fair PLS based on correlation measure

Conclusion

In this chapter, we examined two specific biomedical datasets which first exploratory analysis reflected a bias. To remove the underlying adverse effects, we proposed two fair learning approaches based on standard dimension reduction techniques (PCA and PLS) that enable to select components that restore the variability of the original data while limiting the bias effect. Fairness is imposed as a threshold to control the independence between the PCA and PLS components and the sensitive variable causing bias. In order to select the optimal number of unbiased components, we evaluated the performances of each algorithm in terms of Mean Squared Error. In fact, in such biomedical context, we aim to identify relevant variables or groups of interacted variables involved in the development of a disease. Through numerical experiments and comparisons, we show the effectiveness of the fair PLS strategy based on correlation measure and identify a set of latent variables related to the indicator of the disease but decorrelated with the sensitive variable. From a methodological perspective, an interesting research direction would be to impose fairness as an in-processing step of the standard PLS technique. Our goal would be to directly create fair components by penalizing PLS optimization problem using independence criteria between the conditional distributions such as covariance, Hilbert-Schmidt Independence Criterion (HSIC) [START_REF] Mooij | Regression by dependence minimization and its application to causal inference in additive noise models[END_REF] or Wasserstein-2 distance [START_REF] Villani | Optimal transport: old and new[END_REF] (see Figure 5.10). The primary developments of these strategies are detailed hereafter (Section 5).

Further developments

In this section, we provide in-processing fair approaches based on standard PLS method.

PLS under fair constraint based on covariance measure

The first fair approach aims to directly create fair components by penalizing PLS objective function (Equation (5.2)) using a covariance measure. It is somehow an extension of [START_REF] Zafar | Fairness contraints: A flexinle approach for fair classification[END_REF] which also introduced a fair covariance measure in a constraint-based framework but incorporated into a classifier formulation.

Let S be the sensitive variable and T the latent variable scores matrix output of PLS algorithm. Then, the covariance between T and S can be written as: We recall that we aim to find an optimal subspace where the latent variables are independent from the sensitive feature S, i.e. we are looking for T such that E[T |S = 1] -E[T |S = 0] is very small, and such that each T is a good predictor of the outcome. Hence, we must find a trade-off between fairness, i.e. low covariance coefficients and good predictions. We are therefore looking for principal factors (w 1 , . . . , w p ) defined as: ∀h ∈ {1, . . . , p} w h = arg max w =1

, {f PLS (w) -λ Cov(Xw, S)} , (5.3) where λ is a positive penalty parameter.

5.2 PLS under fair constraint based on Hilbert-Schmidt Independance Criterion

In this subsection, we introduce another fair approach based on Hilbert Schmidt independence criterion. Suppose (X , Y) are drawn from distribution P X ,Y , where X ∈ R p is a p-dimensional explanatory variable and Y ∈ R is a response variable. Denote by ∀i ∈ {1, . . . , n} , X i := (X 1 i , . . . , X p i ) a p-dimensional vector of variables.

Hilbert Schmidt Independence Criterion

Hilbert Schmidt Independence Criterion (HSIC) measures allow to simultaneously take into account many forms of dependence between two random variables X and Y relying on Reproducing Kernel Hilbert Spaces (RKHS). The reader can refer to [START_REF] Aronszajn | Theory of reproducing kernels[END_REF] for a complete bibliography on RKHS spaces.

Definition 7 Let X an arbitrary set and H be a Hilbert space of real-valued functions on X with a scalar product <, > H . The Hilbert space H is said to be a RKHS if, for all X i in X the application h → h(X i )∀h ∈ H is a continuous linear form.

The particularity and interest of RKHS spaces is the Riesz representation theorem. This representation consists in associating to each value of the starting set, a function in the RKHS. Each element is then represented by a random functional variable belonging to a space having good properties.

Proposition 7 Let H be a RKHS space associated to a set X and with a scalar product denoted by <, > H . Then, for all X i in X there is a unique φ in H such that h(X i ) =< h, φ > H , for all h ∈ H.

The HSIC characteristics depend entirely on the choice of this RKHS. This choice consists in associating a mapping function

φ : X -→ H X i -→ φ(X i ).
that assigns to each element of the domain a representative functional φ in the RKHS and a scalar product that defines the nature of the relationships between the representatives and so between the elements of the domain. The application that defines this scalar product is called kernel and is defined as follows:

K : X × X -→ X X i , X j -→ K(X i , X j ).

Let T h be the latent variables created by standard PLS method and defined as T h = Xw h , ∀h ∈ {1, . . . , p}. Fair PLS adresses the problem of finding the subspace Xw such that the dependence between the projected data Xw and the outcome is maximized but the dependence between the projected data and the outcome is minimized. In order to measure this independence, we use the Hilbert-Schmidt Independence Criterion.

Hence, we need to minimize HSIC(P (Xw,S) , H K , H L ) = 1 (n-1) 2 tr(KHLH), where K and L are the kernels associated to Xw and S. Let us consider linear kernels i.e. ∀x, x ∈ X , K(x, x ) = x T x . Thus, we need to find the following principal factors (w 1 , . . . , w p ): ∀h ∈ {1, . . . , p} , w h = arg min w =1 tr(w T X T HLHXw).

The optimal solutions to this optimization problem is U := [v 1 , . . . , v p ] where p denotes the final number of latent variables created, v 1 , . . . , v p are the eigenvectors associated to the k smallest eigenvalues λ 1 , . . . , λ p obtained by eigen decomposing the symmetric and real matrix Q = X T HLHX.

This measure is added as a regularization term to the equation (5.2). Hence, we look for principal factors w 1 , . . . , w p defined as: ∀h ∈ {1, . . . , p} , w h = arg max w =1

Cov(Xw, Y ) -µ tr(w T X T HLHXw),

(5.5)

where µ > 0 is the regularization parameter.

PLS through Wasserstein-2 distance

In this subsection, we introduce a third in-processing fair approach based on Wasserstein-2 distance. We sought to impose independence between PLS features with respect to the protected variable S i.e. that either the distribution of the features, or its conditional distribution does not depend on S. Hence, fairness quantification can be implemented using distance between conditional distributions. The Wasserstein distance appears as an appropriate tool for comparing probability distributions. We refer to [START_REF] Villani | Optimal transport: old and new[END_REF] and references there in for more details on Wasserstein's distance.

Our last regularization strategy consists in ensuring that the Wasserstein-2 distance between the conditional distributions µ Xw|S=0 and µ Xw|S=1 remain close to each other. Hence, we need to find principal factors (w 1 , . . . , w p ) defined as: ∀h ∈ {1, . . . , p} , w h = arg max w =1

Cov(Xw, Y ) -µ W 2 2 (µ Xw|S=0 , µ Xw|S=1 ), (5.6)

where µ > 0 is a regularization parameter.

In [START_REF] Risser | Using wasserstein-2 regularization to ensure fair decisions with neural-network classifiers[END_REF], such a regularization term is approximated by calculating the gradient of the Wasserstein distance between the distributions of the two groups.

Chapter 6 Conclusion 1 Synthesis

This thesis, raised from the MATHOBIOMIX project, is devoted to the study and the development of statistical methods and cluster-based algorithms of specific biomedical datasets. In collaboration with the Institute for the Study of Cardiovascular and Metabolic Diseases (I2MC), this project fits specifically within the framework of liver based diseases which we strive to understand the biological foundations using genomic and metagenomic databases.

Chapter II and Chapter III result in the development of graph clustering algorithms dedicated to the detection of densely connected variables in complex and respectively noisy networks, charcateristic of real biological systems. The associated scripts give rise to the development of two R packages.

In Chapter IV, as part of a medical research project, we conducted a statistical study on metagenomic datasets characterizing patient in their early development of liver fibrosis disease, which confirms the need for specific tools due to the nature of the features under consideration. An adapted strategy combining a variety of exploratory, clustering and predictive methods are thus proposed to reveal the main biological interactions implied in the development of the disease.

While the datasets, submitted to statistical studies, highlighted undesirable effects on standard methods, we developed, in Chapter V, fair learning approaches to limit the bias. This study could be extended in the future in order to incorporate differently the limiting factor and eventually, to find the best model improving standard learning procedures under bias. Algorithms, 14(2), 66,

  FIGURE 1.1 -Progression de la Non-Alcoholic Fatty Liver Disease (NAFLD)

FIGURE 1

 1 FIGURE 1.2 -Structure de la molécule d'ADN et d'ARN

  FIGURE 1.3 -Stratégie globale a) et ciblée b) de séquenc ¸age en métagénomique.

FIGURE 1

 1 FIGURE 1.5 -Analyse en Composantes Principales

  4.1.1 Caractéristiques d'un grapheUn graphe G(V, E) est la donnée d'une ensemble V fini de sommets et d'un ensemble E de couples de sommets (u, v) ∈ V 2 . S'il existe une notion d'ordre au sein de ces couples, on parle de graphe orienté, sinon on parle de graphe non orienté (Figures 1.6 a) et c)). De plus, si un poids est associé à chacune des arêtes composant le graphe, on dit que le graphe est pondéré (Figure 1.6 b)). Définition 4 -Une paire (u, v) ∈ V 2 est appelée une arête, ou un arc, et est représentée graphiquement par un trait (Figure 1.6 a)), ou une flèche (Figure 1.6 b)), -Deux sommets u et v sont adjacents si ils sont reliés par une arête, -Une boucle est une arête reliant un sommet à lui-même, -Le degré d(u) d'un sommet u ∈ V est le nombre d'arêtes incidentes à ce sommet, -Le chemin d'un sommet u vers un sommet v est une séquence (s 0 , . .

  FIGURE 1.6 -Exemple de graphe non orienté a), b) non orienté pondéré et orienté c).

FIGURE 1

 1 FIGURE 1.9 -Exemple d'arbre hiérarchique indicé

  (0, α). Three different type of networks were simulated using this protocol with different parameterizations. (a) The first type of network consisted in simulating n = 100 observations of 40 variables with K = 2 clusters of 20 variables. The additional noise was simulated with α ranging from 0.25 to 1.5. (b) The same protocol was used for the second type of networks, but K = 5 clusters of 7 variables were simulated. (c) The third kind of networks consisted in varying the number of the observations from n = 5 to n = 30, with K = 5 clusters having 50 to 100 variables.

Figure 2

 2 Figure 2.1 -Boxplots of the scores R (see Section 4.1.2) obtained on simulated datasets using the standard and the greedy CORE-clustering algorithms as well as the Louvain and the Spectral clustering algorithms. A score of 1 reflects a purely accurate detection of the simulated clusters and the lower this score, the lower the accuracy. The boxplots were obtained by reproducing 30 times the procedure of Section 4.1.1. (a) Two simulated clusters with noise levels ranging from 0.25 to 1.5. (b) Same as (a) with five simulated clusters. (c) Five clusters simulated using 30, 15, 10 and 5 observations and a noise level of 0.5.

Figure 2

 2 Figure 2.2 -CORE-clusters obtained using Alg. 2 on the yeast dataset of Spellman et al. (1998) and the granularity coefficient τ = 30. CORE-clusters containing 30 to 59 variables are then estimated.

Figure 2

 2 Figure 2.3 -Results obtained using the greedy CORE-clustering algorithm on the U.S. road network (≈ 2.4 × 10 7 nodes) in about one minute. (a and b) CORE-clusters obtained using τ = 5 × 10 4 and ξ = 10 -3 represented in the U.S. and in New-York City urban area. (c) CORE-clusters obtained using τ = 1 × 10 5 and ξ = 10 -3 represented in New-York City urban area only. (d) CORE-clusters obtained using τ = 5 × 10 4 and ξ = 0.5 represented in New-York City urban area only. Background: Google maps. Clusters represented using MI Map Tools: GeoPlotter .

  88) 1. (0.82) 1. (0.79) 1. (0.80) 1. (0.40) 1. (0.75) 1. (0.59) 1. (0.79) µ = 0.4 1. (0.87) 1. (0.79) 1. (0.78) 1. (0.83) 1. (0.46) 1. (0.80) 1. (0.75) 1. (0.84) µ = 0.6 1. (0.88) 1. (0.79) 1. (0.78) 1. (0.84) 1. (0.44) 1. (0.74) 1. (0.88) 1. (0.78) µ = 0.8 1. (0.92) 1. (0.72) 1. (0.75) 1. (0.69) 1. (0.47) 1. (0.62) 1. (0.93) 0.96 (0.71) µ = 0.9 1. (0.92) 1. (0.59) 1. (0.70) 1. (0.48) 1. (0.49) 1. (0.44) 1. (0.93) 0.89 (0.59)

  Propositions 3 and 4 below show that the connected components indicators (1 C i ) i∈{1,...,k} are solutions of 0 -minimization problems. Proposition 3 The minimization problem arg min

Figure 3

 3 Figure 3.2 -Clustering results, in terms of NMI, of the 1 -spectral clustering (in red), SpectACl (in green) and spectral clustering (in blue) algorithms applied on perturbed graphs for perturbations ranging from 0.01 to 0.5 and the associated 50% confidence intervals.

Figure 3

 3 Figure 3.3 -At the top, clustering results, in terms of NMI, of the self-tuning 1 -spectral clustering (in red and in green, where only classified nodes are taken into accounts) and the Self-Tuning Spectral Clustering (in blue) algorithms applied on perturbed graphs for perturbations ranging from 0.01 to 0.5. At the bottom, the associated estimated number of clusters of both methods across the 100 perturbed versions of the graphs.

Figure 3

 3 Figure 3.4 -From left to right: (a) the co-regulatory network representing co-regulations between genes in KIRC, (b) the components discovered by applying the 1 -spectral clustering algorithm and highlighted with different colors, (c) the five components we particularly focus on.

  Figure3.5 -Boxplots representing the association between the cluster expression (averaged expression across all genes from the cluster) and kidney cancer grades, which range from 1 (low grade) to 4 (high grade), grade x indicating that the grade could not be evaluated.

Figure 3

 3 Figure 3.6 -Kaplan-Meier curves representing the association between high/low cluster expression and survival.

vFigure

  Figure 2, Champion et al

  Figure 2, Champion et al

Figure 4

 4 Figure 4.3 -Visualization of liver 16S rDNA sequences by principal component analyses according to countries and fibrosis scores. Barplot depicting the frequencies of liver microbial composition of each patient at the phylum level for the overall cohort (total) or according to the fibrosis scores (F 0 , F 1 , F 2 ).

  Figure 2, Champion et al

Figure 4

 4 Figure 4.5 -Mean frequencies of discriminating taxa. Boxplot representing the frequencies of A Proteobacteria, B Actinobacteria C Firmicutes phyla throughout two groups of liver fibrosis scores (red= F 0 , green= F 1 , blue= F 2 ).

Figure 4 Figure 4 Figure 4 Figure 4 Figure 4 Figure 4 Figure 4

 4444444 Figure 2, Champion et al

  B,C,D).

Figure 4 FigureFigure 4

 44 Figure 3, Champion et al

FigureFigure 4 Figure 4

 44 Figure 4, Champion et al

  Figure 4, Champion et al

Figure 4 .

 4 Figure 4.19 -Discriminant analyses of the 16S rDNA OTUs variables using fairness strategies. Venn diagram depicting the liver microbial taxonomies of common OTUs identified by standard (sPLS-DA) and fair approaches (fairtree, random forest, 1 -spectral clustering) as signatures of low fibrosis scores (green= sPLSDA, red= fair algorithms).

Figure 5 ,

 5 Figure 5, Champion et al Figure 4.20 -Identification of clusters by wordclouds representation with or without TFIDF normalization. Wordclouds representing taxa of all significant bacteria according to A, their frequencies or B, after TF-IDF normalization. The size of the name of bacteria is proportional to the frequency of the cluster in the cohorts.

  Figure 6, Champion et al

FigureFigure 4

 4 Figure 6, Champion et al

FigureFigure 4 Figure 4

 44 Figure 6, Champion et al

  Figure 5.1 -Fair post processing procedure

Figure 5

 5 Figure 5.2 -1) Mean Squared Error of the fair predictions, 2) Total explained variance of the fair principal components and 3) Number of components removed after thresholding, depending on the different threshold values.

Figure 5

 5 Figure 5.3 -Conditional distributions of first unfair and fair principal components with respect to the country

  Figure 5.5 -1) Mean Squared Error of the fair predictions, 2) Total explained variance of the fair latent variables and 3) Number of variables removed after thresholding, depending on the different threshold values.

  Figure 5.6 -1) Mean Squared Error of the fair predictions, 2) Total explained variance of the fair principal components and 3) Number of components removed after thresholding, depending on the different threshold values.

Figure 5 . 7 -

 57 Figure 5.7 -Conditional distributions of the first two unfair and fair principal components with respect to the country

  Figure 5.9 -1) Mean Squared Error of the fair predictions, 2) Total explained variance of the fair latent variables and 3) Number of variables removed after thresholding, depending on the different threshold values.

  Figure 5.10 -Fair in processing procedure

  Cov(T, S) = E[T S] -E[T ]E[S]= pE[T |S = 1] -p E[T ],where E[S] = p andE[T ] = E[E[T |S]] = p E[T |S = 1] + (1 -p) E[T |S = 0]. Hence, Cov(T, S) = p E[T |S = 1] -p 2 E[T |S = 1] -p(1 -p) E[T |S = 0] = p (E[T |S = 1](1 -p) -(1 -p) E[T |S = 0]) = p(1 -p) (E[T |S = 1] -E[T |S = 0]).

  Champion, A.C. Brunet, R. Burcelin, J.M. Loubes, L. Risser (2021). Detection of Representative Variables in Complex Systems with Interpretable Rules Using Core-Clusters.

  p}. Ce chemin (s 0 , . . . , s p ) forme un cycle si s 0 = s p .-Un graphe non-orienté est dit simple s'il est sans boucle et s'il admet au maximum une arête entre deux sommets, -Un sous-graphe induit de G est un graphe

  77 0.85 0.77 0.77 0.18 0.25 0.15 0.19 G 2 0.77 0.85 0.77 0.77 0.18 0.25 0.40 0.47 Table 2.1 -Coherence of the CORE-clusters tested Appendix 6.1 on G 1 and G 2 . Corresponding representative variables are given table 2.2.

Table 2 .

 2

			11	3	11	5	11	6	11
	G 2	4	11	3	11	5	11	3	11

2 -Representative variables of the CORE-clusters tested Appendix 6.1 on G 1 and G 2 . Corresponding coherences are given table 2.1.

  Table 2.1 the coherence of four pairs of tested CORE-clusters on G 1 , G 2 and their maximal spanning trees. Corresponding estimated representative variables are given in Table 2.2. The tested CORE-clusters are: (Test 1) The reference blocks of variables S 1 Ref and S 2 Ref , i.e. S 1 T 1

Table 2 .

 2

3 -Portion of representative variable estimates contained in the proper reference block of variables (main value) and corresponding to the ground truth representative variables (between brack-

Table 3

 3

				1	1	1	0.99	0.99
	0.25	1	0.92	1	1	0.71	0.99	0.88
	0.5	0.99	0.65	1	0.76	0.51	0.96	0.40
	0.1 0.01	1	0.99	1	1	1	0.99	1
	0.1	1	0.98	1	1	0.96	0.99	0.99
	0.25 0.99	0.84	1	0.98	0.63	0.99	0.69
	0.5	0.91	0.60	1	0.51	0.49	0.79	0.35
	0.25 0.01 0.99	0.96	1	1	0.98	0.99	0.98
	0.1	0.99	0.91	1	0.99	0.79	0.99	0.87
	0.25 0.94	0.69	1	0.78	0.54	0.95	0.47
	0.5	0.54	0.54	0.75	0.32	0.46	0.27	0.30
	0.5 0.01 0.95	0.84	0.99	0.93	0.82	0.99	0.86
	0.1	0.83	0.68	0.97	0.73	0.55	0.90	0.51
	0.25 0.52	0.54	0.73	0.34	0.46	0.32	0.30
	0.5	0.22	0.46	0.11	0.21	0.43	0.10	0.25

.1 -NMI scores obtained after clustering perturbed graphs of different sizes using the 1spectral clustering algorithm. All results are averaged over 100 replicates.

  A multicentric observational study was conducted in the Second Department of Surgery, Emergency Mures ¸County Hospital of Romania, the Department of Systems Medicine of the Tor Vergata University of Rome, the Institut d'Investigacio Biomedica de Girona IdibGi, the Endocrinology and Nutrition Department of Dr. Josep Trueta University Hospital, and the University Hospital of Innsbruck. All research procedures performed in this study were in strict accordance with a Table 4.1 -Baseline characteristics of patients with biopsy-proven fibrosis predefined protocol and adhered to the Good Clinical Practice guidelines and the Declaration of Helsinki. The study was approved by the Coordinating Ethics Committee of the Emergency Mures Çounty Hospital, Romania (registration 4065/2014), the Institutional review board & Ethics Committee and the Committee for Clinical Research (CEIC) of Dr. Josep Trueta University Hospital, Girona, Spain; the Policlinico Tor Vergata Ethics Committee, Rome, Italy as part of the FLORINASH Study the Institutional Ethics Commission at the medical University of Innsbruck (amendment to AN20170016369

	Characteristics	All patients	Stage F 0	Stage F 1	Stage F 2	p-value p-value p-value
		N = 82	N = 34	N = 37	N = 11	F 0 vs F 1 F 2 vs F 0 F 2 vs F 1
	Age (years)	41.50 ± 11.52 39.5 ± 12.77	39 ± 9.53	50 ± 9.15	0.99	10.16	0.03 *
	Female (n)	47(57%)	15(18%)	26(32%)	6(7.3%)	0.65	0.99	0.97
	Height (m)	1.67 ± 0.08	1.67 ± 0.08	1.7 ± 0.08	1.62 ± 0.07	0.99	0.61	0.14
	Smoker (n)	22(27%)	10(12%)	10(12%)	2(3%)	0.99	0.99	0.97
	Weight (kg)	118.5 ± 23.99 120 ± 22.55	118 ± 21.59 115.8 ± 35.77	0.99	0.99	0.97
	BMI (kg/m 2 )	42.65 ± 7.73	43.25 ± 6.9	41.6 ± 7.2	41.52 ± 11.41	0.99	0.99	0.97
	Waist (cm)	121 ± 18.37	124.5 ± 19.4	120 ± 15.61	120 ± 24.12	0.99	0.99	0.99
	Blood Glucose (mg/dl)	95.7 ± 25.76	95 ± 27.46	99 ± 21.22	95 ± 34.63	0.99	0.99	0.97
	Treated Diabetes (n)	7(8.5%)	1(1.2%)	2(2%)	4(4.7%)	0.99	0.02 *	0.027 *
	Systolic (mm Hg)	130 ± 19.47 130.5 ± 20.76 124 ± 17.43	134 ± 18.45	0.65	0.99	0.73
	Diastolic (mm Hg)	80.0 ± 11.59	80.5 ± 11.4	75 ± 10.31	90 ± 15.3	0.88	0.99	0.97
	Treated hypertension (n)	20(24%)	8(9.7%)	5(6%)	7(8.2%)	0.99	0.15	0.027 *
	Treated dyslipidemia (n)	6(7.3%)	2(2%)	3(3.6%)	1(1.2%)	0.99	0.99	0.99
	Total Cholesterol (mg/dL) 189.1 ± 39.78 190.0 ± 36.93 200.0 ± 43.11 167.0 ± 38.71	0.99	0.99	0.97
	HDL Cholesterol (mg/dL) 43.91 ± 13.38	47 ± 11.73	43 ± 13.48	42 ± 16.62	0.99	0.61	0.97
	GOT (U/l)	20.85 ± 17.56 18.50 ± 18.14	22 ± 18.97	22 ± 7.54	0.99	0.99	0.97
	GPT (U/l)	27.50 ± 25.23 23.50 ± 17.50	29 ± 31.91	30 ± 14.16	0.65	0.99	0.97
	GGT (U/l)	29 ± 23.04	27.50 ± 18.04	30 ± 25.84	32 ± 23.3	0.65	0.61	0.99
	HCT (%)	41 ± 4.03	40 ± 4.09	41.1 ± 3.05	40.5 ± 6.13	0.99	0.99	0.94
	Leukocytes (G/L)	7.84 ± 2.63	7.48 ± 2.4	8.1 ± 2.39	7.8 ± 3.7	0.99	0.61	0.97
	Neutrophils (G/L)	5 ± 2.44	4.8 ± 2.36	5.15 ± 2.28	5.3 ± 3.2	0.99	0.99	0.99

Table 4

 4 ).

	OTU name	Phylum	Family	Genus
	Cluster 30 Proteobacteria	Moraxellaceae	Acinetobacter
	Cluster 74	Firmicutes	Ruminococcaceae Faecalibacterium
	Cluster 35 Actinobacteria Microbacteriaceae	Rhodoluna
	Cluster 28 Actinobacteria Micrococcaceae	Kocuria
	Cluster 21 Proteobacteria Caulobacteraceae	Caulobacter
	Table 4.2 -Identification of specific bacterial signatures (unfair analyses) in F 0 patients
	OTU name	Phylum	Family	Genus
	Cluster 43 Proteobacteria Pseudomonadaceae Pseudomonas
	Cluster 31 Proteobacteria Pseudomonadaceae Pseudomonas
	Cluster 25 Proteobacteria Enterobacteriaceae Multi-affiliation
	Cluster 37 Proteobacteria Rhodobacteraceae	Paracoccus
	Cluster 40	Bacteroidetes	Chitinophagaceae Ferruginibacter

.3 -Identification of specific bacterial signatures (unfair analyses) in F 1 patients

Table 4

 4 ).

	OTU name	Phylum	Family	Genus
	Cluster 122	Firmicutes	Lachnospiraceae	Multi-affiliation
	Cluster 59 Actinobacteria	Corynebacteriaceae	Corynebacterium 1
	Cluster 34	Bacteroidetes	Weeksellaceae	Cloacibacterium
	Cluster 48	Firmicutes	Peptostreptococcaceae	Romboutsia
	Cluster 42 Proteobacteria	Burkholderiaceae	Tepidimonas
	Cluster 26	Firmicutes	Streptococcaceae	Lactococcus
	Cluster 53	Bacteroidetes	Weeksellaceae	Cloacibacterium
	Cluster 18 Proteobacteria	Enterobacteriaceae	Pantoea
	Cluster 115 Proteobacteria	Burkholderiaceae	Delftia
	Cluster 50 Actinobacteria	Microbacteriaceae	Clavibacter
	Cluster 91 Actinobacteria	Corynebacteriaceae	Corynebacterium
	Cluster 54 Proteobacteria	Burkholderiaceae	Ralstonia

.4 -Identification of specific bacterial signatures (unfair analyses) in F 2 patients

  ). VS F 1 and F 0 VS F 2 Cluster 11 Proteobacteria Pseudomonadaceae PseudomonasF 0 VS F 1 Cluster 89 Actinobacteria Corynebacteriaceae Corynebacterium 1 F 0 VS F 1 and F 0 VS F 2 VS F 1 and F 0 VS F 2 Table 4.5 -Identification of clusters of cohort-independent 16S rDNA associated with the different low scores of fibrosis from sPL-SDA clustering)(Figure 4.19). Interestingly, from all selected OTUs eight common OTUs were from the same phylum i.e. Proteobacteria (Table

	OTU name	Phylum	Family	Genus	Significance
	Cluster 20 Proteobacteria Burkholderiaceae	Ralstonia	F 0 VS F 1
	Cluster 15 Proteobacteria Xanthobacteraceae	Bradyrhizobium	F 0 VS F 1
	Cluster 16 Proteobacteria Enterobacteriaceae	Multi-affiliation	F 0 VS F 1
	Cluster 31 Proteobacteria Pseudomonadaceae	Pseudomonas	F 0 VS F 1
	Cluster 24 Proteobacteria Enterobacteriaceae	Kluyvera	F 0 VS F 1
	Cluster 14 Proteobacteria Xanthomonadaceae Stenotrophomonas	F 0 VS F 1
	Cluster 4	Proteobacteria Enterobacteriaceae	Multi-affiliation	F 0 VS F 1
	Cluster 25 Proteobacteria Enterobacteriaceae F 0 Cluster 5 Multi-affiliation Bacteroidetes Flavobacteriaceae Flavobacterium F 0

  ).

	OTU name		Phylum	Family		Genus	Significance
	Cluster 14 Proteobacteria Xanthomonadaceae Stenotrophomonas	F 0 VS F 1
	Cluster 77 Proteobacteria Enterobacteriaceae	Multi-affiliation	F 1 VS F 2 and F 0 VS F 2
	Cluster 35 Actinobacteria Microbacteriaceae		Rhodoluna	F 0 VS F 1
	Cluster 36	Bacteroidetes	Flavobacteriaceae	Flavobacterium	F 0 VS F 1
	Cluster 16 Proteobacteria Enterobacteriaceae	Multi-affiliation	F 0 VS F 1 and F 0 VS F 2
	Cluster 312 Proteobacteria Enterobacteriaceae		Kosakonia	F 0 VS F 1
	Cluster 24 Proteobacteria Enterobacteriaceae		Kluyvera	F 0 VS F 1
	Table 4.7 -Identification of clusters of cohort-independent 16S rDNA associated with the different
	low scores of fibrosis from Random Forest
	OTU name	Phylum	Family		Genus	Significance
	Cluster 31 Proteobacteria Pseudomonadaceae	Pseudomonas	F 0 VS F 1
	Cluster 15 Proteobacteria Xanthobacteraceae	Bradyrhizobium	F 0 VS F 1
	Cluster 341 Proteobacteria Enterobacteriaceae	Enterobacter	F 0 VS F 1
	Cluster 20 Proteobacteria Burkholderiaceae		Ralstonia	F 0 VS F 1
	Cluster 4	Proteobacteria Enterobacteriaceae	Multi-affiliation	F 0 VS F 1
	Cluster 2	Proteobacteria Enterobacteriaceae Escherichia-Shigella	F 0 VS F 1
	Cluster 27 Proteobacteria Enterobacteriaceae	Serratia	F 0 VS F 1
	Cluster 3	Proteobacteria Pseudomonadaceae	Pseudomonas	F 0 VS F 1
	OTU name Cluster 24 Proteobacteria Enterobacteriaceae Phylum Family Cluster 31 Proteobacteria Pseudomonadaceae Cluster 22 Proteobacteria Burkholderiaceae Cluster 66 Proteobacteria Pasteurellaceae Cluster 16 Proteobacteria Enterobacteriaceae	Genus Kluyvera Pseudomonas Multi-affiliation Haemophilus F 1 VS F 2 and F 0 VS F 2 F 0 VS F 1 Significance F 0 VS F 1 F 0 VS F 1 Multi-affiliation F 0 VS F 1
	Cluster 335 Proteobacteria Burkholderiaceae Table 4.8 -Identification of clusters of cohort-independent 16S rDNA associated with the different Janthinobacterium F 0 VS F 1
	Cluster 42 Proteobacteria Burkholderiaceae low scores of fibrosis from Fair l 1 -spectral clustering		Tepidimonas	F 0 VS F 1
	Cluster 341 Proteobacteria Enterobacteriaceae		Enterobacter	F 0 VS F 1
	Cluster 248 Proteobacteria	Reyranellaceae		Reyranella	F 0 VS F 1 and F 0 VS F 2
	Cluster 231 Proteobacteria Enterobacteriaceae	Multi-affiliation	F 0 VS F 1
	Cluster 36 OTU name Cluster 64 Actinobacteria Corynebacteriaceae Bacteroidetes Flavobacteriaceae Phylum Family Cluster 31 Proteobacteria Pseudomonadaceae Pseudomonas Flavobacterium Genus Lawsonella Pseudomonas putida Species F 1 VS F 2 and F 0 VS F 2 Significance F 0 VS F F 1 VS F 2 and F 0 VS F 2 Cluster 339 Proteobacteria Burkholderiaceae Limnohabitans Cluster 15 Proteobacteria Xanthobacteraceae Bradyrhizobium Bradyrhizobium sp. F 0 VS F F 1 VS F 2 Cluster 15 Proteobacteria Xanthobacteraceae Bradyrhizobium Cluster 25 Proteobacteria Enterobacteriaceae Multi-affiliation Multi-affiliation F 0 VS F F 0 VS F 1 Cluster 44 Actinobacteria Intrasporangiaceae Multi-affiliation Cluster 14 Proteobacteria Xanthomonadaceae Stenotrophomonas Multi-affiliation F 0 VS F F 0 VS F 1 and F 0 VS F 2 Cluster 25 Proteobacteria Enterobacteriaceae Multi-affiliation Cluster 16 Proteobacteria Enterobacteriaceae Multi-affiliation Multi-affiliation F 0 VS F F 0 VS F 1 Cluster 45 Proteobacteria Burkholderiaceae Comamonas Cluster 24 Proteobacteria Enterobacteriaceae Kluyvera Multi-affiliation F 0 VS F F 0 VS F 1 Cluster 20 Proteobacteria Burkholderiaceae Ralstonia Multi-affiliation F 0 VS F
	Cluster 14 Proteobacteria Xanthomonadaceae Stenotrophomonas Cluster 4 Proteobacteria Enterobacteriaceae Multi-affiliation Multi-affiliation F 0 VS F 1 F 0 VS F
	Table 4.6 -Identification of clusters of cohort-independent 16S rDNA associated with the different Table 4.9 -Microbial signatures common to all strategies
	low scores of fibrosis from Fair-tree

-pénalisée du spectral clustering afin d'améliorer ses performances sur des graphes bruités. Considérons dans un premier temps un graphe idéal G(V, E) non dirigé, non pondéré et sans boucle modélisant un système non bruité. Le graphe ainsi défini est composé d'un ensemble de p noeuds V = {1, . . . , p} et d'un ensemble d'arêtes E ⊆ V × V connectant chaque noeud. Celui-ci peut être modélisé par une matrice dite d'adjacence symétrique notée A = (A ij ) (i,j)∈E de taille p × p et vérifiant :

https://es.sourceforge.net/projects/core-clustering/

https://archive.ics.uci.edu/ml/datasets/Yeast

http://users.diag.uniroma1.it/challenge9/

https://es.sourceforge.net/projects/core-clustering/

https://es.sourceforge.net/projects/core-clustering/ 2. https://github.com/championcamille/

Consider the matrix X of p predictors previously introduced and a random variable Y = (Y 1 , . . . , Y n ) T representing the target variable we want to predict from the observed X. PLS iteratively creates decorrelated components T 1 , . . . , T p , called latent variables, defined as linear combinations of the p original centred variables (X 1 , . . . , X p ):

where the vectors (w 1 , . . . , w p ), called principal factors, are solutions of the following optimization problems :

(5.2)

Contrary to PCA, the optimal number of PLS latent variables is automatically determined by cross-validation [START_REF] Lachenbruch | Estimation of error rates in discriminant analysis[END_REF]. Note that PLS is commonly combined with a regression step. The resulting method, called Partial Least Square Regression (PLSR) yields two variables B and B 0 such that Ŷ = B 0 + T B is the best approximation of Y , where T is the matrix obtained by concatenation of the centred latent variables.

In what follows, we will note f PLS (w) = Cov(Xw, Y ) the objective function to optimize in PLS (equation 5.2).

The biomedical studies conducted so far on biased dataset led us to develop fair learning approaches which rely on such standard dimensionality reduction techniques. Based on specific independence criteria, these strategies aim to explain the total variablities in the dataset while limiting the bias effect.

satisfying for all X i , X j ∈ X ,

This independence criterion, proposed by [START_REF] Gretton | Measuring statistical dependence with hilbert-schmidt norms[END_REF], is measured by computing the Hilbert-Schmidt norm of the cross covariance operator associated with their RKHSs. In order to make HSIC a practical criterion for testing independence, it can be approximited given finite number of observations as:

where K and L are the kernels of X and Y and H a centering matrix verifying (I n -1 n 1 n,n ) where 1 n,n is a matrix which all coefficients are equal to 1 and I n the identity matrix of size n. Based on this result, in order to maximize the dependence between two kernels we need to increase the value of this empirical estimate.

Kernel based PLS algorithm

In many situations, non linear transformations are required to successfully apply learning algorithms. Consider a non linear transformation of X i ∈ X ⊂ R p into a feature space H. Using the straightfoward link between a RKHS and H, PLS has been extended into a non linear kernel form [START_REF] Rosipal | Kernel partial least squares regression in reproducing kernel hilbert space[END_REF]. This extension implies the construction of a PLS model in H.

Denote Φ as the (n × M ) matrix of mapped X -space data φ(X i ) into an M -dimensional feature space H. Instead of an explicit mapping of the data, we can directly use the following matrix using the kernel trick [START_REF] Schölkopf | Nonlinear component analysis as a kernel eigenvalue problem[END_REF]:

where K represents the (n × n) kernel Gram matrix of the cross dot products between all input data points {φ(X i )} 1≤i≤n : K ij = K(X i , X j ), where K(., .) is a selected kernel function.

Motivated by the theory of RKHS, [START_REF] Rosipal | Kernel partial least squares regression in reproducing kernel hilbert space[END_REF] derived the PLS algorithm for the (nonlinear) kernel PLS model whose purpose is to find principal factors (w 1 , . . . , w p ), solutions of the following optimization problems:

(5.4)

PLS under fair constraint based on Hilbert-Schmidt Independance Criterion

Several authors have considered the framework of penalized regression either to achieve fairness or to provide guarantees for independence [START_REF] Pérez-Suay | Fair kernel learning[END_REF]. In this subsection, we describe how HSIC can be used in Partial Least Square technique as an independence criterion between the latent variables and the sensitive variable S.

-M. Minty, T. Canceill, S. Lê, P. Dubois, O. Amestoy, P. Loubieres, J. Christensen, C.

R packages

-R package CoreClustering, Available on sourceforge 1 , -R package l1spectral, Available on CRAN and GitHub 2 .