- Email: [email protected]
Un estimateur de distributions ajusté pour une aide à la décision de la neutralité génétique des populations
Disponible en ligne sur
www.sciencedirect.com IRBM 33 (2012) 24–28
Article original
Un estimateur de distributions ajusté pour une aide à la décision de la neutralité génétique des populations A distributions adjusted estimator for decision support to neutral population genetics M. Troudi a,∗,b , L. Cherni c , A. Ben Ammar Gaaied c , M. Moalla d a
Institut des Hautes Études Commerciales de Carthage, Carthage Présidence-2016, Tunis, Tunisie b UR traitement du signal, de l’image et reconnaissance des formes, ENIT, Tunisie c Laboratoire de génétique moléculaire, faculté des sciences de Tunis, d’immunologie et de biotechnologie, Tunisie d Laboratoire d’informatique, de parallélisme et de productique (LIP2), faculté des sciences de Tunis, Tunisie Rec¸u le 15 aoˆut 2011 ; rec¸u sous la forme révisée le 12 d´ecembre 2011 ; accepté le 16 d´ecembre 2011 Disponible sur Internet le 1er février 2012
Résumé L’estimation de la neutralité d’une population est souvent réalisée à partir d’une seule mesure, basée sur des simulations de modèles aléatoires. Les cas litigieux posent une problématique quant à la prise de décision concernant la neutralité de la population étudiée. Nous proposons dans cet article d’apporter une aide à la décision de cette problématique en couplant un simulateur de populations génétiquement neutres à un estimateur non paramétrique ajusté et rapide des distributions simulées. Le choix du noyau d’Epanechnikov pour la méthode d’estimation non paramétrique du noyau permet une approche analytique du pas optimal, ce qui entraîne une amélioration de la complexité et par conséquent un gain en temps de calcul. L’estimation rapide et fiable des distributions simulées permet d’évaluer la crédibilité des décisions prises quant à la neutralité de ces populations. Ces résultats sont illustrés par deux populations tunisiennes dont la neutralité génétique n’est pas tranchée. © 2012 Elsevier Masson SAS. Tous droits réservés. Abstract The estimate of the neutrality of a population is often made from a single measurement based on simulations of random patterns. Disputed cases pose a problem regarding the decision-making about the neutrality of the study population. We propose in this paper to provide assistance to the decision of this problem by coupling a genetically neutral populations simulator to a fast adjusted non-parametric estimator of simulated distributions. The choice of the Epanechnikov kernel for the non-parametric method of estimating the kernel allows an analytical approach of the optimal bandwidth, which leads to an improvement in the complexity and therefore a gain in computing time. The fast and reliable estimation of simulated distributions allows to assess the credibility of the decisions about the neutrality of these populations. These results are illustrated by two Tunisian populations whose neutrality is not genetically determined. © 2012 Elsevier Masson SAS. All rights reserved.
1. Introduction La majorité des études réalisées en génétique des populations se basent sur l’hypothèse de la neutralité des populations étudiées. Cette notion est en rapport direct avec la diversité génétique des espèces définie par l’existence de plusieurs versions d’un même gène appelées allèles. Ces allèles assurent généralement les mêmes fonctions tout en présentant des variations
∗
Auteur correspondant. Adresse e-mail : [email protected] (M. Troudi).
1959-0318/$ – see front matter © 2012 Elsevier Masson SAS. Tous droits réservés. doi:10.1016/j.irbm.2011.12.003
structurales dues à des mutations ponctuelles sur certains nucléotides appelés sites polymorphes et notés S. L’apparition d’une mutation sur un site polymorphe est un phénomène naturellement rare. Par conséquent, les populations neutres contiennent un nombre relativement faible d’allèles jeunes caractérisés par leurs faibles fréquences dans la population étudiée. Il existe dans la littérature un nombre important des tests de neutralité parmi lesquels nous pouvons citer le test HKA, le test de Tajima, le test de Fu, etc. Ces tests donnent généralement des résultats quelquefois divergents en raison de leur conception même. Les puissances de ces tests ont fait l’objet de plusieurs études comparatives [1]. Le test statistique de Tajima [2] a été construit afin d’évaluer la neutralité génétique d’une
M. Troudi et al. / IRBM 33 (2012) 24–28
population en se basant sur l’estimation du paramètre mutationnel par deux méthodes différentes. La première méthode développée par Watterson [3] donne autant d’importance aux allèles de faible fréquence ou de fréquence élevée. En revanche, l’impact des allèles de faible fréquence est négligeable pour l’estimateur de Tajima [4]. La neutralité génétique est estimée à partir de la densité de probabilité de la distribution théorique d’une statistique D, déduite à partir de populations simulées construites de manière à garantir une neutralité théorique. L’estimation de cette densité de probabilité par la méthode non paramétrique de l’histogramme a permis d’approcher par la suite cette distribution par une loi Gamma. Nous proposons dans ce papier de remplacer l’approche paramétrique par un algorithme d’estimation des densités de probabilité par la méthode du noyau ajusté rapide [5] dans l’objectif d’apporter une aide à la décision quant à l’interprétation des résultats de neutralité litigieux. En effet, une estimation non paramétrique suffisamment précise et systématique garantit une approche plus fine de la neutralité génétique. Par ailleurs, il est parfois difficile de décider de la neutralité d’une population lorsque le résultat du test est proche du seuil de décision. Sachant que ce résultat est une variable aléatoire, l’estimation non paramétrique de sa distribution par des méthodes rapides permet de quantifier le taux de confiance du résultat obtenu. 2. Évaluation de la neutralité génétique L’évaluation de la neutralité des populations se base sur le modèle d’arbre aléatoire utilisé pour représenter une généalogie de gène [6]. En effet, les allèles sont supposés dériver d’un ancêtre commun suite à des mutations. L’espérance du nombre de mutations ayant eu lieu depuis la divergence de deux séquences est appelée paramètre mutationnel noté θ. Ainsi, θ = 2N × 2 × μ, μ étant le taux de mutations par génération et par site polymorphe et N la taille de la population. Dans la littérature, θ peut être estimé par plusieurs méthodes parmi lesquelles nous citons particulièrement les estimateurs de Watterson [3] et de Tajima [4]. 2.1. Estimateurs de Watterson et de Tajima
n−1 i=1
n’est autre que θ = 4Nμ pour chaque génération et pour chaque brin d’ADN. La moyenne du nombre de différences nucléoˆ tides des séquences d’ADN prises deux à deux, que l’on note π, permet donc d’avoir une estimation directe de θ. π est estimé soit à partir de la moyenne des différences nucléotides entre les séquences d’ADN prises deux à deux, soit en se basant sur les fréquences alléliques. On a alors : kij 2 kij s i
kij étant la différence nucléotide entre les séquences i et j, S est le nombre de sites polymorphes et hi l’estimation de la diversité nucléotide (hétérozygotie) pour le ie site polymorphe. 2.2. Test statistique de neutralité de Tajima 2.2.1. Principe Lorsqu’une population ne vérifie pas la neutralité, les valeurs estimées du paramètre mutationnel sont significativement différentes selon la méthode d’estimation utilisée. En effet, l’estimateur de Watterson, θˆ w , alloue le même poids aux différents allèles présents dans l’échantillon quelles que soient leurs ˆ tient compte fréquences. En revanche, l’estimateur de Tajima, π, des fréquences alléliques. La valeur du paramètre mutationnel est par conséquent peu affectée par la présence d’allèles de faible fréquence. Ainsi, l’écart entre θˆ w et πˆ est d’autant plus élevé que le nombre d’allèles rares de faible fréquence est important dans l’échantillon. Dans le cas d’une population neutre, les valeurs du paramètre mutationnel estimées par les deux méthodes sont sensiblement égales et leur différence fluctue autour de 0. D’où l’idée de développer une statistique basée sur la différence entre les deux estimateurs [2]. Soit d la différence de valeur entre les deux estimateurs de θ. − d=π
S a1
(3)
E (d) = 0 et Vˆ (d) = e1 S + e2 S (S − 1) avec e1
La méthode d’estimation de θ proposée par Watterson [3] se base sur le nombre de sites polymorphes S comptabilisés sur les échantillons d’ADN dans l’échantillon. L’espérance de S correspond au nombre total de mutations ayant pu se réaliser tout au long de la généalogie. Ainsi, θ est estimé à partir du nombre de sites polymorphes de l’échantillon : 1 1 S θw = et a2 = avec a1 = a1 i i2
25
=
N +1 1 − 2 3 (N − 1) a1 a1
et e2 =
2a12 N 2 + N + 3 − 9a1 (N − 1) (N + 2) + 9Na2 (N − 1)
9Na12 (N − 1) a12 + a2
Soit D = √ d
n−1
(1)
i=1
Tajima [4] a proposé d’estimer θ en basant sur l’observation qu’à chaque coalescence due à une mutation, un nouvel allèle apparaît, faisant passer le nombre d’allèles de i à i + 1. La modélisation mutationnelle adoptée étant le modèle à nombre infini de sites, le nombre de nucléotides différents entre deux séquences est la somme des mutations subies par l’ancêtre commun qui
(4)
Vˆ (d)
, une statistique centrée réduite représentant
la différence entre les deux estimateurs. Pour tester l’hypothèse de neutralité, la densité de probabilité de D est estimée pour des populations neutres construites par simulation à partir des paramètres de la population étudiée. 2.2.2. Méthodologie Une population est considérée comme présentant un écart à la neutralité lorsque la valeur de D, notée Dp , déduite à partir d’un
26
M. Troudi et al. / IRBM 33 (2012) 24–28
échantillon de la population étudiée, est considérée comme ayant une faible probabilité de faire partie de la distribution théorique de D issue des populations simulées. Ainsi, la population est supposée neutre lorsque la fonction de répartition de D pour Dp est supérieure à 0,02 pour un risque de première espèce de 5 % (F(Dp ) = P[D < Dp] > 0,02). Les différentes étapes permettant d’évaluer la neutralité pour un échantillon de taille N sont : • l’estimation de θ à partir de l’échantillon par l’estimateur de Tajima πˆ ; • la génération d’un nombre important (généralement 1000) de populations neutres simulées. Ces populations ont la même taille N de l’échantillon réel et le paramètre mutationnel πˆ estimé à partir de l’échantillon ; • la détermination de D pour chaque population neutre générée permet d’obtenir une distribution empirique de la variable aléatoire D ; • l’évaluation de la neutralité par la détermination de la valeur F(Dp ). Dans [2], Tajima estime les densités de probabilité de la statistique Dp par la méthode de l’histogramme et en déduit que cette statistique suit une loi bêta. Par conséquent, il propose de déterminer F(Dp ) de manière paramétrique en supposant que la statistique D suit une loi bêta dont les paramètres sont estimés à partir de l’échantillon. Cette hypothèse, largement utilisée dans les logiciels de base en génétique des populations, mérite d’être approfondie, notamment en évaluant une distance entre les densités de probabilités estimées et la loi bêta correspondante. Cela fera l’objet d’une prochaine publication où les densités de probabilités seront estimées en faisant appel à la méthode du noyau adaptée aux densités à supports bornés ou semi-bornés [8]. Des mesures de l’écart quadratique moyen entre l’estimation paramétrique et non paramétrique pourraient corroborer l’hypothèse émise par Tajima. Cette méthode d’évaluation de la neutralité est largement utilisée dans les logiciels de base en génétique des populations. Il nous paraît cependant intéressant d’avoir recours à des méthodes d’estimation non paramétriques et plus particulièrement à la méthode du noyau avec ajustement du pas pour une estimation plus précise de la densité de probabilité et par conséquent de la neutralité.
3. Estimation des densités de probabilité par l’algorithme du noyau ajusté rapide Plusieurs méthodes permettant l’estimation non paramétrique des densités de probabilité sont citées dans la littérature. Les principales sont la méthode de l’histogramme [7], la méthode du noyau [9], la méthode des fonctions orthogonales [10] et la méthode des k plus proches voisins [7]. Nous proposons d’utiliser dans ce papier la méthode du noyau avec pas ajusté par l’algorithme du plug-in rapide [5].
3.1. Méthode du noyau ajusté L’estimateur à noyau a été développé par Parzen en 1962 [9]. Il est défini par : N X − Xi 1 fˆ N (x) = (6) K NhN hN i=1
avec (X1 , X2 ,. . .. . ., XN ) N réalisations d’un échantillon et K une densité de probabilité appelée noyau. Ainsi, le noyau K est centré sur le point x dont on veut estimer l’image par f : il s’agira de sommer les contributions des différents xi en normalisant par l’entité hN appelée « paramètre de lissage » ou plus simplement « pas ». Les noyaux usuels sont essentiellement le noyau gaussien, le noyau uniforme et le noyau d’Epanechnikov ou noyau optimal. Le choix du pas hN est d’une importance primordiale pour une bonne estimation de f(x). En effet, l’obtention de la convergence en moyenne quadratique exige que NhN tende vers l’infini lorsque N tend vers l’infini. En ce qui concerne la convergence en moyenne quadratique intégrée, la condition exigée est plus restrictive puisque N(hN )2 doit tendre vers l’infini lorsque N tend vers l’infini. L’étude de la convergence de cet estimateur montre que la valeur optimale de hN notée h∗N est : 1
1
1
h∗N = N − 5 .(J (f ))− 5 .(M (K)) 5 avec M (K) =
−∞
=
+∞
K2 (u) du et J (f )
+∞
−∞
(7)
2
f (x)
dx
Ainsi, l’expression théorique du pas optimal au sens de l’EQMI est fonction de la taille de l’échantillon N, de l’intégrale du noyau choisi élevé au carré M(K) ainsi que de l’intégrale de la dérivée seconde élevée au carré de fˆ N notée J(f). Or, cette dernière entité est directement liée à f, fonction inconnue que nous cherchons à estimer. Les principales méthodes permettant de donner une estimation de la quantité J(f) sont les méthodes rule-of-thumb, les méthodes cross-validation et les algorithmes plug-in [7]. Nous avons retenu pour notre système, la méthode plug-in pour laquelle nous avons développé une variante rapide basée sur le noyau optimal [5]. 3.2. Méthode du plug-in rapide J(f) étant l’intégrale de la dérivée seconde élevée au carré de la densité à estimer f, il est possible de l’exprimer analytiquement lorsque le noyau utilisé est deux fois dérivable, ce qui est le cas du noyau optimal dont l’expression est : ⎧ √ ⎪ |x| > 5 si ⎨0 √ (8) K (x) = 3 x2 ⎪ |x| ≤ 5 1− si ⎩ √ 5 4 5
M. Troudi et al. / IRBM 33 (2012) 24–28
Sa dérivée seconde s’exprime par : ⎧ √ si |x| > 5 0 ⎪ ⎪ ⎪ √ ⎨ |x| = 5 K (x) = indé√fini si ⎪ ⎪ √ ⎪ ⎩ − 3 5 si |x| < 5 50
4. Présentation du système développé Le recueil et la lecture des données biologiques sont réalisés par les spécialistes en génétique des populations qui fournissent les entrées du système constituées de la taille de l’échantillon de la population étudiée notée N ainsi que du paramètre mutationnel selon Tajima noté π et de la valeur de Dp déduite de l’échantillon. Le système met en œuvre l’algorithme1 permettant de simuler 1000 populations théoriquement neutres de taille N et de paramètre mutationnel π et de déduire la valeur de D pour chacune des populations simulées. On obtient ainsi une distribution de la statistique D. Cette distribution est estimée par l’algorithme2 du noyau rapide ajusté et la valeur de F(Dp ) calculée numériquement. Si F(Dp ) est proche de 0,02, l’algorithme 1 est répété 100 fois afin d’obtenir une distribution de F(Dp ) et déduire la probabilité d’avoir une valeur de F(Dp ) inférieure à 0,02.
Ainsi, nous proposons d’approximer J(f) par : J (f ) =
+∞
−∞
1 dx = N 2 h6N
La quantité
1 K Nh3N i=1
+∞
−∞
N i=1
N
N
K
i=1
K
x − Xi hN
x − XI hN
x − X1 hN
2
2 dx
peut s’exprimer par
i ∈ AN (X1 .,...,XN )
K
x − Xi hN
avec AN un sous-ensemble d’entiers tel que |x − Xi | √ AN (X1 ., ..., XN ) = 0 ≤ i ≤ N; < 5 hN En remplac¸ant la dérivée seconde du noyau par sa valeur, ⎤ 2
N 2 x − X 9 i ⎣ ⎦ = K IAN (X1 .,...,XN ) (x) hN 500 ⎡
i ∈ AN (x)
i=1
=
27
9 β (x) 500
β(x) est une fonction constante par intervalles et formant une partition sur la droite réelle. Ainsi, J(f) est composé par la somme finie des dérivées secondes de la fonction noyau optimal et s’exprime par : +∞ 1 9 β (x) dx (9) Jˆ (f ) = 500 N 2 h6N −∞ Cette expression analytique de J(f) ne tient pas du tout compte des points indéfinis qui seront traités de la même manière √ que les points dont la valeur absolue est strictement supérieure à 5, c’est-à-dire que leur contribution sera considérée comme nulle. En effet, comme le nombre de points non définis pour β(x) est fini, on peut considérer la contribution de ces points dans J(f) comme négligeable. Cependant, les simulations réalisées nous ont permis de déterminer empiriquement la puissance de hN qui permettra d’avoir une convergence pour une valeur de 4,5. Ainsi, J(f) s’exprimera par : +∞ 1 9 β (x) dx (10) J (f ) = 500 N 2 h4.5 −∞ N
5. Application à des populations tunisiennes Le système ci-dessus présenté est illustré par l’estimation de la neutralité génétique de deux populations tunisiennes présentant des valeurs de F(Dp ) proches de 0,02 d’où la difficulté de décider de leurs neutralité à partir d’une seule observation (Fig. 1 et 2). En optimisant les temps de calcul grâce à l’algorithme du noyau ajusté rapide, non seulement les moyennes et la variance de F(Dp ) sont estimées de manière fiable, mais également la probabilité pour que cette valeur soit supérieure à 0,02. Les deux populations de Tunisie étudiées sont nommées P1 et P2 et leurs paramètres N, π et S ainsi que la valeur Dp mesurée sur l’échantillon sont présentés dans le Tableau 1. Les moyennes et variances estimées pour les deux populations sont présentées dans le Tableau 2. On observe qu’environ 80 % des tests donneront une valeur de F(Dp ) supérieure ou égale à 0,02 pour la population P1, ce qui permettra d’étayer la conclusion que cette population est plutôt neutre. En revanche, on constate qu’il est difficile de trancher
Fig. 1. Estimation de la distribution de F(Dp ) pour la population P1.
28
M. Troudi et al. / IRBM 33 (2012) 24–28
particulièrement pour les cas litigieux. Nous avons proposé dans ce papier d’apporter une amélioration à ces mesures en développant notre propre simulateur de populations génétiquement neutres couplé à un estimateur de densités de probabilité basé sur la méthode du noyau ajusté rapide. Ce système a permis de donner une estimation de la fréquence des mesures ce qui représente une aide à la décision intéressante pour conclure de manière plus fiable à la neutralité ou à l’absence de neutralité d’une population donnée. Les perspectives de ce travail sont nombreuses tant au niveau applicatif que fondamental. Ainsi, les applications directes concernent tous les travaux qui requièrent l’estimation de densités de probabilité univariées. L’étude du cas des densités à support borné ainsi que la généralisation vers l’estimation des densités de probabilité multivariées représentent des pistes intéressantes pouvant apporter des améliorations importantes pour des problématiques actuelles dans les domaines de la reconnaisse de formes et de la reconstitution 3D. Fig. 2. Estimation de la distribution de F(Dp ) pour la population P2.
Déclaration d’intérêts Tableau 1 Caractéristiques des populations P1 et P2.
Les auteurs n’ont pas transmis de déclaration de conflits d’intérêts.
Population
N
π
Dp
P1 P2
55 40
7,60471 6,33846
−1,717 −1,755
Tableau 2 Valeurs moyennes et variance de F(Dp ). Population
P1 P2
F(Dp ) Moyenne
Variance
P[F(Dp ) < 0,02]
0,0237 0,0208
(0,0048)2 (0,004)2
0,2027 0,5081
pour le cas de P2 puisque seuls 50 % des tests donneront une valeur de F(Dp ) supérieure ou égale à 0,02. 6. Discussion et conclusion Afin d’estimer la neutralité d’une population, les généticiens ont souvent recours à des logiciels fermés qui leur permettent de tirer des conclusions à partir d’une seule mesure. Or, cette mesure se basant sur des simulations de modèles aléatoires, il paraît imprudent de se prononcer à partir d’une seule observation
Références [1] Zhai W, Nielsen R, Slatkin M. An investigation of the statistical power of neutrality tests based on comparative and population genetic data. Mol Biol Evol 2009;26:273–83. [2] Tajima F. Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics 1989;123:585–95. [3] Watterson GA. On the number of segregating sites in genetical models without recombination. Theor Popul Biol 1974;7:256–76. [4] Tajima F. Evolutionary relationships of DNA sequences in finite populations. Genetics 1983;105:437–60. [5] Troudi M, Alimi AM, Saoudi S. Analytical plug-in method for kernel density estimator applied to genetic neutrality study. EURASIP J Adv Signal Process. Volume 2008, Article ID 379082, 8 pages, doi:10.1155/2008/379082. [6] Wright S. The distribution of gene frequencies in populations. Genetics 1937;23:307–20. [7] Delaigle A, Gijbels I. Practical bandwidth selection in deconvolution kernel density estimation. Comput Stat Data Anal 2004;45(2):249–67. [8] Saoudi S, Ghorbel F, Hillion A. Some statistical properties of the Kernel-diffeomorphism estimator. Appl Stochastic Models Data Anal 1997;10:39–58. [9] Parzen E. On estimation of a probability density function and mode. Ann Math Stat 1962;33:1065–76. [10] Zribi m, Ghorbel F. An unsupervised and non-parametric Bayesian classifier. Pattern Recogn Lett 2003;24:97–112.
www.sciencedirect.com IRBM 33 (2012) 24–28
Article original
Un estimateur de distributions ajusté pour une aide à la décision de la neutralité génétique des populations A distributions adjusted estimator for decision support to neutral population genetics M. Troudi a,∗,b , L. Cherni c , A. Ben Ammar Gaaied c , M. Moalla d a
Institut des Hautes Études Commerciales de Carthage, Carthage Présidence-2016, Tunis, Tunisie b UR traitement du signal, de l’image et reconnaissance des formes, ENIT, Tunisie c Laboratoire de génétique moléculaire, faculté des sciences de Tunis, d’immunologie et de biotechnologie, Tunisie d Laboratoire d’informatique, de parallélisme et de productique (LIP2), faculté des sciences de Tunis, Tunisie Rec¸u le 15 aoˆut 2011 ; rec¸u sous la forme révisée le 12 d´ecembre 2011 ; accepté le 16 d´ecembre 2011 Disponible sur Internet le 1er février 2012
Résumé L’estimation de la neutralité d’une population est souvent réalisée à partir d’une seule mesure, basée sur des simulations de modèles aléatoires. Les cas litigieux posent une problématique quant à la prise de décision concernant la neutralité de la population étudiée. Nous proposons dans cet article d’apporter une aide à la décision de cette problématique en couplant un simulateur de populations génétiquement neutres à un estimateur non paramétrique ajusté et rapide des distributions simulées. Le choix du noyau d’Epanechnikov pour la méthode d’estimation non paramétrique du noyau permet une approche analytique du pas optimal, ce qui entraîne une amélioration de la complexité et par conséquent un gain en temps de calcul. L’estimation rapide et fiable des distributions simulées permet d’évaluer la crédibilité des décisions prises quant à la neutralité de ces populations. Ces résultats sont illustrés par deux populations tunisiennes dont la neutralité génétique n’est pas tranchée. © 2012 Elsevier Masson SAS. Tous droits réservés. Abstract The estimate of the neutrality of a population is often made from a single measurement based on simulations of random patterns. Disputed cases pose a problem regarding the decision-making about the neutrality of the study population. We propose in this paper to provide assistance to the decision of this problem by coupling a genetically neutral populations simulator to a fast adjusted non-parametric estimator of simulated distributions. The choice of the Epanechnikov kernel for the non-parametric method of estimating the kernel allows an analytical approach of the optimal bandwidth, which leads to an improvement in the complexity and therefore a gain in computing time. The fast and reliable estimation of simulated distributions allows to assess the credibility of the decisions about the neutrality of these populations. These results are illustrated by two Tunisian populations whose neutrality is not genetically determined. © 2012 Elsevier Masson SAS. All rights reserved.
1. Introduction La majorité des études réalisées en génétique des populations se basent sur l’hypothèse de la neutralité des populations étudiées. Cette notion est en rapport direct avec la diversité génétique des espèces définie par l’existence de plusieurs versions d’un même gène appelées allèles. Ces allèles assurent généralement les mêmes fonctions tout en présentant des variations
∗
Auteur correspondant. Adresse e-mail : [email protected] (M. Troudi).
1959-0318/$ – see front matter © 2012 Elsevier Masson SAS. Tous droits réservés. doi:10.1016/j.irbm.2011.12.003
structurales dues à des mutations ponctuelles sur certains nucléotides appelés sites polymorphes et notés S. L’apparition d’une mutation sur un site polymorphe est un phénomène naturellement rare. Par conséquent, les populations neutres contiennent un nombre relativement faible d’allèles jeunes caractérisés par leurs faibles fréquences dans la population étudiée. Il existe dans la littérature un nombre important des tests de neutralité parmi lesquels nous pouvons citer le test HKA, le test de Tajima, le test de Fu, etc. Ces tests donnent généralement des résultats quelquefois divergents en raison de leur conception même. Les puissances de ces tests ont fait l’objet de plusieurs études comparatives [1]. Le test statistique de Tajima [2] a été construit afin d’évaluer la neutralité génétique d’une
M. Troudi et al. / IRBM 33 (2012) 24–28
population en se basant sur l’estimation du paramètre mutationnel par deux méthodes différentes. La première méthode développée par Watterson [3] donne autant d’importance aux allèles de faible fréquence ou de fréquence élevée. En revanche, l’impact des allèles de faible fréquence est négligeable pour l’estimateur de Tajima [4]. La neutralité génétique est estimée à partir de la densité de probabilité de la distribution théorique d’une statistique D, déduite à partir de populations simulées construites de manière à garantir une neutralité théorique. L’estimation de cette densité de probabilité par la méthode non paramétrique de l’histogramme a permis d’approcher par la suite cette distribution par une loi Gamma. Nous proposons dans ce papier de remplacer l’approche paramétrique par un algorithme d’estimation des densités de probabilité par la méthode du noyau ajusté rapide [5] dans l’objectif d’apporter une aide à la décision quant à l’interprétation des résultats de neutralité litigieux. En effet, une estimation non paramétrique suffisamment précise et systématique garantit une approche plus fine de la neutralité génétique. Par ailleurs, il est parfois difficile de décider de la neutralité d’une population lorsque le résultat du test est proche du seuil de décision. Sachant que ce résultat est une variable aléatoire, l’estimation non paramétrique de sa distribution par des méthodes rapides permet de quantifier le taux de confiance du résultat obtenu. 2. Évaluation de la neutralité génétique L’évaluation de la neutralité des populations se base sur le modèle d’arbre aléatoire utilisé pour représenter une généalogie de gène [6]. En effet, les allèles sont supposés dériver d’un ancêtre commun suite à des mutations. L’espérance du nombre de mutations ayant eu lieu depuis la divergence de deux séquences est appelée paramètre mutationnel noté θ. Ainsi, θ = 2N × 2 × μ, μ étant le taux de mutations par génération et par site polymorphe et N la taille de la population. Dans la littérature, θ peut être estimé par plusieurs méthodes parmi lesquelles nous citons particulièrement les estimateurs de Watterson [3] et de Tajima [4]. 2.1. Estimateurs de Watterson et de Tajima
n−1 i=1
n’est autre que θ = 4Nμ pour chaque génération et pour chaque brin d’ADN. La moyenne du nombre de différences nucléoˆ tides des séquences d’ADN prises deux à deux, que l’on note π, permet donc d’avoir une estimation directe de θ. π est estimé soit à partir de la moyenne des différences nucléotides entre les séquences d’ADN prises deux à deux, soit en se basant sur les fréquences alléliques. On a alors : kij 2 kij s i
kij étant la différence nucléotide entre les séquences i et j, S est le nombre de sites polymorphes et hi l’estimation de la diversité nucléotide (hétérozygotie) pour le ie site polymorphe. 2.2. Test statistique de neutralité de Tajima 2.2.1. Principe Lorsqu’une population ne vérifie pas la neutralité, les valeurs estimées du paramètre mutationnel sont significativement différentes selon la méthode d’estimation utilisée. En effet, l’estimateur de Watterson, θˆ w , alloue le même poids aux différents allèles présents dans l’échantillon quelles que soient leurs ˆ tient compte fréquences. En revanche, l’estimateur de Tajima, π, des fréquences alléliques. La valeur du paramètre mutationnel est par conséquent peu affectée par la présence d’allèles de faible fréquence. Ainsi, l’écart entre θˆ w et πˆ est d’autant plus élevé que le nombre d’allèles rares de faible fréquence est important dans l’échantillon. Dans le cas d’une population neutre, les valeurs du paramètre mutationnel estimées par les deux méthodes sont sensiblement égales et leur différence fluctue autour de 0. D’où l’idée de développer une statistique basée sur la différence entre les deux estimateurs [2]. Soit d la différence de valeur entre les deux estimateurs de θ. − d=π
S a1
(3)
E (d) = 0 et Vˆ (d) = e1 S + e2 S (S − 1) avec e1
La méthode d’estimation de θ proposée par Watterson [3] se base sur le nombre de sites polymorphes S comptabilisés sur les échantillons d’ADN dans l’échantillon. L’espérance de S correspond au nombre total de mutations ayant pu se réaliser tout au long de la généalogie. Ainsi, θ est estimé à partir du nombre de sites polymorphes de l’échantillon : 1 1 S θw = et a2 = avec a1 = a1 i i2
25
=
N +1 1 − 2 3 (N − 1) a1 a1
et e2 =
2a12 N 2 + N + 3 − 9a1 (N − 1) (N + 2) + 9Na2 (N − 1)
9Na12 (N − 1) a12 + a2
Soit D = √ d
n−1
(1)
i=1
Tajima [4] a proposé d’estimer θ en basant sur l’observation qu’à chaque coalescence due à une mutation, un nouvel allèle apparaît, faisant passer le nombre d’allèles de i à i + 1. La modélisation mutationnelle adoptée étant le modèle à nombre infini de sites, le nombre de nucléotides différents entre deux séquences est la somme des mutations subies par l’ancêtre commun qui
(4)
Vˆ (d)
, une statistique centrée réduite représentant
la différence entre les deux estimateurs. Pour tester l’hypothèse de neutralité, la densité de probabilité de D est estimée pour des populations neutres construites par simulation à partir des paramètres de la population étudiée. 2.2.2. Méthodologie Une population est considérée comme présentant un écart à la neutralité lorsque la valeur de D, notée Dp , déduite à partir d’un
26
M. Troudi et al. / IRBM 33 (2012) 24–28
échantillon de la population étudiée, est considérée comme ayant une faible probabilité de faire partie de la distribution théorique de D issue des populations simulées. Ainsi, la population est supposée neutre lorsque la fonction de répartition de D pour Dp est supérieure à 0,02 pour un risque de première espèce de 5 % (F(Dp ) = P[D < Dp] > 0,02). Les différentes étapes permettant d’évaluer la neutralité pour un échantillon de taille N sont : • l’estimation de θ à partir de l’échantillon par l’estimateur de Tajima πˆ ; • la génération d’un nombre important (généralement 1000) de populations neutres simulées. Ces populations ont la même taille N de l’échantillon réel et le paramètre mutationnel πˆ estimé à partir de l’échantillon ; • la détermination de D pour chaque population neutre générée permet d’obtenir une distribution empirique de la variable aléatoire D ; • l’évaluation de la neutralité par la détermination de la valeur F(Dp ). Dans [2], Tajima estime les densités de probabilité de la statistique Dp par la méthode de l’histogramme et en déduit que cette statistique suit une loi bêta. Par conséquent, il propose de déterminer F(Dp ) de manière paramétrique en supposant que la statistique D suit une loi bêta dont les paramètres sont estimés à partir de l’échantillon. Cette hypothèse, largement utilisée dans les logiciels de base en génétique des populations, mérite d’être approfondie, notamment en évaluant une distance entre les densités de probabilités estimées et la loi bêta correspondante. Cela fera l’objet d’une prochaine publication où les densités de probabilités seront estimées en faisant appel à la méthode du noyau adaptée aux densités à supports bornés ou semi-bornés [8]. Des mesures de l’écart quadratique moyen entre l’estimation paramétrique et non paramétrique pourraient corroborer l’hypothèse émise par Tajima. Cette méthode d’évaluation de la neutralité est largement utilisée dans les logiciels de base en génétique des populations. Il nous paraît cependant intéressant d’avoir recours à des méthodes d’estimation non paramétriques et plus particulièrement à la méthode du noyau avec ajustement du pas pour une estimation plus précise de la densité de probabilité et par conséquent de la neutralité.
3. Estimation des densités de probabilité par l’algorithme du noyau ajusté rapide Plusieurs méthodes permettant l’estimation non paramétrique des densités de probabilité sont citées dans la littérature. Les principales sont la méthode de l’histogramme [7], la méthode du noyau [9], la méthode des fonctions orthogonales [10] et la méthode des k plus proches voisins [7]. Nous proposons d’utiliser dans ce papier la méthode du noyau avec pas ajusté par l’algorithme du plug-in rapide [5].
3.1. Méthode du noyau ajusté L’estimateur à noyau a été développé par Parzen en 1962 [9]. Il est défini par : N X − Xi 1 fˆ N (x) = (6) K NhN hN i=1
avec (X1 , X2 ,. . .. . ., XN ) N réalisations d’un échantillon et K une densité de probabilité appelée noyau. Ainsi, le noyau K est centré sur le point x dont on veut estimer l’image par f : il s’agira de sommer les contributions des différents xi en normalisant par l’entité hN appelée « paramètre de lissage » ou plus simplement « pas ». Les noyaux usuels sont essentiellement le noyau gaussien, le noyau uniforme et le noyau d’Epanechnikov ou noyau optimal. Le choix du pas hN est d’une importance primordiale pour une bonne estimation de f(x). En effet, l’obtention de la convergence en moyenne quadratique exige que NhN tende vers l’infini lorsque N tend vers l’infini. En ce qui concerne la convergence en moyenne quadratique intégrée, la condition exigée est plus restrictive puisque N(hN )2 doit tendre vers l’infini lorsque N tend vers l’infini. L’étude de la convergence de cet estimateur montre que la valeur optimale de hN notée h∗N est : 1
1
1
h∗N = N − 5 .(J (f ))− 5 .(M (K)) 5 avec M (K) =
−∞
=
+∞
K2 (u) du et J (f )
+∞
−∞
(7)
2
f (x)
dx
Ainsi, l’expression théorique du pas optimal au sens de l’EQMI est fonction de la taille de l’échantillon N, de l’intégrale du noyau choisi élevé au carré M(K) ainsi que de l’intégrale de la dérivée seconde élevée au carré de fˆ N notée J(f). Or, cette dernière entité est directement liée à f, fonction inconnue que nous cherchons à estimer. Les principales méthodes permettant de donner une estimation de la quantité J(f) sont les méthodes rule-of-thumb, les méthodes cross-validation et les algorithmes plug-in [7]. Nous avons retenu pour notre système, la méthode plug-in pour laquelle nous avons développé une variante rapide basée sur le noyau optimal [5]. 3.2. Méthode du plug-in rapide J(f) étant l’intégrale de la dérivée seconde élevée au carré de la densité à estimer f, il est possible de l’exprimer analytiquement lorsque le noyau utilisé est deux fois dérivable, ce qui est le cas du noyau optimal dont l’expression est : ⎧ √ ⎪ |x| > 5 si ⎨0 √ (8) K (x) = 3 x2 ⎪ |x| ≤ 5 1− si ⎩ √ 5 4 5
M. Troudi et al. / IRBM 33 (2012) 24–28
Sa dérivée seconde s’exprime par : ⎧ √ si |x| > 5 0 ⎪ ⎪ ⎪ √ ⎨ |x| = 5 K (x) = indé√fini si ⎪ ⎪ √ ⎪ ⎩ − 3 5 si |x| < 5 50
4. Présentation du système développé Le recueil et la lecture des données biologiques sont réalisés par les spécialistes en génétique des populations qui fournissent les entrées du système constituées de la taille de l’échantillon de la population étudiée notée N ainsi que du paramètre mutationnel selon Tajima noté π et de la valeur de Dp déduite de l’échantillon. Le système met en œuvre l’algorithme1 permettant de simuler 1000 populations théoriquement neutres de taille N et de paramètre mutationnel π et de déduire la valeur de D pour chacune des populations simulées. On obtient ainsi une distribution de la statistique D. Cette distribution est estimée par l’algorithme2 du noyau rapide ajusté et la valeur de F(Dp ) calculée numériquement. Si F(Dp ) est proche de 0,02, l’algorithme 1 est répété 100 fois afin d’obtenir une distribution de F(Dp ) et déduire la probabilité d’avoir une valeur de F(Dp ) inférieure à 0,02.
Ainsi, nous proposons d’approximer J(f) par : J (f ) =
+∞
−∞
1 dx = N 2 h6N
La quantité
1 K Nh3N i=1
+∞
−∞
N i=1
N
N
K
i=1
K
x − Xi hN
x − XI hN
x − X1 hN
2
2 dx
peut s’exprimer par
i ∈ AN (X1 .,...,XN )
K
x − Xi hN
avec AN un sous-ensemble d’entiers tel que |x − Xi | √ AN (X1 ., ..., XN ) = 0 ≤ i ≤ N; < 5 hN En remplac¸ant la dérivée seconde du noyau par sa valeur, ⎤ 2
N 2 x − X 9 i ⎣ ⎦ = K IAN (X1 .,...,XN ) (x) hN 500 ⎡
i ∈ AN (x)
i=1
=
27
9 β (x) 500
β(x) est une fonction constante par intervalles et formant une partition sur la droite réelle. Ainsi, J(f) est composé par la somme finie des dérivées secondes de la fonction noyau optimal et s’exprime par : +∞ 1 9 β (x) dx (9) Jˆ (f ) = 500 N 2 h6N −∞ Cette expression analytique de J(f) ne tient pas du tout compte des points indéfinis qui seront traités de la même manière √ que les points dont la valeur absolue est strictement supérieure à 5, c’est-à-dire que leur contribution sera considérée comme nulle. En effet, comme le nombre de points non définis pour β(x) est fini, on peut considérer la contribution de ces points dans J(f) comme négligeable. Cependant, les simulations réalisées nous ont permis de déterminer empiriquement la puissance de hN qui permettra d’avoir une convergence pour une valeur de 4,5. Ainsi, J(f) s’exprimera par : +∞ 1 9 β (x) dx (10) J (f ) = 500 N 2 h4.5 −∞ N
5. Application à des populations tunisiennes Le système ci-dessus présenté est illustré par l’estimation de la neutralité génétique de deux populations tunisiennes présentant des valeurs de F(Dp ) proches de 0,02 d’où la difficulté de décider de leurs neutralité à partir d’une seule observation (Fig. 1 et 2). En optimisant les temps de calcul grâce à l’algorithme du noyau ajusté rapide, non seulement les moyennes et la variance de F(Dp ) sont estimées de manière fiable, mais également la probabilité pour que cette valeur soit supérieure à 0,02. Les deux populations de Tunisie étudiées sont nommées P1 et P2 et leurs paramètres N, π et S ainsi que la valeur Dp mesurée sur l’échantillon sont présentés dans le Tableau 1. Les moyennes et variances estimées pour les deux populations sont présentées dans le Tableau 2. On observe qu’environ 80 % des tests donneront une valeur de F(Dp ) supérieure ou égale à 0,02 pour la population P1, ce qui permettra d’étayer la conclusion que cette population est plutôt neutre. En revanche, on constate qu’il est difficile de trancher
Fig. 1. Estimation de la distribution de F(Dp ) pour la population P1.
28
M. Troudi et al. / IRBM 33 (2012) 24–28
particulièrement pour les cas litigieux. Nous avons proposé dans ce papier d’apporter une amélioration à ces mesures en développant notre propre simulateur de populations génétiquement neutres couplé à un estimateur de densités de probabilité basé sur la méthode du noyau ajusté rapide. Ce système a permis de donner une estimation de la fréquence des mesures ce qui représente une aide à la décision intéressante pour conclure de manière plus fiable à la neutralité ou à l’absence de neutralité d’une population donnée. Les perspectives de ce travail sont nombreuses tant au niveau applicatif que fondamental. Ainsi, les applications directes concernent tous les travaux qui requièrent l’estimation de densités de probabilité univariées. L’étude du cas des densités à support borné ainsi que la généralisation vers l’estimation des densités de probabilité multivariées représentent des pistes intéressantes pouvant apporter des améliorations importantes pour des problématiques actuelles dans les domaines de la reconnaisse de formes et de la reconstitution 3D. Fig. 2. Estimation de la distribution de F(Dp ) pour la population P2.
Déclaration d’intérêts Tableau 1 Caractéristiques des populations P1 et P2.
Les auteurs n’ont pas transmis de déclaration de conflits d’intérêts.
Population
N
π
Dp
P1 P2
55 40
7,60471 6,33846
−1,717 −1,755
Tableau 2 Valeurs moyennes et variance de F(Dp ). Population
P1 P2
F(Dp ) Moyenne
Variance
P[F(Dp ) < 0,02]
0,0237 0,0208
(0,0048)2 (0,004)2
0,2027 0,5081
pour le cas de P2 puisque seuls 50 % des tests donneront une valeur de F(Dp ) supérieure ou égale à 0,02. 6. Discussion et conclusion Afin d’estimer la neutralité d’une population, les généticiens ont souvent recours à des logiciels fermés qui leur permettent de tirer des conclusions à partir d’une seule mesure. Or, cette mesure se basant sur des simulations de modèles aléatoires, il paraît imprudent de se prononcer à partir d’une seule observation
Références [1] Zhai W, Nielsen R, Slatkin M. An investigation of the statistical power of neutrality tests based on comparative and population genetic data. Mol Biol Evol 2009;26:273–83. [2] Tajima F. Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics 1989;123:585–95. [3] Watterson GA. On the number of segregating sites in genetical models without recombination. Theor Popul Biol 1974;7:256–76. [4] Tajima F. Evolutionary relationships of DNA sequences in finite populations. Genetics 1983;105:437–60. [5] Troudi M, Alimi AM, Saoudi S. Analytical plug-in method for kernel density estimator applied to genetic neutrality study. EURASIP J Adv Signal Process. Volume 2008, Article ID 379082, 8 pages, doi:10.1155/2008/379082. [6] Wright S. The distribution of gene frequencies in populations. Genetics 1937;23:307–20. [7] Delaigle A, Gijbels I. Practical bandwidth selection in deconvolution kernel density estimation. Comput Stat Data Anal 2004;45(2):249–67. [8] Saoudi S, Ghorbel F, Hillion A. Some statistical properties of the Kernel-diffeomorphism estimator. Appl Stochastic Models Data Anal 1997;10:39–58. [9] Parzen E. On estimation of a probability density function and mode. Ann Math Stat 1962;33:1065–76. [10] Zribi m, Ghorbel F. An unsupervised and non-parametric Bayesian classifier. Pattern Recogn Lett 2003;24:97–112.