Je Remercie

Jean-Louis Foulley

Julien Chiquet

Laurent Schibler

Keywords: modèle de mélange, modèle graphique, données RNA-seq, classification, inférence de réseaux, sélection de modèle. 4 Clustering and network inference mixture model, graphical model, RNA-seq data, clustering, network inference, model selection . . . . . . . . . . . . . . . . . . . . . . . . . .

Cette thèse regroupe des contributions méthodologiques à l'analyse statistique des données issues des technologies de séquençage du transcriptome (RNA-seq). Les difficultés de modélisation des données de comptage RNA-seq sont liées à leur caractère discret et au faible nombre d'échantillons disponibles, limité par le coût financier du séquençage.

Une première partie de travaux de cette thèse porte sur la classification à l'aide de modèles de mélange. L'objectif de la classification est la détection de modules de gènes co-exprimés. Un choix naturel de modélisation des données RNA-seq est un modèle de mélange de lois de Poisson. Mais des transformations simples des données permettent de se ramener à un modèle de mélange de lois gaussiennes. Nous proposons de comparer, pour chaque jeu de données RNA-seq, les différentes modélisations à l'aide d'un critère objectif permettant de sélectionner la modélisation la plus adaptée aux données. Par ailleurs, nous présentons un critère de sélection de modèle prenant en compte des informations biologiques externes sur les gènes. Ce critère facilite l'obtention de classes biologiquement interprétables. Il n'est pas spécifique aux données RNA-seq. Il est utile à toute analyse de co-expression à l'aide de modèles de mélange visant à enrichir les bases de données d'annotations fonctionnelles des gènes.

Une seconde partie de travaux de cette thèse porte sur l'inférence de réseau à l'aide d'un modèle graphique. L'objectif de l'inférence de réseau est la détection des relations de dépendance entre les niveaux d'expression des gènes. Nous proposons un modèle d'inférence de réseau basé sur des lois de Poisson, prenant en compte le caractère discret et la grande variabilité inter-échantillons des données RNA-seq. Cependant, les méthodes d'inférence de réseau nécessitent un nombre d'échantillons élevé. Dans le cadre du modèle graphique gaussien, modèle concurrent au précédent, nous présentons une approche non-asymptotique pour sélectionner des sous-ensembles de gènes pertinents, en décomposant la matrice variance en blocs diagonaux. Cette méthode n'est pas spécifique aux données RNA-seq et permet de réduire la dimension de tout problème d'inférence de réseau basé sur le modèle graphique gaussien.

contribué à la bonne ambiance de ce lieu : Denis Laloë, Tatiana Zerjal, Dominique Montagu, Xavier Rognon, Agathe Vieaud, Michèle Boichard, Francis Minvielle, Patricia Huan, Wendy Brand-Williams, André Neau, Gwendal Restoux, Grégoire Leroy, Thomas Heams, Simon Boitard, Eléonore Charvolin, Yvelise Fricot, Marie-Hélène Pinard Vander-Laan, Valérie Bonjean, Nadjet Boushaba, Sonia Eynard, Fréderic Hospital et en particulier Gilles Monneret, Mathieu Tiret, Jean-Noël Hubert et Belén Jimenez pour leur aide et leur soutien, et François Deniau pour sa gentillesse.

Je remercie également les doctorants d'Orsay pour leur présence : Thomas Morzadec, Pierre-Antoine Guihéneuf, Vincent Pecastaing, Emilien Joly, Elodie Vernet, Pierre Gaillard, YueYuan Gao et Loïc Lacouture (merci pour ton accueil au MAP5). Merci en particulier à Solenne Thivin, Céline Abraham, Emilie Devijver et Valérie Robert pour leur aide et les bons moments passés en leur compagnie. Merci à ceux qui ont animé le couloir du 440 : Christine Keribin, Kevin Bleakley, Laura Brocco, Célia Barthélémy, Elodie Maillot, Romain Bar, Patrick Bouvier, Rémi Foucherau, Chú Cng Mp, Mohmoh A. Sedki, Sébastien Miquel, Florence Duc, Magda et Jana Khayal. Merci en particulier à Valérie Robert, Yann Vasseur, Vincent Brault et Jana Kalawoun pour l'entraide "fraternelle".

Enfin, je remercie ma famille pour son soutien et son amour inconditionnel ainsi que toutes les personnes qui m'ont aidée et inspirée au cours de ces années de thèse et toutes les personnes qui comptent beaucoup pour moi : con mucho cariño, je pense à vous tendrement. Le premier génome séquencé est celui d'une bactérie durant les années 1970 [START_REF] Sanger | DNA sequencing with chainterminating inhibitors[END_REF]. Ce génome ne comportait que 5 000 nucléotides et le coût de séquençage était très élevé. Les méthodes ont depuis beaucoup évolué, notamment grâce au développement de la technique de réaction en chaîne par polymérase (PCR pour polymerase chain reaction), qui permet de répliquer artificiellement un brin d'ADN. Cette technique de PCR est toujours très utilisée. On séquence à l'heure actuelle le génome de nombreux individus à un prix toujours décroissant. En septembre 2001, le coût de séquençage d'un génome complet humain s'élevait à presque 100 millions de dollars. En avril 2015, ce coût est d'environ 4000 dollars 1 . Cette baisse du coût de séquençage de l'ADN entraine 1. http://www.genome.gov/sequencingcosts/ 1.1.1 -Mesurer l'expression des gènes la production de très grandes quantités de données : [START_REF] Stephens | Big Data : Astronomical or Genomical ?[END_REF] prédisent que entre 100 millions et 2 milliards de génomes humains seront séquencés d'ici 2025. Cependant, connaître la séquence d'ADN d'un individu ne suffit pas pour comprendre les mécanismes d'expression de ses gènes. Pour les comprendre, le contenu en ARN et protéines des cellules doit être analysé.

Les protéines sont les acteurs principaux des phénomènes biologiques observés et leur présence ou absence peut expliquer une caractéristique précise comme un phénotype ou une pathologie. Cependant, mesurer la quantité des protéines est très difficile du fait de leur instabilité et de leur dégradation rapide dans la cellule. Il est de plus difficile d'associer une protéine donnée à une position du génome. La mesure de la quantité d'ARN matures, stables, est plus facile. Cette quantité de transcrits est donc considérée comme représentative du niveau d'expression du gène (i.e. de l'abondance des protéines produites). Dans cette thèse, l'expression des gènes désigne la quantité de transcrits présents dans la cellule associée à l'ensemble des gènes considérés.

Pour mesurer cette quantité de transcrits, plusieurs méthodes existent et suivent le protocole général suivant. En amont, les conditions expérimentales et le nombre d'échantillons par condition expérimentale doivent être déterminés en fonction du type d'analyses souhaitées. L'ARN est ensuite extrait de la cellule et retranscrit en un brin d'ADN complémentaire (ADNc) à l'aide d'enzymes spécifiques. Le brin complémentaire de l'ADNc est également produit, puis l'ARN présent est retiré du mélange. Seuls les ARN matures sont collectés. Par la suite, les brins d'ADNc sont amplifiés par PCR afin de disposer d'une quantité assez importante pour la suite du protocole. Il convient alors de mesurer la quantité d'ADNc. Pour effectuer cette mesure, plusieurs méthodes existent : L'alignement est rarement parfait, un nombre de bases peut différer entre le read et le génome de référence comme l'illustrent les bases en gras.

Les puces à ADN

Le nombre moyen de reads alignés par position du génome est appelé la couverture de séquençage. Plus elle est grande, plus le séquençage est complet. Le nombre de reads alignés pour un échantillon expérimental est appelé la profondeur de séquençage. L'algorithme d'alignement peut également s'appuyer sur une annotation du génome déjà connue (i.e. la liste des positions de chaque gène sur le génome). L'annotation déjà connue d'un génome donné peut être téléchargée sous la forme de fichiers .GFF à partir des bases de données publiques d'annotation, sur le site NCBI ou Ensembl (Biomart) 4 . L'étape d'alignement ne requiert pas cette annotation du génome et peut être effectuée sans cette information.

De nombreuses méthodes existent pour effectuer cet alignement, par exemple TopHat [START_REF] Trapnell | TopHat : Discovering splice junctions with RNA-Seq[END_REF]. Les programmes d'alignement transforment les fichiers de reads bruts .FASTQ en fichier binaire de format BAM ou texte de format SAM signifiant Sequence Alignment/Map. Les fichiers .SAM et .BAM contiennent l'information de position de chaque read sur le génome. Ils sont volumineux et doivent être manipulés à l'aide d'outils spécifiques, par exemple samtools [START_REF] Li | The Sequence Alignment/Map format and SAMtools[END_REF].

4. Comptage des reads : Une fois la position de chaque read connue, les reads alignés sur chaque région génomique d'intérêt sont dénombrés. L'annotation de référence téléchargée à partir des bases de données publiques permet de connaître la délimitation des régions génomiques d'intérêt (en particulier, la position des gènes). Cependant, les programmes permettent aussi de détecter de nouvelles régions génomiques encore non identifiées. Les régions génomiques nouvellement détectées viennent ainsi enrichir les bases de données existantes. De nombreux programmes existent également pour compter l'abondance des transcrits par région génomique, en particulier Cufflink [START_REF] Trapnell | Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation[END_REF] ou HTSeq-counts [START_REF] Anders | HTSeq A Python framework to work with high-throughput sequencing data[END_REF] . Le fichier produit en sortie de cette étape d'estimation d'abondance des reads par région génomique est un tableau de comptage, comme l'illustre le Tableau 1.1.

La technologie RNA-seq vient révolutionner la mesure de l'expression des gènes. En s'appuyant sur la puissance de calcul à disposition (notamment pour l'étape d'alignement des reads), elle permet la mesure haut-débit de l'expression simultanée de plusieurs milliers de gènes à la fois, contrairement à la technologie SAGE dont le débit est beaucoup moins important. De plus, le RNA-seq mesure l'expression sans limite d'amplitude : un gène peut avoir un alignement de 0 read ou plus de 100 000 reads, contrairement à la technologie des puces à ADN. L'alignement des séquences peut être effectué sur un génome de référence potentiellement incomplet (alignement de novo) ce qui permet la découverte de nouveaux gènes, ce qui n'est pas possible avec les données de puces à ADN. En particulier, cette technologie permet d'examiner l'expression spécifique aux allèles, qui sont les variantes d'un même gène résultant d'une mutation et ayant la même fonction que le gène initial selon ses modalités propres. La technologie RNA-seq offre ainsi un niveau de précision supérieur à celui des données de puces à ADN.

De nombreuses méthodes statistiques sont utiles à différentes étapes de mesure de l'expression des gènes par la technologie RNA-seq : contrôle qualité des reads, alignement des reads et estimation de l'abondance des transcrits. Dans cette thèse, on s'intéressera uniquement aux méthodes statistiques destinées à analyser le tableau de comptages (Table 1.1) après estimation de l'abondance des reads.

Modélisation statistique des données RNA-seq

Modélisation et normalisation

Les notations utilisées tout au long de la thèse sont fixées ici. Par abus de langage, toute région génomique d'intérêt est assimilée à un gène. On utilise également le terme d'échantillon pour désigner les réplicats biologiques ou techniques. Les réplicats techniques désignent les échantillons collectés sur un même individu pour une même condition expérimentale. Les réplicats biologiques désignent les échantillons collectés sur des individus différents. Les jeux de données analysés dans cette thèse ne contiennent pas de réplicat technique : il s'agit toujours de réplicats biologiques d'une même condition ou de conditions expérimentales différentes. On utilise le terme d'échantillon pour désigner les réplicats biologiques et l'on précise les différentes conditions expérimentales si besoin. Le tableau de comptage d'expression de gènes (exemple : tableau 1.1) est représenté par une matrice, notée y. Cette matrice d'expression est de taille (n × p), n est le nombre d'échantillons séquencés et p le nombre de gènes. L'entrée y ij désigne l'expression du gène j pour l'échantillon i. Pour tout échantillon i ∈ 1, . . . , n, le vecteur y i est le vecteur d'expression de l'échantillon i pour les p gènes. Pour tout gène j ∈ 1, . . . , p, le vecteur y j est le vecteur d'expression du gène j pour les n échantillons. Si les échantillons correspondent à différentes conditions expérimentales, cette condition expérimentale est notée c(i) pour chaque échantillon i.

Modélisation des données RNA-seq

Le nombre de reads alignés pour un échantillon donné est fortement dépendant de la profondeur de séquençage qui est une contrainte spécifique de la technologie RNA-seq. Pour un échantillon donné, ce nombre peut être considéré fixe. Sachant ce nombre N i = p j=1 y ij , le vecteur aléatoire Y i modélisant la distribution des reads sur les p régions génomiques du génome pour l'échantillon i peut être modélisé par une loi multinomiale de paramètres N i et (π i1 , . . . , π ip ) :

p(Y i1 = y i1 , . . . , Y ip = y ip ) = N i y i1 ! . . . y ip ! π y i1 i1 . . . π y ip ip .
Chaque composante Y ij de ce vecteur multinomial Y i suit une loi binomiale de paramètres N i et π j . Le nombre de reads alignés par échantillon N i est élevé et π j , la proportion de reads alignés sur le gène j est petite. L'approximation d'une loi binomiale par une loi de Poisson est applicable :

Y ij ∼ P(λ ij ), p(Y ij = y ij ) = λ y ij ij y ij ! e -λ ij , E(Y ij ) = λ ij , V (Y ij ) = λ ij .
Le modèle de Poisson est utile en particulier pour l'analyse de réplicats techniques d'une même condition [START_REF] Marioni | comparison with gene expression arrays RNA-seq : An assessment of technical reproducibility and comparison with gene expression arrays[END_REF]. Cependant, il modélise mal la grande variabilité inter-échantillons observée dans les données RNA-seq lorsque les différents réplicats biologiques sont des réplicats biologiques et non techniques d'une même condition expérimentale. Afin de prendre en compte cette grande variabilité, des modèles de loi de Poisson surdispersée ou de loi binomiale négative ont été proposés comme solutions alternatives au modèle de Poisson : Y ij ∼ NB(λ ij , φ j ). Le paramètre λ ij est modélisé de manière appropriée et φ j ≥ 0 est le paramètre de dispersion du gène j. Plusieurs paramétrisations de la loi binomiale négative sont possibles, mais il est courant d'adopter la convention suivante :

Y ij ∼ NB(λ ij , φ j ), p(Y ij = y ij ) = φ j φ j + λ ij φ j Γ(φ j + y ij ) y ij ! Γ(φ j ) λ ij φ j + λ ij y ij , E(Y ij ) = λ ij , V (Y ij ) = λ ij + λ 2 ij φ j .
Que ce soit pour la loi de Poisson ou la loi binomiale négative, il convient de modéliser le paramètre λ ij de manière appropriée, par exemple en supposant que λ ij = e i g j où g j est un paramètre spécifique au gène j et e i est un paramètre spécifique à l'échantillon i. Ces modélisations sont adaptées en fonction du type d'analyse effectuée : analyse différentielle, classification ou inférence de réseaux. Les paramètres doivent prendre en compte les biais techniques de la technologie RNA-seq que nous détaillons dans la section suivante.

Prise en compte de la profondeur de séquençage par échantillon Le nombre de reads alignés pour un gène et un échantillon donné est une mesure relative et non absolue de l'expression du gène. Ce nombre dépend du nombre total de fragments d'ADN alignés pour l'échantillon donné comme l'illustre l'exemple du tableau 1.2. Pour que l'expression d'un gène soit comparable entre plusieurs réplicats biologiques, il convient de prendre en compte le nombre total de reads alignés pour chaque échantillon, aussi appelé taille de librairie ou profondeur de séquençage par échantillon. Le comptage par million (count per million CPM) du gène j pour l'échantillon i correspond au nombre de reads y ij alignés sur cette région, normalisé par le nombre total de reads alignés pour l'échantillon i (noté

N i = p j=1 y ij ) divisé par un million : cpm(y ij ) = y ij N i /10 6 .
Le log-comptage par million (log-count per million log-cpm) du gène j pour l'échantillon i a été utilisé par [START_REF] Law | Voom : precision weights unlock linear model analysis tools for RNA-seq read counts[END_REF]. Le log-cpm correspond au logarithme du nombre de reads y ij alignés sur cette région, normalisé par le nombre total de lectures alignées pour l'échantillon i exprimé en un million : log-cpm(y ij ) = log 2 y ij + 0.5 (N i + 1)/10 6 .

Prise en compte de la longueur du gène

Pour comparer l'expression de différents gènes au sein d'un échantillon donné, il convient de prendre en compte la longueur de la région génomique sur laquelle les 1.2.2 -Analyse différentielle fragments sont alignés. En effet, la longueur d'un gène peut varier entre 400 et 2 millions de paires de bases. La longueur du gène peut donc expliquer à elle seule les différences d'expression observées. Pour un gène j, on suppose cette longueur connue. La longueur du gène, notée L j et exprimée en nombre de paires de bases, ne peut être estimée à partir du simple tableau de comptages RNA-seq. Elle doit être extraite des bases de données publiques d'annotation du génome. La méthode de normalisation la plus simple pour la comparaison de mesure d'expression entre deux gènes consiste à diviser chaque comptage par la longueur du gène L j correspondante. On obtient ainsi une mesure d'expression du gène par paire de bases. Cependant, cette mesure d'expression ne prend pas en compte les différences de taille de librairie. La prise en compte des différences de taille de librairie est réalisée dans le calcul du nombre de transcripts per million (TPM) [START_REF] Li | The Sequence Alignment/Map format and SAMtools[END_REF] ou dans le calcul du nombre de lectures par kilobase d'exon par million de lectures alignées (reads per kilobase of exon per million reads mapped, RPKM) [START_REF] Mortazavi | Mapping and quantifying mammalian transcriptomes by RNA-Seq[END_REF] :

rpkm(y ij ) = y ij L j 10 3 N i 10 6 .
Autres biais des données RNA-seq et modélisation D'autres biais des données RNA-seq existent. Un exemple est le biais GC. Le contenu en GC d'un gène désigne la proportion de nucléotides G et C dans la séquence du gène. Si le contenu en GC est sous-représenté ou sur-représenté, le séquençage du read est plus difficile, ce qui peut être pris en compte dans le calcul d'unité d'expression des données RNA-seq. [START_REF] Risso | GC-content normalization for RNA-seq data[END_REF]; [START_REF] Hansen | Removing technical variability in RNAseq data using conditional quantile normalization[END_REF] ajustent à l'aide d'une régression loess les log-comptages sur le contenu en GC de chaque gène.

Dans cette thèse, nous ne présentons pas toutes les méthodes possibles de correction des biais des données RNA-seq. De plus, ces méthodes dépendent fortement du type d'analyse statistique effectuée en aval. On peut cependant distinguer deux grandes approches de modélisation des données RNA-seq :

1. Les modèles discrets (loi de Poisson ou loi binomiale négative) qui modélisent directement les comptages des reads et prennent en compte les biais des données RNA-seq directement dans l'écriture des paramètres du modèle, 2. Les modèles gaussiens qui modélisent les comptages normalisés et prennent en compte les biais des données RNA-seq préalablement à l'utilisation du modèle.

Analyse différentielle

L'objectif de l'analyse différentielle est de détecter les gènes différentiellement exprimés entre différentes conditions expérimentales testées. Pour que l'analyse différentielle soit pertinente, il convient de planifier l'expérience avec soin, et de choisir, en fonction du budget disponible, la profondeur de séquençage (le nombre de reads séquencés par échantillon) et le nombre d'échantillons à séquencer. [START_REF] Busby | Scotty : a web tool for designing RNA-Seq experiments to measure differential gene expression[END_REF]

H 0j = {λ 1j = λ 2j },
vs. Puisque le nombre de gènes est élevé, il convient de corriger le seuil de rejet de l'hypothèse nulle afin de contrôler correctement le nombre de gènes détectés différentiellement exprimés à tort, aussi appelé le nombre de faux positifs ou False Discovery Rate (FDR). Différentes méthodes existent pour contrôler le nombre de faux positifs, correspondant à différentes manières de mesurer la quantité de faux positifs. Si l'on souhaite contrôler le nombre de faux positifs en probabilité, on adoptera la procédure de Bonferroni améliorée par [START_REF] Hochberg | A sharper Bonferroni procedure for multiple tests of significance[END_REF]. Si l'on souhaite contrôler le taux de faux positifs en espérance, on adopte la procédure [START_REF] Benjamini | Controlling the false discovery rate : a practical and powerful approach to multiple testing[END_REF]. La procédure de Benjamini-Yekutieli est une adaptation de la procédure de Benjamini et Hochberg dans le cadre du contrôle du taux de faux positifs avec hypothèse de dépendance entre les différents tests.

H 1j = {λ 1j = λ 2j }.
L'étape 2 doit être adaptée aux données RNA-seq. Pour les données microarray qui sont continues, il est possible de construire des modèles d'analyse statistique en se basant sur des hypothèses de normalité des données [START_REF] Smyth | Linear models and empirical Bayes methods for assessing differential expression in microarray experiments[END_REF]. Ces techniques d'analyse, adaptées aux données gaussiennes ne peuvent pas être appliquées directement aux données RNA-seq qui sont des données de comptage, discrètes et positives. De nombreux travaux de recherche ont été effectués sur ce sujet.

Plusieurs modèles sont possibles, comme des lois de Poisson [START_REF] Wang | DEGseq : an R package for idenetifying differentially expressed genes from RNA-seq data[END_REF][START_REF] Marioni | comparison with gene expression arrays RNA-seq : An assessment of technical reproducibility and comparison with gene expression arrays[END_REF], des lois de Poisson sur-dispersées [START_REF] Auer | A two-stage Poisson model for testing RNA-seq data[END_REF], des lois binomiales négatives [START_REF] Anders | Differential expression analysis for sequence count data[END_REF][START_REF] Robinson | A scaling normalization method for differential expression analysis of RNA-seq data[END_REF] ou des lois normales [START_REF] Law | Voom : precision weights unlock linear model analysis tools for RNA-seq read counts[END_REF]. [START_REF] Love | Moderated estimation of fold change and dispersion for RNA-Seq data with DESeq2[END_REF] et utilisé dans le package DESeq2 est un modèle linéaire généralisé :

Le modèle proposé par

Y ij ∼ NB(µ ij , φ j ), µ ij = s i λ ij , 1.2.3 -Détection de gènes co-exprimés log 2 (λ ij ) = D i. β j .
Le comptage Y ij du gène j pour l'échantillon i est modélisé à l'aide d'une loi binomiale négative de moyenne µ ij et de paramètre de dispersion φ j spécifique au gène j. La moyenne ajustée µ ij est composée d'un paramètre spécifique à la taille de la librairie s i et d'un paramètre λ ij proportionnel à l'expression du gène j dans l'échantillon i. Le vecteur coefficient β j modélise les variations de l'expression du gène j en fonction des conditions expérimentales de chaque échantillon résumées dans la matrice de design D.

Cette modélisation prend en compte la variabilité de la profondeur de séquençage pour les différents réplicats biologiques. Il existe différentes façons de estimer les facteurs de normalisation pour les différences de profondeur de séquençage entre échantillons, par exemple en ne considérant qu'un sous-ensemble de gènes non-différentiellement exprimés. En effet, les gènes différentiellement exprimés entre différentes conditions présentent généralement des comptages très déséquilibrés entre échantillons qui viennent fausser l'estimation de la profondeur de séquençage. Les deux méthodes de normalisation correspondantes les plus utilisées sont la méthode Trimmed M-Means (TMM) [START_REF] Robinson | A scaling normalization method for differential expression analysis of RNA-seq data[END_REF] et la méthode median-of-ratios utilisée dans les packages DESeq et DESeq2 [START_REF] Anders | Differential expression analysis for sequence count data[END_REF]. Ces deux méthodes sont particulièrement appréciées pour corriger le biais induit par des tailles de librairie très différentes au sein d'un même jeu de données [START_REF] Dillies | A comprehensive evaluation of normalization methods for Illumina high-throughput RNA sequencing data analysis[END_REF]. Concernant la prise en compte de la longueur du gène, il a été montré que l'utilisation des RPKMs réduisait la puissance de détection des gènes différentiellement exprimés [START_REF] Oshlack | Transcript length bias in RNA-seq data confounds systems biology[END_REF][START_REF] Dillies | A comprehensive evaluation of normalization methods for Illumina high-throughput RNA sequencing data analysis[END_REF].

Une fois le modèle choisi, il convient d'estimer les paramètres des modèles correspondants. Pour le modèle de loi binomiale négative, les deux méthodes les plus utilisées ont été proposées par [START_REF] Robinson | A scaling normalization method for differential expression analysis of RNA-seq data[END_REF] et [START_REF] Love | Moderated estimation of fold change and dispersion for RNA-Seq data with DESeq2[END_REF] et sont disponibles respectivement dans les packages R edgeR et DESeq2. L'originalité de ces méthodes réside dans l'estimation du paramètre de dispersion de la loi binomiale négative. Les données RNA-seq ont généralement peu de réplicats biologiques et le paramètre de dispersion est de ce fait mal estimé. Les méthodes utilisent des techniques de partage d'information entre les gènes afin d'estimer le paramètre de dispersion de chaque gène. Le package DESeq2 modélise la tendance moyenne-variance afin d'estimer ce paramètre de dispersion, tandis que le package edgeR utilise un compromis entre une dispersion commune à tous les gènes et une dispersion spécifique à chaque gène. On parle alors de shrinkage du paramètre de dispersion vers une dispersion commune.

Les contributions méthodologiques à l'analyse différentielle des données RNA-seq sont nombreuses et font l'objet de recherches actives [START_REF] Oshlack | From RNA-seq reads to differential expression results[END_REF][START_REF] Dillies | A comprehensive evaluation of normalization methods for Illumina high-throughput RNA sequencing data analysis[END_REF]. Dans cette thèse, nous ne nous intéressons pas à l'analyse différentielle en tant que telle. Nous avons cependant proposé une méthode de filtrage des données RNA-seq qui est une étape préliminaire à l'analyse différentielle. Cette contribution est motivée par la nécessité de réduire la taille du jeu de données afin d'améliorer la puissance de détection des tests d'analyse différentielle décrits plus haut.

Détection de gènes co-exprimés

Au-delà de l'analyse différentielle des données d'expression, d'autres types d'analyses peuvent être effectués afin de mieux comprendre les processus sous-jacents aux données observées et de formuler des hypothèses biologiques pour des recherches futures. L'une d'entre elles est la détection de gènes co-exprimés. Nous détaillons le principe de cette analyse, les différentes méthodes utilisées en pratique et les problèmes auxquels nous nous intéressons.

L'analyse de co-expression des données d'expression de gènes consiste à grouper les gènes ayant des profils d'expression similaires à travers différentes conditions. Elle est particulièrement intéressante dans le cas où plus de deux conditions expérimentales sont étudiées et permet d'aller plus loin que l'analyse différentielle dans l'étude du transcriptome.

Plusieurs techniques de classification non supervisée permettent de classer les gènes en fonction de leur profil d'expression : la classification hiérarchique [START_REF] Ward | Hierarchical grouping to optimize an objective function[END_REF], la méthode des k-means [START_REF] Macqueen | Some methods for classification and analysis of multivariate observations[END_REF], les réseaux de neurones [START_REF] Tamayo | Interpreting patterns of gene expression with self-organizing maps : methods and application to hematopoietic differentiation[END_REF] ou les modèles de mélange [START_REF] Mclachlan | Finite Mixture Models, Willey Series in Probability and Statistics[END_REF]. Pour l'analyse de co-expression des données d'expression de puces à ADN, [START_REF] Yeung | Modelbased clustering and data transformations for gene expression data[END_REF] ont montré l'intérêt des modèles de mélange. Nous adoptons le même point de vue pour l'analyse de co-expression des données RNA-seq. L'objectif de l'analyse de co-expression est l'enrichissement des bases de données d'annotation fonctionnelle des gènes. Généralement, les gènes co-exprimés sont impliqués dans les mêmes processus biologiques [START_REF] Eisen | Cluster analysis and display of genome-wide expression patterns[END_REF][START_REF] Jiang | Cluster analysis for gene expression data : a survey[END_REF]. Ces groupes de gènes co-exprimés permettent de formuler des hypothèses sur les fonctions possibles des gènes : on peut ainsi inférer les fonctions d'un gène orphelin, dont le rôle n'a pas encore été identifié, en regardant les fonctions déjà connues des gènes classés dans le même groupe de co-expression.

Les méthodes de classification utilisées en pratique sur les données RNA-seq

En classification des données RNA-seq, certains auteurs classent les réplicats biologiques et non les gènes. Plusieurs méthodes ont été utilisées en pratique : clustering hiérarchique avec la distance euclidienne sur les données RNA-seq après transformation pour stabiliser la variance (Variance stabilizing transformation, VST) [START_REF] Anders | Differential expression analysis for sequence count data[END_REF], clustering hiérarchique avec la distance de corrélation de Pearson sur les données RNA-seq exprimées en RPKM [START_REF] Severin | RNA-Seq Atlas of Glycine max : a guide to the soybean transcriptome[END_REF] ou classification hiérarchique avec une distance entre réplicats biologiques calculée à partir d'un modèle de Poisson log linéaire [START_REF] Witten | Classification and clustering of sequencing data using a poisson model[END_REF]. On remarque que deux stratégies existent : soit les auteurs transforment les données pour utiliser les méthodes classiques de classification [START_REF] Anders | Differential expression analysis for sequence count data[END_REF][START_REF] Severin | RNA-Seq Atlas of Glycine max : a guide to the soybean transcriptome[END_REF], soit les auteurs proposent un modèle adapté aux données RNA-seq [START_REF] Witten | Classification and clustering of sequencing data using a poisson model[END_REF].

