Comparaison et évaluation des méthodes de pénalisation du modèle de Cox avec les données mRNA-seq 

Séquençage des banques

Le séquençage consiste à déterminer la séquence de nucléotides de chaque fragment (Fig. 1.2.C). Pour cela, des sondes fluorescentes sont attachés au premier nucléotides de chaque fragment. Ces sondes sont spécifiques : une couleur différente pour chacun des nucléotides A, C, G, T. Une image de fluorescence est enregistrée, et cette opération est répétée séquentiellement, pour chaque nucléotide jusqu’à la fin des fragments. Ainsi, la séquence de chaque fragment présent dans chaque banque est connue. Les séquences ainsi générées de chacun de ces fragments sont appelées « lectures ».
Par abus de langage, le terme de séquençage fait référence à l’ensemble du processus permettant de mesurer le niveau d’expression des gènes (RNA-seq).

Alignement sur un génome de référence

Pour quantifier l’expression des gènes dans un échantillon, les lectures sont alignées sur un génome de référence contenant les séquences de nucléotides de l’ensemble des gènes (Fig. 1.2.D). Le niveau d’expression d’un gène se traduit alors par le nombre de lectures qui se sont alignées sur une partie de sa séquence. On parle de « données de comptage ».

Taille des banques et profondeur de séquençage

Pour un échantillon donné, la « taille de la banque » peut se définir de deux manières :
— le nombre total de lectures séquencées.
— le nombre total de lectures séquencées qui ont été alignés sur le génome de référence.
Dans la suite, nous utiliserons la deuxième définition.
Il est important de noter que la taille des banques peut varier d’un patient à l’autre pour un même protocole. Par exemple, quelques fragments d’ADN ne reçoivent pas d’adaptateurs, et ne vont pas être séquencés. De plus, certaines lectures ne sont pas alignées sur le génome (e.g. problème de séquençage de l’adaptateur qui empêche la lecture d’être associée à un patient). Ensuite, l’amplification PCR de séquences d’ADN riches en bases nucléotidiques guanine (G) et cytosine (C) est difficile [MAMMEDOV et collab., 2008], et implique un biais dans les données RNA-seq [BENJAMINI et S PEED, 2012; RISSO et collab., 2011]. Enfin, le nombre de fragments peut varier d’une grille à l’autre, et le multiplexage (A) Chaque couleur (bleu, rouge, noir) correspond à des séquences d’ARN différentes.
Pour chaque échantillon (biopsie d’une tumeur), les molécules d’ARN sont découpées en fragments, transformées en simple brin, et des adaptateurs sont ajoutés. Ces derniers permettent de les identifier, et les différents échantillons peuvent ainsi être mélangés pour être séquencés en parallèle. (B) Les fragments sont placés sur une grille (flow cell) et amplifiés (PCR) pour former des clusters. Le signal obtenu lors de l’étape suivante seraainsi plus important. (C) Des sondes fluorescentes spécifiques d’un nucléotide (A -jaune, C – vert, G – orange, ou T – violet) sont ajoutées et vont venir s’attacher au premier nucléotide de chaque fragment. Une image est alors prise, et cette procédure est répétée jusqu’à ce que les séquences nucléotidiques de chaque fragment soient déterminées. (D)
Les séquences ainsi obtenues, les lectures, sont alors alignées sur un génome de référence pour déterminer de quel gène ils sont issus.

Biais de la normalisation CPM et ajustement

Dans le cas où plusieurs gènes seraient surexprimés (resp. sous-exprimés) entre deux patients, la normalisation CPM introduit un biais. Prenons l’exemple de deux patients A et B pour lesquels trois gènes sont séquencés avec la même profondeur (Fig. 1.3). Imaginons que les trois gènes s’expriment de manière équivalente chez le patient A (i.e. même quantité d’ARN transcrite par chacun des gènes), et que les gènes du patients B ont le même niveau d’expression que ceux du patient A, sauf pour le troisième gène qui ne s’exprime pas. Dans ce cas, les données CPM-normalisées seront de 3,3 ×10 5 pour les trois gènes du patient A, de 5,0 × 10 5 pour les deux premiers gènes du patient B, et de 0 pour le dernier. Ainsi, la présence d’un gène différentiellement exprimé biaise les résultats de la normalisation CPM : les deux premiers gènes ont des niveaux d’expression CPM-normalisés différents, alors que nous avons supposé que leur niveau d’expression réel était les mêmes chez les deux patients.

