---
title: "Description de la variabilité pharmacogénétique dans la population française à l'aide de méthodes bioinformatiques"
author: [Sacha Schutz]
date: "2020-03-17"
bibliography: export.bib
csl: citation.csl
link-citations: true
subject: ""
keywords: [Markdown, Example]
subtitle: "Sous la direction d'Emmanuelle Genin et Tristan Montier"
lang: "fr"
titlepage: true,
titlepage-rule-color: "360049"
titlepage-background: "background1.pdf"
table-use-row-colors: true
...

\pagebreak

# Introduction
Mon projet de thèse vise à étudier l'impact de la variabilité génétique dans la réponse aux médicaments. L'objectif est de décrire la diversité génétique qui existe dans les gènes impliqués dans cette réponse dans la population française. Pour ce faire, je dispose de données de séquençage, exomes et génomes provenant du projet Frex et FranceGenRef.
Différentes études ont déjà été réalisées pour décrire cette diversité à partir de panels d'individus provenant des différentes populations à travers le monde. Peu d'études se sont intéressées à cette diversité à des échelles géographiques plus fines comme nous l'envisageons dans ce projet.

## Les gènes impliqués dans la réponse aux médicaments
De nombreux gènes sont impliqués dans la réponse aux médicaments. Ces gènes peuvent agir de différentes manières, en influençant la métabolisation du médicament (pharmacinétique) ou étant la cible des principes actifs (pharmacodynamique).
Parmi les gènes influençant la cinétique, les principaux étudiés sont ceux impliqués dans la biotransformation du médicament, le plus souvent de nature hépatique. Il s'agit en l’occurrence :

* **des réactions enzymatiques de phase 1**
qui sont des réactions d'oxydation réalisées par une famille d'enzymes appelées *Cytochrome P450* et dont les principaux représentants sont le *CPY3A4* et le *CYP2D6*. Ces pharmacogènes sont impliqués dans le métabolisme de nombreux médicaments comme les psychotropes et le tamoxifène.
* **Les réactions enzymatiques de phase 2** sont des réactions de conjugaison de différents radicaux. Par exemple, la famille des *UGT* permet l'ajout d'un acide glucuronique à des composés liposolubles augmentant ainsi la solubilité du médicament.

La connaissance de la variabilité de ces gènes dans la population permet d'ajuster les stratégies thérapeutiques.

Pour cette première année de thèse, j'ai orienté mes recherches sur l'identification des variations du gène codant pour le cytochrome CYP2D6 à partir des données de séquençage haut débit en utilisant et en développant des méthodes bioinformatiques.
Ces méthodes ont été utilisées sur les données de séquençage de la population française (FREX et FranceGenRef) afin de décrire la diversité génétique.