Le constat est identique pour la classification des gènes, [START_REF] Li | The Sequence Alignment/Map format and SAMtools[END_REF] utilisent un algorithme de type k-means sur les données RNA-seq log-transformées, mais des techniques adaptées au caractère discret des données RNA-seq existent. Une méthode de classification de type k-means couplée à une distance entre gènes basée sur un modèle loglinéaire de Poisson a été proposée pour la classification des tags de la technologie SAGE, ancêtre du RNA-seq [START_REF] Huang | Clustering analysis of SAGE transcription profiles using a Poisson approach[END_REF]. Plus récemment, deux modèles spécifiques aux données RNA-seq ont été proposés : un modèle de mélange de loi de Poisson [START_REF] Rau | Coexpression analysis of high-throughput transcriptome sequencing data with Poisson mixture models[END_REF] et un modèle de mélange de lois binomiales négatives [START_REF] Si | Model-based clustering for RNA-seq data[END_REF].

Problématique abordée : Le choix de modélisation est une question importante et aucune méthode de classification n'est clairement établie en analyse de co-expression des données RNA-seq. Dans le cadre de la classification par modèle de mélange, nous soulevons la question du choix de modélisation (discrète ou continue) des données RNAseq.

Les bases de données d'annotation fonctionnelle des gènes

Le terme annotation fonctionnelle désigne l'ensemble des fonctions que l'on peut associer aux produits du gène (protéines ou ARN). Plus généralement, le terme d'annotation du génome se réfère à l'ensemble de méta-données associées aux régions codantes du génome, en particulier, la position sur le génome de la région codante des gènes. Dans cette thèse, le terme annotation désignera l'annotation fonctionnelle et non l'annotation au sens général si aucune précision n'est donnée. Plusieurs bases de données ont été construites à partir des références bibliographiques collectées. La plus connue est la base de données Gene Ontology (GO) [START_REF] Ashburner | Gene ontology : tool for the unification of biology. The Gene Ontology Consortium[END_REF], construite sous la forme de graphe acyclique orienté hiérarchisant les termes d'annotations. D'autres bases de données ont été créées, comme la base de données Kyoto Encyclopedia of Genes and Genomes (KEGG) [START_REF] Kanehisa | KEGG : Kyoto encyclopedia of genes and genomes[END_REF]. Plus récemment, la base de données MSigDB (Molecular Signatures databases) [START_REF] Liberzon | Molecular signatures database (MSigDB) 3.0[END_REF] a regroupé et facilité l'utilisation de cette information en proposant notamment des listes de fonctions filtrées manuellement. Ces bases de données contiennent des informations riches sur un large spectre d'espèces (e.g., humain, souris, mouche, Arabidopsis thaliana). Cependant, notre connaissance des annotations fonctionnelles des gènes est loin d'être complète [START_REF] Tipney | An introduction to effective use of enrichment analysis software[END_REF]. Les analyses de co-expression permettent de compléter ces bases de données.

Problématique abordée : Nous verrons comment les informations déjà existantes dans les bases de données peuvent contribuer à améliorer la qualité de l'analyse de co-expression et la pertinence des groupes de gènes détectés.

Inférence de réseaux de régulation de gènes

L'objectif de l'analyse différentielle des données d'expression de gènes ou la classification de ces données par modèle de mélange n'est pas de modéliser la dépendance entre les gènes. Or, la biologie des systèmes s'intéresse aux interactions entre les différents acteurs du système biologique. Les gènes interagissent entre eux, le plus souvent par l'intermédiaire de produits facteurs de transcription. Les facteurs de transcription sont des protéines produites par des gènes qui viennent réguler ou initier la transcription d'autres gènes. Cependant, mesurer la quantité de protéines présentes dans la cellule est très difficile du fait de l'instabilité des protéines qui se dégradent très rapidement. On préfère mesurer la quantité de transcrits et considérer que cette quantité est représentative du niveau d'expression du gène (i.e. du niveau de protéines produites). Des variations simultanées de ces qualités d'expression pour deux gènes peuvent ainsi être synonymes d'interaction entre ces deux gènes. Un graphe, ensemble de noeuds et ensemble d'arêtes joignant ces noeuds, formalise ces interactions et permet de les visualiser. Dans le cas particulier des réseaux de régula-tion de gènes, les noeuds représentent les gènes et les arêtes représentent les interactions entre les gènes. Les interactions entre les gènes peuvent prendre de multiples formes et ne peuvent se résumer à un seul type d'interaction. Le formalisme de graphe n'est donc qu'une simplification de la réalité utile pour formuler des hypothèses sur le rôle des gènes. Par exemple, la détection de modules de gènes très connectés dans le réseau peut éventuellement permettre de découvrir de nouveaux gènes d'intérêt ou de nouvelles voies biologiques (biological pathways). D'un point de vue méthodologique, il existe de nombreuses manières de reconstituer un graphe. La méthode la plus simple pour inférer un réseau de régulation de gènes à partir de données d'expression consiste à calculer les corrélations simples paire à paire entre les gènes et à fixer un seuil en dessous duquel les corrélations sont considérées comme insignifiantes afin de déterminer l'absence d'arête dans le graphe. Cette méthode d'inférence de réseaux est très populaire en analyses des données d'expression [START_REF] Butte | Discovering functional relationships between RNA expression and chemotherapeutic susceptibility using relevance networks[END_REF][START_REF] De La Fuente | Discovery of meaningful associations in genomic data using partial correlation coefficients[END_REF]. La méthode Weighted Correlation Network Analysis (WGCNA), proposé par [START_REF] Langfelder | WGCNA : an R package for weighted correlation network analysis[END_REF] est également très utilisée, et implémentée dans le package R WGCNA. Cette méthode d'inférence est basée sur le feuillage de la matrice de covariance empirique. Des méthodes basées sur le calcul de l'information mutuelle ont également été proposée [START_REF] Basso | Reverse engineering of regulatory networks in human B cells[END_REF][START_REF] Margolin | ARACNE : an algorithm for the reconstruction of gene regulatory networks in a mammalian cellular context[END_REF]. Une implémentation de méthode d'inférence de réseaux basée sur l'information mutuelle est disponible dans le package R/Bioconductor minet [START_REF] Meyer | minet : A R/Bioconductor package for inferring large transcriptional networks using mutual information[END_REF]. Cependant, une corrélation très élevée entre le niveau d'expression de deux gènes n'implique pas un lien direct entre les deux gènes. [START_REF] Reverter | Combining partial correlation and an information theory approach to the reversed engineering of gene co-expression networks[END_REF] ont proposé une méthode combinant le calcul de corrélation partielle et d'information mutuelle (Partial Correlation Information Theory (PCTI) algorithm) pour supprimer les liens non directs entre les gènes du graphe. Afin de ne détecter que les liens directs entre les gènes, de nombreuses méthodes basées sur l'utilisation du modèle graphique gaussien ont été utilisées [START_REF] Schäfer | A shrinkage approach to large-scale covariance matrix estimation and implications for functional genomics[END_REF][START_REF] Peng | Partial Correlation Estimation by Joint Sparse Regression Models[END_REF]. Toutes les méthodes décrites précédemment permettent de reconstituer des réseaux de dépendance non orientés. Les réseaux bayésiens permettent d'inférer des graphes orientés [START_REF] Chen | BNArray : An R package for constructing gene regulatory networks from microarray data by using Bayesian network[END_REF] et sont particulièrement utiles sur les données d'intervention (knock-out). La modélisation des données d'expression par des équations différentielles ordinaires ou stochastiques permet également de reconstituer des réseaux de gènes. Ces méthodes sont particulièrement intéressantes en analyse des données cinétiques.

Des comparaisons entre les différentes méthodes d'inférence de réseaux sont proposées par [START_REF] Werhli | Comparative evaluation of reverse engineering gene regulatory networks with relevance networks, graphical gaussian models and bayesian networks[END_REF] et Allen et al. (2012). Il serait cependant difficile de hiérarchiser les méthodes d'inférence de réseaux. La dépendance existant entre deux gènes peut prendre de multiples formes et ne peut être résumée à un simple indice statistique. Combiner les résultats des différentes méthodes semble être l'approche la plus informative [START_REF] Novère | Quantitative and logic modelling of molecular and gene networks[END_REF]. Dans cette thèse, nous focalisons notre attention sur les modèles graphiques non orientés et l'inférence de réseaux de dépendance directe entre les gènes.

Les méthodes d'inférence utilisées en pratique sur les données RNA-seq

Pour inférer des réseaux de gènes à partir des données RNA-seq, [START_REF] Giorgi | Comparative study of RNA-seq and Microarray-derived coexpression networks in Arabidopsis thaliana[END_REF] ont utilisé le coefficient de corrélation de Pearson calculé sur les données RNA-seq transformées pour stabiliser la variance (Variance Stabilizing Transformation, VST). [START_REF] Hong | Canonical correlation analysis for RNA-seq co-expression networks[END_REF] ont proposé une méthode basée sur le calcul d'un coefficient de corrélation canonique. [START_REF] Iancu | Utilizing RNA-Seq data for de novo coexpression network inference[END_REF] ont utilisé la méthode WGCNA sur les données RNA-seq normalisées et transformées pour stabiliser la variance (VST). Ces méthodes appliquées sur les données RNA-seq sont donc principalement des méthodes d'inférence de réseaux de dépendances marginales : le réseau est reconstitué à partir d'un indice calculé marginalement sur toutes les paires de gènes possibles. Certaines méthodes visent cependant à reconstituer des réseaux de dépendances conditionnelles, où la dépendance entre deux gènes est calculée conditionnellement aux valeurs prises par les autres gènes. Un exemple est l'utilisation du modèle graphique gaussien sur les données RNA-seq exprimées en RPKM par [START_REF] Cai | Utilizing RNA-Seq Data for Cancer Network Inference[END_REF]. Plus récemment, [START_REF] Allen | A local Poisson graphical model for inferring networks from sequencing data[END_REF] ont proposé un modèle log-linéaire de Poisson, adapté aux données RNA-seq après transformation. Notons également que ces méthodes d'inférence de réseaux nécessitent un nombre d'échantillons relativement élevé.

Problématiques abordées : Comme en analyse différentielle et en analyse de co-expression, le choix de modélisation des données RNA-seq (discrète ou continue) est une question importante en inférence de réseaux. Comme l'analyse différentielle, l'inférence de réseaux souffre également du faible nombre d'échantillons à disposition. Nous proposerons également une solution palliative à ce problème.

Chapitre 2

Cadre statistique

Pour chaque modèle statistique, nous distinguons les unités statistiques, indexées par l'indice i = 1, . . . , n, et les variables, indexées par l'indice j = 1, . . . , p.

Classification

Les méthodes de classification

On considère une matrice y de données d'observations de taille n × p. On souhaite classer les n observations, y i pour i = 1, . . . , n, en fonction des p variables y j pour j = 1, . . . , p. Plusieurs techniques de classification non supervisée permettent de classer des observations en fonction de variables.

La classification hiérarchique

La Les modèles de mélange Ils offrent une définition rigoureuse de la notion de classes, les observations appartenant à une même classe sont des échantillons issus de variables aléatoires suivant la même loi. Les données observées sont donc la réalisation de variables aléatoires de différentes lois. Une présentation détaillée des modèles de mélange est proposée par [START_REF] Mclachlan | Finite Mixture Models, Willey Series in Probability and Statistics[END_REF].

Cette liste n'est pas exhaustive. D'autres méthodes existent pour classer des observations comme les réseaux de neurones [START_REF] Tamayo | Interpreting patterns of gene expression with self-organizing maps : methods and application to hematopoietic differentiation[END_REF], les méthodes basées sur la théorie des graphes ou les machines à vecteurs de support. Nous concentrons notre attention sur les modèles de mélange car ils offrent une définition rigoureuse de la notion de classes et un cadre probabiliste facilitant le choix de modèle.

Les modèles de mélange

On suppose qu'il existe une partition des observations y 1 , . . . , y n en K classes C 1 , . . . , C K . Cette partition est notée z = (z 1 , . . . , z n ) où z i = (z i1 , . . . , z iK ) est un vecteur binaire indiquant l'appartenance de l'observation i aux K classes : z ik = 1 si l'observation y i appartient à la classe k, et z ik = 0 sinon. Cette partition z est une donnée non observée et l'objectif de l'analyse est l'estimation de cette partition z à partir des données observées y. Pour tout i, on suppose que l'observation y i appartient à la classe C k si elle est la réalisation d'une variable aléatoire Y i conditionnellement à l'événement Z ik = 1 où Z i = (Z i1 , . . . , Z iK ) est le vecteur indiquant l'appartenance de la variable aléatoire Y i aux classes : Z ik = 1 si Y i appartient à la classe k, Z ik = 0 sinon. La distribution de la variable Y i conditionnellement à l'appartenance à la classe k suit une loi de densité f k (y i ). La variable Z i suit une loi multinomiale de paramètre (1, p), où p = (p 1 , . . . , p K ) est le vecteur des proportions de chaque composante du modèle de mélange tel que p k ∈]0, 1[ pour tout k et K k=1 p k = 1. Chaque paramètre p k de ce vecteur est égal à la probabilité que l'observation i provienne de la composante k : p(Z ik = 1). La densité jointe du couple (Y i , Z i ) s'écrit :

f (y i , z i ) = K k=1 [p k f k (y i )] z ik .
La densité marginale de la variable aléatoire Y i est un mélange de densité :

f (y i ) = K k=1 p k f k (y i ).
La distribution de la variable Z i conditionnellement à l'évenement {Y i = y i } est une loi multinomiale de paramètre (1, t i ) où t = (t i1 , . . . , t iK ) est le vecteur des probabilités conditionnelles t ik = p(Z ik = 1 | Y i = y i ). Le théorème de Bayes conduit à la formulation suivante :

t ik = p k f k (y i ) K t=1 p t f t (y i )
.

Ces probabilités conditionnelles permettent d'attribuer une classe d'appartenance à l'observation i à l'aide de la règle du maximum a posteriori (MAP) :

z ik = 1 si arg max ℓ t iℓ = k, 0 sinon.
En pratique, ces probabilités conditionnelles ne sont pas connues car les proportions (p 1 , . . . , p K ) et les densités de chaque composante du mélange f 1 (y i ), . . . , f K (y i ) ne sont pas connues. Elles doivent être estimées à partir des données observées y. Afin de réaliser cette estimation, on suppose que les couples (y 1 , z 1 ), . . . , (y n , z n ) sont indépendants et identiquement distribués et qu'ils sont la réalisation des couples de variables aléatoires (Y 1 , Z 1 ), . . . , (Y n , Z n ). On suppose également que les densités f k (y i ) s'écrivent sous la forme f k (y i ; a k ) et que les paramètres de chaque variable aléatoire (a 1 , . . . , a K ) sont tous distincts. La densité de l'observation y i s'écrit :

f (y i ; K, θ K ) = K k=1 p k f k (y i ; a k ).
Le vecteur des paramètres est θ K = (p 1 , . . . , p K , a 1 , . . . , a K ). Le choix de la distribution des variables aléatoires qui composent le mélange est important. Nous distinguons ici des modèles qui seront utilisés dans cette thèse.

Le modèle de mélange de lois gaussiennes

Le modèle de mélange de lois gaussiennes suppose que la densité des données pour une observation s'écrit :

f (y i ; K, θ K ) = K k=1 p k Φ(y i ; ν k , Σ k ),
où Φ(•) est la densité d'une loi normale multivariée de dimension p, de moyenne le vecteur ν k et de variance la matrice Σ k . Les modèles de mélange de lois gaussiennes sont identifiables à une permutation des composantes près : deux modèles de mélange ayant la même densité ont exactement les mêmes paramètres à une réindexation des composantes près. Ce n'est pas le cas des modèles de lois de Bernoulli, uniformes ou binomiales. Il existe plusieurs manières de décomposer la matrice de variance Σ k [START_REF] Celeux | Gaussian parsimonious clustering models[END_REF] :

Σ k = L k D k A k D -1 k . Le paramètre L k = Σ 1
p caractérise le volume de la composante k. D k est la matrice orthogonale des vecteurs propres de Σ k , elle indique l'orientation de la composante. Cette matrice modélise la dépendance entre les différentes variables y 1 , . . . , y p : si elle est diagonale, les variables sont supposées indépendantes. La matrice diagonale A k contient les valeurs propres de Σ k , normalisées et ordonnées par ordre décroissant. Cette matrice A k indique la forme de la composante k. En imposant que l'un ou plusieurs de ces paramètres L k , A k ou D k soient communs à toutes les composantes, on obtient une collection de 8 modèles de mélange gaussiens. Dans cette thèse, on utilise généralement le modèle p k L k ADA, noté p k L k C, où les proportions et les volumes des composantes sont libres, spécifiques à chaque composante et l'orientation et la forme des composantes sont communes à toutes les composantes du mélange.

Le modèle de mélange de lois de Poisson

Le modèle de mélange de lois de Poisson suppose que la densité des données pour une observation y i s'écrit :

f (y i ; K, θ K ) = K k=1 p k p j=1 P(y ij ; λ jk ),
où P(•; λ jk ) dénote la loi de Poisson univariée de moyenne λ jk . Ce modèle suppose l'indépendance conditionnelle des variables y 1 , . . . , y p sachant les classes. Cette hypothèse d'indépendance conditionnelle est forte et n'est généralement pas utilisée dans le cas des modèles de mélange de lois gaussiennes, pour lesquels la matrice Σ k modélise la dépendance conditionnelle entre les variables. Cependant, l'implémentation d'une loi de Poisson mutivariée n'est pas simple [START_REF] Karlis | An EM algorithm for multivariate Poisson distribution and related models[END_REF], et l'hypothèse d'indépendance facilite l'implémentation de l'estimation des paramètres.

Estimation des paramètres du modèles de mélange

De manière générale, la vraisemblance d'un modèle de mélange s'écrit :

Classification ℓ(θ K ; y) = n i=1 K k=1 p k f k (y i ; a k ) .
(2.1)

Il n'existe pas de forme explicite des estimateurs du maximum de vraisemblance des paramètres du modèle θ K . Pour effectuer l'estimation des paramètres, on considère la vraisemblance complétée des données (y, z) :

ℓ(θ K ; y, z) = n i=1 K k=1 (p k f k (y i ; a k )) z ik .
La partition z n'est pas observée. L'algorithme Espérance Maximisation, introduit par [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF], permet l'estimation des paramètres d'un modèle comportant des variables latentes non observées. Après une étape d'initialisation des paramètres du modèle θ (0) K , cet algorithme alterne deux étapes successives à chaque itération (b) :

Espérance : La première étape consiste à effectuer la maximisation de l'espérance de la log vraisemblance complétée conditionnellement aux observations y et aux paramètres courants estimés notés θ

(b) K : Q(θ K |θ (b) K ) = E log ℓ(θ K ; y, z)|y, θ (b)
K . Pour les modèles de mélange, cette étape revient à calculer les probabilités conditionnelles, notées t

K : t (b) ik = E Z ik |y, θ (b) K , = p (b) k f k (y i ; a (b) k ) K t=1 p (b) t f t (y i ; a (b) t )
.

Maximisation : La deuxième étape de l'algorithme consiste à trouver θ

(b+1) k , les pa- ramètres maximisant la quantité Q(θ K |θ (b) K ) en θ K . La maximisation de la quantité Q(θ K |θ (b) K ) garantit d'augmenter la log vraisemblance L(θ K |y).
Pour les proportions du modèle, cette étape de maximisation repose sur la formule suivante :

p (b+1) k = 1 n n i=1 t (b) ik .
Dans le cas du modèle de mélange de lois gaussiennes, les paramètres de moyennes et de variances sont mis à jour de la manière suivante pour la forme p k L k C :

ν (b+1) k = n i=1 t (b) ik y i n i=1 t (b) ik , Σ (b+1) k = n i=1 t (b) ik (y i -ν (b) k )(y i -ν (b) k ) T n i=1 t (b) ik .
Dans le cas du modèle de mélange de lois de Poisson, les paramètres de moyenne λ jk pour tout j ∈ 1, . . . , p se mettent à jour de la manière suivante :

λ (b) jk = n i=1 t (b) ik y ij n i=1 t (b) ik ( 1 p p j=1 y ij )
.

L'algorithme EM alterne les étapes d'Espérance et de Maximisation et met à jour les paramètres du modèle à chaque itération (b) jusqu'à convergence, parfois lente, vers un maximum local, qui n'est d'ailleurs pas forcement le maximum global de la fonction de vraisemblance. Diverses stratégies d'initialisation de l'algorithme ont été proposées pour améliorer la convergence de l'algorithme et la qualité de l'estimation, comme la stratégie du petit EM [START_REF] Biernacki | Choosing starting values for the EM Algorithm for getting the highest likelihood in multivariate Gaussian mixture models[END_REF]. Cette stratégie consiste à initialiser l'algorithme EM à l'aide d'un algorithme de classification de type k-means [START_REF] Macqueen | Some methods for classification and analysis of multivariate observations[END_REF], puis à effectuer un petit nombre (par exemple 10) d'itérations de l'algorithme EM. Les paramètres estimés à l'issue de ce premier lancement de l'algorithme sont utilisés pour initialiser l'algorithme EM complet. Une fois l'estimation des paramètres effectuée θK , on peut déduire la partition des observations ẑ à l'aide de la règle du MAP :

ẑik = 1 si arg max ℓ t iℓ ( θK ) = k, 0 sinon, t iℓ ( θK ) = pℓ f ℓ (y i ; âℓ ) K t=1 pt f t (y i ; ât )
.

Sélection de modèle

Jusqu'à présent, nous avons considéré que le nombre de classes K était connu. Or, sur un jeu de données réels, ce nombres de classes n'estpas connu et nous devons le déduire à partir des données observées y. Les modèles de mélange fournissent un cadre rigoureux pour choisir ce nombre de classes.

Le critère BIC

Une première solution au choix du nombre de classes K consiste à choisir le modèle qui maximise la vraisemblance intégrée :

f (y; K) = θ K f (y; K, θ K )π(θ K )dθ K ,
où π(θ K ) est une distribution a priori non informative sur les paramètres θ K . Le calcul de cette vraisemblance intégrée est rarement possible. Le critère Bayesian Information Criterion (BIC), proposé par [START_REF] Schwarz | Estimating the Dimension of a Model[END_REF], est une approximation asymptotique du logarithme de la vraisemblance intégrée, réalisée à l'aide d'une approximation de Laplace [START_REF] Lebarbier | Le critère BIC : fondements théoriques et interprétation[END_REF]). Ce critère s'écrit :

BIC(K; y) = log f (y; K, θK ) - ν K 2 log(n),
où θK est l'estimateur du maximum de vraisemblance des paramètres du modele de mélange et ν K le degré de liberté du modèle à K composantes. Ce critère BIC fournit une approximation correcte de la vraisemblance lorsque le nombre d'observations n tend vers l'infini. Le critère BIC sélectionne le modèle le plus proche de la loi des données au sens de la pseudo-distance de Kullback-Leibler KL(f, g) entre deux lois de probabilité de densités respectives f et g :

KL(f, g) = log f (x) g(x) f (x)dx.

Le critère ICL

Une alternative au critère BIC est le critère Integrated Completed Likelihood criterion (ICL) proposé par [START_REF] Biernacki | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF] :

ICL(K; y) = BIC(K; y) -Ent(K; y),
(2.2) où Ent(K; y) est un terme d'entropie qui mesure la confiance que l'on peut accorder à la classification :

Ent(K; y) = - n i=1 K k=1 t ik ( θK ) log t ik ( θK ) ≥ 0. (2.3)
En réalité, le critère ICL est une approximation du logarithme de la log-vraisemblance intégrée :

f (y, z; K) = θ K f (y, z; K, θ K )π(θ K )dθ K .
Le terme d'entropie additionnel présent dans l'écriture du critère ICL (2.3) favorise les modèles qui présentent la partition des données la plus classifiante possible.

Les données représentées en Figure 2.1 sont simulées à partir d'un mélange de quatre lois gaussiennes bivariées représentées par les quatre différentes couleurs. Le critère ICL sélectionne une partition à trois classes. Les deux composantes en croix en haut à droite de la figure sont rassemblées en une seule classe. On obtient ainsi trois classes bien distinctes. À l'inverse, le critère BIC sélectionne le modèle fournissant la meilleure approximation de la densité des données. Il identifie correctement les composantes du modèle de mélange, et détecte bien quatre composantes. Cependant, les éléments appartenant aux deux composantes en croix (symbolisées par les triangles verts et les croix violettes) sont classés dans l'une ou l'autre des composantes avec incertitude et la confiance accordée à la classification en deux classes est faible.

Inférence de réseaux

Les méthodes d'inférence de réseaux

On considère une matrice y de taille n×p correspondant à n observations y i pour i = 1, . . . , n et p variables y j pour j = 1, . . . , p. On suppose que les vecteurs y 1 , . . . , y p sont n

E = {(j, j ′ ); j ∈ V, j ′ ∈ V, j = j ′ tels que les variables Y j et Y j ′ sont dépendantes }.
Il existe plusieurs méthodes d'inférence de réseaux, dont nous exposons brièvement le principe.

Réseaux de dépendance marginale

Ces réseaux de corrélation simples sont également appelés graphes de covariance (covariance graph ou relevance network ) [START_REF] Hastie | The Elements of Statistical Learning[END_REF]. Ils modélisent des dépendances marginales en s'intéressant à la distribution jointe de chaque paire de variables. Deux variables aléatoires Y j et Y j ′ sont indépendantes marginalement si l'on peut écrire la densité jointe des deux variables aléatoires comme le produit de deux densités p(y j ; y j ′ ) = p(y j )p(y j ′ ) où p(y j ; y j ′ ), p(y j ) et p(y j ′ ) sont respectivement les densités des lois de (y j ; y j ′ ), y j et y j ′ . Les réseaux de dépendance marginales ne sont pas orientés. Dans le cas particulier où les variables aléatoires Y 1 , . . . , Y p suivent des lois normales, le vecteur

Y i = (Y 1 i , . . . , Y p i ) suit une loi normale multivariée N p (µ, Σ).
Dans un contexte de grande dimension, inférer le réseau de dépendance marginale revient à calculer l'estimateur régularisé de la matrice Σ. Le problème d'estimation de cette matrice Σ est non convexe mais des méthodes d'estimation ont été proposées [START_REF] Chaudhuri | Estimation of a covariance matrix with zeros[END_REF][START_REF] Bien | Sparse estimation of a covariance matrix[END_REF].

Réseaux basés sur l'information mutuelle

La corrélation simple entre deux variables aléatoires capture uniquement les dé-pendances linéaires entre des variables. Les réseaux basés sur l'information mutuelle I(Y j , Y j ′ ) prennent en compte des dépendances non linéaires.

I(Y j , Y j ′ ) = p(y j , y j ′ ) log p(y j , y j ′ ) p(y j ) p(y j ′ ) dy j dy j ′ .
Ces réseaux ont été étudié par [START_REF] Meyer | minet : A R/Bioconductor package for inferring large transcriptional networks using mutual information[END_REF] et [START_REF] Butte | Discovering functional relationships between RNA expression and chemotherapeutic susceptibility using relevance networks[END_REF]. Dans cette thèse, nous ne nous intéressons pas aux méthodes basées sur l'information mutuelle.

Réseaux bayésiens

Ces modèles sont aussi appelés modèles graphiques orientés. Les arêtes entre les différents noeuds sont orientées, en considérant les ordres de conditionnement entre les variables aléatoires. Par exemple, on considère trois variables aléatoires Y 1 , Y 2 , Y 3 de réalisations y 1 , y 2 , y 3 et la factorisation suivante des densités de probabilités :

p(y 1 , y 2 , y 3 ) = p(y 1 | y 2 , y 3 )p(y 2 | y 3 )p(y 3 )
Dans le réseau bayésien correspondant, il existe une arête orientée du noeud 3 vers le noeud 2 due au facteur p(y 2 | y 3 ). On dit que le noeud 3 est un parent du noeud 2 et que 2 est un noeud enfant du noeud 3. Les parents du noeud 1 sont {2, 3} : le réseau comporte donc une arête du noeud 2 vers le noeud 1 et une arête du noeud 3 vers le noeud 1. En revanche, aucune arête n'est dirigée vers le noeud 3. La structure d'un réseau bayésien est un graphe acyclique orienté. Cette structure est induite par la factorisation de la densité sous la forme suivante, où pa(y j ) désigne les parents du noeud j :

p(y 1 , . . . , y p ) = p j=1 p(y k | pa(y j )).
Ces réseaux ont été introduits par [START_REF] Pearl | Causality[END_REF] et sont l'objet de nombreux développements méthodologiques [START_REF] Friedman | Using Bayesian Networks to Analyze Expression Data[END_REF]. De manière générale, les graphes orientés sont plus compliqués à manipuler et à inférer que les graphes non orientés. On peut cependant remarquer que certaines méthodes d'inférence de réseaux non orientés peuvent servir à l'inférence de réseaux orientés [START_REF] Hastie | Statistical Learning with Sparsity : The Lasso and Generalizations[END_REF]. Nous ne nous intéressons pas à ce type de réseau dans cette thèse et concentrons notre attention sur les réseaux non orientés.

Réseaux de dépendance conditionnelle

Les modèles graphiques décrivant les relations de dépendance conditionnelle entre des variables aléatoires sont appelés des champs de Markov ou des réseaux de Markov, ou plus simplement des modèles graphiques non orientés [START_REF] Hastie | The Elements of Statistical Learning[END_REF]. Une définition de la notion de dépendance conditionnelle sousjacente aux modèles graphiques non orientés est détaillée par [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF]. On note Y A = {Y j ; j ∈ A}. Le vecteur Y = (Y 1 , . . . , Y p ) vérifie la propriété de Markov par rapport au graphe G = (V, E) s'il vérifie la proposition suivante : pour tout ensemble de noeuds S séparant le graphe en deux sous ensembles de noeuds disjoints