La base de données TCGA

Description de la base de données

The Cancer Genome Atlas (TCGA) (https://www.cancer.gov/tcga) est une base de données publiques multicentriques qui provient de la collaboration entre deux instituts américains, le National Cancer Institute (NCI) et le National Human Genome Research Institute (NHGRI). Ce projet a été lancé en 2005 et a pour but d’obtenir et d’analyser le transcriptome et le génome de tissus tumoraux et de quelques tissus sains environnant la tumeur obtenus à partir de biopsies. A ce jour, plus de 20 000 échantillons provenant de plus de 11 000 patients et de 33 types de cancer ont été collectés, caractérisés et analysés.
Les données génomiques et transcriptomiques sont associées à des données cliniques et de survie (e.g. âge, sexe, taille de la tumeur, présence de métastases ou non, temps de suivi du patient après le diagnostic, décès observé au cours du suivi ou non). Différents évènements et temps de survie associée peuvent être étudiés à partir des données cliniques. LIU et collab. [2018b] ont analysé les données de survie TCGA et donnent des recommandations sur l’utilisation de différents évènements, dont la survie globale et la survie sans progression. Pour cela, une procédure permettant d’évaluer la qualité des données de survie a été mise en place, et des recommandations sont faites pour chaque cancer. L’analyse du transcriptome des tumeurs et des données de survie a permis l’émergence de nouveaux biomarqueurs et de nouvelles thérapies pour des utilisations cliniques [DUMBRAVA et M ERIC-B ERNSTAM, 2018].

Choix des cancers étudiés

Les analyses que nous allons présenter dans les chapitres suivants sont menées sur différents cancers. Nous avons décrit les noms des cancers et quelques caractéristiques des données dans le tableau 1.2. Dans la suite du manuscrit, nous utiliserons les acronymes de TCGA pour nommer les cancers, et nous nous référerons à ce tableau pour leur signification.
Le premier critère retenu pour le choix des cancers est le nombre de patients. Nous avons retiré les cancers pour lesquels les ARNm (resp. miRNA) ont été séquencés pour moins de 80 patients dans les études menées par TCGA. Ainsi, les cancers dont les acronymes sont DLBC, KICH, CHOL et UCS ne seront pas étudiés pour les ARNm (resp. DLBC, KICH, CHOL, GBM et UCS pour les miARN). Ensuite, LIU et collab. [2018b] recommandent de ne pas utiliser le mélanome cutané (SKCM – Skin Cutaneous Melanoma) car les données transcriptomiques de tumeurs primaires et de métastases sont fournies sans distinction dans TCGA. De plus, les auteurs suggèrent que les phéochromocytomes et les paragangliomes (PCPG – Pheochromocytoma and Paraganglioma) ne possèdent pas assez d’évènements et nécessitent des temps de suivi plus importants à la fois pour la survie globale et la survie sans progression. Toujours selon cette même étude, parmi les cancers restants, la survie globale est utilisée, sauf pour ceux dont l’acronyme est BRCA, LGG, PRAD, READ, TGCT, THCA et THYM pour lesquels la survie sans progression est préférée [LIU et collab., 2018b]. Après cette première étape de choix, 26 cancers sont retenus pour les ARNm et 25 pour les miARN.
TABLEAU 1.1 – Cancers retenus et données de survie utilisées.
Le cancer GBM n’est pas utilisé pour les données miRNA-seq, et est indiqué en gras dans le tableau. Les données de séquençage des miARN ne sont pas disponible pour ce cancer. « OS » : Overall Survival (i.e la survie globale est utilisée comme temps de survie). « PFI » : Progression Free Interval (i.e la survie sans récidive est utilisée comme temps de survie).

Modélisation des données de survie

Dans cette partie, nous nous attachons à décrire les principaux outils mathématiques permettant de modéliser les données de survie. Pour plus de détails, nous référons le lecteur au livre « The Statistical Analysis of Failure Time Data » [KALBFLEISCH et P RENTICE, 2011].

Les censures, singularité des données de survie

On parle de « censure » (ou donnée censurée) lorsque l’évènement étudié (i.e. décès du patient ou nouvel évènement associé à la tumeur) n’a pas lieu au cours du suivi. Une censure est typiquement observée lorsque le patient guéri, quitte l’étude, ou lorsqu’il change d’hôpital et sort de la cohorte étudiée. Le « temps de censure » est alors défini comme le temps entre le diagnostic et la censure. Les données de survie sont donc bidimensionnelles et contiennent une variable de temps (le temps de suivi) et une variable binaire indiquant si l’évènement est observé ou non (par convention, 1 si l’évènement est observé, 0 sinon).
Dans le cas d’une censure, le temps entre le diagnostic et l’évènement n’est pas connu, mais le fait que l’évènement n’a pas été observé durant le suivi du patient est une information importante qui doit être prise en compte dans les modèles. L’enjeu des modèles de survie réside donc dans la prise en compte de toute l’information présente dans ces données hétérogènes censurées et non censurées. L’estimateur de Kaplan-Meier et le modèle de Cox sont deux approches largement utilisées pour étudier des données de survie et prendre en compte les censures, et nous nous attacherons à les décrire dans les paragraphes suivant.

La p-valeur du modèle de Cox univarié

La p-valeur du modèle de Cox univarié permet de tester l’association entre les indices pronostiques du jeux de données test et les temps de survie. Le calcul de cette métrique consiste à apprendre un modèle de Cox univarié sur le jeu de données test où l’unique variable explicative est l’indice pronostique. L’hypothèse nulle du test est « H 0.
Le test du rapport de vraisemblance [BUSE, 1982] est préféré au test de Wald [KALBFLEISCH et P RENTICE, 2011] car il possède de meilleures propriétés lorsque la taille de l’échantillon est faible [KENDALL, 1946]. La statistique du test du rapport de vraisemblance s’exprime comme la différence de la log-vraisemblance du modèle après estimation du coefficient β et la log-vraisemblance du modèle nul : λ . Asymptotiquement, cette statistique suit une loi du χ2 à un degré de liberté. Ainsi, la p-valeur du modèle de Cox univarié permet de tester la plausibilité d’obtenir les indices pronostiques du jeu de données test sous l’hypothèse nulle qu’aucun lien n’existe entre les indices pronostiques (et donc le niveau d’expression des gènes) et la survie (i.e. β = 0). Elle peut donc être interprétée comme une mesure du niveau de corrélation entre les indices pronostiques et la survie, mesuré à travers un modèle linéaire généralisé.
La fonction coxph du package survival [THERNEAU, 2020] permet d’apprendre un modèle de Cox avec l’indice pronostique comme seule variable explicative de la survie, et de tester son association avec la survie par un test du rapport de vraisemblance.

La « malédiction de la grande dimension »

Définition

La « malédiction de la grande dimension » fait référence à l’apparition de phénomènes qui ont lieu lorsque le nombre p de variables excède le nombre n de patients. La dimension de l’espace dans lequel « vivent » les variables devient alors si grand que les données paraissent éparses et éloignées. Il devient alors difficile de dégager des généralités, et de nombreux algorithmes statistiques classiques sont mis à mal par un tel passage à l’échelle.
Dans le cadre des données d’expression génétiques, pour les sous-types de cancer les plus étudiés, la base de données TCGA contient typiquement 500 patients pour plus de 20 000 gènes. Dans ce contexte, les conséquences classiques de la malédiction de la grande dimension sont le manque de stabilité [J ARDILLIER et collab., 2018] et le surajustement [PAVLOU et collab., 2015]. Dans ce cas, le modèle « colle » trop aux données d’apprentissage (i.e. les estimateurs du modèle changent radicalement lorsque les données d’apprentissage changent), et il est difficile de généraliser les résultats sur un nouveau jeu de données.
Deux stratégies distinctes permettent de remédier à ce fléau de la grande dimension, et elles visent toutes deux à réduire le nombre de variables considérées :
— le pré-filtrage univarié des gènes.
— les méthodes de pénalisation.

Pré-filtrage univarié des gènes

Le pré-filtrage univarié des gènes consiste à assigner un score à chaque gène individuellement (e.g. variabilité du gène, test d’association entre l’expression du gène et la survie), et à sélectionner uniquement les gènes dont le score est inférieur (ou supérieur) à un certain seuil. Cette approche a l’avantage d’être intuitive et simple mathématiquement, mais ne tient pas compte des corrélations entre les gènes [BØVELSTAD et collab., 2007]. Des variables très corrélées ont des scores proches et donc une forte probabilité d’être toutes sélectionnées, ce qui entraîne une redondance d’information. En revanche, de telles approches peuvent être intéressantes à considérer en tant que première étape de filtrage afin de supprimer les gènes a priori non-informatifs (e.g. très faible variabilité, faible association individuelle avec la survie), avant d’inclure uniquement les variables retenues dans un modèle multivarié. BENNER et collab. [2010] ont monté sur données simulées que le pré-filtrage peut permettre d’obtenir de meilleures prédictions avec le modèle de Cox, et l’étude de l’impact du pré-filtrage sur la prédiction fera l’objet du chapitre3.

Choix du poids α accordé à la norme l 1 dans la pénalisation elastic net

Dans le paragraphe précédent consacré aux formes classiques de pénalisation, nous avons vu que les pénalisations elastic net et adaptive elastic net possède un hyperparamètre supplémentaire, α. Ce paramètre permet de répartir les poids assignés aux normes l 1 et l 2 dans la pénalisation. Plus α est élevé, plus le poids accordé à la norme l 1 est élevé, et plus celui accordé à la norme l 2 est faible. Ainsi, plus α est élevé, plus le nombre de gènes sélectionnés est faible.
Ce paramètre α peut-être fixé par l’utilisateur [J IANG et collab., 2016], ou fixé par validation croisée sur un vecteur de valeurs comprises entre 0 et 1. Dans ce dernier cas, tout comme pour le choix du poids λ accordé à la pénalisation, la déviance du modèle de Cox est classiquement retenue [MILANEZ -A LMEIDA et collab., 2020; OJEDA et collab., 2016].
Une double validation croisée menée à la fois sur λ et α est alors effectuée, et le processus de calcul est augmenté d’un facteur correspondant à la taille du vecteur des α testés.
Pour KIRC et avec la pénalisation elastic net, cette déviance reste relativement stable suivant les valeurs de α, et le nombre de gènes sélectionnés commence à se stabiliser à partir de 0,3 (Fig. 1.9). Cette tendance se retrouve pour l’ensemble des jeux de données étudiés (données non montrées). Ainsi, pour que les avantages de la pénalisation ridge soient pris significativement en compte, nous avons choisi α = 0,3 pour les pénalisations elastic net et adaptive elastic net pour l’ensemble des cancers.

Résumé

Ainsi, malgré les efforts de recherche importants mis en œuvre pour une meilleure compréhension biologique du cancer, la recherche de nouveaux traitements, et une meilleure prise en charge des patients, de nombreux enjeux demeurent.
Tout d’abord, le « principe d’unicité du cancer » (partie 1.1.3) rend difficile la prise en charge et le suivi des patients. En effet, la diversité du microenvironnement tumoral, des caractéristiques génomiques des cellules tumorales, et des types cellulaires obligent à un suivi individualisé.
Ensuite, nous avons vu que les molécules d’ARN caractérisent le niveau d’expression d’un gène, et jouent un rôle important dans l’agressivité des cancers (partie 1.2.4). Certains gènes s’expriment anormalement (i.e. sur-expression ou sous expression par rapport aux gènes de tissus sains), et le niveau d’expression des gènes peut permettre de classer les patients en terme de risque associé à la survie. La technologie « RNA-seq » permet de mesurer ces niveaux d’expression (partie 1.3). Le coût du séquençage a fortement diminué au cours des dix dernières années, mais il demeure encore trop élevé pour une utilisation routinière en clinique. Le coût et la qualité du séquençage sont croissants de la profondeur de séquençage.
Le modèle de Cox permet de modéliser les données de survie et de les lier à des covariables explicatives (e.g. cliniques, génomiques). Les données RNA-seq comportent environ 20 000 gènes et typiquement 500 patients (partie 1.5.1). Le nombre de variables explicatives (∼ 20 000) dépasse donc largement le nombre d’échantillons (∼ 500), et ce cas est appelé la « malédiction de la grande dimension ». Les données paraissent alors éparses et éloignées, et les propriétés mathématiques classiques sont mises à mal. Deux méthodes permettent de remédier à ce passage à l’échelle : le pré-filtrage univarié des gènes, et les méthodes de pénalisation de la pseudo-vraisemblance.

Table des matières
Liste des figures 
Liste des tableaux 
Contributions 
1 Introduction 
1.1 Cancer, recherche, et prédiction de la survie
1.2 Les molécules d’ARN
1.3 Le séquençage : mesure du niveau d’expression des gènes
1.4 Normalisation des données de comptage RNA-seq
1.5 La base de données TCGA
1.6 Modélisation des données de survie
1.7 Métriques d’évaluation de la qualité de prédiction
1.8 La « malédiction de la grande dimension »
1.9 Conclusions
1.10 Références
2 Comparaison et évaluation des méthodes de pénalisation du modèle de Cox avec les données mRNA-seq 
2.1 État de l’art et objectifs du chapitre
2.2 Comparaison des capacités de prédiction des méthodes de pénalisation du modèle de Cox sur données réelles
2.3 Combinaison des données cliniques et mRNA-seq pour prédire la survie
2.4 Procédure de simulation des données de survie et d’évaluation des pénalisations
2.5 Capacités de prédiction sur données simulées
2.6 Capacités de sélection des méthodes de pénalisation du modèle de Cox
2.7 Influence des performances de sélection sur la prédiction
2.8 Conclusions
2.9 Références
3 Pré-filtrage univarié des données d’expression génétiques et prédiction avec le modèle de Cox multivarié 
3.1 Pré-filtrage univarié des données d’expression génétique
3.2 Impact du pré-filtrage sur les capacités de prédiction du modèle de Cox pénalisé
3.3 Comparaison des méthodes de pénalisation du modèle de Cox après préfiltrage
3.4 Comparaison du pré-filtrage bi-dimensionnel à l’algorithme Iterative Sure Independance Screening
3.5 Conclusions et perspectives
3.6 Références
4 Impact de la profondeur de séquençage des données miRNA-seq sur la prédiction 
4.1 Contexte de l’étude
4.2 Choix des cancers étudiés
4.3 Comparaison des prédictions obtenues avec les ARNm et les miARN
4.4 Prédiction de la survie avec les variables cliniques et les données miRNA-seq
4.5 Dégradation des données de séquençage
4.6 Impact de la taille des banques et du nombre de patients sur la prédiction
4.7 Cause de la dégradation de la qualité de prédiction
4.8 Conclusions
4.9 Références
5 Conclusions et perspectives 
5.1 Conclusions générales
5.2 Résultats préliminaires et perspectives
5.3 Conclusions pratiques
5.4 Références
A Figures et Tableaux Annexes
A.1 Chapitre 2
A.2 Chapitre 3
A.3 Chapitre 4
B Liste des acronymes
C Glossaire
Résumé

Lire le rapport complet