## Le gène CYP2D6
Le gène CYP2D6 est un gène complexe sur bien des aspects rendant difficile son analyse en routine. Il est accompagné de deux pseudogènes paralogues (CYP2D7 et CYP2D8) présentant un degré de similarité élevé. Cette structure rend l'alignement des reads et l'interprétation des variants particulièrement ardus. De plus, ce gène est polymorphe avec de nombreuses variations décrites, notamment des variations du nombre de copies ainsi que des réarrangements complexes avec ses pseudo-gènes. La base de données [PharmVar](https://www.pharmvar.org/) répertorie l'ensemble des allèles connus de ce gène. Elle adopte une nomenclature spécifique pour chaque haplotype (star allèle) en considérant aussi bien les SNPs, les variations du nombre de copie et les réarangements complexes. Par exemple, l'allèle CYP2D6\*14 est défini par 3 SNPs majeurs (1759G>A, 2851C>T, 4181G>C) dont l'effet est une diminution de la fonction du gène. L'allèle CYP2D6\*5 correspond à la délétion totale du gène et donc une perte de fonction. L'allèle CYP2D6\*2x2 correspond à une duplication de l'haplotype \*2.
A ce jour, plus de 145 allèles ont été répertoriés dans PharmVar.
Une fois ces haplotypes identifiés, il faut encore prédire le phénotype. Des scores empiriques existent permettant de classer le patient en 4 groupes : métaboliseur lent (PM), intermédiaire (IM), normal (EM) et ultra rapide (UM). Il reste cependant de nombreux variants du gène dont l'effet n'est pas connu.

# Outils et données à disposition
## Etats de l'art
Plusieurs outils bioinformatiques permettent l'identification des haplotypes à partir des données de séquençage haut-débit.

- **Pharmcat** est un outil prenant en entrée un fichier VCF et génère un rapport HTML avec les différents star allèles détéctés. Cependant, il ne tient pas compte des variations du nombre de copies.

- **Stargazer** est un pipeline qui infère les haplotypes avec une approche populationnelle en utilisant *Beagle* et 2500 références du projet 1000Genome. Les CNVs sont détéctés avec *GATK depthCoverage* qui est aujourd'hui obsolète. Par ailleurs l'approche populationnelle pour l'inférence des haplotypes est discutable avec des variants rares.

- **Aldy** infère les haplotypes par une approche combinatoire en utilisant de la programmation par contrainte avec la librarie OR-tools. Les CNVs sont détectés en comparant les variations de profondeurs de séquençage à celle du CYP2D8. L'outil depuis sa version 3.0 supporte les génomes hg19 et hg38.

- **Cyrius**, publié par illumina en Janvier 2021, propose une approche similaire en améliorant la détection des réarangements complexes entre le CYP2D6 et CYP2D7.

Ce sont ces deux derniers outils que nous avons retenus pour l'analyse.

## Données françaises à disposition
Au sein de l'unité, deux jeux de données sont à disposition.

- Le projet FREX (French Exome Project) est constitué de 574 exomes provenant de 6 régions françaises.
- Le Projet France GenRef est un projet pilote constitué de 852 génomes provenant de la cohorte GAZEL et PREGO.

# Distribution des haplotypes du CYP2D6 dans la population française

La distribution des haplotypes a été calculée sur les données du projet FranceGenRef avec Aldy. Notre laboratoire n'étant pas hébergeur des données brutes, je n'ai pas pu récupérer l'intégralité des génomes. Mais seulement un sous ensemble des reads mappant sur des gènes d'intérêt en pharmacogénétique, dont le CYP2D6 et ses pseudo-gènes. Ces données étaient suffisantes pour Aldy, mais pas pour Cyrius qui a besoin des génomes complets. Aldy n'a pas fonctionné sur les exomes car la référence utilisée pour détécter les CNVs est un pseudo gène qui n'a pas été séquencé.

## Distribution des allèles


\pagebreak
Les 5 principaux allèles sont:

- Allele *1 : Allèle de Référence - 32%
- Allele *4 : 1 snp délétère avec un effet sur l’épissage - 18%
- Allele *2 : 1 snp avec effet Normal - 15%
- Allele *41 : Combinaison de 3 snps avec diminution de fonction - 8%
- Allele *68 : gène hybride CYP2D6-CYP2D7 - 5%

La présence de l'allèle \*68 est très intéressante. En effet il s'agit d'un réarrangement complexe avec diminution de l'activité. Ce variant n'est probablement pas détecté en diagnostic clinique par les méthodes usuels. La prédiction du phénotype peut donc en être altéré.

## Comparaison des fréquences selon l’origine géographique

Le CPIC fournit une base de données avec les fréquences des haplotypes obtenus dans différentes populations. J'ai comparé les fréquences observées dans FranceGenRef avec cette base.
Dans l'ensemble, les fréquences des haplotypes observées se rapprochent de celles de la population européenne sauf pour l'haplotype \*1 et \*2.


## Distribution des phénotypes

Le phénotype a été prédit à partir des diplotypes obtenus à l'aide d'un score défini par le [CPIC](https://cpicpgx.org/) (Clinical Pharmacogenetics Implementation Consortium.)
La majorité de la population est prédite métaboliseur normal ou intermédiaire. Le nombre de métaboliseurs lents et ultra rapides n'est pas négligeable. Pour ces derniers, l'adaptation thérapeutique est essentielle.


# Iscard : Un outil pour la détéction des CNVs
Aldy et Cyrius détectent les CNVs en comparant les profondeurs de séquençage à une référence interne. Si un grand réarrangement emporte la référence, le résultat peut être faussé. Par ailleurs aucun de ces outils ne détaillent la localisation précise du CNV. Pour y remedier, j'ai crée un outil en python appelé "Iscard" disponible à cette adresse: https://iscard.readthedocs.io/en/latest/quickstart.html

## Apprentissage suppervisé
A partir d'un set de témoins portant 2 copies du gène, Iscard calcul la profondeur de séquençage sur des intervalles fixes puis les normalise.
2 modèles sont ensuite calculés :

- **Le modèle inter-sample**, Pour chaque intervalle, la moyenne et l’écart-type est calculé. Nous pouvons alors mesurer une déviation de la profondeur sur un nouvel échantillon avec z-score.
- **Le modèle intra-sample**, cherche sur le jeu d'entraînement des régions situées sur différents chromosomes qui sont linéairement corrélées entre elles. Nous pouvons alors calculer sur un nouvel échantillon les profondeurs des deux régions afin de tester si elles s'écartent du modèle linéaire.


Iscard s'utilise en ligne de commande de la façon suivante. La documentation officielle est disponible à cette adresse [https://iscard.readthedocs.io/en/latest/quickstart.html](https://iscard.readthedocs.io/en/latest/quickstart.html)

# Création du modèle
iscard learn -i training/*.bam -r region.bed -o model.h5

# Test d'un échantillon
iscard test -i sample.bam -m model.h5 -o sample.result

# Affichage des résultats
iscard plot -i sample.result -m model.h5 -o sample.png

#### Résultats

Iscard a besoin de faire un apprentissage sur plusieurs fichiers BAM sans duplication ni délétion. N'ayant pas ces données pour le CYP2D6, j'ai réalisé mes tests avec les données du séquençage du gène PKD1 effectué dans mon laboratoire.
Le modèle a été construit à partir de 120 fichiers BAMs issus d'un séquençage en capture Roche (13 gènes) sur Illumina NextSeq. J'ai ensuite testé 10 échantillons connus pour porter des CNVs identifiés par MLPA (Multiplex Ligation-dependent Probe Amplification).
Sur l'ensemble des échantillons, tous les CNVs détectés par MLPA ont pu être identifiés visuellement par Iscard. 3 exemples sont représentés sur la figure suivante:


J'ai tout de même réalisé une analyse sur des données CYP2D6 provenant de 3 contrôles qualités séquencés en amplicon sur Séquenceur IonTorrent. L'entraînement a été réalisé sur seulement 5 contrôles à deux copies. Les profondeurs des 3 échantillons sont représentées sur la figure suivante:


La délétion attendue est visible sur le graphique. Cependant le signal pour la duplication est difficilement interprétable. La technologie de séquençage utilisée et le faible nombre de témoins dans le jeu d'apprentissage ne permettent pas de conclure de façon fiable.

## Borne d'un CNV par programmation probabiliste
Aldy et Cyrius ne donnent pas d'information sur la localisation exacte du CNV. Pour identifier les bornes du CNVs, je me suis intéressé à un nouveau paradigme de programmation basé sur l'inférence bayésienne: la programmation probabiliste.
En deux mots, cette méthode permet d'estimer l'incertitude sur une variable à l'aide d'une distribution de probabilité ajustée à partir des données observées.
Cela consiste à décrire un modèle probabiliste et faire tourner un échantillonneur MCMC pour trouver la distribution a posteriori de la variable à estimer.
Dans le cas du bornage du CNVs, j'ai cherché la distribution de probabilité de la position de début et de fin du CNV. J'ai réalisé cette analyse à l'aide de la librairie Python [PyMC3](https://docs.pymc.io/) en modélisant de la façon suivante:

#### Formalisation du problème
Posons pour l'exemple le bornage d'une duplication. Ce problème correspond à l'analyse d'un signal unidimensionnel où il faut trouver le début $x_{1}$ et la fin $x_{2}$ de la déviation.
Appelons X ce signal défini sur l'intervalle [0,1000]. Cette variable suit une loi normale de moyenne $mu_{1}$ à l'extérieur de la déviation et $mu_{2}$ à l'intérieur de la déviation.

Le problème peut être formalisé de la façon suivante :

$$ x_{1} \sim uniform(0,1000)$$
$$x_{2} \sim uniform(0,1000)\\$$
$$\mu_{1} \sim Normal(0,1)\\$$
$$\mu_{2} \sim Normal(0,1)\\ $$
$$x \sim Normal(\mu,0.1) \\$$

$$
\mu = s, avec : s = \left\{\begin{matrix}
\mu_{1} \; si \; x < x1 \; ou \; x > x2
\\
\mu_{2} \; si \; x \geq x1 \; et \; x \leq x2
\end{matrix}\right.
$$

\pagebreak

Le modèle est implémenté en python comme suit :

with pm.Model() as model:
    
    mu1 = pm.Normal("mu_normal")
    mu2 = pm.Normal("mu_error")
    
    x1 =  pm.DiscretUniform("x1",lower=0, upper=1000)
    x2 =  pm.DiscretUniform("x2",lower=0, upper=1000)
    
    mu = pm.math.switch(pm.math.or_(x1 >= x, x2 <= x), mu1, mu2)
    
    x = pm.Normal("observed",mu = mu, signal = 0.1, observed=)

#### Résultats

La programmation bayésienne est très puissante lorsqu'il faut estimer de façon probabiliste des variables à partir de donnée observée parfois incomplètes. Ces méthodes sont par exemple utilisées en astrophysique pour estimer les caractéristiques d'une exoplanète à partir du peu de données collectées. Elle reste cependant lourde en temps de calcul et n'est pas justifiée si un calcul direct plus rapide est possible. Il serait interessant d'utiliser cette méthode pour inférer la probabilité des haplotypes et des phénotypes à partir des données de séquençages haut débit.

# Cutevariant: Logiciel d'interprétation des variants

Je développe également un logiciel d'interprétation des variants de séquençage haut débit avec un étudiant en Master 2. Une publication est en soumission. Il est prévu d'intégrer des extensions dans ce logiciel spécialisées aux questions de pharmacogénétique.


# Perspectives

- **Décrire les autres pharmacogènes**
Seul le CYP2D6 a été présenté ici. D'autres gènes sont en cours d'analyse dont les principaux gènes actionnables. Des premiers résultats nous ont montré par exemple une grande variabilité du gène DPYD impliqué dans la toxicité du 5-FU.

- **Tester l'apport du séquençage de 3ème génération** ( PacBio / Nanopore)
Les technologies de séquençage en long read permettraient de résoudre les haplotypes sans avoir à les inférer depuis les génotypes. Les réanragements complexes pourraient également être mieux détectés.

- **Tester les prédictions phénotypes par machine learning**
Une trentaine d'individus de France GenRef n'as pas été associé à un phénotype de métaboliseur alors que l'haplotype était connu. Il s'agit d'haplotype rare dont l'effet nous est inconnu. Des méthodes d'intelligence artificielle ont été publiées récemment s'aidant de réseau de neurones à convolution. (https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1008399)