A et B, Y A et Y B sont indépendants conditionnellement à Y S : p(y A ; y B | y S ) = p(y A | y S )p(y B | y S ).
Par ailleurs, on définit la factorisation de la densité p du vecteur (Y 1 , . . . , Y p ) par rapport au graphe G de la manière suivante. Une clique C est un ensemble de noeuds du graphe totalement connectés entre eux : pour tout j, j ′ ∈ C, (j, j ′ ) ∈ E.

On note C l'ensemble des cliques du graphe G. Pour chaque clique C ∈ C, on définit une fonction de potentiel (ou potential function) notée ψ C et associant à y C = {y j ; j ∈ C} un nombre réel strictement positif. La distribution p se factorise sur le graphe G si l'on peut écrire :

p(y 1 , . . . , y p ) = 1 Z C∈C ψ C (y C ), (2.4

Les modèles graphiques non-orientés

Le modèle graphique gaussien

Dans le cadre du modèle graphique gaussien, les données sont un échantillon de taille n issu d'une loi normale mutivariée de dimension p, de moyenne nulle 0 et de variance Σ, matrice définie positive de taille p × p. Les observations (y 1 , . . . , y n ) sont indépendantes et identiquement distribuées :

y i ∼ N p (0, Σ) pour tout i ∈ 1, . . . , n.
Dans ce modèle, les dépendances conditionnelles entre deux variables Y j et Y j ′ conditionnellement aux autres variables sont directement liées aux coefficients de la matrice la matrice de covariance inverse, notée Θ = Σ -1 . En effet, le coefficient de corrélation partielle ρ jj ′ entre les variables j et j ′ vérifie la relation suivante :

ρ jj ′ = θ jj ′ √ θ jj θ j ′ j ′
. Un coefficient nul θ jj ′ = 0 indique l'indépendance des variables j et j ′ conditionnellement à toutes les autres variables du jeu de données [START_REF] Whittaker | Graphical Models in Applied Multivariate Statistics[END_REF][START_REF] Lauritzen | Graphical Models[END_REF]. Inférer le graphe de dépendance entre les p variables revient donc à détecter les coefficients non nuls de la matrice Θ. Ce problème est connu sous le nom de sélection de covariance, en anglais covariance selection model [START_REF] Dempster | Covariance Selection[END_REF] ou Gaussian concentration graph model. La matrice Θ est aussi appelée matrice de concentration.

Notons que l'on peut présenter cette distribution gaussienne multivariée sous la forme d'un modèle graphique à l'aide de la factorisation de la densité suivante, comme présentée à la formule (2.4) :

p(y 1 , . . . , y p ) = exp - 1 2 j,j ′ ∈V θ jj ′ y j y j ′ -A(Θ) , où A(Θ) = -1 2 log det Θ 2π .
Cette fonction A(Θ) correspond au logarithme de la fonction de partition Z(Θ). Cette écriture permet de voir que la factorisation du vecteur aléatoire Y par rapport au graphe G est déterminée par les paramètres θ jj ′ . La structure du graphe à inférer est donc totalement déterminée par Θ.

Une première approche pour obtenir la structure du graphe consiste à calculer S la matrice de variance covariance empirique et d'inverser cette matrice afin d'obtenir une estimation de la matrice Θ :

S = 1 n n i=1 (y i -ȳ)(y i -ȳ) T ,
où ȳ est le vecteur de moyenne empirique. Lorsque le nombre d'échantillons est inférieur au nombre de variables, cette approche n'est pas possible puisque la matrice S n'est pas inversible. Une autre approche consiste à tester l'inclusion ou l'exclusion des arêtes du graphe à l'aide d'un algorithme de type forward ou backward [START_REF] Dempster | Covariance Selection[END_REF] mais cette approche est coûteuse et peu efficace. En grande dimension (cas n < p), on distingue deux grands types de méthodes d'inférence du graphe : certaines méthodes sont basées sur des tests multiples [START_REF] Schäfer | A shrinkage approach to large-scale covariance matrix estimation and implications for functional genomics[END_REF][START_REF] Wille | Low-order conditional independence graphs for inferring genetic networks low-order conditional independence graphs for inferring genetic networks[END_REF] tandis que d'autres méthodes sont basées sur des méthodes de régularisation ℓ 1 [START_REF] Meinshausen | High-dimensional graphs and variable selection with the lasso[END_REF][START_REF] Huang | Combining gene annotations and gene expression data in model-based clustering : weighted method[END_REF][START_REF] Banerjee | Model selection through sparse maximum likelihood estimation for multivariate gaussian or binary data[END_REF][START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF]. Nous concentrons notre attention sur ces méthodes d'inférence de graphe à l'aide de pénalisation ℓ 1 , en détaillant d'abord le lasso graphique (ou graphical lasso) basé sur le calcul de la vraisemblance pénalisée exacte, puis la sélection au voisinage d'un noeud (ou neighbourhood selection) basée sur le calcul de régressions pénalisées.

Le lasso graphique

On considère la log-vraisemblance pénalisée suivante, où λ est le paramètre de régularisation de la pénalité ℓ 1 imposée sur les éléments de la matrice Θ à estimer.

L λ (Θ) = log det(Θ) -trace(SΘ) -λ || Θ || 1 .
(2.5)

Le problème d'estimation de la matrice Θ à partir de cette vraisemblance pénalisée est connu sous le nom du lasso graphique ou graphical lasso. Il s'agit d'un problème d'estimation convexe. Pour le résoudre de manière efficace sur des problèmes de grande dimension, [START_REF] Banerjee | Model selection through sparse maximum likelihood estimation for multivariate gaussian or binary data[END_REF][START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF] ont proposé un algorithme de descente du gradient par bloc. Pour le choix du paramètre de régularisation λ, [START_REF] Banerjee | Model selection through sparse maximum likelihood estimation for multivariate gaussian or binary data[END_REF] ont proposé un critère de sélection de modèle. Pour un seuil fixé α ∈]0, 1[, le choix de paramètre de régularisation λ(α) garantit que la probabilité de joindre par erreur deux composantes connexes disjointes du réseau soit bornée par α :

λ(α) = (max j ′ >j σj ′ σj ) t n-2 (α/2p 2 ) n -2 + t 2 n-2 (α/2p 2 )
, où t n-2 (α) désigne le quantile d'ordre (100α) de la distribution de student à n -2 degrés de liberté, et σj est la variance empirique de la variable j.

D'autres choix de paramètres de régularisation ont été proposés (Stability Approach to Regularization Selection criterion (StARS), extendedBIC ...) et les développements méthodologiques autour de ce choix sont nombreux.

Sélection de voisinage (ou neighbourhood selection)

Une méthode alternative au lasso graphique est la sélection d'arêtes au voisinage de chaque noeud du graphe à l'aide de régression pénalisée. [START_REF] Meinshausen | High-dimensional graphs and variable selection with the lasso[END_REF] ont utilisé le lien existant entre les coefficients dans la régression de la variable j sur les autres variables du réseau et les coefficients de la matrice Θ :

y j = j ′ =j β jj ′ y j ′ + ǫ j , β j k = θ jk θ jj .
Les paramètres β jj ′ sont les coefficients dans la régression de y j sur les autres variables et ǫ j ∼ N (0, σ j ). Dans le contexte de grande dimension, [START_REF] Meinshausen | High-dimensional graphs and variable selection with the lasso[END_REF] proposent d'effectuer des régressions pénalisées de chaque variable du réseau sur les autres variables afin d'identifier les coefficients θ jj ′ non nuls. Les coefficients θ jj ′ sont estimées non nuls si au moins un des coefficients β jj′ et β j′j est non nul (OR rule), ou, de manière alternative, si ces deux coefficients sont non nuls (AND rule).

Pour le choix du paramètre de régularisation λ dans chaque régression pénalisée, le paramètre λ(α), pour α ∈]0, 1[, garantit que la probabilité de joindre par erreur deux composantes connexes disjointes du réseau soit borné par α :

λ(α) = 2σ j √ n Φ-1 α 2pn 2 ,
où Φ = 1 -Φ, Φ est la distribution cumulée de la la normale N (0, 1) et σj est la variance empirique de la variable j. Cette procédure est asymptotiquement consistante pour l'estimation des éléments non nuls de la matrice Θ [START_REF] Meinshausen | High-dimensional graphs and variable selection with the lasso[END_REF].

Les modèles graphiques discrets

On considère un ensemble de variables aléatoires Y 1 , . . . , Y p prenant ses valeurs dans un ensemble discret {0, 1}. Ce modèle graphique discret, aussi appelé modèle d'Ising, est un exemple de modèle graphique discret. La factorisation de la densité suivante par rapport au graphe G s'écrit de la manière suivante :

p(y 1 , . . . , y p ) = exp   j∈V θ jj y j + (j,j ′ )∈E θ jj ′ y j y j ′ -A(Θ)   , A(Θ) = log   y∈{0,1} p exp   j∈V θ jj y j + (j,j ′ )∈E θ jj ′ y j y j ′     .
Ce modèle graphique discret est aussi appelé modèle d'Ising. Contrairement au modèle graphique gaussien, pour lequel la fonction de partition A(Θ) se calcule bien, le calcul de A(Θ) dans le cas discret implique une sommation sur 2 p termes. Une estimation des paramètres Θ similaire à celle du lasso graphique est donc impossible. La méthode d'inférence de réseaux basée sur la sélection de voisinage est donc particulièrement intéressante dans le cas discret. Pour le modèle d'Ising, cette estimation revient à effectuer, pour chaque variable y j la régression logistique suivante, en ajoutant une pénalité ℓ 1 pour forcer la nullité de certains coefficients Θ :

log p(y j = 1 | y V \j ) p(y j = 0 | y V \j ) = θ jj + j ′ ∈V,j ′ =j θ jj ′ y j ′ .
Chapitre 3

Contributions

Notations

Dans cette thèse, l'indexation des données y dépend du modèle statistique considéré, en conservant les conventions suivantes : l'indice i = 1, . . . , n indexe toujours les unités statistiques ou observations, l'indice j = 1, . . . , p indexe toujours les variables. En fonction du type de modèle considéré, l'indice j indexe les gènes (analyse différentielle, inférence de réseau) ou les échantillons (classification par modèle de mélange).

Filtrage des données RNA-seq

On rappelle que l'on dispose d'une matrice y de taille n× p correspondant à n échantillons y i pour i = 1, . . . , n pour p gènes y j pour j = 1, . . . , p. La quantité y ij représente l'expression du gène j pour l'échantillon i. Dans le cadre de l'analyse différentielle des données d'expression de gènes (section 1.2.2), l'étape de correction des p-values pour les tests multiples dans la détection de gènes différentiellement exprimés est une étape importante afin de garder un taux de détection de faux positifs relativement bas. Cependant, ce contrôle du nombre de faux positifs se fait au détriment de la capacité du test à détecter les gènes réellement différentiellement exprimés, aussi appelée puissance de détection. Pour cette raison, il est essentiel de retirer, avant toute analyse différentielle, les gènes dont le profil d'expression ne parait pas informatif. On considère l'exemple suivant (tableau 3.2) : le gène 4 comporte un alignement de 48 reads pour l'échantillon 3 de la condition 1, et uniquement pour cet échantillon. Ce gène ne semble pas apporter d'information et pourra difficilement être considéré comme différentiellement exprimé entre la condition 1 et la condition 2. Il serait donc raisonnable de supprimer le gène 4 avant d'effectuer l'analyse différentielle. gène 1 gène 2 gène 3 gène 4 gène 5 gène 6 . . . Les méthodes de filtrage utilisées en pratique consistent à supprimer les gènes :

-dont la somme totale des reads est inférieure à un certain seuil fixé arbitrairement [START_REF] Sultan | of the Human Transcriptome[END_REF]) -ayant au moins un comptage nul dans chaque condition expérimentale [START_REF] Bottomly | Evaluating gene expression in C57BL/GJ and DBA/2J mouse striatum using RNA-seq and microarrays[END_REF], -dont la valeur moyenne exprimée en RPKM est inférieure à une certaine valeur fixée arbitrairement [START_REF] Mortazavi | Mapping and quantifying mammalian transcriptomes by RNA-Seq[END_REF], -dont la valeur maximale exprimée en CPM est inférieure à une certaine valeur fixée arbitrairement [START_REF] Robinson | A scaling normalization method for differential expression analysis of RNA-seq data[END_REF] . Toutes ces méthodes sont basées sur le choix d'un seuil à fixer arbitrairement. Nous avons proposé un niveau de seuil des données qui ne soit pas arbitrairement choisi. Ce seuil est ajusté aux données et correspond à la valeur maximisant un indice de similarité entre les réplicats biologiques d'une même condition. On considère deux réplicats biologiques i et i ′ appartenant à une même condition expérimentale, un seuil fixé s et le tableau de contingence tableau 3.3. Pour un seuil donné, on calcule un indice de Jaccard, coefficient de similarité entre les deux réplicats i et i ′ à partir du tableau de contingence tableau 3.3 :

J s (y i , y i ′ ) = a a + b + c .
Réplicat i Nombre de gènes j tels que y ij > s Nombre de gènes j tels que y ij ≤ s Réplicat i ′ Nombre de gènes j tels que y i ′ j > s a b

Nombre de gènes j tels que y i ′ j ≤ s c d

Table 3.3 -Tableau de contingence du nombre de gènes dont le comptage normalisé est inférieur/ supérieur au seuil s pour deux échantillons i et i ′ .

On note T le cardinal de l'ensemble des couples d'échantillons possibles appartenant à une même condition. On calcule la moyenne de ces coefficients sur l'ensemble des couples d'échantillons appartenant à une même condition possible :

J ⋆ s (y) = 1 T i<i ′ c(i)=c(i ′ ) J s (y i , y i ′ ).
Le seuil à choisir, guidé par les données, est celui qui maximise la quantité J ⋆ s (y) :

s ⋆ = argmax s J ⋆ s (y).
Le choix du seuil selon cette méthode garantit que la similarité entre les profils d'expression des échantillons d'une même condition expérimentale soit la plus forte possible. Le seuil fixé est spécifique à chaque jeu de données. Les caractéristiques statistiques et pratiques de ce seuil sont détaillés dans le chapitre 4 qui correspond à l'article publié suivant [START_REF] Rau | Data-based filtering for replicated high-throughput transcriptome sequencing experiments[END_REF].

Classification des données d'expression par modèle de mélange

On rappelle que l'on dispose d'une matrice de mesures d'expression de n gènes y = (y 1 , . . . , y n ). Pour chaque gène i = 1, . . . , n, le vecteur y i indique l'expression du gène i pour les p conditions expérimentales j = 1, . . . , p. On souhaite classer les gènes en fonction de leur niveau d'expression à travers les différentes conditions expérimentales à l'aide d'un modèle de mélange. Dans ce contexte, les unités statistiques sont les gènes y 1 , . . . , y n tandis que les échantillons y 1 , . . . , y p correspondent aux variables du modèle (voir tableau 3.1). Contrairement à l'analyse différentielle ou à l'inférence de réseaux, le faible nombre d'échantillons disponibles (p petit) ne pose pas de problème statistique puisque nous classons les gènes et non les échantillons.

Transformation des données et comparaison de modèles pour la classification des données RNA-seq

Comme détaillé en section 2.1.2, on suppose que les données y sont la réalisation d'un mélange de K variables aléatoires :

f (y i ; K, θ K ) = K k=1 p k f k (y i ; a k ).
Le choix de la densité f k est une question importante. Le tableau 3.4 récapitule les différentes lois et méthodes utilisées sur les données de puces à ADN et les données RNA-seq. Pour l'analyse différentielle, [START_REF] Robinson | A scaling normalization method for differential expression analysis of RNA-seq data[END_REF] et [START_REF] Anders | Differential expression analysis for sequence count data[END_REF] ont d'abord proposé des modèles basés sur des lois discrètes, puis [START_REF] Law | Voom : precision weights unlock linear model analysis tools for RNA-seq read counts[END_REF] ont proposé une transformation des données RNA-seq et l'utilisation du modèle basé sur des lois gaussiennes initialement développé pour l'analyse différentielle des données de puces à ADN. Le choix entre l'utilisation de modèle discret (edgeR, DESeq2) ou continu (lima + "voom") dépend principalement du design de l'expérience et du nombre de réplicats disponibles.

Pour l'analyse de co-expression des données RNA-seq, [START_REF] Rau | Coexpression analysis of high-throughput transcriptome sequencing data with Poisson mixture models[END_REF] ont proposé un modèle de mélange de lois de Poisson. Nous proposons une transformation simple des données RNA-seq permettant l'utilisation du modèle de mélange gaussien, modèle bien établi pour l'analyse des données de puces à ADN [START_REF] Yeung | Modelbased clustering and data transformations for gene expression data[END_REF]. Afin de choisir entre l'utilisation du modèle gaussien sur données transformées et l'utilisation du modèle de Poisson, nous proposons de calculer le critère BIC de chaque modèle et de comparer ces BIC. Le calcul du critère BIC est effectué en prenant en compte la transformation des données. La prise en compte de la transformation des données dans le calcul du BIC a déjà été proposé par [START_REF] Thomas | The Fractal Morphology of the Built-Up Landscape[END_REF] en analyse de données d'urbanisme. La prise en compte de la transformation dans le calcul a un intérêt très pratique pour les données RNA-seq. Nous illustrons l'utilité de cette comparaison de modélisation discrète ou continue dans le chapitre 5. Ce chapitre a fait l'objet d'une communication orale aux 47èmes Journées de Statistique de la SFdS. de mélange (Yeung & al, 2000, Biernacki & al, 2006) ou Gaussienne sur données transformées ? Table 3.4 -Tableau récapitulatif des lois utilisées et des packages R implémentant les méthodes correspondantes pour l'analyse différentielle et la classification par modèle de mélange.

3.2.2 Un critère de sélection de modèle pour la classification par modèle de mélange de données annotées d'expression de gènes.

En classification des données RNA-seq par modèle de mélange, le problème statistique est inversé par rapport aux problèmes d'analyse différentielle ou d'inférence de réseaux : les observations à classer y 1 , . . . , y n sont les gènes et les variables y 1 , . . . , y p sont les échantillons. Dans ce contexte, le manque d'échantillons, caractéristique des jeux de données RNA-seq à notre disposition, ne pose donc pas de problème statistique. Ce faible nombre d'échantillons (i.e. de variables) nous incite cependant à chercher des variables informatives externes sur la nature des gènes afin d'améliorer la qualité de la classification. Ces variables informatives externes peuvent être extraites des bases de données d'annotations de gènes, décrites en section 1.2.3. Une première approche consiste à inclure ces données directement dans le modèle en les considérant comme des variables au même titre que les échantillons. Une autre approche consiste à utiliser ces annotations pour améliorer l'estimation des paramètres du modèle de mélange [START_REF] Pan | Incorporating gene functions as priors in model-based clustering of microarray gene expression data[END_REF][START_REF] Huang | Combining gene annotations and gene expression data in model-based clustering : weighted method[END_REF]. Cependant, ces annotations sont généralement utilisées pour valider a posteriori la pertinence des modules détectés. Nous préférons donc utiliser ces variables externes uniquement dans l'étape de sélection de modèle, et non dans l'étape d'estimation des paramètres du modèle de mélange. Une première solution pour intégrer ces annotations externes dans l'étape de sélection de modèle est l'utilisation du critère SICL. Nous détaillons le principe du critère SICL ci-dessous.

Le critère SICL [START_REF] Baudry | Enhancing the selection of a model-based clustering with external categorical variables[END_REF] a proposé un critère de type ICL prenant en compte plusieurs variables catégorielles externes pouvant éventuellement contribuer à expliquer la partition des données inférée. On note ces variables externes u = (u 1 , . . . , u R ) avec r = 1, . . . , R 3.2.2 -Un critère de sélection de modèle pour la classification par modèle de mélange de données annotées d'expression de gènes. 45 où u r iℓ = 1 indique que l'observation i appartient à la catégorie ℓ de la r ième variable externe catégorielle, u r iℓ = 0 sinon. On souhaite obtenir une classification z à partir des données y telle que la classification z soit la plus cohérente possible avec les variables catégorielles externes u. En supposant que y et u sont conditionnellement indépendants sachant la classification z, le critère de classification supervisée Supervised Integrated Completed Likelihood (SICL) est une approximation asymptotique du logarithme de la log-vraisemblance complétée :

f (y, u, z; K) = f (y, u, z; K, θ K )π(θ K )dθ K .
Le critère SICL est défini de la manière suivante : rassemblées une seule classe. Le critère SICL, prenant en compte cette information externe, sélectionne une classification à 4 classes, isolant la composante en bas à gauche de la figure (cercles) des points en haut à gauche de la classification (triangles). La partition sélectionnée par le critère SICL est ainsi plus facilement interprétable en fonction de la classification externe symbolisée par les formes des points.

SICL(K; y) = ICL(K) + R r=1 Ur ℓ=1 K k=1 n r kℓ log n r kℓ n k• , ( 3

Le critère ICAL

Même si les termes d'annotations sont des variables binaires (i.e. un gène est annoté ou n'est pas annoté), les bases de données d'annotation fonctionnelle sont souvent incomplètes ce qui rend l'utilisation du critère SICL peu approprié. On considère un ensemble de G termes d'annotation extrait de ces bases de données d'annotation. Ces termes sont indexés par l'indice g, l'information extraite des bases de données pour chaque terme g est notée u g où :

u g i = 1 si
b g ik = 1 avec probabilité p g k if u g i = 1, 0 si u g i = 0.
Le critère ICAL est une approximation de la log-vraisemblance complétée annotée des données :

log f (y, z, b 1 , . . . , b G ; K) = log θ K f (y, z, b 1 , . . . , b G ; K, θ K )π(θ K )dθ K . On définit n g = card{i : u g i = 1} et n g k = card{i : ẑik = 1 et u g i = 1}. Ce critère s'écrit : ICAL(K; y) = ICL(K) + G g=1 K k=1 n g k log n g k n g .
Ce critère est présenté au chapitre 6 et fait l'objet d'une publication acceptée [START_REF] Gallopin | A model selection criterion for model-based clustering of annotated gene expression data[END_REF].

3.3.1 -Un modèle pour l'inférence de réseaux RNA-seq 47

Inférence de réseaux à l'aide de modèle graphique

On rappelle que l'on dispose d'une matrice y de taille n × p correspondant à n échantillons y i pour i = 1, . . . , n pour p gènes y j pour j = 1, . . . , p. On suppose que les vecteurs y 1 , . . . , y p sont n réalisations issus de variables aléatoires Y 1 , . . . , Y p représentant le niveau d'expression des gènes. On souhaite inférer les dépendances entre les p variables aléatoires à partir des n échantillons y 1 , . . . , y n . Pour représenter et visualiser ces dépendances, on considère le graphe G = (V, E) où V = {1, . . . , p} est l'ensemble des noeuds représentant les variables aléatoires Y 1 , . . . , Y p et E ⊂ V × V est l'ensemble des arêtes du graphe représentant les dépendances entre les variables aléatoires :

E = {(j, j ′ ); j ∈ V, j ′ ∈ V, j = j ′ tels que les variables Y j et Y j ′ sont dépendantes }.
3.3.1 Un modèle pour l'inférence de réseaux RNA-seq L'utilisation du modèle graphique gaussien est bien établi en inférence de réseaux à partir de données de puces à ADN. Comme pour l'analyse différentielle et la classification, les modèles d'inférence de réseaux doivent être adaptés au caractère discret des données RNA-seq. Dans le cas du modèle graphique de Poisson, la factorisation de la distribution du vecteur aléatoire (Y 1 , . . . , Y p ) par rapport au graphe G est la suivante :

p(y 1 , . . . , y p ) = exp   j∈V (β j y j -log y j ! ) + (j,j ′ )∈E β jj ′ y j y j ′ -A(β)  
où A(β) est le logarithme de la fonction de partition Z(β) garantissant que p soit une distribution (voir (2.4)). Comme dans le cas du modèle d'Ising, le calcul de A(β) est impossible ne permettant pas l'estimation des paramètres β à partir de la vraisemblance pénalisée exacte, comme détaillé en section 2.2.2.

Une approche alternative est la sélection des arêtes aux voisinage de chaque noeud du réseau (neighbourhood selection). [START_REF] Allen | A local Poisson graphical model for inferring networks from sequencing data[END_REF] ont proposé l'utilisation de régression linéaire généralisée de Poisson. Ce modèle prend en compte le caractère discret des données RNA-seq, mais il ne prend pas en compte la grande dispersion inter-échantillons. Afin de prendre en compte cette dispersion, nous avons proposé de remplacer la régression linéaire généralisée de Poisson par une régression linéaire généralisée mixte, basée sur une loi de Poisson hiérarchique Y j ∼ P(m j ) où le paramètre de moyenne m j est également une variable aléatoire : m j = µ j exp(ε j ) avec ε j ∼ N n (0, σ 2 j I n ) et µ j est un vecteur de taille n dont chaque entrée est exp( j ′ =j β jj ′ y ij ′ ) pour i ∈ 1, . . . , n. Contrairement au modèle graphique gaussien ou au modèle graphique log-linéaire de Poisson, le modèle hiérarchique log-linéaire de Poisson proposé ne requiert pas de transformation préalable des données pour son utilisation sur les données RNA-seq.

Ce travail est exposé dans le chapitre 7 et a fait l'objet d'une publication [START_REF] Gallopin | A Hierarchical Poisson Log-Normal Model for Network Inference from RNA Sequencing Data[END_REF].

Sélection de covariance diagonale par bloc pour le modèle

graphique gaussien en grande dimension.

Le faible nombre d'échantillons disponibles pour les jeux de données RNA-seq limite fortement la performance de tout type de méthode d'inférence de réseaux (modèle graphique gaussien ou modèles graphiques discrets). Ce constat correspond à un résultat déjà connu pour le cas de la régression linéaire gaussienne. Dans ce cadre, Verzelen (2012) a défini les cas d'ultra-grande dimension de la manière suivante, où p est le nombre de gènes, n le nombre d'échantillons et d le degré maximal du réseau, i.e. le nombre maximal de connexions entre un gène et les autres gènes du réseau :

d log( p d ) n ≥ 1 2 .
Dans le contexte d'ultra-grande dimension, il convient donc de réduire le plus possible le nombre de variables à inclure dans le réseau afin de réduire le nombre de paramètres à estimer. Ce constat a été à l'origine de la contribution suivante.

Dans le cadre du spécifique du modèle graphique gaussien, on estime la matrice Θ à l'aide du graphical lasso [START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF], afin d'inférer la structure du graphe de dépendances conditionnelles correspondant, comme détaillé en section 2.2.2 :

ℓ λ (Θ) = log det(Θ) -trace(SΘ) -λ || Θ || 1 (3.2) Θ(λ) = argmin Θ {ℓ λ (Θ)} .
Afin de réduire le nombre de paramètres à estimer, Mazumder and Hastie (2012) et [START_REF] Witten | New insights and faster computations for the Graphical Lasso[END_REF] ont remarqué la propriété suivante : si la solution Θ(λ) du problème 3.2 est diagonale par bloc à une permutation des variables près, alors le problème d'inférence peut être résolu dans chaque bloc indépendamment. De plus, il existe un moyen simple de savoir si la solution Θ(λ) a une structure diagonale par bloc : Θ(λ) et la matrice de covariance empirique seuillée S (λ) = (σ jj ′ 1 {|σ jj ′ |>λ} ) jj ′ ont la même structure diagonale par bloc. Ce constat est à l'origine de la règle du screening diagonal par bloc (block diagonal screening rule) pour une valeur fixée du paramètre de régularisation λ :

Etape 1 Identification des blocs, correspondant aux composantes connexes du graphe inféré à partir de la matrice d'adjacence (1

{|σ jj ′ |>λ} ) jj ′ , Etape 2 Estimation des paramètres Θ(λ) 1 , . . . , Θ (λ)
K dans chacun des blocs détectés à l'étape 1.

Le paramètre de régularisation λ utilisé à l'étape 1 est le même que les paramètres de régularisation utilisés pour l'estimation des paramètres dans chaque bloc à l'étape 2. Afin d'améliorer les estimations, [START_REF] Tan | The Cluster Graphical Lasso for improved estimation of Gaussian graphical models[END_REF] ont proposé le Cluster Graphical Lasso : à l'étape 1, ils identifient K composantes connexes à l'aide d'un clustering hiérarchique réalisé à partir de la matrice de similarité S =| S |. Puis ils infèrent les réseaux dans chaque bloc en utilisant différents paramètres de régularisation ρ 1 , . . . , ρ K . Ils montrent ainsi que l'utilisation de différents paramètres de régularisation λ à l'étape 1 et dans les sous-problèmes de l'étape 2 permet d'améliorer les performances de l'inférence.

Dans cet esprit, nous proposons d'inférer le réseau en deux étapes en utilisant des paramètres de régularisation spécifiques à chaque étape. Nous concentrons notre attention 3.3.2 -Sélection de covariance diagonale par bloc pour le modèle graphique gaussien en grande dimension. 49

sur le choix du paramètre de régularisation dans l'étape 1. Nous notons B l'ensemble des partitions possibles des p variables en différents blocs. Nous reformulons le choix du seuil à appliquer à la matrice de covariance S en un problème de sélection de modèle parmi la collection de modèle suivante F = (F B ) B∈B :

F B =      φ(0,Σ B ) où Σ B =       Σ 1 0 0 0 . . . 0 0 0 Σ K       ,Σ k matrice
fB = φ(0, ΣB ) où ΣB =    Σ1 0 0 0 . . . 0 0 0 ΣK    et Σk = S k pour k ∈ {1, . . . , K}.
Une exploration de toutes les partitions de B n'est pas possible. Comme suggéré par la règle du screening diagonale par bloc, nous concentrons notre attention sur le sous-ensemble de partitions suivant, où Λ = (σ j,j ′ ) j>j ′ :

B Λ = B(λ), λ ∈ Λ où B(λ) sont les composantes connexes du graphe (1 {|σ j,j ′ |>λ} ) jj ′ .
On sélectionne la partition B suivant un critère de vraisemblance pénalisée, où 1 n n i=1 log(ŝ B (y i )) est la log-vraisemblance du modèle et pen un terme de pénalité à définir :

B = argmin B∈B Λ - 1 n n i=1 log(ŝ B (y i )) + pen(B) . (3.3)
Le critère BIC [START_REF] Schwarz | Estimating the Dimension of a Model[END_REF] correspond au terme de pénalité suivant : pen(B) = log n 2 D B où D B est le nombre de paramètres du modèle F B . Ce critère est asymptotique consistant : si les données sont issues d'une variable aléatoire de densité f ⋆ = φ(0, Σ B ⋆ ) alors BIC sélectionne la vraie partition des données B ⋆ lorsque le nombre d'observations n tend vers l'infini. Dans notre contexte, le nombre d'observations n est limité. Pour cette raison, nous considérons un critère de sélection non-asymptotique basé sur l'heuristique de pente proposé par [START_REF] Birgé | Minimal penalties for Gaussian model selection[END_REF]. Cette heuristique consiste à prendre le choix de pénalité suivant : pen(B) = κD B , où κ est un coefficient à déterminer. Ce terme de pénalité dépend non seulement du nombre de paramètres définissant chaque modèle mais également de la complexité de la collection de modèles considérée. Idéalement, on aimerait choisir le modèle qui minimise le risque de Kullback-Leiber avec la vraie distribution de densité f ⋆ . On note fB(f ⋆ ) cet estimateur. Il est inconnu car f ⋆ est inconnue. C'est le modèle oracle que l'on ne peut atteindre en pratique. Mais on peut avoir une procédure qui est telle que l'on contrôle le risque de Kullback-Leibler entre f ⋆ et le modèle choisi f B . Le modèle sélectionné f B grâce à cette pénalité vérifie la borne oracle :

R(f ⋆ , f B ) ≤ CR(f ⋆ , fB(f ⋆ ) ),
où C est une constante. Nous disposons également d'une borne minimax sur la qualité de l'estimation qui minore le risque R(f ⋆ , f B ) par un terme du même ordre que la borne obtenue dans la borne oracle, à un terme logarithmique près.

Le calcul du coefficient κ peut être effectué explicitement en fonction des conditions imposées sur les paramètres des modèles considérés. En pratique, il est plus simple de calibrer ce paramètre à partir des données à l'aide des deux méthodes détaillées par [START_REF] Baudry | Slope heuristics : overview and implementation[END_REF]. Ces deux méthodes consistent à détecter le coefficient κ à l'aide d'une régression robuste réalisée entre la log-vraisemblance et la dimension du modèle pour les modèles complexes où à l'aide du saut de dimension détecté sur la fonction représentant la dimension du modèle en fonction du paramètre κ à calibrer.

Cette méthode d'inférence de réseaux, ainsi que les résultats théoriques et illustrations sur des données simulées et réelles, sont détaillés au chapitre 8. Ce chapitre fait l'objet d'un article en cours d'écriture avec Emilie Devijver.

Chapitre 4

Data-based filtering for replicated high-throughput transcriptome experiments Résumé. Le séquençage haut-débit RNA-seq est très utilisé en analyse différentielle des données d'expression des gènes à travers différentes conditions. Les tests sont réalisés sur un grand nombre de gènes ce qui implique un contrôle très stringent du nombre de faux positifs au détriment de la puissance de détection des tests. Afin de modérer l'impact de cette correction pour les tests multiples, les données sont d'abord filtrées afin de retirer de l'analyse les gènes peu exprimés ou présentant peu de signal. Cependant, l'impact de ces filtres sur l'analyse différentielle n'a pas été beaucoup étudié. Nous proposons une méthode de filtrage dont le seuil est calibré à l'aide d'un indice de Jaccard mesurant la similarité entre les différents réplicats d'une même condition, calculé sur les données. Par comparaison avec d'autres méthodes de filtrage régulièrement utilisées sur les données RNA-seq, notre méthode filtre correctement les gènes peu exprimés et contribue à augmenter la puissance de détection des gènes moyennement ou fortement exprimés. De plus, le seuil de filtre, calibré à partir des données, varie en fonction du jeu de données considéré, ce qui souligne l'intérêt de notre méthode. La méthode proposée est implémentée dans le package R package HTSFilter disponible sur Bioconductor.

Abstract. RNA sequencing is now widely performed to study differential expression among experimental conditions. As tests are performed on a large number of genes, very stringent false discovery rate control is required at the expense of detection power. Ad hoc filtering techniques are regularly used to moderate this correction by removing genes with low signal, with little attention paid to their impact on downstream analyses. We propose a data-driven method based on the Jaccard similarity index to calculate a filtering threshold for replicated RNA-seq data. In comparisons with alternative data filters regularly used in practice, we demonstrate the effectiveness of our proposed method to correctly filter lowly expressed genes, leading to increased detection power for moderately to highly expressed genes. Interestingly, this data-driven threshold varies among experiments, highlighting the interest of the method proposed here. The proposed filtering method is implemented in the R package HTSFilter available on Bioconductor.

Introduction

Over the past five years, next-generation high-throughput sequencing (HTS) technology has become an essential tool for genomic and transcriptomic studies. In particular, the use of HTS technology to directly sequence the transcriptome, known as RNA sequencing (RNA-seq), has revolutionized the study of gene expression by opening the door to a wide range of novel applications. Unlike microarray data, which are continuous, RNA-seq data represent highly heterogeneous counts for genomic regions of interest (typically genes), and often exhibit zero-inflation and a large amount of overdispersion among biological replicates ; as such, a great deal of methodological research (e.g., [START_REF] Anders | Differential expression analysis for sequence count data[END_REF]Robinson et al., 2010;[START_REF] Dillies | A comprehensive evaluation of normalization methods for Illumina high-throughput RNA sequencing data analysis[END_REF] has recently focused on appropriate normalization and analysis techniques that are adapted to the characteristics of RNA-seq data ; see [START_REF] Oshlack | From RNA-seq reads to differential expression results[END_REF] for a review of RNA-seq technology and analysis procedures.

As with data arising from previous technologies, such as microarrays or serial analysis of gene expression (SAGE), HTS data are often used to conduct differential analyses. In recent years, several approaches for gene-by-gene tests using gene-level HTS data have been proposed, with the most popular making use of Poisson [START_REF] Wang | DEGseq : an R package for idenetifying differentially expressed genes from RNA-seq data[END_REF], overdispersed Poisson [START_REF] Auer | A two-stage Poisson model for testing RNA-seq data[END_REF], or negative binomial distributions [START_REF] Anders | Differential expression analysis for sequence count data[END_REF]Robinson et al., 2010). Because a large number of hypothesis tests are performed for gene-by-gene differential analyses, the obtained p-values must be adjusted to address the fact that many truly null hypotheses will produce small p-values simply by chance ; to address this multiple testing problem, several well-established procedures have been proposed to adjust p-values in order to control various measures of experiment-wide false positives, such as the false discovery rate (FDR). Although such procedures may be used to control the number of false positives that are detected, they are often at the expense of the power of an experiment to detect truly differentially expressed (DE) genes, particularly as the number of genes in a typical HTS dataset may be in the thousands or tens of thousands. To reduce this impact, several authors in the microarray literature have suggested the use of data filters in order to identify and remove genes which appear to generate an uninformative signal [START_REF] Bourgon | Independent filtering increases detection power for high-throughput experiments[END_REF] and have no or little chance of showing significant evidence of differential expression ; only hypotheses corresponding to genes that pass the filter are subsequently tested, which in turn tempers the correction needed to adjust for multiple testing.

In recent work, [START_REF] Bourgon | Independent filtering increases detection power for high-throughput experiments[END_REF] advocate for the use of independent data filtering, in which the filter and subsequent test statistic pairs are marginally independent under the null hypothesis and the dependence structure among tests remains largely unchanged pre-and post-filter, ensuring that post-filter p-values are indeed true p-values. For such an independent filter to be effective, it must be positively correlated with the test statistic under the alternative hypothesis ; indeed, it is this correlation that leads to an increase in detection power after filtering. In addition, Bourgon et al. demonstrate that non-independent filters for which dependence exists between the filter and test statistic (e.g., making use of condition labels to filter genes with average expression in at least one condition less than a given threshold), can in some cases lead to a loss of control of experiment-wide error rates.

Several ad hoc data filters for RNA-seq data have been used in recent years, including 4.2.1 -Types of filters used for RNA-seq data 57 filtering genes with a total read count smaller than a given threshold [START_REF] Sultan | of the Human Transcriptome[END_REF] and filtering genes with at least one zero count in each experimental condition [START_REF] Bottomly | Evaluating gene expression in C57BL/GJ and DBA/2J mouse striatum using RNA-seq and microarrays[END_REF] ; however, selecting an arbitrary threshold value to filter genes in this way does not account for the overall sequencing depth or variability of a given experiment. One exception to these ad hoc filters is the work of [START_REF] Ramsköld | An abundance of ubiquitously expressed genes revealed by tissue transcriptome sequence data[END_REF], in which a comparison between expression levels of exonic and intergenic regions was used to find a threshold for detectable expression above background in various human and mouse tissues, where expression was estimated as Reads Per Kilobase per Million mapped reads (RPKM) [START_REF] Mortazavi | Mapping and quantifying mammalian transcriptomes by RNA-Seq[END_REF]. The threshold of 0.3 RPKM identified in this work has in turn been applied to several other studies (e.g., [START_REF] Łabaj | RNA-seq precision in quantitative expression profiling[END_REF][START_REF] Cánovas | SNP discovery in the bovine milk transcriptome using RNA-seq technology[END_REF][START_REF] Sam | A comparison of single molecule and amplification based sequencing of cancer transcriptomes[END_REF]. However, to our knowledge, although filters for read counts are routinely used in practice, little attention has been paid to the choice of the type of filter or threshold used or its impact on the downstream analysis.

In this paper, we propose a novel data-based procedure to choose an appropriate filtering threshold based on the calculation of a similarity index among biological replicates for read counts arising from replicated high-throughput transcriptome sequencing data. This technique provides an intuitive data-driven way to filter RNA-seq data and to effectively remove genes with low, constant expression levels. Our proposed filtering threshold may be useful in a variety of applications for RNA-seq data, including differential expression analyses, clustering and co-expression analyses, and network inference.

Methods

Types of filters used for RNA-seq data

Data filters are routinely used in practice for differential analyses of RNA-seq data. Most such filters are applied to data that have been normalized in some way, rather than directly to the raw counts, in order to account for systematic inter-sample biases typical of RNA-seq data, e.g., differences in library size [START_REF] Robinson | A scaling normalization method for differential expression analysis of RNA-seq data[END_REF][START_REF] Anders | Differential expression analysis for sequence count data[END_REF] or GC content [START_REF] Risso | GC-content normalization for RNA-seq data[END_REF][START_REF] Hansen | Removing technical variability in RNAseq data using conditional quantile normalization[END_REF]. In particular, the Trimmed M-Means (TMM) library size normalization [START_REF] Robinson | A scaling normalization method for differential expression analysis of RNA-seq data[END_REF] and the normalization included in the DESeq Bioconductor package [START_REF] Anders | Differential expression analysis for sequence count data[END_REF] have been found to be robust methods to correct for library size biases, even in the presence of widely different library compositions [START_REF] Dillies | A comprehensive evaluation of normalization methods for Illumina high-throughput RNA sequencing data analysis[END_REF].

We consider two broad categories of filters for RNA-seq data, based on the filtering criterion used : mean-based filters and maximum-based filters. We note that although variance-based filters are routinely used for microarray data [START_REF] Bourgon | Independent filtering increases detection power for high-throughput experiments[END_REF], they have not been applied to RNA-seq data ; this is likely due to the small number of replicates available in most RNA-seq datasets (and thus, the difficulty in obtaining accurate estimates of per-gene variances) and the fact that the variance is assumed to be a function of the mean under a negative binomial model.

Methods

Mean-based filters

In mean-based filters, genes with mean normalized counts across all samples less than or equal to a pre-specified cutoff are filtered from the analysis. Some authors [START_REF] Sultan | of the Human Transcriptome[END_REF] have also proposed filtering genes with a total read count less than or equal to a given threshold s ; we note that this is equivalent to mean-based filters for threshold s divided by the number of samples.

In addition to normalized counts, we also consider mean-based filters for the RPKM [START_REF] Mortazavi | Mapping and quantifying mammalian transcriptomes by RNA-Seq[END_REF] measure, which was initially proposed to simultaneously normalize RNA-seq data for biases due to library size and gene length. However, we note that it has been shown that counts, rather than RPKM values, are preferable for the differential analysis of RNA-seq data [START_REF] Oshlack | Transcript length bias in RNA-seq data confounds systems biology[END_REF]. For this reason, after filtering genes with a RPKM mean filter, raw counts are used for the subsequent differential analysis. A comparison of differential analysis methods developed for counts and RPKM values is beyond the scope of this work.

Maximum-based filters

In maximum-based filters, genes with maximum normalized counts across all samples less than or equal to a pre-specified threshold are filtered from the analysis. As above, in addition to normalized counts we also consider maximum-based filters for RPKM values, which we refer to as a RPKM maximum filter.

A generalization of the maximum-based filter has also been proposed in the edgeR analysis pipeline (Robinson et al., 2010) based on counts per million (CPM), calculated as the raw counts divided by the library sizes and multiplied by one million. Genes with a CPM value less than a given cutoff (e.g., 1 or 100) in more samples (ignoring condition labels) than the size of the smallest group are subsequently filtered from the analysis. To distinguish this approach from the other maximum-based filters, we refer to this strategy as a CPM filter.

We note that maximum-based filters are not independent filters as described by [START_REF] Bourgon | Independent filtering increases detection power for high-throughput experiments[END_REF] ; in particular, for extremely large filtering thresholds, maximumbased filters do not guarantee control of the Type I error rate if p-values are computed using the pre-filter null distribution. For the threshold values typically used in practice (e.g., based on a quantile, or the data-based threshold proposed below), this is usually not a concern (see Supplementary 1 Figure 28). Although it may be difficult to verify that conditional and unconditional p-value distributions coincide for real data, it may be useful to examine histograms of each (for example, as shown in Figure 4.3).

A data-based threshold for maximum-based filters

For each of the filter types previously defined, a biologically pertinent cutoff (or alternatively, number of genes to be filtered) must be chosen ; in practice, arbitrary thresholds are routinely used with little or no discussion of their impact on the downstream analysis. To address this issue, we propose a data-based choice for the threshold 1. http ://bioinformatics.oxfordjournals.org/content/suppl/2013/06/24/btt350.DC1/Full S uppM at.pdf .1 -Definition of the constants used to calculate the Jaccard similarity index for a pair of samples j and j ′ and a given threshold s. The constant a represents the number of genes with normalized counts greater than s in both samples j and j ′ , and so on.

to be used in maximum-based filters. The main idea underlying this choice is to identify the threshold that maximizes the filtering similarity among replicates, that is, one where most genes tend to either have normalized counts less than or equal to the cutoff in all samples (i.e., filtered genes) or greater than the cutoff in all samples (i.e., non-filtered genes).

To define this filtering similarity, we begin with some notation. Let y gj represent the observed normalized read count (e.g., after scaling raw counts by library size) for gene g in sample j and let C(j) represent the experimental condition of sample j, with g ∈ {1, . . . , G} and j ∈ {1, . . . , J}. Typically, in the context of differential analyses, the number of conditions is equal to 2. We denote the full vector of read counts in a given sample as y j . We now wish to define a similarity index between a pair of replicates within the same condition {(y j , y j ′ ) : C(j) = C(j ′ )} after binarizing the data for a fixed cutoff s (1 if y gj > s and 0 otherwise). We note that a variety of similarity indices have been proposed since the early 1900s ; however, in a comparison among a set of similarity indices (see the Supplementary Materials2 ) we found the Jaccard index [START_REF] Jaccard | Étude comparative de la distribution florale dans une portion des Alpes et des Jura[END_REF] to be simple, natural, and easy to interpret for the analysis of high-throughput sequencing data. This index is defined as follows :

J s (y j , y j ′ ) = a a + b + c (4.1)
where a, b, and c are defined in Table 4.1. We note that J s (y j , y j ′ ) takes on values from 0 (dissimilar) to 1 (similar). Because multiple replicates and/or conditions are typically available in HTS experiments, we extend the definition of the pairwise Jaccard index in Equation ( 4.1) to a global Jaccard index by averaging the indices calculated over all pairs in each condition :

J ⋆ s (y) = mean {J s (y j , y j ′ ) : j < j ′ and C(j) = C(j ′ )} . (4.2)
Using the global Jaccard index defined in Equation ( 4.2) as a measure of similarity, we now wish to identify the cutoff s ⋆ for normalized counts that corresponds to the Methods greatest similarity possible among replicates, that is, the value of s corresponding to the maximum value of the global Jaccard index :

s ⋆ = argmax s J ⋆ s (y). (4.3)
In practice, for the calculation of the data-based global filtering threshold in Equation (4.3), we calculate the value of the global Jaccard index in Equation ( 4.2) for a fixed set of threshold values and fit a loess curve [START_REF] Cleveland | Robust locally weighted regression and smoothing scatterplots[END_REF] through the set of points ; the value of s ⋆ is subsequently set to be the maximum of these fitted values.

Once the data-driven filter threshold for normalized counts s ⋆ has been identified, the subsequent steps to be taken may change for different applications. To perform an analysis of differential expression between two experimental conditions, we propose using this threshold s ⋆ in a maximum-based filter, as defined in Section 4.2.1 ; in the following, we refer to this technique as the Jaccard filter.

The Bioconductor package HTSFilter

The proposed filtering method is implemented in the HTSFilter package, currently available as part of the Bioconductor project [START_REF] Gentleman | Bioconductor : Open software development for comptuational biology and bioinformatics[END_REF] within the statistical environment R (R Development Core Team, 2009). The HTSFilter package is compatible with a variety of data classes and analysis pipelines, including matrix and data.frame objects, the S4 class CountDataSet in the DESeq pipeline [START_REF] Anders | Differential expression analysis for sequence count data[END_REF], and the S3 class DGEList in the edgeR pipeline (Robinson et al., 2010). A package vignette describes the use of the HTSFilter package within each of these pipelines.

Results

In the following, we apply the normalization approach proposed by [START_REF] Anders | Differential expression analysis for sequence count data[END_REF] for mean-and maximum-based filters, although other types of normalization may be appropriate for some data. For gene-by-gene comparisons between two conditions, we illustrate the use of the proposed filter in conjunction with the model proposed in the DESeq Bioconductor package (version 1.8.3), which has been developed in order to model count data with a small number of replicates in the presence of overdispersion [START_REF] Anders | Differential expression analysis for sequence count data[END_REF] ; p-values are adjusted for multiple testing using the Benjamini-Hochberg procedure [START_REF] Benjamini | Controlling the false discovery rate : a practical and powerful approach to multiple testing[END_REF] to control the false discovery rate (FDR). We note that the filtering method proposed here may also be used in conjunction with other popular methods, e.g., edgeR (Robinson et al., 2010). See the Supplementary Materials for additional discussion of the normalization and statistical testing methods used in this work.

In practice, there may be some question about the appropriate point in the analysis pipeline to apply data filters : Should normalized data first be filtered, then normalization factors re-estimated and the model fit (i.e., mean and dispersion parameters estimated) ? Should normalization factors and model parameters be estimated based on the full data, and the data filtered only at the end of the analysis pipeline ? The difference between the two options is nontrivial, particularly as the differential analysis approaches implemented in the DESeq and edgeR packages both borrow information across genes (whether all or only those passing the filter) to obtain per-gene parameter estimates. In this work, we present results based on the application of filters applied as late in the pipeline as possible, i.e., after library size and dispersion parameter estimation ; a more detailed discussion of this issue is included in the Supplementary Materials.

Description of data

We applied our proposed Jaccard index filter, in addition to the alternative filter types described above, on the following data :

-Sex-specific expression of liver cells in human. [START_REF] Sultan | of the Human Transcriptome[END_REF] obtained high-throughput transcriptome sequencing data from a human embryonic kidney and a B cell line, with two biological replicates each. The raw read counts and phenotype tables were obtained from the ReCount online resource [START_REF] Frazee | ReCount : a multi-experiment resource of analysis-ready RNA-seq gene count datasets[END_REF]. -Differential striatal expression between inbred mouse strains. [START_REF] Bottomly | Evaluating gene expression in C57BL/GJ and DBA/2J mouse striatum using RNA-seq and microarrays[END_REF] performed RNA-seq experiments for ten biological replicates of the C57BL/6J inbred mouse strain and eleven for the DBA/2J strain, and the results were compared with those arising from two different microarray platforms. The raw read counts and phenotype tables were obtained from the ReCount online resource [START_REF] Frazee | ReCount : a multi-experiment resource of analysis-ready RNA-seq gene count datasets[END_REF]). -MiTF repression in a human melanoma cell line. [START_REF] Strub | Essential role of microphthalmia transcription factor for DNA replication, mitosis and genomic stability in melanoma[END_REF] obtained high-throughput sequencing data to compare gene expression in a melanoma cell line expressing the Microphtalmia Transcription Factor (MiTF) to one in which small interfering RNAs (siRNAs) were used to repress MiTF, with three biological replicates in each group. The raw read counts and phenotype tables are available in the Supplementary Materials of [START_REF] Dillies | A comprehensive evaluation of normalization methods for Illumina high-throughput RNA sequencing data analysis[END_REF]. -Simulated data. To investigate the effect of the various filtering methods on downstream results, we developed a simulation framework as described in Section 4.3.3. For the three real datasets (described in further detail in Supplementary Table 1), gene annotations for Mus musculus (NCBIM37) and Homo sapiens (GRCh37.p7) were obtained from Ensembl version 67 [START_REF] Birney | An Overview of Ensembl[END_REF] using the Biomart tool [START_REF] Kasprzyk | EnsMart : A generic system for fast and flexible access to biological data[END_REF], and the length of each gene in base pairs was calculated. In the Bottomly, Sultan, and Strub data, about 4%, 2%, and 5% (respectively) of the genes had been retired from Ensembl ; for these genes, RPKM-based filters were not used. Readers wishing to examine the complete analyses in detail may find a Sweave document containing commented R code for the analyses of each of the datasets in the Supplementary Materials.

Comparison of filters on real data

For both real datasets, differential analyses were performed using a negative binomial model [START_REF] Anders | Differential expression analysis for sequence count data[END_REF] for unfiltered data and after filtering the data using the techniques described above. For the Jaccard index filter, the global Jaccard index in Equation ( 4.2) was calculated for a range of threshold values s, yielding a largely unimodal distribution of values for all datasets considered (Figure 4.1 and Supplementary Figures 11 and18). We also note that the data-driven threshold values s ⋆ identified for the Bottomly, Sultan, and Strub data were not equal ; in the case of the Bottomly data, the threshold for normalized counts was found to be 14.7, while this threshold was found to be 11.5 for the Sultan data and 103.5 for the Strub data. These differences in filtering threshold among experiments are due to both sequencing depth and variability within the data ; in particular, experiments with greater sequencing depth will tend to have higher filtering thresholds, and those with greater variability will tend to have lower filtering thresholds. Among the data considered in our study, the Strub data have the highest sequencing depth (1.5 × 10 8 ) coupled with low intra-condition variability (minimum correlation among replicates equal to 0.98), and they also have the highest threshold considered here. On the other hand, the Sultan data have a much lower sequencing depth (1.8 × 10 6 ), and thus have a much lower threshold.

For each dataset, the Jaccard filter was applied with the corresponding data-based threshold calculated above. For the alternative mean-and maximum-based filters for normalized counts and RPKM values, cutoffs were chosen based on the 15% quantile of the respective criterion. For the CPM filter, as suggested in the edgeR pipeline Robinson et al. (2010), genes with a CPM value less than 1 in more samples than the size of the smallest group are subsequently filtered from the analysis.

Among the filters considered, it may immediately be seen in Figure 4.2 that in the Strub data, with the exception of the CPM and Jaccard filters, most of the filters considered here appear to be ineffective as they are largely unable to filter genes with very low levels of expression and small log-fold changes ; a similar phenomenon may be observed in the Bottomly and Sultan data (Supplementary Figures 7 and14). Although these techniques are thus unlikely to (incorrectly) filter truly DE genes, the number of statistical tests is not markedly reduced, and as such the power to detect differential expression will remain largely unchanged as compared to the unfiltered data. With the exception of the RPKM-based filters, all are able, to some extent, to identify and remove genes contributing to a peak of raw p-values close to one (Figure 4.3 and Supplementary Figures 13 and20), a phenomenon due to the discretization of p-values for small counts ; indeed, histograms of raw p-values following the application of most filters appear to be roughly uniformly distributed under the null hypothesis. In addition, we note that the proposed Jaccard filter is able to more effectively remove genes with moderate log base means and small log fold changes (see Figure 4.2). Similar conclusions may be drawn from the volcano plots shown in Supplementary Figures 8, 15, and 22.

It is also of interest to consider the effect of each filter on the number of DE genes identified at various levels of expression ; in Figure 4.4 and Supplementary Figures 5 and12, we note that in all data sets the Jaccard filter leads to more discoveries at all but very weak levels of expression (i.e., mean expression less than 10), with this difference being particularly marked for moderate levels of expression (i.e., mean expression greater than 50). We note that a large number of the missed discoveries for the Jaccard filter at very low levels of expression correspond to genes with zero read counts in one condition and a small number of read counts in the other ; for example, in the Sultan data 50.3% of the 449 discoveries among genes with mean normalized read counts less than 10 had zero read counts in one of the two conditions. Among the other filter types, the CPM filter appears to come closest to the Jaccard filter, with the remaining filters performing similarly.

In considering the overlap of genes filtered using each method (Supplementary Figures 9, 16, and 23), it is clear that a large number of genes may be filtered regardless of the technique used. However, the Jaccard index is better able to filter a large number of weakly expressed genes in all three of the datasets considered here, leading to a more moderate correction for multiple testing ; the direct consequence of this is a larger number of discoveries at moderate to high levels of expression. In order to determine whether this advantage is due to the filtering type (i.e., maximum) or the threshold used for each (i.e., using the 15% quantile or the data-based threshold identified by the global Jaccard index), we consider a set of simulation studies in the next section.

Simulated data

Data were simulated using a negative binomial model, with parameters chosen based on the Bottomly, Sultan, and Strub datasets. Briefly, for each dataset, genes with zero counts in one of the two groups and mean less than 5 were removed (see Supplementary Table 1). Differential analyses were performed on unfiltered data using the DESeq Bioconductor package [START_REF] Anders | Differential expression analysis for sequence count data[END_REF] as described in the Supplementary Materials. After adjusting raw p-values for multiple testing using the Benjamini-Hochberg procedure [START_REF] Benjamini | Controlling the false discovery rate : a practical and powerful approach to multiple testing[END_REF], we identified 570, 2515, and 2485 differentially expressed genes, respectively, at the 5% significance level in the Bottomly, Sultan, and Strub data. In addition, the parametric gamma regressions fitted to the per-gene dispersion estimates α and gene expression means µ (Supplementary Figure 25) were identified for each dataset : Subsequently, simulation parameters for each dataset were fixed as follows. For genes identified as being differentially expressed, means for each condition were set to be the empirically calculated means from each condition from the normalized data ; for genes not identified as being differentially expressed, means for each condition were both set to be the global mean (across both conditions) from the normalized data. Note that this allows genes to be simulated as differentially expressed across the full range of mean expression values (Supplementary Figure 27). Per-gene dispersion parameters were set to be the fitted values from the regression equations defined in Equations ( 4 .5 -Sensitivity (over simulated 300 datasets) to detect differentially expressed genes for a variety of filter types and cutoffs, with simulation parameters based on the Strub data. None : no filter. CPM : genes with a CPM less than one in more than half the samples are filtered. Max : maximum-based filter, using the threshold based on the Jaccard index (J), quantiles (5%, 15%, 30%), or values (5,10,25). Mean : mean-based filter, using the threshold based on quantiles (5%, 15%, 30%), or values (5,10,25). RPKM.max : maximum RPKM filter, using the threshold based on quantiles (5%, 15%) or the value 0.3. RPKM.mean : maximum RPKM filter, using the threshold based on quantiles (5%, 15%) or the value 0.3.

Bottomly and Sultan data, dispersion parameters for genes with overall mean expression less than 20 were fixed to be equal to 10 -10 to simulate negligible overdispersion, as shot noise appears to dominate biological noise at low expression levels in these data (Supplementary Figure 26). Once these parameters were fixed for each gene, a negative binomial model was used to simulate 300 individual datasets each for the parameters based on the Bottomly, Sultan and Strub data, with 21 samples (10 in one condition and 11 in the other), 4 samples (2 in each condition), and 6 samples (3 in each condition), respectively. Real lengths corresponding to the genes in each dataset were used for the calculation of RPKM values.

Comparison of filters on simulated data

In order to assess performance on simulated data, we focus on the sensitivity of detecting differentially expressed genes after each data filter, defined as the proportion of truly DE genes detected among all truly DE genes. In addition, we construct Receiver Operating Characteristic (ROC) curves of each filter, based on the filtering sensitivity, defined as the proportion of correctly unfiltered genes (i.e., DE and unfiltered) among all truly DE genes, and the filtering specificity, defined as the proportion of correctly filtered genes (i.e., non-DE and filtered) among all non-DE genes.

In Figure 4.5, we note that the sensitivity to detect differentially expressed genes greatly varies among the filtering types for simulations based on the Strub data, as well as among different thresholds within each filtering type ; in addition, the RPKM maximum and RPKM mean filters actually lead to lower detection power than unfiltered data. Similar results may be seen for the simulations based on the Bottomly and Sultan data in Supplementary Figure 30. For the simulation setting shown in Figure 4.5, larger thresholds appear to yield the highest detection sensitivity (i.e., maximum or meanbased filters for normalized counts using the 30% quantile as a threshold). However, this trend is reversed for the other simulation settings, where smaller cutoffs (i.e., cutoffs of 5% or 15% for the maximum or mean-based filters) lead to higher detection sensitivity. This highlights the difficulty in pre-selecting a fixed threshold for a given filtering method, as well as the advantage of our proposed Jaccard filter. Indeed, the Jaccard filter appears to lead to high detection sensitivity for all simulations with the exception of those based on the Bottomly data ; for these data, because many weakly expressed genes were simulated to be differentially expressed (see Supplementary Figure 27), unfiltered data had the highest sensitivity.

To assess the role of the choice of threshold for each filter type, we constructed ROC curves of the filtering sensitivity and specificity over varying cutoffs in Figure 4.6 and Supplementary Figure 29. We note that the mean and maximum RPKM-based filters tend to have lower filtering sensitivity than the others, and the maximum filters (for both normalized counts and CPM values) tend to be quite similar ; however, for all simulation settings, the maximum filter for normalized counts has a slight advantage over that for CPM values. In other words, regardless of the threshold used for filtering, maximum-based filters appear to more effectively filter non-DE genes than the remaining methods. Finally, we note also that the data-based threshold using the global Jaccard index appears to find a good compromise between filtering sensitivity and filtering specificity.

Conclusions and discussion

Data filtering has proven to be of great practical importance for the differential analysis of high-throughput microarray and RNA-seq data by identifying and removing genes with uninformative signal prior to testing. In recent years, many ad hoc procedures have been used to filter RNA-seq data, such as filtering genes with a total or mean normalized read count less than a specified threshold. However, despite its impact on the downstream analyses, no clear recommendations have yet been provided concerning the choice of filtering technique.

Among the filter types considered here, we have found that filters using the maximum normalized count appear to be best able to correctly filter genes with low levels of expression and little evidence of differential expression. In addition, we have proposed a method to calculate a data-driven and non pre-fixed filtering threshold value for normalized counts from replicated RNA-seq data, based on the global Jaccard similarity index. In particular, our proposed filtering technique was found to remove from the analysis a large number of genes with little or no chance of showing evidence of differential expression, and therefore to increase detection power at moderate to high levels of expression through a moderation of the correction for multiple testing. As such, we recommend that genes with a normalized count value less than this data-driven threshold in all samples be filtered from subsequent differential analyses. We emphasize that the datadriven threshold value may vary greatly among RNA-seq experiments due to differences in sequencing depth and intra-condition variability (see Supplementary Figure 31 for the data-based thresholds calculated on three additional RNA-seq datasets) ; as such, the threshold value must be recalculated for each data set of interest.

The impact of the proposed filtering method has been investigated here in the context of differential analyses. We anticipate that it will also be useful in a variety of other applications, for example detecting genes that are specifically expressed in one condition or ubiquitously expressed across several conditions, which is often a crucial biological question. In addition, we anticipate that such filtering will be useful, for example, in co-expression or network reconstruction analyses to remove genes with low, constant levels of expression. Finally, we note that although this filter was presented here for the analysis of RNA-seq data, it can readily be applied to other types of replicated high-throughput sequencing data, such as CHIP-seq data. Abstract. High-throughput transcriptome sequencing data (RNA-seq) are made up of highly heterogeneous counts. Although they are often modeled with discrete distributions, including the Poisson and negative binomial distributions, Gaussian models on transformed data could alternatively be considered. We show how the likelihood of these different models can be objectively compared. We focus attention on the problem of clustering gene profiles, where Poisson mixtures on count data are compared with Gaussian mixtures on transformed data. On dispose d'une matrice de mesures d'expression de n gènes y = (y 1 , . . . , y n ). Pour chaque gène i (i = 1, . . . , n), le vecteur y i indique l'expression du gène i pour les p conditions expérimentales j (j = 1, . . . , p). On suppose que les données y sont la réalisation d'un mélange de K variables aléatoires de lois de Poisson de densité :

f (y i ; K, θ K ) = K k=1 p k p j=1 P(y ij ; µ ijk ).
(5.1)

Suivant [START_REF] Rau | Coexpression analysis of high-throughput transcriptome sequencing data with Poisson mixture models[END_REF], les paramètres (p 1 , . . . , p K ) sont les proportions de chaque composante du mélange et p j=1 P(y ij ; µ ijk ) est la densité d'un vecteur de p variables aléatoires indépendantes de lois de Poisson de moyennes respectives µ ijk = w i s j λ jk pour k = 1, . . . , K. On note

y i• = p j=1 y ij , y •j = n i=1 y ij et y •• = n i=1 p j ′ =1 y ij ′ .
Les facteurs w i = y i• prennent en compte le niveau d'expression de chaque gène, ils sont calculés avant tout estimation des paramètres du modèle. Les facteurs s j corrigent le biais technique spécifique aux données RNA-seq lié aux différences de profondeur de séquençage entre échantillons (voir Chapitre 1). L'expression s j = y •j y•• correspond à un estimateur du maximum de vraisemblance. Les paramètres λ k = (λ 1k , . . . , λ pk ) correspondent aux profils d'expression des gènes de la composante k. Ainsi, le modèle classe les gènes en fonction de leur dynamique d'expression (λ k ) et non en fonction de leur niveau d'expression absolu (w i ). Les paramètres p k et λ jk sont estimés par l'algorithme EM sous les contraintes K k=1 p k = 1 et p j=1 λ jk s j = 1 pour tout k. L'implémentation de l'estimation des paramètres de ce modèle est proposée dans le package HTSCluster.

Transformation des données RNA-seq et modèle de mélange de lois gaussiennes

Une alternative à ce modèle de mélange de lois de Poisson est un modèle de mélange de lois gaussiennes, classiquement utilisé pour les données de puces à ADN. Dans un cadre différent, celui de l'analyse différentielle d'expression de gènes, [START_REF] Law | Voom : precision weights unlock linear model analysis tools for RNA-seq read counts[END_REF] ont proposé une transformation logarithmique des données RNA-seq afin d'utiliser les modèles linéaires gaussiens développés initialement pour l'analyse des données de puces. Dans le même esprit, on propose ici le même type de transformation pour l'utilisation d'un modèle de mélange de lois gaussiennes.

Les données y sont transformées de sorte que l'objectif de classification reste le plus proche de celui du modèle de mélange de Poisson précédent (modélisation de la dy-namique d'expression entre conditions). Chaque comptage y ij est divisé par le facteur N j = y •j 10 6 afin de corriger le biais technique spécifique aux données RNA-seq lié aux différences de profondeur de séquençage entre échantillons (voir Chapitre 1). Le facteur N j est le nombre comptages pour condition j, exprimé en million. Il correspond globalement au facteur s j du modèle de mélange de lois de Poisson. Afin de modéliser la variation d'expression du gène, on compare le comptage normalisé y ij /N j à m i = 1 p p j ′ =1 y ij ′ N j ′ , l'expression moyenne du gène i à travers les p conditions expérimentales. Le facteur m i correspond globalement au facteur w i dans le modèle de mélange de lois de Poisson. Une constante est ajoutée au numérateur et au dénominateur pour éviter les problèmes numériques. On nomme cette transformation des données t :

t(y ij ) = log y ij /N j + 1 m i + 1 .
On modélise le vecteur des données transformées ỹi = t(y i ) par un mélange de K lois gaussiennes de densité :

g(ỹ i ; K, η k ) = K k=1 p k φ(ỹ i ; ν k , Σ k ).
(5.2)

Les paramètres (p 1 , . . . , p K ) sont les proportions de chaque composante du mélange et φ(ỹ i ; ν k , Σ k ) est la densité d'une loi normale de dimension p de moyenne ν k et de variance-covariance Σ k . Une implémentation de l'estimation des paramètres de ce modèle est proposée dans le package Rmixmod [START_REF] Lebret | Rmixmod : The R Package of the Model-Based Unsupervised, Supervised and Semi-Supervised Classification Mixmod Library[END_REF].

Transformation des données et comparaison de modèles

La vraisemblance du modèle de mélange de lois de Poisson sur les données brutes par rapport à la mesure de comptage s'écrit :

ℓ f (y 1 , . . . , y n ; K, θ K ) = n i=1 f (y i ; K, θ K ).
La vraisemblance du modèle de mélange gaussien sur les données transformées par rapport à la mesure de Lebesgue s'écrit :

ℓ g (ỹ 1 , . . . , ỹn ; K, η K ) = n i=1 g(ỹ i ; K, η k ).
La comparaison de ces deux modèles est a priori difficile car les vraisemblances des modèles s'écrivent par rapport à des mesures différentes. Cependant, la densité du modèle de Poisson peut être définie comme la limite d'une densité par rapport à la mesure de Lebesgue. Pour tout x ∈ R et 0 < ∆ < 1 2 , on considère :

f ∆ (x) =      exp(-λ) λ t t! 1 ∆ , si x ∈]t -∆ 2 ; t + ∆ 2 [, t ∈ N 0, sinon.
Illustration sur des données simulées

Cette densité est bien définie par rapport à la mesure de Lebesgue et sa limite, lorsque ∆ tend vers 0, est la densité d'une loi de Poisson. Dans ce contexte, on peut écrire l'égalité des densités dỹ i = t ′ (y i )dy i à partir de la relation ỹi = t(y i ). On obtient la relation suivante :

g(ỹ i ; K, η k )dỹ i = g(t(y i ); K, η k )t ′ (y i )dy i .
Cette relation permet de réécrire la vraisemblance du modèle de mélange sur données transformées en fonction des données initiales y = (y 1 , . . . , y n ) :

ℓ g (y 1 , . . . , y n ; K, η K ) = n i=1 g(t(y i ); K, η k )t ′ (y i ).
Les deux modèles peuvent alors être comparés par un critère de vraisemblance pénalisée comme le BIC :

BIC f (K; y) = n i=1 log f (y i ; K, θK ) - ν f 2 log(n), BIC g (K; y) = n i=1 log g(t(y i ); K, ηk ) + n i=1 log t ′ (y i ) - ν g 2 log(n).
Les quantités θK et ηK sont les estimateurs du maximum de vraisemblance des paramètres des modèles respectifs, ν f et ν g sont les nombres de paramètres des modèles respectifs. Le modèle s'ajustant le mieux aux données est le modèle maximisant le critère BIC associé. Cette prise en compte de la transformation appliquée aux données dans le calcul du BIC a été utilisée auparavant dans un autre domaine [START_REF] Thomas | The Fractal Morphology of the Built-Up Landscape[END_REF].

Illustration sur des données simulées

Paramètres de simulation

Afin d'illustrer la comparaison de modèles proposée, nous simulons des données sous le modèle de mélange de lois de Poisson détaillé à l'équation (5.1), en fixant le nombre de conditions expérimentales p à 3, le nombre de gènes n à 5000, les facteurs de normalisation (s 1 , s 2 , s 3 ) = ( 13 , 1 3 , 1 3 ) et les facteurs w i à partir du jeu de données réelles de Mach et al. (2014) (jeu de données décrit dans la section suivante) en sélectionnant aléatoirement trois conditions expérimentales et n gènes parmi les gènes du jeu de données ayant au moins 20 comptages par gènes. Les valeurs des w i varient ainsi de 20 à 1 800 000 comptages. Nous fixons ensuite le nombre de classes K = 4, les proportions de chaque classe (p 1 , p 2 , p 3 , p 4 ) = ( 14 , 1 4 , 1 4 , 1 4 ) et les paramètres λ jk pour j = 1, 2, 3 et k = 1, 2, 3, 4 tels que j λ jk s j = 1 : λ =   1.5 1 0.5 1.5 0.5 1.5 1 1 1 0.5 1.5 0.5   .

Résultats

La figure 5.1 (gauche) illustre les comptages simulées transformés pour la condition 1 versus la condition 2. On constate que l'hypothèse de modélisation des données transformées par un modèle de mélange de lois gaussiennes est cohérente : les différentes composantes simulées peuvent être représentées par des ellipses.

Conformément au résultat attendu, la figure 5.1 (droite) montre que le BIC du modèle de mélange gaussien, ajusté pour la t-transformation des données (5.2) est inférieur au BIC du modèle de mélange de lois de Poisson pour un nombre de classes supérieur ou égal à 4. On remarque également que le modèle de mélange de lois gaussiennes surestime fortement le nombre de composantes du mélange et sélectionne 14 classes au lieu des 4 classes simulées.

• • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • •• • • • • • • • • • • •• • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • -1.

Illustration sur des données réelles

Pour plusieurs jeux de données RNA-seq, nous effectuons la classification des gènes à l'aide du modèle de mélange de Poisson sur les données de comptage brutes, et à l'aide du modèle de mélange gaussien sur les données transformées.

Illustration sur des données réelles

Description des données

Cellules hépatiques chez l'humain : [START_REF] Sultan | of the Human Transcriptome[END_REF] ont analysé l'expression des gènes dans les cellules humaines embryonnaires du rein (HEK293T) et dans les cellules de la lignée Ramos B en effectuant le séquençage de deux réplicats biologiques dans chaque type de cellule par la technologie RNA-seq. Après avoir supprimé les gènes peu exprimés, nous effectuons la classification des 4959 gènes restants.

Evolution embryonnaire chez la drosophile : Dans le cadre du projet modEN-CODE visant à completer l'annotation fonctionnel du génome de la mouche Drosophila melanogaster, [START_REF] Graveley | The developmental transcriptome of Drosophila melanogaster[END_REF] ont analysé les différences d'expression géniques à différents stades du développement de la mouche. On s'intéresse à 12 cellules d'embryons collectées toutes les deux heures pendant 24 heures. Ces 12 réplicats biologiques sont séquencés à l'aide de la technologie RNA-seq. Le tableau de comptage des données a été extrait de la base de données ReCount [START_REF] Frazee | ReCount : a multi-experiment resource of analysis-ready RNA-seq gene count datasets[END_REF]. Les 13164 gènes exprimés sont classés. Tissus de l'intestin grêle chez le porc : [START_REF] Mach | Extensive expression differences along porcine small intestine evidenced by transcriptome sequencing[END_REF] ont analysé les différences d'expression entre trois tissus (le duodenum, le jejunum et l'ileum) de l'intestin grêle de quatre porcelets sains. Après avoir sélectionné les gènes différentiellement exprimés entre ces trois tissus à l'aide d'un modèle linéaire généralisé basé sur une loi négative binomiale [START_REF] Robinson | A scaling normalization method for differential expression analysis of RNA-seq data[END_REF], on effectue la classification des 4021 gènes restants.

Résultats

Sur les Figures 5.2, nous constatons que deux des trois jeux de données de Graveley et Mach sont mieux ajustés par le modèle de mélange de lois gaussiennes sur les données transformées que par le modèle de Poisson. Le mélange de Poisson s'ajuste mieux au jeu de données de Sultan. Cependant, la différence entre les valeurs des deux BIC n'est pas comparable entre les différents jeux de données : sur les données de Sultan, les valeurs des BIC pour le mélange de Poisson et le mélange gaussien sont relativement proches alors que le BIC du mélange gaussien est très supérieur au BIC du mélange de Poisson sur les données de Graveley et Mach. On constate également que le choix de modèle réalisé par le critère BIC est plus parcimonieux pour le modèle ajustant le mieux les données : le BIC du mélange de Poisson sélectionne 27 classes sur les données Sultan alors que le BIC du mélange gaussien sélectionne plus de 30 classes. Sur les données Graveley et Mach, le critère BIC du mélange de Poisson ne parvient à sélectionner un modèle (ce BIC augmente toujours, même pour des très grandes valeurs de K), alors que le BIC du mélange gaussien sélectionne moins de 50 classes pour les deux jeux de données.

La technologie RNA-seq évolue et les jeux de données générés au moment de l'émergence de la technologie ne ressemblent pas forcement aux jeux de données plus récents. Les plans d'expérience sont de plus en plus complexes et le nombre de réplicats biologiques à disposition augmente. Le modèle de mélange de lois gaussiennes semble mieux

Conclusion

prendre en compte ces plans d'expérience complexes. Cependant, une généralisation de ce résultat nécessiterait de comparer plus de trois jeux de données.

Conclusion

À l'aide du critère BIC, on peut déterminer si une transformation des données fournit un meilleur ajustement des données et comparer différentes stratégies de modélisation.

Cette comparaison est pertinente pour la classification des données RNA-seq par modèle de mélange, pour lesquels deux alternatives existent : modélisation des données à l'aide de modèle de mélange discret ou à l'aide de modèle de mélange de lois gaussiennes. Nous avons constaté qu'il n'existe pas de modèle le mieux adapté aux données RNA-seq. La comparaison de critère BIC permet ainsi de sélectionner le modèle le plus adapté aux données. Sur les trois jeux de données étudiés, la transformation logarithmique des données que nous avons proposée fournit souvent des modèles beaucoup plus convaincants que le modèle de mélange de Poisson.

Dans ce chapitre, nous avons utilisé le critère BIC pour comparer les différentes stratégies. On peut cependant utiliser d'autres critères basés sur la log-vraisemblance pénalisée des données, comme le critère ICL ou le critère SICL .

Chapitre 6 A model selection criterion for model-based clustering of annotated gene expression data

Note. La version intégrale de ce chapitre est disponible en version publiée [START_REF] Gallopin | A model selection criterion for model-based clustering of annotated gene expression data[END_REF]. La pré-publication est disponible à l'adresse https://hal.inria.fr/ hal-01088870.

Résumé. En analyse de co-expression des données de gènes, on souhaite souvent interpréter les modules détectés à l'aide d'informations externes, telles que la liste potentiellement incomplète des propriétés fonctionnelles attribuées à un ensemble de gènes. Dans le cadre des modèles de mélange, nous proposons un critère de sélection de modèle prenant en compte ces informations externes fournissant ainsi un outil pertinent pour sélectionner un modèle et un nombre de classes. Ce critère Integrated Completed Annotated Likelihood (ICAL) est défini en ajoutant un terme d'entropie à la vraisemblance pénalisée du modèle afin de mesurer la concordance entre la partition inférée et les annotations externes. Le critère ICAL conduit à un choix de modèle plus facilement interprétable au regard des annotations fonctionnelles des gènes déjà répertoriées. On illustre l'intérêt de ce critère de sélection de modèle pour des mélanges gaussiens sur des données simulées et sur un jeu de données réel d'expression de gènes RNA-seq.

Abstract. In co-expression analyses of gene expression data, it is often of interest to interpret clusters of co-expressed genes with respect to a set of external information, such as a potentially incomplete list of functional properties for which a subset of genes may be annotated. Based on the framework of finite mixture models, we propose a model selection criterion that takes into account such external gene annotations, providing an efficient tool for selecting a relevant number of clusters and clustering model. This criterion, called the Integrated Completed Annotated Likelihood (ICAL), is defined by adding an entropy term to a penalized likelihood to measure the concordance between a clustering partition and the external annotation information. The ICAL leads to the choice of a model that is more easily interpretable with respect to the known functional gene annotations. We illustrate the interest of this model selection criterion in conjunction with Gaussian mixture models on simulated data and on real gene expression RNA-seq data.

Conclusion et perspectives

Chapitre 7

A hierarchical Poisson log-normal model for network inference from RNA sequencing data Résumé. L'inférence de réseau de gènes à partir des données transcriptomiques est un défi méthodologique important et un aspect clef de la biologie des systèmes. Des méthodes ont été proposées pour l'inférence de réseau à partir des données de puces à ADN, mais ces méthodes ne sont pas directement applicables aux données RNA-seq qui sont des données discrètes et très hétérogènes. Dans ce chapitre, nous proposons un modèle hiérarchique log-normal de Poisson pénalisé afin d'inférer des réseaux de gènes sur les données RNA-seq. Basé sur des lois de Poisson, ce modèle prend en compte le caractère discret des données RNA-seq. Grâce à sa structure hiérarchique, il modèlise également la grande variabilité des données RNA-seq pour lesquelles la variance empirique est supérieure à la moyenne empirique. Nous comparons la méthode proposée à deux alternatives : un modèle graphique gaussien pénalisé appliqué sur les données log-tranformées, et un modèle graphique de Poisson log linéaire appliqué sur les données transformées par une transformation de puissance. Ces méthodes sont comparées à l'aide de simulations et à l'aide d'un jeu de données réels microRNA-seq extrait de cellules cancéreuses du sein. Sur les données simulées avec une large dispersion inter échantillons, le modèle proposé est plus performant que les autres méthodes appliquées sur données transformées en terme de sensibilité, de spécificité et d'aire sous la courbe ROC. Ces résultats montrent la nécessité de proposer des méthodes d'inférence de réseaux spécifiques aux données RNA-seq.

Abstract. Gene network inference from transcriptomic data is an important methodological challenge and a key aspect of systems biology. Although several methods have been proposed to infer networks from microarray data, there is a need for inference methods able to model RNA-seq data, which are count-based and highly variable. In this work we propose a hierarchical Poisson log-normal model with a Lasso penalty to infer gene networks from RNA-seq data ; this model has the advantage of directly modelling discrete data and accounting for inter-sample variance larger than the sample mean. Using real microRNA-seq data from breast cancer tumors and simulations, we compare this method to a regularized Gaussian graphical model on log-transformed data, and a Poisson log-linear graphical model with a Lasso penalty on power-transformed data.

Materials and Methods

We first define the notation that will be used throughout this paper. Let Y ij be the random variable corresponding to the gene expression measure for the sample i (i = 1, . . . , n) for the gene j (j = 1, . . . , p), with y ij being the corresponding observed value of Y ij . Note that i always indexes samples and j always indexes genes with n the number of samples and p the number of genes. A network represents gene interactions. The nodes are random variables modelling the gene expression levels and the edges indicate the dependencies between those variables. In this section we provide a short description of the models that will be compared for gene network inference from RNAseq data.

Gaussian graphical model

The underlying assumption of this model is that the data are normally distributed. In the case of untransformed RNA-seq data, this assumption is not valid since data counts cannot take negative values. We investigated a variety of Box-Cox transformations to lead to approximately normal data [START_REF] Box | An analysis of transformations[END_REF], where the δ value was chosen to maximize the log-likelihood of the transformed data :

y ij → f (y ij ) =    y δ ij -1 δ , if δ = 0, log (y ij ), if δ = 0.
Since gene expression data may contain zero counts, we usually use (y + 1) instead of y in the Box-Cox formula above. Let z i = (f (y i1 ), ..., f (y ip )) be the transformed vector of expression values for p genes for the ith biological sample (i = 1, ..., n). We assume that z i ∼ N (µ, Σ). The edges of the inferred network correspond to non-zero partial correlations, i.e. the non-zero elements of matrix Σ -1 [START_REF] Whittaker | Graphical Models in Applied Multivariate Statistics[END_REF][START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF].

Let S be the empirical covariance matrix. The log-likelihood of the model is :

L(Σ -1 ) = log(det(Σ -1 )) -trace(SΣ -1 ).
(7.1)

A common assumption in the context of gene networks is that the matrix Σ -1 is sparse.

We add an ℓ 1 penalty to the log-likelihood (8.1) so that some coefficients in the estimated Σ -1 matrix are precisely equal to 0 :

log(det(Σ -1 )) -trace(SΣ -1 ) -λ Σ -1 ℓ 1 . (7.2)
Network inference using a Gaussian graphical model has been extensively studied and used over the past years. Many methods exist to compute the penalized maximum likelihood estimate of the Σ matrix above. We use the method implemented in the glasso R package [START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF] which makes use of a coordinate descent algorithm.

The choice of the regularization parameter λ has also been extensively studied [START_REF] Giraud | Graph selection with GGMselect[END_REF]. We choose to perform model selection by maximizing the Bayesian

Hierarchical log-normal Poisson graphical model

We note that the Poisson model presented above requires a transformation of the data to account for the high dispersion. Here we propose to deal with it directly with a hierarchical log-normal Poisson model. The count expression of gene j for sample i ∈ 1, . . . , n is modeled as :

Y ij ∼ P(θ ij ) with log(θ ij ) = j ′ =j β jj ′ ỹij ′ + ε ij , ε j = (ε 1j , ..., ε nj ) ∼ N (0, σ 2 j I n ).
As before, the notation ỹ corresponds to a standardization of the log-transformed data.

Here, the vector Y j ∼ P(θ j ) and θ j is itself a random variable :

θ j = µ j exp(ε j ) with ε j ∼ N n (0, σ 2 j I n ) and µ j = exp( j ′ =j β jj ′ ỹij ′ ).
Note that the variance of the random variable P(θ j ) is larger than its mean if σ 2 j is positive. As previously, an edge is present in the graph between genes j and j ′ if one or both parameters β jj ′ and β j ′ j are different from zero.

In this model, the likelihood for gene j can be written as :

L(β j , σ j ) = R n i=1 exp(-µ ij + y ij log(µ ij ) -log(y ij !) 1 (2π) n/2 σ n j exp - 1 2σ 2 j ||ε j || 2 2 dε j . (7.5)
Similar to the previous model, we assume that the vector β j is sparse. We add an ℓ 1 penalty to a function of the log-likelihood (7.5) so that some coefficients in the estimated β j vector are set to 0 :

-2L(β j , σ j ) + λ β j ℓ 1 .
Estimation of parameters β j and σ j was done using the R function glmmixedlasso [START_REF] Schelldorfer | Glmmlasso : An algorithm for high-dimensional generalized linear mixed models using 1-penalization[END_REF], based on a Laplace approximation of the penalized likelihood and a coordinate descent algorithm.

An important aspect of this method is the choice of the regularization parameter λ. To choose a common λ parameter for all the gene-by-gene regressions, we propose to use a two stage approach for this parameter. First, for each gene j, a λ j parameter is chosen by maximizing the BIC criterion defined as BIC = L(β j , σ j )ν log (n)/2, where L(β j , σ j ) is the unpenalized log-likelihood and ν is the number of free parameters in the model. Then the mean of the λ j parameters is taken as the regularization parameter and used for all the regressions : λ = p j=1 λ j /p. Since BIC is an asymptotic criterion, taking the average of the regularization parameters over all the regressions helps to improve network inference performance.

Results

Results

Simulation study

Multivariate Poisson data simulation

In order to simulate multivariate Poisson data, we use a method described by Karlis [START_REF] Karlis | Multivariate Poisson regression with covariance structure[END_REF]. As an illustration, for a two dimensional multivariate Poisson distribution, we simulate three independent Poisson variables (X 1 , X 2 , X 12 ) and sum them up (Y 1 = X 1 + X 12 and Y 2 = X 2 + X 12 ) so that the resulting variables are not independent : cov(Y 1 , Y 2 ) = 0 if E(X 12 ) = 0. In the general case, a sample y of dimension (n × p) where p is the number of nodes in the network, n the number of samples is obtained by summing samples from (p + p(p -1)/2) independent Poisson random variables. The adjacency matrix A ∈ {0, 1} p×p encodes the underlying graph structure : A ij = 0 means that the expression level of genes i ∈ 1, . . . , p and j ∈ 1, . . . , p are conditionally independent given the other gene expression levels. In order to sum the (p + p(p -1)/2) terms accordingly, we fix the matrix B of dimension (p×(p+p(p-1)/2) : B = [I p ; P ⊙(I p tri(A) ′ )] ′ where P is a permutation matrix of dimension (p×(p(p-1)/2)) of vector (1,1,0,...,0), ⊙ denotes the matrix multiplication element by element and tri(A) is the vector of dimension (p(p-1)/2)×1 containing the elements of the upper triangular adjacency matrix. The matrix product y = BX gives a count data table of size n × p : n samples from a p-dimensional Poisson random variable whose underlying dependency structure is encoded in the known A matrix.

RNA-seq data are known to be overdispersed relative to a Poisson distribution with the sample variance of a gene expression vector larger than the sample mean. In our simulation study, we also consider the possibility of inflating the variance of the independent Poisson random variables used in the X matrix of the formula above by simulating independent variables according to a log-normal Poisson model. For gene j and sample i, we sample

X ij ∼ P(µ ij ) with log(µ ij ) = θ j + ε ij , ε ij ∼ N (0, σ 2 j ).
We use this log-normal Poisson distribution only for the first p columns of the matrix, the other columns being sampled from a simple Poisson distribution.

Simulation settings

The three methods were compared on two sets of simulations : multivariate Poisson data and overdispersed multivariate Poisson. For each type of data, we simulated 50 different adjacency matrices A with a scale-free structure. This implies that degrees of the edges are assumed to follow a power law distribution, i.e. few nodes in the network are well connected and most of the nodes have only one or two neighbours. The number of nodes p was set to 50. With a scale-free structure, the maximum degree of a node is k max = 35 and the average degree is less than 2. To avoid the ultra-high dimensional setting, defined as k log( p k )/n ≥ 1 2 for Gaussian linear regression [START_REF] Verzelen | Minimax risks for sparse regressions : Ultra-high dimensional phenomenons[END_REF], we set the number of biological samples to n = 100. For each of the 50 different adjacency matrices, 1225 samples of size n were simulated from Poisson random variables (adding extra inter-sample variance or not) and summed up as explained above to obtain the final data set of size 100 × 50. We chose to use Poisson distributions of mean µ = 100 to build the X data matrix, resulting in data counts ranging from around 100 to 2500. In the case of Poisson data with inflated variance, the parameter σ j was set to 0.25, which is slightly smaller than the amount of dispersion observed in the real data presented below.

To evaluate the different methods, we tried to infer the adjacency matrix A from the simulated dataset y (100×50) and compared the inferred matrix A pred with the real adjacency matrix A used to simulate the data. For each type of data (with and without extra inter-sample variance) and for each network inference method (Gaussian, loglinear Poisson, and the proposed hierarchical log-normal Poisson graphical models), Receiver Operating Characteristic (ROC) curves were constructed by varying values of the regularization parameter from an empty network (sensitivity equal to 0) to a full network (specificity equal to 0). The sensitivity and specificity values were also compared for the different methods using the chosen regularization parameter (with the BIC criterion for the Gaussian graphical model, StARS criterion for the log-linear Poisson graphical model and the mean-BIC criterion presented above for the hierarchical log-normal Poisson model). Note that in the case of the Poisson graphical model, a power transformation is applied only in the simulation setting inducing inflated variance.

Results

ROC curves, averaged over the 50 simulated datasets, are presented in Figure 7.1 for the two simulation settings (multivariate Poisson data with or without inflated variance). It can be noticed that in the first setting, with no over-dispersion, the log-linear Poisson model outperforms the Gaussian graphical model applied to transformed data. This result was already observed [START_REF] Allen | A log-linear graphical model for inferring genetic networks from high-throughput sequencing data[END_REF]. As expected, in this case the performance of the log-linear Poisson model and the proposed hierarchical model are very similar. When adding extra variability to the data, we are compelled to use a powertransformation of the data to apply the log-linear Poisson model [START_REF] Allen | A log-linear graphical model for inferring genetic networks from high-throughput sequencing data[END_REF], since the data no longer respect the Poisson assumption of equal mean and variance. The performance of the log-linear Poisson model in this case is considerably deteriorated, and is now comparable to the poor performance of the Gaussian graphical model on log-transformed data. The proposed hierarchical log-normal Poisson model therefore outperforms the two other methods in this case, keeping in mind that the data were simulated under a closely related model that was deemed to be a reasonable choice to approximate the dynamics of RNA-seq data. It has to be pointed out that for the overdispersed data, performances of the three methods are considerably worse compared to the simple case of multivariate Poisson data due to the presence of additional variability.

Sensitivity and specificity obtained by each method for the chosen regularization parameters are represented in diamond-shape squares on the ROC curves (Figure 7.1) and are summarized in Table 7.1.

The regularization parameter chosen with the mean-BIC criterion for the proposed hierarchical log-normal Poisson model offers a higher sensitivity than the Poisson or Gaussian graphical models, even when no over-dispersion was simulated (0.84 compared to 0.71 and 0.57, respectively), while keeping a high specificity (0.97 compared to 0.99 and 0.98, respectively). The number of correctly detected edges is therefore larger for the proposed model compared to the other two methods, even in the case of multivariate [START_REF] Allen | A log-linear graphical model for inferring genetic networks from high-throughput sequencing data[END_REF] on power transformed data (y α ), Hierarchical model : proposed model as detailed in the Methods section and applied on the raw data.

Poisson data with no over-dispersion. When adding extra inter-sample variability, the differences between the three methods are even larger, even if the performances deteriorate for all methods (sensitivity equal to 0.4 for the proposed model compared to 0.1 for the Gaussian graphical model and 0.05 for the Poisson graphical model). These very low sensitivity values can partly be explained by the fact that scale-free structures were considered for the simulated graphs, therefore generating only a small number of edges compared to a random graph structure that are difficult to correctly detect. This also explains, on the other hand, the high specificity values. In fact, as the models infer very few edges for low numbers of biological replicates, they have less chance to detect incorrect edges. Both the ROC curves and the sensitivity/specificity for the chosen regularization [START_REF] Allen | A log-linear graphical model for inferring genetic networks from high-throughput sequencing data[END_REF] on power transformed data (y α ), Hierarchical model : proposed model as detailed in the Methods section and applied on the raw data.

parameter therefore show much better performances for the proposed hierarchical model than the Gaussian graphical model on log-transformed data or the Poisson graphical model on power-transformed data, especially in the case of overdispersed multivariate Poisson data.

Figure 7.2 represents the relationships between the degree of the nodes in the estimated network and in the simulated structure for both the Poisson graphical model and the proposed hierarchical model. It can be observed that, as expected, in the case of no over-dispersion, both methods perform quite similarly, as already seen in the ROC curves above. In the case of over-dispersion, however, even if the sensitivity was quite poor for all methods (Table 7.1), the structure of the graph was much better preserved with the proposed model than with the Poisson graphical model on power transformed data.

To ensure that these results do not depend on the scale-free structure of the graphs, we have drawn ROC curves and performed similar model selection on data simulated with an Erdos-Rényi structure [START_REF] Erdos | On Random Graphs[END_REF]) (Figure 7.3 and Table 7.2). For Erdos-Rényi graphs, each pair of nodes are connected with the same probability, independently of the other pairs of nodes. As previously observed [START_REF] Allen | A log-linear graphical model for inferring genetic networks from high-throughput sequencing data[END_REF], although the differences among the three methods are less pronounced for Erdos-Rényi structures than for scale-free structures, the same general conclusions hold.

Real data analysis

Data description

The three methods were applied to a publicly available microRNA-seq data set available at The Cancer Genome Atlas (TCGA) Data Portal (http://cancergenome.nih. gov/). We selected 100 samples from breast invasive carcinoma (BRCA) tumors. To avoid being in an ultra high-dimensionality setting [START_REF] Verzelen | Minimax risks for sparse regressions : Ultra-high dimensional phenomenons[END_REF], we reduced the number of microRNAs used for network inference to 50 (among 863). To do so, we first removed all microRNAs that had at least one null count. Among the remaining 207, we selected the microRNAs with the largest inter-sample variance (as suggested by (Allen and Liu, 2012)). These microRNAs are the most likely to be linked to breast cancer development since they are selected among the most highly variable microRNAs. Note that we did not perfom any normalization for differences in library sizes on this data set, as contrary to differential analyses [START_REF] Anders | Differential expression analysis for sequence count data[END_REF][START_REF] Robinson | A scaling normalization method for differential expression analysis of RNA-seq data[END_REF], differences in library sizes have no impact on the network inference results since we do not compare two different biological samples, but relate the expression of genes within each biological sample. Since each miRNA has an equal number of nucleotides, there is no need for a gene length correction either.

Modelling the data

Shapiro-Wilk tests on miRNA expression vectors showed that the data, even for highly expressed miRNAs, could not be directly modelled as a normal distribution [START_REF] Shapiro | An analysis of variance test for normality[END_REF]. We therefore used a Box-Cox transformation [START_REF] Box | An analysis of transformations[END_REF] prior to applying a Gaussian graphical model to these data. The optimal Box-Cox parameter to make the data as normally distributed as possible was found to be close to zero, which corresponds to a log-transformation of the data (Figure 7.4).

For these data, the Poisson assumption is not verified either, as shown in Figure 7.5, since the sample variance is considerably larger than the sample mean for all miRNAs. As suggested in [START_REF] Allen | A log-linear graphical model for inferring genetic networks from high-throughput sequencing data[END_REF], we therefore applied the power-transformation implemented in the PoiClaClu package prior to applying the log-linear Poisson graphical model.

The Gaussian graphical model with the BIC criterion detected 48 edges, the loglinear Poisson graphical model with the StARS criterion [START_REF] Allen | A log-linear graphical model for inferring genetic networks from high-throughput sequencing data[END_REF] detected 74 edges, and the proposed hierarchical log-normal Poisson graphical model detected 369 edges among the 50 miRNAs considered here. As shown in Figure 7.5, these data exhibit significant over-dispersion with respect to the Poisson assumption. We are therefore close to the second simulation setting presented above. In this case, the sensitivity of the proposed hierarchical model is expected to be much higher than for the other two methods, which explains the much larger number of detected edges. Figure 7.6 presents the network inferred by the hierarchical model. Table 7.3 presents the biological functions of the most highly connected nodes found with the proposed hierarchical model. It can be noticed that a large majority of these miRNAs are already known to be related to breast cancer. Further biological validation would be interesting for the remaining ones that could be new potential therapeutic targets.

Discussion

Network inference from RNA-seq data is an important methodological challenge. This work is a pioneer study to provide some guidelines on the best methods to achieve this goal. There are two main approaches. The first and simplest idea is to perform a transformation of the data and apply previously proposed methods for microarray studies based on Gaussian graphical models, for example using a Box-Cox transformation. Another possibility is to apply methods specifically developed for the analysis of count data using Poisson graphical models, either with a power transformation of the data or by accounting for over-dispersion directly in the model using for example a hierarchical log-normal Poisson graphical model as proposed here. We found in both simulation study and real data application that the power transformation did not work well to correct for over-dispersion. It has to be noted that the same α parameter was used here for all the genes. It might be possible to improve the performance of this method if a different coefficient was estimated for each gene. This is, however, not possible with the method proposed by [START_REF] Witten | Classification and clustering of sequencing data using a poisson model[END_REF], which finds the optimal value by maximizing the adequacy criterion for a group of genes. In this work the best suited methodology for network inference from RNA-seq data currently appears to be the proposed hierarchical Poisson log-normal model, which seems to be able to appropriately deal with highly dispersed count data. However, the implementation of this approach based on the R package glmmixedlasso [START_REF] Schelldorfer | Glmmlasso : An algorithm for high-dimensional generalized linear mixed models using 1-penalization[END_REF] is quite slow for a large number of biological samples and more research is needed to optimize this function.

It has to be pointed out that in high-dimensional settings (number of genes much larger than the number of biological samples), all methods were unsurprisingly found to perform very poorly, despite the ℓ 1 regularization. As for microarray studies, the limited number of biological replicates available in RNA-seq experiments considerably restrains the number of genes that can be included in the network. Future research is needed to tackle this issue. A first possibility may be to try to reduce the number of parameters to be estimated. In fact, in a first step we aim at finding the regulatory relationships between genes without necessarily estimating their strength precisely. Therefore, in the regression models presented above, instead of trying to estimate one parameter for each gene we could infer parameters for groups of genes. Alternatively, to face the problem of small numbers of biological replicates, instead of inferring regulatory networks within each experimental condition, it would be interesting to use joint graphical model approaches [START_REF] Guo | Joint estimation of multiple graphical models[END_REF] to jointly infer a network in multiple conditions, thus highlighting the common or differing patterns across conditions. miRNA reference hsa-mir-451 BC [START_REF] Kovalchuk | Involvement of microRNA-451 in resistance of the MCF-7 breast cancer cells to chemotherapeutic drug doxorubicin[END_REF] hsa-let-7b BC [START_REF] Peter | Let-7 and miR-200 microRNAs : Guardians against pluripotency and cancer progression[END_REF]) hsa-mir-486 BC [START_REF] Dalmay | MicroRNAs and the hallmarks of cancer[END_REF] hsa-let-7f-2 cancer [START_REF] Peter | Let-7 and miR-200 microRNAs : Guardians against pluripotency and cancer progression[END_REF]) hsa-mir-150 no reference hsa-mir-145 BC [START_REF] Zou | miR-145 inhibits tumor angiogenesis and growth by N-RAS and VEGF[END_REF] hsa-mir-24-2 BC [START_REF] Srivastava | miR-24-2 controls H2AFX expression regardless of gene copy number alteration and induces apoptosis by targeting antiapoptotic gene BCL-2 : a potential for therapeutic intervention[END_REF]) hsa-mir-200c BC (Gregory et al., 2008), (Peter, 2009) hsa-mir-143 BC (Stahlhut Espinosa and[START_REF] Stahlhut Espinosa | The role of microRNAs in cancer[END_REF] hsa-mir-142 no reference Table 7.3 -Ten most highly connected genes in the network inferred by the proposed hierarchical model. BC corresponds to miRNAs known to be linked to Breast Cancer, with the corresponding references. Chapitre 8

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 0 5

Block diagonal covariance selection for gaussian graphical model in high dimension

Résumé. Les modèles graphiques gaussiens permettent d'inférer et de visualiser les dépendances entre des variables. Ces modèles sont difficiles à estimer lorsque la taille de l'échantillon est plus petite que le nombre de variables. Afin de réduire la dimension du problème, nous proposons une procédure de sélection de modèle non-asymptotique : nous approchons la matrice de covariance par une matrice diagonale par blocs. Pour détecter la structure de cette matrice, nous seuillons la matrice de covariance empirique, le seuil étant choisi à l'aide de l'heuristique de pente. Grâce à cette structure, le problème d'estimation est décomposé en plusieurs sous-problèmes indépendants : dans chaque bloc, nous inférons les dépendances entre variables à l'aide du graphical Lasso. Nous justifions cette procédure par des résultats théoriques. Nous illustrons aussi cette méthode sur des données simulées et des données réelles.

Abstract. Gaussian graphical models provide networks to recover and visualize dependencies between continuous variables. Inferring the graph is difficult when the sample size is smaller than the number of variables. To reduce the number of parameters to estimate in the model, we propose a non-asymptotic model selection procedure : we approximate the covariance matrix of the model by a block diagonal matrix. To detect the structure of this matrix, we threshold the sample covariance matrix, choosing the threshold using the slope heuristic. Subsequently, the estimation problem is divided into several independent problems based on the block diagonal structure : in each block, we estimate the network of dependancies using the Graphical lasso algorithm. We justify our procedure from a theoretical point of view and illustrate it on simulated and real datasets.

Introduction

We are interested in designing a Gaussian graphical model (GGM) from n observations described with p variables in a situation where p is large and n small in comparison to p. Let y be the centered and scaled data matrix of size n × p, y ij denoting the value of observation i for variable j. The vectors (y 1 , . . . , y n ) are assumed to be independent and to follow a p-multivariate normal distribution φ(0, Σ) with Σ j,j = 1 for each j ∈ {1, . . . , p}. The edges of the corresponding graph are the non-zero coefficients of the matrix Θ = Σ -1 , estimated by the graphical lasso defined by (8.1), where S is the empirical covariance matrix. Null coefficients in the inverse covariance matrix Θ indicate conditionally independent variables [START_REF] Whittaker | Graphical Models in Applied Multivariate Statistics[END_REF]. To estimate Θ and recover the graph, we seek to maximise the following penalized log-likelihood :

ℓ λ (Θ) = log det(Θ) -trace(SΘ) -λ || Θ || 1 (8.1) Θ(λ) = argmin Θ∈S ++ p {ℓ λ (Θ)} , for λ ≥ 0,
where S ++ p (R) is the set of symmetric and positive matrices on R. A standard algorithm to solve this problem is the graphical lasso implemented in the R glasso package [START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF]. Even if the graphical lasso is designed to solve high-dimensional problems (n < p), the graphical lasso performs poorly when the number of observations n is too small [START_REF] Giraud | Estimation of Gaussian graphs by model selection[END_REF][START_REF] Verzelen | Minimax risks for sparse regressions : Ultra-high dimensional phenomenons[END_REF]. Besides, [START_REF] Mazumder | Exact covariance thresholding into connected components for large-scale Graphical Lasso[END_REF] and [START_REF] Witten | New insights and faster computations for the Graphical Lasso[END_REF] have shown the following property of the graphical lasso solution : if the solution Θ(λ) of the problem (8.1) has a block diagonal structure for some ordering of the variables, the inference problem can be solved in each block independently. We always consider block diagonal structure for some ordering of the variables. Based on the sub gradient equations used in the graphical lasso algorithm [START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF], there is an easy way to check if the solution Θ(λ) will have a block diagonal structure : the solution Θ(λ) and the thresholded sample covariance matrix S (λ) = ( Σjj ′ 1 {| Σjj ′ |>λ} ) jj ′ have the same block diagonal structure. This result induces the following block diagonal screening rule to speed up computations and improve estimation in the graphical lasso. For a given λ, solving the optimization problem in (8.1) amounts to :

Step 1 : identify of the connected components of the graph encoded in the adjacency matrix

(1 {| Σjj ′ |>λ} ) jj ′ ,
Step 2 : solve the graphical lasso problem with regularization parameter λ in each connected component, independently. This block diagonal screening rule has been included by default in the last version of the glasso package [START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF]. It substantially reduces the algorithm complexity from O(p 3 ) to O(p 2 + K k=1 p 3 k ) where p k is the number of variables in block k. This reduction is of great interest and has been used for sparse quadratic discriminant analysis [START_REF] Le | Sparse Quadratic Discriminant Analysis and Community Bayes[END_REF]. We note that the parameter λ used for thresholding the matrix | S | in step 1 is the same parameter λ used for regularization in each graphical lasso sub-problem in step 2. [START_REF] Tan | The Cluster Graphical Lasso for improved estimation of Gaussian graphical models[END_REF] proposed a variation of this two-step procedure, called the Cluster Graphical Model : they detect the connected components by clustering the variables into K groups using a complete linkage clustering on the sample covariance matrix S =| S |, and then, infer a network on each group independently, each with its own regularization parameter λ. Their numerical results suggest that choosing different parameters λ for step 1 and for each subproblem in step 2 improves network inference performance. They choose the number of blocks K (which is equivalent to the selection of threshold λ for S (λ) ) using a leave-one-out algorithm to recast the unsupervised clustering into a supervised one and select the number K giving the smallest mean square error. They also argue that the method is not sensitive to the number of clusters.

In the spirit of the Cluster Graphical Model, we propose to solve the network inference problem using different parameters λ for step 1 and sub-problems in each block for step 2. We pay particular attention to the detection of the connected components in step 1 of the procedure and turn this detection into a model selection problem. We construct a collection of models with different block diagonal structures for the covariance matrix Σ, obtained by thresholding the sample covariance matrix S. Then, we choose a model among this collection using a model selection criterion. [START_REF] Akaike | A model selection criterion for model-based clustering of annotated gene expression data Bibliographie[END_REF] and [START_REF] Schwarz | Estimating the Dimension of a Model[END_REF] have introduced penalized model selection criteria. These criteria are based on asymptotic approximations which are not ensured to hold for small values of n. More recently, [START_REF] Birgé | Minimal penalties for Gaussian model selection[END_REF] have introduced a nonasymptotic criterion : the slope heuristic. The slope heuristic is especially useful when n is small with respect to the model dimension and when the model family is somewhat biased with respect to the sampling distribution. [START_REF] Baudry | Slope heuristics : overview and implementation[END_REF] have provided some practical tools to implement the slope heuristic from [START_REF] Birgé | Minimal penalties for Gaussian model selection[END_REF]. It has proven to be effective in a variety of practical situations : for exemple, [START_REF] Rau | Coexpression analysis of high-throughput transcriptome sequencing data with Poisson mixture models[END_REF] select the number of components in Poisson mixture models on gene expression data using the slope heuristic. [START_REF] Bouveyron | The discriminative functional mixture model for a comparative analysis of bike sharing systems[END_REF] select the number of components in discriminative functional mixture model on data describing bike sharing systems using the slope heuristic. However, the theoretical justification of the slope heuristic has met with several technical difficulties. [START_REF] Birgé | Minimal penalties for Gaussian model selection[END_REF] and [START_REF] Baraud | Gaussian model selection with an unknown variance[END_REF] prove the existence of minimal penalties in heteroscedastic regression with fixed design, [START_REF] Arlot | Data-driven calibration of penalties for least-squares regression[END_REF] extends it for homoscedastic regression with fixed design. Few papers provide theoretical guarantees : [START_REF] Lebarbier | Detecting multiple change-points in the mean of gaussian process by model selection[END_REF] for multiple change point detection, [START_REF] Castellan | Sélection d'histogrammes ou de modèles exponentiels de polynômes par morceaux à l'aide d'un critère de type Akaike[END_REF] for selecting the best partition for histogram construction and [START_REF] Maugis | A non asymptotic penalized criterion for Gaussian mixture model selection[END_REF] for variable selection in mixture models.

In the context of graphical models with block diagonal structure, we provide an oracle inequality to provide an upper bound for the risk between the true model and the model selected among the collection of models with different block diagonal structures. We also provide a lower bound of this same risk. This oracle inequality theoretically supports our model selection procedure based on the slope heuristic with a penalty proportional to the dimension (up to a logarithm term). The main difficulty to prove these theoretical results arises from the fact that the family of models to be considered is random. More precisely, the collection of the block diagonal structures obtained by thresholding is datadriven, and therefore random. Our oracle inequality is derived from a theorem developed in [START_REF] Massart | Concentration inequalities and model selection[END_REF]. For the lower bound of the risk, some results have been previously obtained by [START_REF] Bickel | Regularized estimation of large covariance matrices[END_REF] or [START_REF] Cai | Optimal rates of convergence for covariance matrix estimation[END_REF]. To obtain our lower bound, we use the Birgé lemma developed in [START_REF] Birgé | A new lower bound for multiple hypothesis testing[END_REF] in conjunction with a discretization of the model collection space.

The chapter is organized as follows. In Section 8.2, after providing basic notations Detecting the block-diagonal structure by model selection and definitions, the procedure we propose to detect the block-diagonal structure by nonasymptotic model selection is described. Section 8.3 describes theorical results supporting our model selection procedure. Section 8.4 investigates the numerical performance of our procedure in a simulation study. It is also illustrated through an application on a real gene expression RNA-seq dataset in Section 8.5. After a short discussion, all the proofs are gathered in Section 8.7.

Detecting the block-diagonal structure by model selection

Let y = (y 1 , . . . , y n ) be iid random vectors in R p from a multivariate distribution with density φ p (0, Σ) where Σ j,j = 1 for all j ∈ {1, . . . , p}. The matrix Σ is assumed to have a block diagonal structure for some ordering of the variables. We note K the number of blocks, B = (B 1 , . . . , B K ) the partition of variables into K blocks, B k the subset of variables in block k, p k the number of variables in block k. The matrix Σ can be written as a block diagonal matrix with Σ k ∈ S ++ p k (R) for k ∈ {1, . . . , K} :

Σ B =    Σ 1 0 0 0 . . . 0 0 0 Σ K    . (8.2)
We define B the set of all possible partitions of variables. This set is large : there are p k=1 S(p, k) possible partitions where S(p, k) denotes the Stirling number of the second kind. Considering each element in the set B would not be possible. Following the block diagonal screening rule described in [START_REF] Mazumder | Exact covariance thresholding into connected components for large-scale Graphical Lasso[END_REF], we consider the grid of threshold parameters Λ = ( Σj,j ′ ) j>j ′ derived from the empirical covariance matrix S. We restrict our attention to the sub-collection B Λ of B :

B Λ = B(λ), λ ∈ Λ where B(λ) are the connected components from E (λ) = (1 {|S j,j ′ |>λ} ) jj ′ . (8.3)
We define the model collection F = (F B ) B∈B where F B is defined as follows for every B ∈ B :

F B = {f B with Σ B ∈ S B } (8.4) S B =      Σ B ∈ S ++ p (R) Σ B = P σ    Σ 1 0 0 0 . . . 0 0 0 Σ K    P -1 σ , P σ a permutation matrix, Σ k ∈ S ++ p k (R) for k ∈ {1, . . . , K}      . (8.5)
The dimension of the model

F B is D B = K k=1 p k (p k -1) 2
. The F B contains all densities in the form φ(0, Σ B ). For a given partition B, we define the maximum likelihood where S k is the maximum likelihood estimator of the covariance matrix for the dataset restricted to the block k. We select the best partition using the following criterion :

B = argmin B∈B Λ - 1 n n i=1 log( fB (y i )) + pen(B) , (8.7)
where 1 n n i=1 log( fB (y i )) is the log-likehood of the model and pen a penalty term to be defined.

A first approach to select B would be to consider the classical BIC criterion proposed by [START_REF] Schwarz | Estimating the Dimension of a Model[END_REF], which leads to choose the following penalty :

pen(B) = log n 2 D B .
Since the Gaussian graphical model fulfills regularity conditions, the BIC criterion asymptotically select the true partition B = B ⋆ where B ⋆ is the true partition of variables into blocks, i.e. y is an iid random vector from a multivariate distribution with parameters 0 and Σ B ⋆ [START_REF] Haughton | On the choice of a model to fit data from an exponential family[END_REF].

In our context, the number of observations n cannot be considered to be as large. For this reason, we consider a non-asymptotic model selection based on the slope heuristic, developed by [START_REF] Birgé | Minimal penalties for Gaussian model selection[END_REF]. This heuristic leads to the following penalty :

pen(B) = κD B , (8.8)
where κ is a coefficient to calibrate. The model selected using this coefficient satisfies an oracle inequality. Theoretical guarantees are described in the following Section 8.3. Although the coefficient κ can be computed based on theory, we prefer to calibrate it from the data, using one of the two methods described in [START_REF] Baudry | Slope heuristics : overview and implementation[END_REF]. One method is the Slope Heuristic Dimension Jump (SHDJ) : we approximate the coefficient κ by twice the coefficient corresponding to the biggest dimension jump on the graph representing the model dimension as a function of the coefficient κ. Another method is the Slope Heuristic Robust Regression (SHRR) : we approximate the coefficient κ by twice the slope of a robust regression between the log-likehood and the model dimension.

The two methods are derived from the same heuristic and give similar results. They are implemented in the R package capushe [START_REF] Baudry | Slope heuristics : overview and implementation[END_REF].

Theoretical results for non-asymptotic model selection

Model selection based on the slope heuristic in conjunction with the calibration of the κ coefficients by dimension jump (SHDJ) or robust regression (SHRR) have been proven

Theoretical results for non-asymptotic model selection to be effective in a variety of practical situations. For example, [START_REF] Rau | Coexpression analysis of high-throughput transcriptome sequencing data with Poisson mixture models[END_REF] select the number of components in Poisson mixture models on RNA-seq gene expression data using the slope heuristic. [START_REF] Bouveyron | The discriminative functional mixture model for a comparative analysis of bike sharing systems[END_REF] select the number of components in discriminative functional mixture models on data describing bike sharing systems using the slope heuristic. However, they did not provide any theoretical justification for their procedures.

In contrast, we do provide theoretical justification for our criterion based on an oracle inequality. [START_REF] Lebarbier | Detecting multiple change-points in the mean of gaussian process by model selection[END_REF] have provided theoretical justification based on an oracle inequality for model selection in multiple change point detection, and [START_REF] Maugis | A non asymptotic penalized criterion for Gaussian mixture model selection[END_REF] for variable selection in mixture models. However, few papers provide a minimax lower bound, which we do have. Remark that the theoretical justification of the slope heuristic has encountered several technical difficulties. The existence of minimal penalties is proved in heteroscedastic regression with fixed design [START_REF] Birgé | Minimal penalties for Gaussian model selection[END_REF][START_REF] Baraud | Gaussian model selection with an unknown variance[END_REF], and for homoscedastic regression with fixed design [START_REF] Arlot | Data-driven calibration of penalties for least-squares regression[END_REF].

For our block-diagonal structure detection procedure, we prove an oracle inequality for a penalty proportional to the dimension (up to a logarithm term) and a lower bound of the risk between the true model and the model selected among the model collection. This ensures that the selected model is close to the oracle, the best one in estimation among our collection. Both inequalities guarantee that our model selection procedure has an optimal rate of convergence, which is a strong theoretical result. Note that these results are non-asymptotical, which means that they hold for a fixed sample size n.

To state the theorem, we recall the definition of the Hellinger distance between two densities f and g defined on R p ,

d 2 H (f, g) = 1 2 R p ( f (x) -g(x)) 2 dx = 1 - R p f (x)g(x)dx,
and the Kullback-Leibler divergence between two densities f and g defined on

R p , KL(f, g) = R p log f (x) g(x) f (x)dx.
In order to properly define the penalty term used in equation (8.8) to select the best partition of variables B, we work with the following model collection :

F bound = (F bound B ) B∈B (8.9) F bound B = φ(0, Σ B ) ∈ F B , Σ B ∈ S bound B (8.10) S bound B = {Σ B ∈ S B |e m ≤ min(Σ B ) ≤ max(Σ B ) ≤ e M , λ m ≤ Λ min (Σ B ) ≤ Λ max (Σ B ) ≤ λ M } ,
where Λ min (A) and Λ max (A) are the smallest and the largest eigenvalues of the matrix A.

The model collection (8.9) is defined such that covariance matrices have bounded coefficients, which is useful for constructing a discretization of this space. If the matrix has bounded coefficients, we can prove that it has bounded eigenvalues. Nevertheless, to simplify the reading, we denote by λ m and λ M bounds on eigenvalues. In the following, we denote by Adj(Σ) the adjacency matrix associated to the covariance matrix Σ.

Theorem 8.3.1 Let y = (y 1 , . . . , y n ) be the observations, arising from a density f * . Consider the model collection F bound defined in (8.9). Suppose that there exists an absolute constant κ ′ > 0 such that for every partition B in the set of all possible partitions of variables B,

pen(B) ≥ κ ′ D B n 2c 2 + ρ log 1 D B ( D B n c 2 ∧ 1) + (1 ∨ τ ) log 0.792p log(p + 1)
,

where c is an absolute constant, depending only on the model collection. Let fB be the maximum likelihood estimator, B Λ ⊂ B as defined in (8.3), and B selected as follows :

B = argmin B∈B Λ - 1 n n i=1 log( fB (y i )) + pen(B) .
Then, f B satisfies :

E(d 2 H (f * , f B )) ≤ CE inf B∈B Λ inf t∈F bound B KL(f * , t) + pen(B) + (1 ∨ τ ) 1 n (8.11)
for some absolute constant C.

This non-asymptotic result is consistent with the point of view adopted in this work where the number of observations n is limited. The proof is presented in Appendix 8.7.2. This theorem is deduced from an adaptation for a random sub-collection of the whole model collection of a general model selection theorem for maximum likelihood estimator developed by [START_REF] Massart | Concentration inequalities and model selection[END_REF]. This adaptation is proved in Appendix 8.7.2. To apply our theorem, the main assumptions to satisfy are the control of the bracketing entropy of each model in the whole model collection and the construction of weights for each model to control the model collection complexity. Remark that the control of the bracketing entropy is a classical tool to bound the Hellinger risk of the maximum likelihood estimator, and has already been done for Gaussian densities in [START_REF] Maugis | A non asymptotic penalized criterion for Gaussian mixture model selection[END_REF] and [START_REF] Genovese | Rates of convergence for the Gaussian mixture sieve[END_REF]. Theorem 8.3.1 provides a lower bound for the penalty, which ensures a good model selection by penalized criterion : the model selected is as good as possible among the model collection. The only assumption made to state Theorem 8.3.1 is a classical one : we work with bounded parameters for each model as detailed in (8.9). Every constant involved in (8.11) depends on those bounds. Even if the bounds are not tractable in practice, this assumption is plausible. To guarantee a good model selection procedure, we need to assume that the true density of the data is not too far from the constructed model collection. Since a covariance matrix can always be considered to be a block-diagonal matrix, with possibly a single block, the block-diagonal covariance matrix assumption is not a strong one.

To complete this analysis, we provide a minimax lower bound for the risk between the true model and the model selected among the model collection. For the lower bound of the risk, some results have been previously obtained by [START_REF] Bickel | Regularized estimation of large covariance matrices[END_REF] and [START_REF] Cai | Optimal rates of convergence for covariance matrix estimation[END_REF]. To obtain our lower bound, we use the lemma developed in [START_REF] Birgé | A new lower bound for multiple hypothesis testing[END_REF] in conjunction with a discretization of the model collection space, already constructed for the oracle inequality.

inf fB sup f ∈F bound B E(d 2 H ( fB , f )) ≥ C 1 D B n (1 + log C 2 D 2 B
).

(8.12)

This theorem is proved in Appendix 8.7.3. Again, this result does not rely on strong assumptions, and the constants involved are explicit. It is also a non-asymptotic result.

This minimax lower bound obviously shows that since the estimator satisfies to (8.11) it is simultaneously approximately minimax on each set F bound B for every B ∈ B. Theorem 8.3.2 and Theorem 8.3.1 lead to the use of the slope heuristic with a penalty proportional to the dimension to select a model among the collection.

Nevertheless, as typically the case, constants are higher in theory than needed (and not always tractable), and we prefer to compute constants from the dataset in practice using the capushe package developed in [START_REF] Baudry | Slope heuristics : overview and implementation[END_REF].

Simulation study

Simulation setting

We simulate n observations from a p-multivariate normal distribution with a null mean and a block diagonal covariance matrix Σ B as defined in (8.2). We fix the number of variables p = 100, the sample size n = 90 and the partition on variable B ⋆ : we vary the number of blocks among K ⋆ ∈ {1, 15}. For each block indexed by k, we design the Σ k matrix as follows, as done in [START_REF] Giraud | Graph selection with GGMselect[END_REF] : Σ k = T T T + D where T is a random lower triangular matrix with values drawn from a uniform distribution between -1 and 1, and D is a diagonal matrix designed to prevent Σ k to have eigenvalues that are too small.

To perform network inference, we use the graphical lasso algorithm proposed in [START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF] and implemented in the R package glasso, version 1.7. We do not use the same thresholding parameter as in the block structure selection, and we denote by ρ the parameter of the graphical lasso algorithm for inferring the network. We compare the following strategies :

1. Glasso : graphical lasso on the set of all variables, with regularization parameter ρ chosen using the following BIC net criterion :

BIC net (ρ) = n 2 log det Θ(ρ) -trace S Θ(ρ) - log(n) 2 df Θ(ρ) ; (8.13)
where S is the sample covariance matrix, and df means the degree of freedom.

2. CGL : cluster graphical lasso as proposed in [START_REF] Tan | The Cluster Graphical Lasso for improved estimation of Gaussian graphical models[END_REF] : first, the partition on variables is detected using an average linkage hierarchical clustering with 8.4.1 -Simulation setting 109 K = 15 clusters. Note that we set the number of clusters to the optimal number K ⋆ . Subsequently, the regularization parameters in each graphical lasso problem ρ 1 , . . . , ρ K ⋆ are chosen from the Corollary 3 from [START_REF] Tan | The Cluster Graphical Lasso for improved estimation of Gaussian graphical models[END_REF] : the inferred network in each block must be as sparse as possible while still remaining a single connected component.

3. Partitions based on model selection : depending on the procedure, the first step is based on asymptotic (BIC) our non asymptotic model selection (SHRR ou SHDJ).

(a) BIC : The partition BBIC is selected using the BIC criterion [START_REF] Schwarz | Estimating the Dimension of a Model[END_REF].

(b) SHRR : The partition BSHRR is selected using the Slope Heuristic Robust Regression.

(c) SHDJ : The partition BSHDJ is detected using the Slope Heuristic Dimension Jump.

Subsequently, the regularization parameters in each graphical lasso problem ρ 1 , . . . , ρ K are chosen using the BIC net criterion :

BIC net (ρ k ) = n 2 log det Θ(ρ k ) -trace S |k Θ(ρ k ) - log(n) 2 df Θ(ρ k ) , (8.14)
where S |k is the sample covariance matrix on variables belonging to the block k and df means the degree of freedom.

4. True partition (truePart) : first, we set the partition on variables to the true partition B ⋆ . Then, the regularization parameters in each graphical lasso problem ρ 1 , . . . , ρ K ⋆ are chosen using the BIC net criterion (8.14).

We compare the performance of the five methods using the sensitivity (SENS), (SENS = TP/(TP + FN)), the specificity (SPEC) SPEC = TN/N = TN/(TN + FP) and the False Discovery Rate (FDR) (FDR = FP/(TP + FP)) where TN, TP, FN, FP are respectively the number of true negative, true positive, false negative, false positive dependencies detected. A network inference procedure is a compromise between sensitivity and specificity : we are looking for a high sensitivity, which measures the proportion of dependencies (presence of edges) that are correctly identified, and a high specificity, which measures the proportion of independencies (absence of edges) that are correctly identified. The accuracy is a global measure of quality of the network inference performance. The False Discovery Rate is the proportion of dependencies wrongly detected. As a rule of thumb, we expect that the FDR remains below 0.05 for common statistical analysis. We compute the Mean Squared Error (MSE) between the true and the estimated covariance matrices. We also investigate the ability to recover the simulated partition on variables B ⋆ using the hierarchical clustering from [START_REF] Tan | The Cluster Graphical Lasso for improved estimation of Gaussian graphical models[END_REF], the BIC, SHRR and SHDJ partition selection methods. Selected partitions for each method are compared with the simulated partition B ⋆ using the Adjusted Rand Index [START_REF] Hubert | Comparing partitions[END_REF].

Results

Block diagonal covariance matrix Σ with K ⋆ = 15 blocks We simulate data from a multivariate normal distribution with a null mean and a block diagonal covariance matrix Σ with K ⋆ = 15 blocks of approximate equal sizes (6 or 7 variables). The performance of the different strategies is illustrated on Figure 8.1. As expected, the true partition strategy (truePart) performs the best : based on the true partition of variables, network inference problem is easier because we solve problems with smaller dimension. The proposed strategies, based on the BIC, SHRR and SHDJ partitions, improve network inference compared to a simple graphical lasso on the set of all variables (glasso) or compared to the cluster graphical lasso (CGL). Note that the strategy based on the BIC partition performs well. This result is not surprising because data has been simulated under the model. Illustrations of the calibration of coefficient κ are presented in Figures 8.2. In addition, we compare the partition selection methods with an average linkage hierarchical clustering with K = 15 as proposed in the cluster graphical lasso [START_REF] Tan | The Cluster Graphical Lasso for improved estimation of Gaussian graphical models[END_REF]. The Figure 8.3 displays the ARI computed over 100 replicated datasets. Despite the fact that the partition with the hierarchical clustering has the true number of clusters (K = 15), the ARI for the hierarchical clustering is lower than the ARI for the three methods (BIC, SHRR and SHDJ) which do not need to specify the number of clusters K in advance, contrary to the hierarchical clustering. The partition of variables selected by BIC, SHRR and SHDJ agreed for n = 90. When n is smaller (n = 30), the robust regression (SHRR) and dimension jump (SHDJ) give slightly better results than the asymptotic criterion BIC.

Full covariance matrix Σ with K ⋆ = 1 blocks We simulate n = 90 observations from a multivariate normal distribution with a null mean and full covariance matrix Σ. The corresponding network of conditional dependencies is almost a clique. As we suspected, solving the graphical lasso problem in this context is too ambitious as proved by [START_REF] Verzelen | Minimax risks for sparse regressions : Ultra-high dimensional phenomenons[END_REF] : detecting the true network requires the estimation of D = 4950 parameters, with only n × p = 900 data points. The log-likelihood of the model displayed in Figure 8.4. In contrast with the Figure 8.2, no linear tendency is observed between the model dimension and the log-likelihood function for complex models. In this context, no relevant block diagonal structure is detected. Calibration by robust regression (left) and by dimension jump (right) : in this extreme setting, no clear linear tendency (red line) between the log-likelihood and the model complexity for complex models (red points) is observed, the biggest largest jump (dotted blue line) is unclear.

Real data analysis

Pickrell et al. analyzed transcriptome expression variation from 69 lymphoblastoid cell lines derived from unrelated Nigerian individuals [START_REF] Pickrell | Understanding mechanisms underlying human gene expression variation with RNA sequencing[END_REF]. The expression of 52580 genes across 69 observations was measured using RNA-seq. The data is extracted from the Recount database [START_REF] Frazee | ReCount : a multi-experiment resource of analysis-ready RNA-seq gene count datasets[END_REF]. First, we filter weakly expressed genes using the HTSFilter package [START_REF] Rau | Data-based filtering for replicated high-throughput transcriptome sequencing experiments[END_REF]. Among the 9191 remaining genes, we identify the 200 most variable genes and restrict our attention to this set of genes for the following network inference analysis.

First, we select the partition B using model selection as described in equation (8.7). The log-likelihood increases with the number of parameters to be estimated in the model as displayed in Figure (8.5). We notice a linear tendency in the relationship between the log-likelihood and the model dimension for complex models (points corresponding to a model dimension higher than 500). This suggests that the use of the slope heuristic is appropriate for selecting a partition B. The model selected by slope heuristic robust regression and by slope heuristic dimension jump described in section 8.2 are the same. The number of blocks detected is KSH = 150 and the corresponding model dimension is D BSH = 283. The partition BSH yields 4 blocks of size 18, 13, 8 and 5, 4 blocks of size 3, 2 blocks of size 2 and 140 blocks of size 1. The partition selected by the Slope Heuristic offers a drastic reduction of the number of parameters to infer, as compared with the graphical lasso performed on the full set of variables which corresponds to a total of D = 19900 parameters to estimate. We also compare the partition BSH selected by the slope heuristic with that selected by BIC BBIC : KBIC = 99 and D BBIC = 2213 with 6 [START_REF] Pickrell | Understanding mechanisms underlying human gene expression variation with RNA sequencing[END_REF] dataset. Calibration by robust regression (left) : the loglikehood of the model is represented as a function of the model dimension. Based on the slope heuristic, the slope of the regression (red line) between the log-likehood and the model dimension for complex models corresponds to the minimal coefficient κ min . The optimal penalty is twice the minimal penalty. Calibration by dimension jump (right) : the dimension of the model is represented as a function of the κ coefficient. Based on the slope heuristic, the biggest jump corresponds to the minimal coefficient κ min . In both cases, the optimal penalty is twice the minimal penalty.

The networks within each cluster of variables are inferred using the graphical lasso algorithm of Friedman [START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF] implemented in the glasso package, version 1.7. The regularization parameter for the graphical lasso on the set of all variables is chosen using the BIC net criterion (8.13). The model inferred based on partition BSH is more parsimonious and easier to interpret than the model inferred on the full set of variables or the model based on the partition BBIC . An illustration of inferred networks in the four largest connected components of the partition BSH are displayed on Figure 8.6. These three networks might be good candidates for further study.

Discussion

In this chapter, we propose a non-asymptotic procedure to detect block diagonal structure for covariance matrices in Gaussian graphical models. Subsequently, we infer the network in each block using the graphical lasso algorithm. The procedure substantially reduces the number of parameters to estimate in the model. Although Gaussian graphical models are widely used in practice, the limited sample size forces the user to restrict the number of variables to include in the model. Usually, this restriction is performed manually, for instance, based on prior knowledge on the role of variables. Second case : p(p -1)/2 ≤ 4D B

We use the Varshamov-Gilbert lemma (see for example [START_REF] Massart | Concentration inequalities and model selection[END_REF], Lemma 4.7) to construct a discretized space for the covariance, and construct FB (r) as previously. Then, Birgé's Lemma (Lemma 8.7.4) involves :

Conclusion générale et perspectives

Cette thèse regroupe des contributions méthodologiques utiles à l'analyse des données RNA-seq qui sont discrètes, hétérogènes et présentent une disproportion entre le faible nombre de réplicats biologiques (n ∼ 100) et le grand nombre de gènes (p ∼ 5000). Nous avons concentré notre attention sur deux analyses différentes : l'inférence de réseau à l'aide de modèle graphique et l'analyse de co-expression à l'aide de modèle de mélange. Dans un premier temps, nous avons proposé une méthode de filtrage des données RNAseq supprimant les gènes peu exprimés à l'aide d'un seuil calibré sur les données. Cette méthode permet de réduire le nombre de gènes préliminairement à l'analyse différentielle, mais elle peut aussi s'avérer utile pour l'analyse de co-expression ou l'inférence de réseaux.

Pour l'inférence de réseau, nous avons proposé un modèle graphique hiérarchique log-linéaire de Poisson modélisant le caractère discret et l'importante variabilité interéchantillon des données RNA-seq. Ce modèle est une alternative à l'utilisation d'un modèle graphique gaussien ou d'un modèle graphique log-linéaire de Poisson sur données transformées. Cependant, les développements autour des modèles graphiques de Poisson sont difficiles car il n'existe pas de manière d'écrire la vraisemblance d'un modèle graphique de Poisson sans imposer des contraintes sévères sur les dépendances modélisées. Dans les modèles graphiques de Poisson, l'estimation des paramètres est effectuée localement, par sélection de voisinage ou neighborhood selection en anglais, et les estimateurs ne correspondent pas au maximum de vraisemblance. Des développements théoriques et méthodologiques sont nécessaires pour pouvoir travailler sur un modèle graphique de Poisson joint proprement défini et n'ayant pas les contraintes actuelles du modèle actuellement proposé.

Dans le cadre de l'analyse de co-expression, nous avons proposé de mettre en concurrence deux modélisations des données RNA-seq (poissonienne ou gaussiennes sur données transformées) en les comparant à l'aide d'un critère de sélection de modèle mesurant l'adéquation de ces modèles aux données, calculés à partir des estimateurs du maximum de vraisemblance des modèles respectifs. Des études préliminaires réalisées sur plusieurs jeux de données réelles montrent que le modèle de mélange gaussien semble bien adapté aux données RNA-seq mais des comparaisons sur d'autres jeux de données sont nécessaires pour valider cette conclusion. Plus généralement, la question du choix de transformation est souvent considérée comme une étape de pré-traitement des données alors qu'elle peut s'interpréter comme un problème de choix de modèle. Nous pourrions ainsi déterminer la meilleure transformation possible des données RNA-seq parmi une liste de transformations candidates : transformations Box-Cox, transformations de puissance ou arc sinus hyberbolique. Ce type de comparaison est utile pour l'utilisation du modèle de mélange de lois gaussiennes en analyse de co-expression ou l'utilisation du modèle graphique gaussien en inférence de réseau des données RNA-seq, mais également pour tout type de données et tout type de modèles probabilistes dont les paramètres sont calculés par maximum de vraisemblance.

Peu importe le type de modélisation, gaussienne ou poissonienne, la performance des méthodes d'inférence de réseaux est limitée par le faible nombre de réplicats à disposition. En pratique, les utilisateurs de ces méthodes doivent sélectionner un sousensemble de gènes pour réduire la dimension du problème d'inférence, mais le choix de ces gènes est souvent arbitraire. Dans le cadre du modèle graphique gaussien, nous avons proposé une méthode pour décomposer le problème d'inférence en plusieurs problèmes de taille plus petite à l'aide d'un argument de sélection de modèle non asymptotique. Cette méthode permet de travailler sur un nombre de gènes de départ plus élevé et réduit la part d'arbitraire dans la sélection des variables à inclure dans le réseau. Afin de réduire la dimension du problème, nous aurions également pu adopter le point de vue inverse : plutôt que de décomposer le problème d'inférence en plusieurs groupes et d'inférer le réseau dans chaque groupe (inférence intra-groupe), nous aurions pu agréger les variables entre elles et inférer un réseau de dépendance entre ces groupes (inférence inter-groupe). Il est cependant difficile de définir ces groupes et de les interpréter (modules de gènes co-exprimés ou co-régulés), de définir un critère d'agrégation adéquate pour représenter ces groupes (gènes moyens ou gènes médians), ni d'interpréter le réseau entre ces groupes (régulation d'un module de gènes sur un autre). Dans tous les cas, nous devons rappeler que la dépendance existant entre deux gènes peut prendre de multiples formes et ne peut être résumée à un simple indice statistique.

La méthode de réduction de dimension proposée permet d'inférer des réseaux sur des sous-ensembles de gènes indépendants à partir d'un ensemble de gènes plus large. Cependant, avec un faible nombre de réplicats (e.g. inférieur à une centaine), il est toujours nécessaire d'effectuer une pré-sélection, ce qui nous empêche de travailler sur la totalité des gènes du jeu de données. Afin de pouvoir travailler sur la totalité du jeu de données, sans omettre certains gènes exprimés de l'analyse, nous nous sommes intéressés aux méthodes de classification qui ne souffrent pas du manque de données puisque le problème statistique est inversé par rapport au problème d'inférence de réseaux : les unités statistiques, nombreuses, correspondent aux gènes et les variables, peu nombreuses, correspondent aux réplicats biologiques. Les hypothèses de départ et les objectifs des deux analyses sont différents, mais il serait intéressant, d'un point de vue biologique, de développer des méthodes de détection de gènes co-régulés à partir de modules de gènes co-exprimés. Afin d'améliorer l'interprétabilité des analyses de co-expression des données RNA-seq par modèle de mélange, nous avons pensé à inclure des variables externes extraites des bases de données d'annotation fonctionnelle de gène dans l'analyse. Nous avons proposé le critère ICAL permettant de sélectionner la partition des données la plus pertinente au regard des données et des annotations externes. Ce critère facilite ainsi l'interprétation de l'analyse de co-expression. Il contribue à améliorer la boucle de rétroaction suivante : les analyses de co-expression successives viennent enrichir les bases de données d'annotation fonctionnelle, et les informations externes de ces bases de données sont à leur tour utilisées pour améliorer la pertinence des nouvelles analyses de co-expression. Cependant, pour que cette boucle de rétroaction soit efficace, le choix des annotations externes à inclure dans l'étape de sélection de modèle doit être effectué par un expert. De manière générale, toute analyse statistique pertinente ne peut être effectuée de manière

4. 1

 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Conclusions and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . II Classification par modèle de mélange 5 Transformation des données et comparaison de modèles pour la classification des données RNA-seq 5.1 Les modèles de mélange pour la classification des données RNA-seq . . .10TABLE DES MATIÈRES 5.2 Transformation des données et comparaison de modèles . . . . . . . . . . 5.3 Illustration sur des données simulées . . . . . . . . . . . . . . . . . . . . 5.4 Illustration sur des données réelles . . . . . . . . . . . . . . . . . . . . . . 5.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 A model selection criterion for model-based clustering of annotated gene expression data III Modèle graphique pour l'inférence de réseaux 7 A hierarchical Poisson log-normal model for network inference from RNA sequencing data 7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2 Materials and Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 Block diagonal covariance selection for gaussian graphical model in high dimension 8.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.2 Detecting the block-diagonal structure by model selection . . . . . . . . . 8.3 Theoretical results for non-asymptotic model selection . . . . . . . . . . . 8.4 Simulation study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.5 Real data analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.7 Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . séquençage haut-débit RNA-seq 1.1 La technologie de séquençage haut-débit RNA-seq 1.1.1 Mesurer l'expression des gènesL'acide désoxyribonucléique, ou ADN, est une macromolécule constituée d'un enchainement de quatre types d'acides aminés : adénine (A), cytosine (C), guanine (G) ou thymine (T). Cette molécule est contenue dans chaque cellule des êtres vivants. Elle contient l'information génétique, appelé le génome, pour chaque individu. Le génome désigne l'ensemble du matériel génétique d'un individu codé dans son ADN. Dans la cellule, l'ADN est structuré sous forme de chromosomes et certaines zones de cette molécule d'ADN, appelé des gènes, représentent des régions codantes qui traduisent en ARN puis protéines l'information contenue dans l'ADN. Ces zones sont appelées les gènes. Le dogme central de la biologie moléculaire schématisé en Figure 1.1 décrit le mécanisme d'expression des gènes. De manière schématique, une molécule d'acide ribonucléique (ARN) est transcrite à partir d'un gène, puis une protéine est traduite à partir de cet ARN. Le passage de la molécule d'ADN à la molécule d'ARN est la transcription, réalisé à l'aide de l'ARN polymérase. Le passage de la molécule d'ARN à la protéine est la traduction. Le transcriptome désigne l'ensemble des molécules d'ARN présentes dans la cellule.

Figure 1

 1 Figure 1.1 -Schéma du dogme central de la biologie moléculaire. Figure extraite de Rau (2010).

  Figure 1.1 -Schéma du dogme central de la biologie moléculaire. Figure extraite de Rau (2010).

  1. Lecture des brins d'ADN : L'ADNc est coupé en petit morceaux de longueur de 200 à 300 paires de bases. Cette longueur dépend du séquenceur utilisé (e.g. Illumina, 454) et évolue progressivement vers des reads de plus en plus longs. Les petits morceaux d'ADNc sont lus. Les lectures correspondantes sont appelées des reads. Les séquences codantes pour chaque read sont répertoriées dans un fichier texte au format FASTQ. Ce fichier contient également des informations sur la qualité du séquençage, comme décrit dans la Figure 1.2.

Figure 1

 1 Figure 1.2 -Extrait du descriptif d'un read dans un fichier .FASTQ. La première ligne contient l'identifiant du read précédé du signe @, la seconde contient la séquence d'ADN, les lignes 3 et 4 permettent l'identification et le contrôle de qualité du read.

Figure 1 . 3 -

 13 Figure 1.3 -Schéma d'alignement de deux reads (en rouge et bleu) sur un génome de référence. En réalité, les reads séquencés sont longs d'au moins 50 bases (A,C, T ou G). L'alignement est rarement parfait, un nombre de bases peut différer entre le read et le génome de référence comme l'illustrent les bases en gras.

  l'observation i appartient à la composante k sachant les données y et les paramètres courant du modèle θ (b)

Figure 2

 2 Figure 2.1 -Données simulées sous un mélange de quatre lois normales bivariées. Les couleurs et formes des points correspondent aux différentes composantes du modèle de mélange.

Figure 3

 3 Figure 3.1 -Indice global J ⋆ s (y) calculé sur les données normalisées Bottomly (Bottomly et al., 2011) pour un seuil variant de 1 à 200. La courbe bleue est une courbe LOESS (locally weighted scatterplot smoothing). La croix rouge indique le maximum de l'indice J ⋆ s (y) et la ligne rouge indique le seuil correspondant et sélectionné. Cette figure est extraite de Rau et al. (2013).

  3.2.1 -Transformation des données et comparaison de modèles pour la classification des données RNA-seq 43

  Données de puces Données RNA-seq données continues données de comptage Analyse Binomiale Négative : DESeq ; EdgeR différentielle Gaussienne : limma (Anders & al, 2010 ; Robinson & al, 2010), (Smyth & al, 2005) ou Gaussienne sur données transformées : limma + "voom" (Law & al, 2014 ) Classification Poisson : HTSCluster par modèleGaussienne : Rmixmod(Rau & al, 2015)

  .1) où U r est le nombre de catégories de la variable catégorielle externe u r , n r kℓ = card i : z ik = 1 and u r iℓ = 1 , et n k. = Ur ℓ=1 n r kℓ . Le terme additional dans l'Equation (3.1) mesure la force du lien entre les variables catégorielles externes u et la classification z.

Figure 3

 3 Figure 3.2 -Données simulées sous un mélange de 4 lois normales bivariées. Les couleurs des points correspondent aux différentes composantes du modèle de mélange tandis que les formes indiquent l'appartenance aux catégories d'une classification externe en deux catégories.

Introduction

Figure 4 . 1 -

 41 Figure 4.1 -Global Jaccard index for the Bottomly data calculated for a variety of threshold values for normalized counts, with a loess curve (blue line) superposed and data-driven threshold value (red cross and red dotted line) equal to s ⋆ = 14.7.

Figure 4 .Figure 4 . 3 -

 443 Figure 4.2 -Log mean expression versus log fold change values for the Strub data.For each filter, genes identified as non-differentially and differentially expressed are drawn in black and colors, respectively, and those filtered from the differential analysis are omitted from the plot. From left to right, the filters are as follows : none, CPM, maximum, mean, RPKM maximum, RPKM mean, and maximum using the global Jaccard index threshold.

Figure 4 . 4 -

 44 Figure 4.4 -Number of DE genes detected in the Sultan data, categorized by normalized base mean for each filter type. Figure made using the ggplot2 package[START_REF] Wickham | ggplot2 : elegant graphics for data analysis[END_REF].

Figure 4

 4 Figure4.5 -Sensitivity (over simulated 300 datasets) to detect differentially expressed genes for a variety of filter types and cutoffs, with simulation parameters based on the Strub data. None : no filter. CPM : genes with a CPM less than one in more than half the samples are filtered. Max : maximum-based filter, using the threshold based on the Jaccard index (J), quantiles (5%, 15%, 30%), or values(5, 10, 25). Mean : mean-based filter, using the threshold based on quantiles (5%, 15%, 30%), or values(5, 10, 25). RPKM.max : maximum RPKM filter, using the threshold based on quantiles (5%, 15%) or the value 0.3. RPKM.mean : maximum RPKM filter, using the threshold based on quantiles (5%, 15%) or the value 0.3.

Figure 4 . 6 -

 46  for the filtering performance on simulated data based on the Sultan data for the CPM, maximum, mean, maximum RPKM and mean RPKM filters over a range of cutoffs. The yellow cross labeled with a "J" corresponds to the filtering sensitivity and specificity for the data-based threshold chosen via the global Jaccard index.

  et comparaison de modèles pour la classification des données RNA-seqRésumé. Les données d'expression issues du séquençage haut-débit (RNA-seq) sont des données de comptage très hétérogènes. Il est naturel de les représenter par des modèles basés sur des lois discrètes comme la loi de Poisson ou la loi binomiale négative. Mais des transformations simples des données peuvent permettre de se ramener à des modèles plus répandus fondés sur des lois gaussiennes. Nous montrons comment comparer objectivement les vraisemblances de ces modèles travaillant sur des données différentes. Nous nous focalisons pour mener ces comparaisons sur des problèmes de classification où les mélanges de Poisson et gaussiens peuvent être mis en compétition.

74

  Les modèles de mélange pour la classification des données RNA-seq 5.1 Les modèles de mélange pour la classification des données RNA-seq 5.1.1 Modèle de mélange de lois de Poisson Rau et al. (2015) ont proposé une paramétrisation du modèle de mélange de lois de Poisson spécifique aux données RNA-seq. Cette paramétrisation prend en compte les différentes conditions expérimentales avec des réplicats biologiques. Nous considérons une version simplifiée de cette modélisation sans prendre en compte les réplicats par condition afin d'alléger les notations.

Figure 5 . 1 -

 51 Figure 5.1 -A gauche, comptages simulés sous un modèle de mélange de lois de Poisson et t-transformés pour la conditions 1 versus la condition 2. Les différentes couleurs et symboles correspondent aux quatre classes simulées. A droite, BIC du modèle de mélange de Poisson et BIC du modèle de mélange gaussien ajusté pour la transformation de données t pour un nombre de classes variant de 3 à 15. Le point sur chaque courbe BIC indique le nombre de classes sélectionné par chacun des modèles.

Figure 5

 5 Figure 5.2 -BIC du modèle de mélange de Poisson sur données brutes et du modèle de mélange gaussien sur données transformées pour les données de Sultan (haut, gauche), de Graveley (haut, droite) et de Mach (bas).

  Poisson log-normal model for network inference from RNA sequencing data 7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2 Materials and Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2.1 Gaussian graphical model . . . . . . . . . . . . . . . . . . . . . . 7.2.2 Log-linear Poisson graphical model . . . . . . . . . . . . . . . . . 7.2.3 Hierarchical log-normal Poisson graphical model . . . . . . . . . . 7.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.3.1 Simulation study . . . . . . . . . . . . . . . . . . . . . . . . . . . Multivariate Poisson data simulation . . . . . . . . . . . . . . . . Simulation settings . . . . . . . . . . . . . . . . . . . . . . . . . . Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.3.2 Real data analysis . . . . . . . . . . . . . . . . . . . . . . . . . . Data description . . . . . . . . . . . . . . . . . . . . . . . . . . . Modelling the data . . . . . . . . . . . . . . . . . . . . . . . . . . 7.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 Block diagonal covariance selection for gaussian graphical model in high dimension 8.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.2 Detecting the block-diagonal structure by model selection . . . . . . . . . 8.3 Theoretical results for non-asymptotic model selection . . . . . . . . . . . 8.4 Simulation study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.4.1 Simulation setting . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.4.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Block diagonal covariance matrix Σ with K ⋆ = 15 blocks . . . . . Full covariance matrix Σ with K ⋆ = 1 blocks . . . . . . . . . . . . 8.5 Real data analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.7 Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.7.1 Model collection and discretization . . . . . . . . . . . . . . . . . Discretization for the adjacency matrices . . . . . . . . . . . . . . SOMMAIRE Discretization for the set of covariance matrices . . . . . . . . . . 8.7.2 Oracle inequality : proof of Theorem 8.3.1 . . . . . . . . . . . . . Model selection theorem for MLE among a random sub-collection Bracketing entropy . . . . . . . . . . . . . . . . . . . . . . . . . . Construction of the weights . . . . . . . . . . . . . . . . . . . . . 8.7.3 Lower bound for the minimax risk : Proof of Theorem 8.3.2 . . . .

Figure 7

 7 Figure 7.1 -ROC curves, averaged over 50 simulated data sets on scale-free graphs. Results are presented for the Gaussian graphical model on log-transformed data (blue), the log-linear Poisson graphical model on power-transformed data (red) and the hierarchical log-normal Poisson model on raw data (black) on multivariate Poisson data (A) and multivariate Poisson data with inflated variance (B). The dotted black lines represent the diagonals.

Figure 7 . 2 -

 72 Figure 7.2 -Relationship between the degree of the nodes in the estimated network and in the simulated network. Results are presented for the log-linear Poisson graphical model without over-dispersion (A) and with over-dispersion (B), for the proposed hierarchical log-normal Poisson graphical model without over-dispersion (C) and with over-dispersion (D). Black dotted lines represent the diagonal, and red lines represent loess curves.

Figure 7

 7 Figure 7.3 -ROC curves, averaged over 30 simulated data sets on Erdos-Rényi graphs. Results are presented for the Gaussian graphical model on logtransformed data (blue), the log-linear Poisson graphical model on power-transformed data (red) and the hierarchical log-normal Poisson model on raw data (black) on multivariate Poisson data (A) and multivariate Poisson data with inflated variance (B). The dotted black lines represent the diagonals.

Figure 7 . 4 -

 74 Figure 7.4 -Optimal parameter for the box-cox transformation of data. Curve obtained with the R package MASS.

Figure 7 . 5 -

 75 Figure 7.5 -Sample mean-variance relationship for the 207 microRNAs.

Figure 7 . 6 -

 76 Figure 7.6 -Network inferred with the hierarchical model. The representation was obtained using the software Gephi (Bastian et al., 2009). The size of nodes represents the number of edges associated with the corresponding gene in the network.

  estimator fB as follows : fB = φ(0, ΣB ) with ΣB = Σk = S k for k ∈ {1, . . . ,K}, (8.6)

  Theorem 8.3.2 Let y = (y 1 , . . . , y n ) be the observations, coming from a density f * . Consider the model collection F bound defined in (8.9), and D B the dimension of the model F bound B for each B ∈ B. Let fB being the maximum likelihood estimator for the model indexed by B. Then, for all B ∈ B, there exists absolute constants C 1 > 0 and C 2 > 0 such that :

Figure 8 . 1 -

 81 Figure8.1 -Performance of network inference methods (glasso : graphical lasso on the set of all variables, CGL : cluster graphical lasso, BIC : network inference based on the partition of variables BBIC , SSHR : network inference based on the partition of variables BSHRR , SHDJ : network inference based on the partition of variables BSHDJ and truePart : network inference based on the partition of variables B ⋆ ) measured by the sensitivity (SENS), the specificity (SPEC), the False Discovery Rate (FDR) of the inferred graph, and the Mean Squared Error (MSE) between the true and estimated covariance matrices over 100 replicated datasets simulated under a p-multivariable normal distribution with a null mean 0 and a block diagonal covariance matrix Σ B ⋆ with p = 100, K = 15, n = 90 and clusters of approximate equal sizes.

Figure 8 . 4 -

 84 Figure 8.4 -Calibration of the κ coefficient on a dataset simulated under a multivariate normal distribution with a full covariance matrix with one K ⋆ = 1 and p = 100, n = 90.Calibration by robust regression (left) and by dimension jump (right) : in this extreme setting, no clear linear tendency (red line) between the log-likelihood and the model complexity for complex models (red points) is observed, the biggest largest jump (dotted blue line) is unclear.

Figure 8 . 5 -

 85 Figure8.5 -Calibration of the κ coefficient on the 200 most variable genes extracted from the[START_REF] Pickrell | Understanding mechanisms underlying human gene expression variation with RNA sequencing[END_REF] dataset. Calibration by robust regression (left) : the loglikehood of the model is represented as a function of the model dimension. Based on the slope heuristic, the slope of the regression (red line) between the log-likehood and the model dimension for complex models corresponds to the minimal coefficient κ min . The optimal penalty is twice the minimal penalty. Calibration by dimension jump (right) : the dimension of the model is represented as a function of the κ coefficient. Based on the slope heuristic, the biggest jump corresponds to the minimal coefficient κ min . In both cases, the optimal penalty is twice the minimal penalty.

Figure 8 . 6 -

 86 Figure 8.6 -Networks inferred on the four largest components detected by slope heuristic. Regularization parameters in each set of variables are chosen using the BIC net criterion (8.14). Numbers indicate genes labels. Red edges indicate positive partial correlation whereas blue edges indicates negative partial correlation between genes.

  (r) ⊂ F bound B , and FB (r) ⊂ F bound B , choosing r = (1 + log(C 2 /D 2 B )) 1/2 , we get that : log(e M /e m ) ρ, and with 0.233 ≤ ρ ≤ 0.234, and C 2 = p(p -1)(e Me m ).

  surée correspond donc à un niveau de fluorescence différentiellement mesuré entre deux conditions expérimentales pour les puces à deux couleurs, ou à l'abondance de chaque ADN ciblé pour les puces à une couleur. Cette mesure d'abondance est continue et relative. Le niveau de précision de la méthode dépend du nombre de fragments d'ADN attachés sur la puce à ADN. Elle nécessite de connaître par avance la séquence d'ADN que l'on souhaite hybrider. Cette technologie ne permet pas de découvrir de nouvelles régions codantes du génome.La technologie SAGE (SerialAnalysis of Gene Expression) : Cette technique a été mise au point dans les années 1990. Contrairement à la technologie des puces à ADN, la technologie SAGE ne requiert pas la connaissance du génome

Cette technique de mesure, aussi appelée technologie microarrays, est apparue dans les années 1990. Les brins d'ADNc sont marqués par une teinture fluorescente ou radioactive. Parallèlement, une puce à ADN est préparée contenant les portions d'ADN de l'individu dont on souhaite mesurer l'expression. Les brins d'ADNc s'hybrident aux ADN présents sur la puce : les nucléotides des brins d'ADNc viennent s'unir aux nucléotides des brins d'ADN de la puce à ADN par complémentarité. Le niveau d'hybridation est alors détecté par radioactivité ou fluorescence pour les différentes conditions expérimentales. L'expression mede l'individu pour être utilisée. Des séquences codantes d'ADNc (de 9 à 14 paires de bases nucléiques), appelées étiquettes ou tags en anglais, sont séquencées à partir des brins d'ADNc. Ces tags sont ensuite recollés bout à bout pour former une molécule d'ADN synthétique, et ainsi révéler la séquence d'ADN (et donc le gène) dont elle est issue. En effectuant cette opération un grand nombre de fois,

  Table1.1 -Tableau de comptages résumant l'abondance des reads par région génomique d'intérêt pour trois réplicats biologiques.

				Modélisation statistique des données RNA-seq
		gène 1 gène 2 gène 3 gène 4 gène 5 gène 6 . . .
	échantillon 1	4	19	2987	0	65	1905 . . .
	échantillon 2	0	18	1206	0	29	121	. . .
	échantillon 3	6	20	12	48	299	169	. . .
	4. http ://www.ensembl.org

Table 1

 1

			Comptages bruts
		gène 1 gène 2 gène 3 . . . nombre total de reads
	échantillon 1	10	10	10	. . .	1000
	échantillon 2	100	100	100	. . .	10000
	Comptages divisés par le nombre total de reads de l'échantillon
		gène 1 gène 2 gène 3 . . . nombre total de reads
	échantillon 1	0.01	0.01	0.01 . . .	1000
	échantillon 2	0.01	0.01	0.01 . . .	10000

.2 -Tableau de comptages résumant l'abondance des reads par région génomique d'intérêt pour deux échantillons (haut). Tableau de comptages normalisés par le nombre total de reads alignés pour chaque échantillon (bas).

  ont formulé quelques recommandations à ce sujet. L'analyse différentielle se décompose en plusieurs étapes : 1. Définition des hypothèses à tester : Pour chaque gène j, l'analyse différentielle détermine si une différence d'expression est observée entre deux ou plusieurs conditions expérimentales à l'aide de tests d'hypothèses. Pour simplifier les notations, on considère la comparaison de deux conditions expérimentales uniquement : condition 1 et condition 2. Le test d'hypothèse détermine s'il existe une différence entre la moyenne λ 1j d'expression du gène j dans la condition 1 et la moyenne λ 2j d'expression du gène j dans la condition 2 :

  ) où la quantité Z est appelée fonction de partition. Il s'agit d'un facteur de normalisation garantissant que la fonction p soit bien une densité de probabilité. Dans cette expression, les fonctions ψ C ne sont pas des densités de probabilité. Le théorème de Hammersley-Clifford montre l'équivalence entre cette factorisa-

tion de la densité p du vecteur aléatoire Y par rapport au graphe G et la propriété de Markov du vecteur aléatoire Y par rapport au graphe G. La notion d'indépendance conditionnelle entre les variables aléatoires est ainsi directement liée à la structure du graphe. Les modèles graphiques non orientés sont plus informatifs que les réseaux de dépendance marginale présentés précédemment. Nous concentrons notre attention sur ces modèles graphiques non orientés. On distingue les modèles graphiques en fonction de la nature, discrète ou continue, des variables aléatoires considérées : modèle graphique gaussien ou modèle d'Ising.

Table 3

 3 Le tableau 3.1 résume les changements d'indices des gènes et échantillons en fonction du modèle considéré.

		Unités statistiques	Variables
		ou observations
	Analyse différentielle	i ∈ 1, . . . , n échantillons	j ∈ 1, . . . , p gènes
	Classification	gènes	échantillons
	Inférence de réseaux	échantillons	gènes

.1 -Indices des gènes et échantillons en fonction du type d'analyse effectuées 40 Filtrage des données RNA-seq

  de taille p k ×p k pour k∈{1,...,K} est le nombre de blocs, B = (B 1 , . . . , B K ) est la partition des variables en K blocs, B k les indices des variables appartenant au bloc k et p k le nombre de variables dans le bloc k. Parmi les modèles F B , on distingue l'estimateur du maximum de vraisemblance fB :

	
	
	
	,
	
	
	où K

http ://bioinformatics.oxfordjournals.org/content/suppl/2013/06/24/btt350.DC1/Full S uppM at.pdf

Remerciements

For data simulated with large inter-sample dispersion, the proposed model performs better than the other methods in terms of sensitivity, specificity and area under the ROC curve. These results show the necessity of methods specifically designed for gene network inference from RNA-seq data.

Introduction

In recent years, high-throughput sequencing technology has become an essential tool for genomic studies. In particular, it allows the transcriptome to be directly sequenced (RNA sequencing), which provides count-based measures of gene expression. Typically, the first biological question arising from these data is to identify genes differently expressed across biological conditions. Because RNA-seq data are known to exhibit a large amount of variability among biological replicates, most methods for differential analysis are based either on overdispersed Poisson [START_REF] Auer | A two-stage Poisson model for testing RNA-seq data[END_REF] or negative binomial models [START_REF] Anders | Differential expression analysis for sequence count data[END_REF]Robinson et al., 2010).

In order to study the relationships between these large numbers of genes, several authors have worked on co-expression networks and used methods based on Pearson correlation [START_REF] Giorgi | Comparative study of RNA-seq and Microarray-derived coexpression networks in Arabidopsis thaliana[END_REF] or canonical correlation [START_REF] Hong | Canonical correlation analysis for RNA-seq co-expression networks[END_REF][START_REF] Iancu | Utilizing RNA-Seq data for de novo coexpression network inference[END_REF], but no specific models have been designed for RNA-seq data. A further question is how these genes interact with each other. Inference of gene networks from transcriptomic data is indeed a key aspect of systems biology that may help unravel and better understand the underlying biological regulatory mechanisms. Various models have been proposed for network inference from microarray data, mainly based on Gaussian graphical models [START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF][START_REF] Meinshausen | High-dimensional graphs and variable selection with the lasso[END_REF]. Until now, very few authors have addressed the question of network inference from RNA-seq data. Some authors simply use methods based on a Gaussian assumption for RNA-seq data [START_REF] Cai | Utilizing RNA-Seq Data for Cancer Network Inference[END_REF]. We propose in this paper to compare various approaches to tackle this issue.

The simplest idea is to perform an appropriate transformation of the data, using for example a Box-Cox transformation [START_REF] Box | An analysis of transformations[END_REF] and apply methods that rely on an assumption of normality. Another possibility is to use models specifically designed for count data with large variability. [START_REF] Allen | A log-linear graphical model for inferring genetic networks from high-throughput sequencing data[END_REF] recently proposed a Poisson log-linear graphical model adapted to count data. This model requires a power transformation of the data [START_REF] Li | Normalization, testing, and false discovery rate estimation for RNA sequencing data[END_REF] when the inter-sample variance is greater than the sample mean. We propose in this paper a hierarchical log-normal Poisson model with a Lasso penalty, which has the advantage of directly modelling inter-sample variability and can therefore be readily applied to the raw data. Performance of these different methods for gene network inference are compared on data simulated under a multivariate Poisson distribution [START_REF] Karlis | Multivariate Poisson regression with covariance structure[END_REF] with various amounts of additional inter-sample variability, as well as on publicly available microRNA-seq data collected on breast invasive carcinoma (BRCA) tumors, downloaded from The Cancer Genome Atlas (TCGA) Data Portal.

Materials and Methods

Information Criterion (BIC) [START_REF] Schwarz | Estimating the Dimension of a Model[END_REF] defined below, where ν represents the number of free parameters in the model :

Note that a single parameter λ is chosen for the entire network.

Log-linear Poisson graphical model

A log-linear Poisson graphical model specifically designed for network inference from count data has been recently proposed [START_REF] Allen | A log-linear graphical model for inferring genetic networks from high-throughput sequencing data[END_REF]. This model is based on a Poisson distribution which assumes the mean and variance to be equal. Therefore, the model does not account for the high dispersion of the data, also called over-dispersion with respect to the Poisson distribution, when the sample variance is higher than the sample mean. To apply it to RNA-seq data, the authors propose to use a power transformation of the data y ij → g(y ij ) = y α ij , with α ∈]0, 1] implemented in the R package PoiClaClu [START_REF] Li | Normalization, testing, and false discovery rate estimation for RNA sequencing data[END_REF]. The coefficient α is chosen to maximize an adequacy criterion between the transformed data y α and a Poisson distribution.

Let z j = (g(y 1j ), ..., g(y nj )) be the transformed vector of expression values for gene j in the n biological samples. It is assumed that the conditional distribution of Z ij given all the other genes z i(-j) = (z i,1 , . . . , z i(j-1) , z i(j+1) , . . . , z ip ) is a Poisson distribution P(µ j ), with log(µ j ) modelled as a linear regression on all the other genes :

The notation z corresponds to a standardization of the log-transformed data. This standardization is a necessity since we model the mean of the gene j and not the random variable itself. An edge is present in the inferred graph if one or both parameters β jj ′ and β j ′ j are different from zero. The log-likelihood for gene j can be written in this case as :

Similar to the previous model, we assume that the vector β j is sparse. We add an ℓ 1 penalty to the log-likelihood (7.4) so that some coefficients in the estimated β j vector are set to 0. Estimation of parameters β j can be obtained by a coordinate gradient algorithm as implemented in the R package glmnet [START_REF] Friedman | Regularization paths for generalized linear models via coordinate descent[END_REF]. We propose to perform the model selection with the Stability Approach to Regularization Selection criterion (StARS), as suggested by [START_REF] Allen | A log-linear graphical model for inferring genetic networks from high-throughput sequencing data[END_REF]. This stability-based method selects the network with the smallest amount of regularization that simultaneously makes the network sparse and replicable under random sampling. Note that we select only one regularization parameter for all the regressions in the network problem.

Simulation study Based on the slope heuristic, the slope of the regression (red line) between the log-likehood and the model dimension for complex models (red points) corresponds to the minimal coefficient κ min . The optimal penalty is twice the minimal penalty. Calibration by dimension jump (right) : the dimension of the model is represented as a function of the κ coefficient. Based on the slope heuristic, the biggest jump (dotted blue line) corresponds to the minimal coefficient κ min . The optimal penalty (blue cross) is twice the minimal penalty.

113

Appendix

Here, we propose a procedure to select relevant subsets of variables based on the data. Therefore, our procedure is of great practical interest to estimate parameters in Gaussian graphical models when the sample size is smaller than the number of parameters to estimate.

The procedure we propose is easy to implement in practice and fast to compute. The calibration of the κ coefficient by robust regression and dimension jump encounter no particular difficulty. Moreover, graphical representation of the log-likelihood or the model dimension can indicate if the block diagonal assumption is wrong, for instance in extreme cases where the covariance matrix is not sparse and no block diagonal structure can be detected.

Our method uses a model selection criterion to detect block diagonal structure. We propose a non-asymptotic criterion supported by strong theoretical results. It is important to highlight that, although BIC can be very efficient when the model is appropriate, it can provide unsatisfactory results in more difficult inference situations, for example, in data with large bias (as is often seen in real data). In such cases, the slope heuristic appears to provide satisfactory results and yield more interpretable models, especially when the model family is somewhat biased with respect to the sampling distribution. The BIC criterion does not deal with the fact that the block diagonal structure assumption is an approximation of the reality, not the truth. The BIC criterion performs well on simulated data because the data has been simulated under the model, with a block diagonal covariance matrix. By taking into account both the model complexity and the size of the model collection considered into the penalty term, the non-asymptotic model selection provides a parsimonious model choice even in situations where the BIC tends to overestimate the number of parameters.

Appendix

In this Appendix, we detail the proof of Theorems 8.3.1 and 8.3.2. First, we describe a discretization of the model collection used, which is useful in the two proofs. Then, in Section 8.7.2, we prove Theorem 8.3.1. We first generalize a model selection theorem for MLE, introduced by Massart, to random model selection. Subsequently, we prove that our model collection satisfies all the assumptions of this Theorem, and deduce the oracle inequality. In Section 8.7.3, we prove Theorem 8.3.2 using Birgé's Lemma (lemma 8.7.4) with the discretization of the model collection obtained in Section 8.7.1.

Model collection and discretization

Discretization for the adjacency matrices

B

, we may identify a corresponding adjacency matrix A B . This matrix of size p 2 could be summarized by the vector of concatenated upper triangular vectors. Then, we construct a discrete space for {0, 1} p(p-1)/2 which is in bijection with : First, we focus on the set {0, 1} p(p-1)/2 . Lemma 8.7.1 Let {0, 1} p(p-1)/2 be equipped with Hamming distance δ. Let {0, 1} p(p-1)/2 B be the subset of {0, 1} p(p-1)/2 of vectors for which the corresponding graph has structure B.

For every α ∈ (0, 1), let β ∈ (0, 1) such that D B ≤ αβp(p -1)/2. There exists some subset R(α) of {0, 1} p(p-1)/2 B with the following properties

where ρ = -α(log(β) + β -1)/ log(αβ) and

Proof. Let R be a maximal subset of {0, 1}

for every x 0 ∈ R, where B A (x, r) = {y ∈ A|δ(x, y) ≤ r}.

Our proof is similar to the proof of Lemma 4.10 in Massart (2007). We consider :

Let α ∈ (0, 1), β ∈ (0, 1) such that D ≤ αβp(p -1)/2. According to [START_REF] Massart | Concentration inequalities and model selection[END_REF], we know that :

Nevertheless, as {0, 1}

As {0, 1}

p(p-1)/2 B corresponds to the stabilizer of x 0 , |{0, 1}

Note that we divide by K! because there are at worst K clusters with the same size.

Using Stirling's approximation, we obtain :

Discretization for the set of covariance matrices

Let α ∈ (0, 1) and β ∈ (0, 1) such that D B ≤ αβp(p -1)/2. Let R(α) as constructed in Lemma 8.7.1, and its equivalent A disc B (α) for adjacency matrices. Let ǫ > 0. Let :

Then,

Proof. Let (Σ, Σ ′ ) ∈ (S disc B (ǫ, α)) 2 with Σ = Σ ′ . If Σ and Σ ′ are close, either they have the same adjacency matrix and they differ only on a coefficient or they differ in their adjacency matrices. In the first case, ||Σ -Σ ′ || 2 2 ≥ ǫ. In the second case, ||Σ -

ǫ, this minimum depending on α and ǫ.

Oracle inequality : proof of Theorem 8.3.1

First, we state the general theorem we use to get the oracle inequality, and its proof. Then, we deduce the oracle inequality by proving that our model collection satisfies all the assumptions.

Model selection theorem for MLE among a random sub-collection

We denote by H [.] (ǫ, S, d H ) the bracketing entropy of the set S with ǫ-brackets according to the Hellinger distance d H . Theorem 8.7.2 Let f * be an unknown density to be estimated from a sample of size n (y 1 , . . . , y n ). Consider {F m } m∈M some at most countable deterministic model collection. Let {w m } m∈M be some family of nonnegative numbers such that :

(8.17)

We assume that for every m ∈ M, H

Moreover, for all m ∈ M, we assume that there exists ψ m on R + such that ψ m is nondecreasing, ξ → ψ m (ξ)/ξ is non-increasing on (0, +∞), and for all ξ ∈ R + , for all u ∈ F m , denoting by

(8.18)

Let η ≥ 0 and consider the collection of η-maximum likelihood estimators { fm } m∈ M. Let pen : M → R + . Suppose that there exists an absolute constant κ > 0 such that for all m ∈ M :

for some absolute positive constant C.

This theorem is a generalization of Theorem 7.11 in [START_REF] Massart | Concentration inequalities and model selection[END_REF] to a random model subcollection of the whole collection. As the proof is adapted from the proof of this theorem, we detail here only differences and we refer the interested reader to [START_REF] Massart | Concentration inequalities and model selection[END_REF].

We denote by γ n the empirical process and by γn the centered empirical process. Following the proof of the Massart's theorem, easy computations lead to :

Appendix where

To bound γn (ŝ m ′ ), we use Massart's arguments. The main difference stands in the control of γn (g m ). As M ⊂ M is random, E(γ n (g m )) = 0. Nevertheless, thanks to the Bernstein inequality, which we may use thanks to the inequality in (8.19), we obtain, for all u > 0, with probability smaller than exp(-u),

where α τ is a constant depending on τ . Then, choosing u = w m for all m ∈ M, where w m is defined in (8.17), some fastidious but straightforward computations similar to those of Massart's lead to Theorem 8.7.2.

We remark that this is a theoretically easy extension, but quite useful in practice, e.g. for controlling large model collections.

Bracketing entropy

Let ǫ > 0 and α > 0. According to Corollary 8.7.1, there exists S ∈ S disc B (ǫ, α) such that :

If we take α = 1ǫ/2D B , we obtain ||Σ B -S|| 2 2 ≤ ǫ. Then we consider :

According to the Proposition 4 in [START_REF] Maugis | A non asymptotic penalized criterion for Gaussian mixture model selection[END_REF]

, d H ) the logarithm of this number, which corresponds to the bracketing entropy, we obtain from Corollary 8.7.1 that :

. We then construct ψ B satisfying Equation (8.18).

For all ξ > 0,

According to [START_REF] Maugis | A non asymptotic penalized criterion for Gaussian mixture model selection[END_REF],

Then, denoting by c = √ log C + √ π, we can define ψ B by :

As we want ξ B such that ψ B (ξ B ) = √ nξ 2 B , we could take :

.

Construction of the weights

We need to control the Bell number, which is the cardinal of B. For this, we use a result of [START_REF] Berend | Improved bounds on bell numbers and on moments of sums of random variables[END_REF], which guarantees the following inequality for p ∈ N :

For every B ∈ B, we know that p ≤ D B ≤ p(p -1)/2. Then, we obtain the following result.

Appendix

Let f * = φ(0, Σ * ) be the true density. Let f be the considered estimator. We define f = argmin

First, we have :

(8.20)

Secondly, we have :

Then, by combining (8.20) and (8.21) we obtain :

We need to design a lower bound for :

For this purpose, we use the Birgé Lemma (Lemma 8.7.4) :

Lemma 8.7.4 Let (P f ) f ∈F a probability family, and (A f ) f ∈F some event pairwise disjoints. Let a 0 = P 0 (A 0 ) and a = min f ∈F P f (A f ). Then :

Then, if use Birgé's Lemma (Lemma 8.7.4) to control max f ∈F B (r) P f (f = f ) in (8.22), we obtain :

if the following inequality is satisfied :

The inequality (8.24) is satisfied if the inequality (8.25) is fulfilled, with :

isolée, sans une connaissance approfondie des aspects biologiques et bio-informatiques du problème.

Plus généralement, nous avons concentré notre attention sur les données RNA-seq et les bases d'annotation de gènes mais d'autres sources d'information pourraient être incluses. Les sources de données se multiplient et se diversifient : l'intégration de données hétérogènes est un défi important. Par ailleurs, le rôle des données n'est plus de valider une hypothèse déjà formulée en amont. Au contraire, les hypothèses sont formulées à partir des données. La part des décisions arbitraires laissées à l'utilisateur des méthodes (choix de modélisation, choix de modèle, sélection de variables pertinentes) doit être réduite le plus possible afin de ne pas conduire à la formulation de fausses hypothèses et de fausses découvertes. Dans ce contexte, les modèles probabilistes et les outils de sélection de modèle associés offrent un cadre rigoureux d'analyse de données qu'il convient d'exploiter dans les années à venir.

Titre : Classification et inférence de réseau pour les données RNA-seq.

Mots clés : modèle de mélange, modèle graphique, données RNA-seq, classification, inférence de réseaux, sélection de modèle Résumé : Cette thèse regroupe des contributions méthodologiques à l'analyse statistique des données issues des technologies de séquençage du transcriptome (RNA-seq). Les difficultés de modélisation des données de comptage RNA-seq sont liées à leur caractère discret et au faible nombre d'échantillons disponibles, limité par le coût financier du séquençage. Une première partie de travaux de cette thèse porte sur la classification à l'aide de modèles de mélange. L'objectif de la classification est la détection de modules de gènes co-exprimés. Un choix naturel de modélisation des données RNA-seq est un modèle de mélange de lois de Poisson. Mais des transformations simples des données permettent de se ramener à un modèle de mélange de lois gaussiennes. Nous proposons de comparer, pour chaque jeu de données RNA-seq, les différentes modélisations à l'aide d'un critère objectif permettant de sélectionner la modélisation la plus adaptée aux données. Par ailleurs, nous présentons un critère de sélection de modèle prenant en compte des informations biologiques externes sur les gènes. Ce critère facilite l'obtention de classes biologiquement interprétables. Il n'est pas spécifique aux données RNA-seq. Il est utile à toute analyse de co-expression à l'aide de modèles de mélange visant à enrichir les bases de données d'annotations fonctionnelles des gènes. Une seconde partie de travaux de cette thèse porte sur l'inférence de réseau à l'aide d'un modèle graphique. L'objectif de l'inférence de réseau est la détection des relations de dépendance entre les niveaux d'expression des gènes. Nous proposons un modèle d'inférence de réseau basé sur des lois de Poisson, prenant en compte le caractère discret et la grande variabilité interéchantillons des données RNA-seq. Cependant, les méthodes d'inférence de réseau nécessitent un nombre d'échantillons élevé. Dans le cadre du modèle graphique gaussien, modèle concurrent au précédent, nous présentons une approche non-asymptotique pour sélectionner des sous-ensembles de gènes pertinents, en décomposant la matrice variance en blocs diagonaux. Cette méthode n'est pas spécifique aux données RNA-seq et permet de réduire la dimension de tout problème d'inférence de réseau basé sur le modèle graphique gaussien.

Title : Clustering and network inference for RNA-seq data.

Keywords : mixture model, graphical model, RNA-seq data, clustering, network inference, model selection Abstract : This thesis gathers methodologicals contributions to the statistical analysis of nextgeneration high-throughput transcriptome sequencing data (RNA-seq). RNA-seq data are discrete and the number of samples sequenced is usually small due to the cost of the technology. These two points are the main statistical challenges for modelling RNA-seq data. The first part of the thesis is dedicated to the co-expression analysis of RNA-seq data using model-based clustering. A natural model for discrete RNA-seq data is a Poisson mixture model. However, a Gaussian mixture model in conjunction with a simple transformation applied to the data is a reasonable alternative. We propose to compare the two alternatives using a data-driven criterion to select the model that best fits each dataset. In addition, we present a model selection criterion to take into account external gene annotations. This model selection criterion is not specific to RNA-seq data. It is useful in any co-expression analysis using model-based clustering designed to enrich functional annotation databases. The second part of the thesis is dedicated to network inference using graphical models. The aim of network inference is to detect relationships among genes based on their expression. We propose a network inference model based on a Poisson distribution taking into account the discrete nature and high inter sample variability of RNA-seq data. However, network inference methods require a large number of samples. For Gaussian graphical models, we propose a nonasymptotic approach to detect relevant subsets of genes based on a block-diagonale decomposition of the covariance matrix. This method is not specific to RNA-seq data and reduces the dimension of any network inference problem based on the Gaussian graphical model.