Bases

1Intuitivement, on sait qu’un échantillon ne sera jamais une représentation absolument fidèle de ce qui existe dans la population entière. Il en découle que l’échantillonnage provoque une déviation des estimateurs que l’on peut calculer (variance d’échantillonnage). On aura donc rarement, même dans une population échantillonnée strictement conforme à Hardy-Weinberg, un F_IS estimé exactement égal à 0, pareil pour le F_ST, pour les déséquilibres de liaison, l’isolement par la distance ou n’importe quel autre paramètre. Le test statistique est là pour nous aider à prendre une décision quant à la disparité observée entre les données et l’attendu. La différence observée peut-elle être expliquée par le hasard et avec quelle probabilité ? Le but d’un test statistique sera donc de fournir une réponse à cette question, en donnant un critère, la valeur P ou P-value1 du test, ou risque de première espèce ou encore probabilité de se tromper en répondant par la négative à cette question (appelée hypothèse nulle). On peut ajouter ici qu’en génétique des populations, la variance d’échantillonnage est d’autant plus importante à prendre en compte que l’échantillonneur lui-même n’arrive qu’à la fin d’un processus d’échantillonnage qui a lieu naturellement : échantillonnage parmi les gamètes disponibles pour fabriquer les zygotes ; échantillonnage des zygotes qui participeront à la reproduction suivante. À cela s’ajoute également l’échantillon de marqueurs génétiques, sensés représenter le génome neutre, pour caractériser la variabilité génétique conditionnée par les événements démographiques de la population cible.

L’hypothèse nulle

2Comme son nom l’indique, c’est une hypothèse qui stipule qu’il ne se passe rien, ou que la population est conforme à une norme, un modèle préétabli (ou modèle nul), par exemple les fréquences génotypiques sont conformes à Hardy-Weinberg, ou les deux populations ont les mêmes fréquences alléliques, ou encore le F_IS (ou le F_ST) n’est pas différent de 0. On nomme cette hypothèse sous le diminutif H₀. L’hypothèse alternative, ou H₁, peut être indéfinie (par exemple, le F_IS est différent de 0) ou au contraire définie (ou orientée) (exemple, le F_IS est plus grand que 0). Dans ce dernier cas, on parle de test unilatéral qui, comme nous le verrons, est en général plus puissant que le premier (ou test bilatéral), sauf si on se trompe de direction (voir plus loin).

Qu’est-ce qu’un test statistique ?

3Un test statistique consiste en un calcul plus ou moins compliqué de la probabilité avec laquelle le hasard (et seulement lui) nous permet d’expliquer la déviation observée dans un échantillon par rapport à ce qui est attendu sous H₀.

4Prenons un exemple très simple. Je tire deux fois à pile ou face. Je peux soit obtenir deux piles avec la probabilité (½)², soit un pile et une face avec la probabilité ½ et deux faces avec la probabilité (½)². Je joue et obtiens deux faces. Nous allons procéder à trois tests statistiques.

• Test unilatéral 1 :
  H₀ : la pièce est bien équilibrée, ce que j’ai observé n’est pas significativement différent de l’attendu ½ / ½.
  H₁ : la pièce n’est pas bien équilibrée, j’obtiens plus de faces qu’attendu.
• Test unilatéral 2 :
  H₀ : la pièce est bien équilibrée, ce que j’ai observé n’est pas significativement différent de l’attendu ½ / ½.
  H₁ : la pièce n’est pas bien équilibrée, j’obtiens moins de faces qu’attendu.
• Test bilatéral :
  H₀ : la pièce est bien équilibrée, ce que j’ai observé n’est pas significativement différent de l’attendu ½ / ½.
  H₁ : la pièce n’est pas bien équilibrée et j’obtiens un résultat significativement différent de l’attendu.
• Pour le test unilatéral 1, la probabilité d’obtenir par hasard autant ou plus de faces est égale à P_u1 = (Somme des probabilités d’obtenir autant ou plus que deux faces)/(Somme des probabilités totales obtenues) = (½)²/1. Donc la P-value du test est P_u1 = 0,25.
• Pour le test unilatéral 2, la probabilité d’obtenir par hasard plus ou autant de piles est égale à P_u2 = (Probabilité de zéro pile + Probabilité de un pile et une face + Probabilité de deux piles)/(Somme totale) = ((1/2)² + 1/2 + (1/2)²)/1. Donc P_u2 = 1.
• Pour le test bilatéral P_b = (probabilité d’avoir deux faces ou deux piles) = ((1/2)² + (1/2)²)/1 = 0,5.

5Plusieurs choses peuvent ici être signalées. Tout d’abord, la plus basse des probabilités obtenues est 0,25. Ce qui illustre la faible puissance du test due à la faiblesse de l’échantillon. Il est difficile de prouver quelque chose avec de trop petits échantillons. Ensuite, on voit bien qu’on est beaucoup plus puissant en unilatéral si on teste dans la bonne direction, et beaucoup moins quand on teste dans la mauvaise. Il faut décider du test que l’on fait, bilatéral ou unilatéral et dans quelle direction, avant de faire le test. Il faut donc bien se poser la question avant, pas après. Si aucune information ne permet de définir dans quelle direction le signal doit avoir lieu, il faut systématiquement procéder à un test bilatéral. Par contre, si on est certain de la direction que le signal est censé prendre, alors le test unilatéral s’impose.

6Par exemple, je mesure le F_IS à partir d’un échantillon quelconque et j’obtiens une valeur légèrement plus grande que 0. Je pose tout d’abord mon hypothèse nulle :

7H₀ : F_IS n’est pas significativement différent de 0. Le test statistique va donc consister à calculer, compte tenu du nombre de loci sur lequel la mesure a été faite, le degré de polymorphisme de ces différents loci (nombre d’allèles, leur distribution) et le nombre d’individus génotypés, la probabilité d’avoir obtenu un F_IS aussi extrême2 ou plus extrême que celui observé, sous l’hypothèse d’une rencontre au hasard des gamètes dans la population d’où ont été tirés les individus génotypés (panmixie). Le test par défaut est unilatéral et le plus souvent pour les valeurs positives (H₁ : F_IS > 0), car la plupart des facteurs influençant ce paramètre génèrent des déficits en hétérozygotes (autofécondation, effet Wahlund…). Cependant, dans certains cas, comme celui d’une reproduction clonale partielle, on s’attend à des déviations dans les deux directions (Balloux et al., 2003 ; De Meeûs et al., 2006). Dans ce cas, et comme les logiciels disponibles ne donnent pas de tests bilatéraux, il faut cumuler les résultats des deux tests unilatéraux (F_IS > 0 et F_IS < 0) en sommant P_min + 1 – P_max, où P_min et P_max correspondent à la plus petite (test unilatéral le plus puissant) et la plus grande des deux P-values des deux tests ou, si P_max = 1 comme c’est le cas ici, en doublant P_min.

Risques de première et de seconde espèce

8En règle générale, on considère (arbitrairement) qu’un test est significatif quand la P-value à laquelle il est associé est inférieure ou égale à 0,05. Mais dans certains cas (que nous verrons plus loin), il peut s’avérer nécessaire d’être plus sévère et de baisser ce seuil. Personnellement, je me sens plus à l’aise avec une P-value < 0,01 pour rejeter H₀ et une P-value > 0,1 pour l’accepter. Le seuil à partir duquel on décide qu’une statistique est significative (rejet de H₀) est appelé risque de première espèce ou erreur de type I et noté α. Il s’agit du risque de se tromper en rejetant H₀ quand elle est vraie. Le risque de seconde espèce, ou erreur de type II, noté β, correspond au risque de se tromper en acceptant l’hypothèse nulle quand elle est fausse. Ce risque, qui est fonction de la puissance du test, est très rarement connu mais peut être appréhendé dans certaines circonstances. L’exemple du pile ou face ci-dessus est typiquement un cas où β est nécessairement très grand puisque, même si la pièce est truquée, on ne pourra jamais le détecter en ne faisant que deux essais.

Le principe des randomisations

9Dans la plupart des situations rencontrées en génétique des populations naturelles (si ce n’est toutes), il ne sera pas possible de procéder au calcul des probabilités exactes telles que dans l’exemple du pile ou face. Cependant, l’utilisation de programmes informatiques va nous permettre, sans beaucoup d’effort, d’estimer avec une excellente approximation, ces P-values. Il s’agit de procédures de ré-échantillonnage ou randomisations. Ces procédures se regroupent en deux grands types. Celles du premier type visent à obtenir un intervalle de confiance de l’estimateur étudié (par exemple, le F_IS), l’autre vise à simuler des populations suivant l’hypothèse nulle afin de pouvoir comparer la valeur observée à celles qu’on peut attendre sous H₀ (obtenues par simulation).

10La plupart des tests décrits dans ce manuel sont disponibles dans le logiciel Fstat 2.9.3. (Goudet, 2002, mise à jour de Goudet, 1995), qui est très convivial. D’autres logiciels sont aussi utiles :

• Genepop 3.4. (Raymond et Rousset, 2003, mise à jour de Raymond et Rousset, 1995b), Genepop 4 (Rousset, 2008), moins convivial, mais qui est le seul à proposer certaines procédures très utiles (comme celles testant des isolements par la distance entre individus) et leur version web ;
• Genetix 4.03, très convivial, en français qui propose des AFC (analyses factorielles des correspondances) ;
• MSA (Dieringer et Schlötterer, 2002), pas très convivial, mais qui propose différents calculs de distances génétiques.
• FreeNA (Chapuis et Estoup, 2007), qui calcule les F-statistiques, leur bootstrap ainsi que la distance de corde (Cavalli-Sforza et Edwards, 1967), avec ou sans correction pour les allèles nuls.

11Il en existe bien sûr bien d’autres que nous utiliserons dans la 2^e partie de ce manuel « Applications à des exemples concrets », mais avec ces trois-ci on peut déjà faire énormément de choses. Ajoutons que ces logiciels sont téléchargeables gratuitement (voir en annexe les liens), chose à ajouter au crédit de leurs auteurs. Nous reviendrons sur d’autres logiciels au moment où nous en aurons besoin.

12Mais avant tout, il y a Create (Coombs et al., 2008) qui permet, à partir d’un fichier texte ou Excel avec toutes les données brutes, de convertir ces données dans un format adéquat pour la plupart des logiciels de génétique des populations. Au moment où je corrige mon manuscrit, Tatiana Giraud m’apprend qu’il en existe un autre PGD-Spider (Lischer et Excoffier, 2012), apparemment assez convivial, mais que je n’ai encore jamais utilisé.

Intervalles de confiance de bootstrap et jackknife

Le bootstrap

13Il s’agit d’un rééchantillonnage répété avec remise. On sélectionne au hasard un des réplicats et, après avoir noté sa valeur, on le remet et ainsi de suite jusqu’à obtention d’autant de mesures qu’il y a de réplicats dans l’échantillon. La procédure est répétée un grand nombre de fois (5 000 pour Fstat). On obtient ainsi une distribution de 5 000 valeurs possibles. En excluant les 2,5 % (0,5 %) plus petites et les 2,5 % (0,5 %) plus grandes de ces valeurs, on obtient l’intervalle de confiance à 95 % (99 %). Ces notions seront plus claires avec la description du bootstrap sur les loci et sur les populations.

Bootstrap sur les loci

14On rééchantillonne au hasard et avec remise les k différents loci disponibles, jusqu’à en avoir k sur lesquels on recalcule l’estimation du paramètre (F_IS ou F_ST). On recommence l’opération un très grand nombre de fois (5 000 fois). Notons que, puisqu’il s’agit d’un échantillonnage des loci avec remise, on peut obtenir plusieurs fois le même locus. On obtient ainsi une distribution des valeurs obtenues sur les 5 000 bootstraps. Il ne reste plus ensuite qu’à regarder les valeurs obtenues de part et d’autre de cette distribution pour obtenir un intervalle de confiance. Par exemple pour 5 000 bootstraps, la valeur obtenue avant les 2,5 % les plus fortes et après les 2,5 % les plus faibles nous donne l’intervalle de confiance à 95 % (voir la figure 10).

Figure 10

15Dans la figure 10, on voit que la valeur observée n’est pas centrée, car le bootstrap génère des distributions décalées. Cette procédure sert à comparer des statistiques F entre différents échantillons ou groupes. Elle n’est pas très puissante, mais fournit la possibilité de faire des graphiques élégants. En général, on l’accompagne d’un autre test plus puissant, par exemple un test de Wilcoxon pour données appariées (par loci si les loci étudiés sont les mêmes) ou un test de Kruskal-Wallis si les loci ne sont pas les mêmes (les deux tests sont implémentés dans tous les logiciels de statistiques). Attention, si les loci ne sont pas les mêmes, la différence observée entre groupes pourra provenir des loci et non des groupes. D’une manière générale, il vaut mieux s’assurer de travailler avec les mêmes loci et que ces derniers soient en nombre suffisant (au moins sept). De toutes les façons, un bootstrap sur les loci ne commence à avoir du sens qu’à partir de quatre loci, et est vraiment puissant bien au-delà (voir Raymond et Rousset, 1995a pour discussion).

Bootstrap sur les populations

16C’est exactement le même principe que le précédent sauf que ce sont les populations (ou ce que l’on considère comme telles, les sous-échantillons) qui sont ici rééchantillonnées. Attention, on ne peut pas faire cela pour le F_ST3. Seul donc le F_IS est concerné. Cette procédure permet de comparer les loci entre eux. Il est en effet important de vérifier si les différents loci convergent vers le même signal, car sinon il sera utile de déterminer les causes responsables des discordances entre loci.

Le jackknife

Jackknife sur les loci

17Ici, il s’agit de prendre chaque locus un à un et de calculer la valeur du F sur ceux qui restent. On obtient ainsi k valeurs sur lesquelles on peut calculer une moyenne et une variance et donc une erreur standard. L’erreur standard du jackknife d’une statistique x évaluée sur n mesures se calcule selon l’équation (à comparer avec l’équation 28) (McIntosh, 2016) :

18Il est ensuite facile de calculer à partir de là l’intervalle de confiance voulu (IC), en faisant l’hypothèse que la distribution des jackknives suit une distribution normale (ce qui n’est probablement pas tout à fait exact, mais passons).

19où F est le F de Wright étudié, et t_n-1,α le paramètre de la loi normale pour n – 1 degré de liberté (n correspond ici au nombre de loci) et au seuil α (α = 0,05 pour un IC de 95 %). Cette valeur du t peut être retrouvée à partir d’une table (tabl. 3) ou d’un programme informatique.

20Pour n très grand et α = 0,05, on a classiquement t = 1,96. Sous Excel, la commande est = LOI.STUDENT.INVERSE(A1;A2-1) où A1 correspond à la case de coordonnées de la colonne A, 1^re ligne où vous avez tapé la valeur pour α et A2 est la case où vous avez tapé la valeur du nombre de réplicas, d’où on retranche 1 pour avoir le degré de liberté.

Jackknife sur populations

21Même chose que pour les loci, mais avec les sous-échantillons. Notons que le F_ST peut se prêter à cette procédure ici, puisqu’on a toujours n – 1 sous-échantillons sur lesquels calculer un F_ST, ce qui n’est pas garanti par le bootstrap (le même sous-échantillon peut être échantillonné n fois par la procédure).

Applications numériques pour le jackknife

22Supposons un jeu de données sur huit échantillons avec cinq loci. Sur l’ensemble des loci F_ST = 0,004, et pour le Locus 1 F_ST1 = 0,002. Le jackknife sur loci (cinq valeurs) donne une erreur standard de StdErrLoci(F_ST) = 0,003. Le jackknife sur populations (huit valeurs) donne une erreur standard de StdErrPop(F_ST1) = 0,001 pour le Locus 1. Alors, les intervalles de confiance à 95 % de F_ST et F_ST1 seront :

23avec γL = 5 – 1 = 4 et γP = 8 – 1 = 7 correspondant aux degrés de liberté des procédures de jackknife sur loci et populations respectivement. En utilisant le tableau 3, nous obtenons alors t_0,05,γL = 2,776 et t_0,05,γP = 2,365, d’où l’on tire que CI(F_ST) = 0,004 ± 0,008 et CI(F_ST1) = 0,002 ± 0,002.

24En règle générale, je préfère utiliser le bootstrap car il ne fait pas d’hypothèse, mais si je souhaite obtenir des intervalles de confiance du F_ST pour les différents loci, je suis bien obligé de le faire avec un jackknife sur populations.

25Comme pour le bootstrap, il faut au minimum cinq réplicats (loci ou populations) pour qu’un jackknife soit effectué par Fstat.

Mise en garde

26Comme nous l’avons vu, le nombre de réplicats à rééchantillonner doit respecter une valeur minimum. Il est nécessaire qu’il y ait au moins cinq loci et/ou sous-échantillons pour que ces procédures soient effectuées. Il est également nécessaire que ces réplicats soient suffisamment variables, un locus monomorphe ou presque pas variable ne pourra pas offrir un réplicat digne de ce nom, même si Fstat effectue la procédure sans problème (c’est-à-dire sans vous prévenir qu’il y a potentiellement un souci). Le calcul d’un intervalle de confiance de jackknife fait l’hypothèse d’une distribution normale, ce qui est probablement faux. Cet intervalle sera donc plus illustratif que statistiquement utile.

Les permutations

27Il s’agit ici de simuler l’hypothèse nulle un grand nombre de fois avec les données. Le programme informatique va utiliser les données (c’est-à-dire les allèles ou les individus des différents sous-échantillons) pour simuler H₀, mesurer la valeur obtenue sous H₀, recommencer un très grand nombre de fois afin d’obtenir une distribution des valeurs possibles sous H₀. La P-value du test correspond donc simplement à la proportion des cas où une valeur aussi grande ou plus grande (unilatéral 1), aussi petite ou plus petite (unilatéral 2), aussi extrême (bilatéral) que la valeur observée a été obtenue dans cette distribution.

28Il existe deux grands types de randomisations : les permutations et les chaînes de Markhov. Les deux sont basées sur le principe de Monte Carlo. Le nom provient d’un clin d’œil de Metropolis à son collègue Stan Ulam et n’est pas sans rapport avec les jeux de hasard ayant cours dans la capitale de Monaco (voir Metropolis, 1987)4. Il s’agit de rééchantillonnages aléatoires (ou pseudo-aléatoires) des données.

29La procédure de permutation correspond à la méthode utilisée dans Fstat (Goudet, 1995). Il s’agit en fait de simuler l’hypothèse nulle un grand nombre de fois (par exemple, 10 000 fois) à partir des données existantes (l’échantillon). Par exemple, dans chaque sous-échantillon et pour chaque locus, les allèles de tous les individus sont réassociés deux à deux au hasard pour reformer des individus fictifs obtenus selon l’hypothèse de panmixie locale. Un F_IS, obtenu donc par hasard sous H₀, est calculé et on recommence la même chose un très grand nombre de fois. La proportion de fois qu’un F_IS aussi grand ou plus grand que celui observé dans les données est apparu au cours du processus correspond à la P-value du test (H₁ étant ici F_IS > 0). Selon un principe analogue, la différenciation (H₀ : les individus se distribuent de façon aléatoire dans les différents sous-échantillons) est testée en assignant chaque individu aléatoirement dans les différents sous-échantillons, un F_ST obtenu sous H₀ est calculé et le processus répété. La proportion de fois qu’un F_ST (ou une autre statistique reflétant la distribution des fréquences alléliques entre sous-échantillons, comme nous le verrons plus loin) simulé sous H₀ a été aussi grand ou plus grand que l’observé procure la P-value du test.

30La chaîne de Markhov correspond à la procédure utilisée dans Genepop (Raymond et Rousset, 1995b). Le principe en est le suivant. Il s’agit de définir une trajectoire aléatoire entre différents tableaux de contingences possibles et dont les sommes marginales sont identiques au tableau des données réelles. La probabilité d’apparition (sous H₀) de chacun des tableaux ainsi obtenus est comparée à celle correspondant au tableau de contingence observé. La probabilité du test est obtenue en comptant le nombre de fois qu’une probabilité s’est montrée inférieure ou égale à celle du tableau observé et en divisant cette valeur par le nombre total de tableaux générés durant le processus. Une description mieux détaillée est présentée dans Rousset et Raymond (1997).

31Les P-values obtenues par ces méthodes constituent d’excellentes approximations des probabilités exactes, pour peu que l’on ait pris garde de mettre en œuvre un nombre suffisant de ces randomisations (un minimum de 1 000 à 10 000 pour les permutations et 10⁶ à 10⁷ pour les chaînes de Markhov), ce qui ne devrait pas représenter un problème avec les micro-ordinateurs d’aujourd’hui.

Tester la panmixie locale

Tester le F_IS

32La panmixie locale peut se tester en prenant les allèles présents dans chaque sous-échantillon et en les réassociant au hasard à l’intérieur de ces sous-populations et ce dans toutes les sous-populations. On mesure alors le F_IS global (moyenne sur l’ensemble des sous-échantillons et des loci) (estimation par f de Weir et Cockerham, 1984). Ce processus est répété un très grand nombre de fois, ce qui permet d’obtenir la distribution des F_IS générés sous l’hypothèse de panmixie locale (H₀). Trois tests sont ensuite possibles (en toute rigueur, il faut choisir lequel avant).

Tester s’il existe un déficit en hétérozygotes

33Il s’agit donc d’un test unilatéral avec H₁ : le F_IS de la population échantillonnée est plus grand que 0. On teste bien ici si les sous-populations échantillonnées sont panmictiques (H₀) et n’utilisent pas un mode fermé de reproduction (du type autofécondation ou croisements entre apparentés) qui doit donner une homozygotie supérieure à celle attendue sous panmixie à tous les loci. La proportion de fois que l’on obtient, au cours de la randomisation, une valeur aussi grande ou plus grande que celle observée nous donne la P-value du test. Si le test est significatif, on peut ensuite essayer d’estimer le taux d’autofécondation ou de croisements frère-sœur qui permet d’expliquer le F_IS observé, ou encore s’il peut être expliqué par un effet Wahlund (voir plus loin).

Tester s’il existe un excès d’hétérozygotes

34C’est le test unilatéral dans l’autre sens avec H₁ : le F_IS de la population est inférieur à 0. La proportion de fois que l’on obtient, dans les randomisations, une valeur aussi faible ou plus petite encore que celle observée nous donne la P-value du test. Ici, ce qui est recherché c’est si les sous-populations se reproduisent de façon asexuée (clonalité) stricte, auquel cas on attend un F_IS < 0 pour tous les loci.

Tester un écart dans n’importe quelle direction (excès ou déficit)

35Il se peut également que l’on s’attende à ce que les différents sous-échantillons ou les différents loci répondent dans toutes les directions (configurations de tests bilatéraux). Dans ce cas, comme nous allons le voir ci-dessous, le problème est assez simple à résoudre, que ce soit par locus (test pour chaque locus) ou sur l’ensemble.

36La P-value bilatérale s’obtient comme suit. Il faut faire les deux tests unilatéraux, ce qui fournit deux P-values. Soit P_min la plus petite de ces deux probabilités (test unilatéral le plus puissant) et P_max la plus grande des deux (test unilatéral le moins puissant des deux). Alors, la P-value bilatérale est simplement égale à P_min + (1 – P_max). Comme mentionné plus haut, il se peut que P_max soit inutilisable ou bien qu’elle ne soit pas calculable. Dans ce cas, on peut aussi multiplier P_min par deux.

Autres méthodes pour tester l’écart à la panmixie

Tests exacts

37Tester la conformité à la panmixie est synonyme de tester une conformité à une rencontre aléatoire des gamètes et donc aux proportions génotypiques attendues sous les hypothèses de Hardy-Weinberg (de la forme p_i², 2p_i(1 – p_i) pour un allèle i quelconque). Ceci peut être également testé avec le test exact de Haldane (1954) (souvent improprement appelé test exact de Fisher qui correspond en fait à autre chose) ou par la généralisation de ce test pour plus de deux allèles (Guo et Thompson, 1992), comme cela est proposé dans Genepop (Raymond et Rousset, 1995b). Je n’encourage cependant pas l’utilisation de ce test pour plusieurs raisons. D’abord, ce test analyse la distribution de tous les génotypes possibles et pas forcément ce qui est affecté par un système de reproduction particulier (telle que l’autofécondation). Pour un locus à plus de deux allèles, l’hypothèse nulle peut être rejetée parce que certaines classes génotypiques manquent au profit de certaines autres, alors que le reste est conforme à l’attendu sous panmixie. Un tel résultat sera difficile à interpréter biologiquement. Ensuite, ce test n’est réalisable que dans chaque sous-échantillon et pour chaque locus séparément. Il n’y a pas de test global possible et on se retrouve obligé de gérer une multitude de P-values, exactes certes, mais, en addition du premier problème, bien difficile à exploiter en termes d’inférence biologique.

Méthode de Rousset et Raymond (1995)

38Dans le logiciel Genepop, ce n’est pas f de Weir et Cockerham (1984) mais c’est un analogue de l’estimateur de Robertson et Hill (1984) qui est utilisé comme statistique pour effectuer le test. Comme nous l’avons déjà évoqué, cet estimateur est biaisé mais montre des variances en général plus faibles (pour des valeurs faibles). Les deux techniques convergent dans la plupart des situations et les différences obtenues ne concernent en général que les résultats à un locus dans un ou quelques sous-échantillons et n’ont donc que très rarement une conséquence interprétative (ou inférentielle) importante. Par conséquent, les deux procédures donnent presque toujours des résultats comparables.

Tester la pangamie

39Tous les tests décrits précédemment ne sont en fait que des approches indirectes, puisque ce n’est que la conséquence de la rencontre aléatoire des gamètes qui est testée, mais pas cette rencontre à proprement parler. Dans certaines circonstances, il est possible de tester la panmixie plus directement, si on a accès aux couples naturellement formés dans le milieu. En génotypant les adultes accouplés (en copulation), on peut tester si ces adultes se sont associés indépendamment de leurs génotypes, c’est-à-dire on peut tester la pangamie. Pour ce faire, il suffit d’adapter un test de Mantel (1967), test que nous détaillerons davantage plus loin pour les tests de corrélations entre matrices de distances, pour tester la corrélation entre la matrice des distances génétiques (apparentement) entre les individus possibles (entre les différents individus accouplés) et la matrice d’accouplement (en codant 0 pour les paires d’individus non accouplés et 1 pour les paires d’individus effectivement trouvés accouplés). Pour des organismes à sexes séparés, les matrices concernent les femelles d’un côté contre les mâles de l’autre. Attention, ce test de Mantel ne peut être effectué par Genepop qui ne gère que des demi-matrices en excluant les valeurs diagonales (dont on a besoin ici). Il faut donc effectuer le test avec un logiciel qui utilise des données en colonne (comme Fstat, ou RT de Manly, 1997). En permutant les cases d’une des deux matrices et en calculant un coefficient de corrélation à chaque fois, on obtient ensuite la probabilité d’observer une valeur aussi extrême ou plus extrême que celle observée. Cette procédure, malgré son intérêt évident, n’a à notre connaissance été utilisée qu’à deux reprises : chez le trématode Schistosoma mansoni chez les rats de Guadeloupe (Prugnolle et al., 2004b) et chez la tique du bétail Rhipicephalus (Boophilus) microplus en Nouvelle-Calédonie (Chevillon et al., 2007a). Dans le premier cas, l’apparentement entre les paires d’individus a été effectué à l’aide du logiciel Kinship V.1.2. (module Relatedness) développé par K. F. Goodnight (http://gsoft.smu.edu/GSoft.html) qui calcule un estimateur d’apparentement non biaisé équivalent de celui décrit dans Queller et Goodnight (1989). Pour les tiques, c’est l’estimateur de Wang (2002) qui a été préféré, car particulièrement robuste aux petits échantillons. Ce dernier fut calculé par le logiciel MER V3 (http://www.zoo.cam.ac.uk/ioz/software.htm#MER). Nous reverrons ce dernier exemple dans la seconde partie de ce manuel.

40Dans tous les cas, l’information apportée par ce test peut s’avérer précieuse pour discuter des hypothèses possibles en vue d’expliquer une déviation du F_IS par rapport aux attendus sous panmixie.

Tester la structuration

Tester le F_ST

41Il s’agit de simuler la migration libre des individus entre sous-échantillons (H₀) en redistribuant au hasard les individus dans ces différents sous-échantillons. On mesure alors le F_ST obtenu avec θ (sur l’ensemble des loci). La répétition de ce processus un très grand nombre de fois (10 000) nous permet d’obtenir une distribution des F_ST possibles sous H₀. L’hypothèse alternative H₁ correspond nécessairement à : « Il y a structuration », ou autrement dit, « le F_ST de l’échantillon est plus grand que 0 ». La P-value est donc donnée par la proportion de cas où le F_ST simulé a été aussi grand ou plus grand que le F_ST observé. Ce test est disponible dans Genetix qui propose également la même procédure avec le θ_RH de Robertson et Hill (1984) et le θ_RH’ de Raufaste et Bonhomme (2000).

La méthode basée sur le G de Goudet et al. (1996)

42En fait, certains travaux ont montré que le calcul d’une autre statistique (G) permettait d’avoir une plus grande puissance du test dans la plupart des situations (voir Goudet et al., 1996). La procédure est rigoureusement identique sauf que l’on mesure un G (logarithme népérien de la vraisemblance du tableau de contingence observé) au lieu d’un F_ST. Cette statistique est calculée à partir d’effectifs alléliques, mais ce sont bien les individus diploïdes qui sont permutés au cours des randomisations (d’où le qualificatif de test génotypique). Une description de la formule du G peut être trouvée dans n’importe quel ouvrage de statistiques (Sokal et Rohlf, 1981) (voir aussi la réponse 7). Un avantage supplémentaire de cette statistique concerne ses propriétés additives, ce qui autorise la mise en place d’un test global sur l’ensemble des loci, comme on le retrouve dans Fstat, et qui se montre nettement plus puissant que toutes les autres méthodes pour combiner des tests indépendants (De Meeûs et al., 2009).

43où l indique le locus et nl est le nombre total de loci, k les sous-échantillons et np le nombre total de sous-échantillons, i l’allèle et na le nombre total d’allèles au locus l dans la population k, N_ikl est le nombre de fois que l’allèle i du locus k est rencontré dans la population l, N_kl est le nombre d’allèles (deux fois la taille du sous-échantillon chez des diploïdes) du locus l dans le sous-échantillon k et est la fréquence moyenne de l’allèle i du locus l dans tout l’échantillon. C’est donc cette statistique qui est calculée sur les données observées et pour chaque randomisation des individus entre sous-échantillons.

Test exact allélique de Rousset et Raymond (1995)

44Il existe une autre solution pour tester la différenciation entre dèmes, mise au point par Raymond et Rousset (1995a). Il s’agit d’un test purement allélique qui fait donc l’hypothèse d’une indépendance totale des allèles dans les individus (panmixie parfaite). Pour que ce test soit valide, il est donc indispensable que les génotypes soient en parfaite conformité avec les attendus sous Hardy-Weinberg, car ce sont les allèles qui sont ici randomisés entre sous-échantillons. Ce test est proposé comme test allélique dans Genepop. Une procédure équivalente, le test « assuming HW » est proposé dans Fstat (qui utilise un test basé sur le G). C’est le test le plus puissant qui existe, mais, parce qu’il est probable qu’aucune population ne soit en conformité avec une panmixie parfaite, je conseillerai de ne jamais appliquer ces procédures et de leur préférer celles utilisant les génotypes (ne supposant donc pas la panmixie). Par ailleurs, le test exact ne peut être effectué que locus par locus, ce qui impose une procédure supplémentaire pour obtenir un test global (voir plus loin le paragraphe sur les tests multiples). Ajoutons enfin que les logiciels cités traitent les données haploïdes en dédoublant chaque allèle (homozygotie artificielle totale). Dans ce cas, on ne peut utiliser le test allélique.

Tester la panmixie globale

45Ceci est fait en réassociant au hasard les allèles des individus de l’ensemble de l’échantillon un très grand nombre de fois. On mesure le F_IT sur l’ensemble des loci. Pour le reste, la procédure est identique à celle présentée pour tester la significativité du F_IS.

46Il peut sembler redondant de tester le F_IT après avoir testé le F_IS et le F_ST, mais dans certains cas cela peut s’avérer utile. En particulier, un F_IT nul associé à d’autres critères (voir plus loin) peut être diagnostique d’une espèce strictement clonale et fortement structurée en de nombreux dèmes (voir De Meeûs et Balloux, 2005 ; Nébavi et al., 2006).

Tester les déséquilibres de liaison

47Ici, plusieurs méthodes sont possibles. Globalement, elles consistent à recombiner au hasard les loci entre eux à l’intérieur de chaque sous-échantillon un très grand nombre de fois et de mesurer (différentes méthodes) une statistique. La statistique observée dans chaque sous-échantillon est ensuite comparée à la distribution obtenue lors des randomisations sous l’hypothèse nulle d’absence d’association statistique entre loci. La statistique peut être une mesure de déséquilibre de liaison par paire de loci (le plus fréquent), ou une mesure multiloci (utilisée par les chercheurs travaillant sur des organismes clonaux). L’avantage des mesures multiloci est qu’elles fournissent une mesure sur l’ensemble des loci, alors qu’il y a autant de mesures (et donc de tests) qu’il y a de paires de loci (potentiellement L(L – 1)/2 où L est le nombre de loci) pour les mesures par paire. Le défaut des mesures multiloci est que leur comportement n’est pas encore bien connu dans toutes les conditions (voir De Meeûs et Balloux, 2004) et qu’il n’existe pas de mesure (et donc de test) multi-échantillons. Dans les tests par paire de loci, on peut utiliser comme statistique la probabilité d’apparition du tableau des génotypes pour les deux loci du sous-échantillon, compte tenu des fréquences génotypiques observées. Dans ce cas, la P-value du test sera simplement la somme des probabilités aussi faibles ou plus faibles que celle observée dans le sous-échantillon (voir le système d’équations 25), divisée par la somme de toutes les probabilités obtenues lors de la procédure de randomisation. Autrement dit, si P_obs est la probabilité du tableau des génotypes observés pour la paire de loci L1_L2 dans le sous-échantillon S1, P_i la probabilité d’occurrence d’un tableau randomisé et Rand le nombre total de randomisations (nombre de fois que les génotypes ont été recombinés librement), alors la P-value du test de déséquilibre de liaison sera :

48C’est ce qui est fait dans Genepop 3.4. (Raymond et Rousset, 2003, mis à jour de Raymond et Rousset, 1995b). On peut aussi calculer une autre statistique, telle qu’un G comme dans le logiciel Fstat 2.9.3. (Goudet, 2002, mise à jour de Goudet, 1995) et Genepop 4 (Rousset, 2008), ou sur un coefficient de corrélation comme dans Genetix 4.03 (Belkhir et al., 2004) ou encore sur un estimateur multilocus comme dans Multilocus 1.3b (Agapow et Burt, 2003, mis à jour d’Agapow et Burt, 2001).

49À partir d’ici, plusieurs points importants doivent être précisés.

Nombre de randomisations

50Certaines procédures de randomisations peuvent être très gourmandes en nombre de randomisations. Ce nombre sera fonction du nombre de combinaisons de génotypes possibles entre les deux loci étudiés. Dans le doute, il faut donc bien veiller à vérifier que deux procédures de randomisations faites indépendamment sur les mêmes données donnent le même résultat. Ceci est particulièrement important pour la procédure (chaîne de Markhov) utilisée dans Genepop où le nombre d’itérations devra atteindre au moins 10⁶, voire 10⁷.

Correction du seuil

51Comme nous l’avons vu, les tests par paire de loci génèrent un grand nombre de tests (autant que de paires de loci). Pour sept loci, par exemple, on a 21 paires de loci possibles. Cette répétition de tests va poser un problème statistique important que nous traiterons dans la section suivante. Ces tests sont par ailleurs non indépendants puisque chaque locus est comparé à chacun des autres loci restants, ce qui signifie que l’information contenue dans chaque locus est utilisée de façon redondante, ce qui pose un problème supplémentaire. Dans le paragraphe qui suit, nous verrons comment corriger le seuil de décision statistique afin de prendre en compte ces difficultés.

Remarques sur les tests de déséquilibres de liaison et leur interprétation

52Comme nous l’avons déjà vu, certaines des procédures que nous utilisons en génétique des populations empiriques requièrent l’utilisation d’un nombre important de loci (au moins cinq) qui devraient être indépendants statistiquement. C’est-à-dire que l’information portée par chacun de ces loci est supposée indépendante. Un déséquilibre de liaison fort risquerait d’apporter une redondance forte conduisant à un risque d’erreur de décision. En fait, l’indépendance des loci ne peut être certaine que si les populations échantillonnées sont de tailles infinies, panmictiques et non structurées et ce depuis un grand nombre de générations, ce qui n’est évidemment jamais le cas. Il y a donc toujours liaison. Le principal est que cette liaison ne nuise pas trop à la détection du signal recherché. Le reste est laissé à l’appréciation de chacun, mais fort heureusement ces tests sont individuellement peu puissants et les procédures qui y sont le plus souvent associées (Bonferroni) rendent la détection de tels déséquilibres peu fréquente. De ma propre expérience sur les populations clonales (déséquilibres de liaison forts à totaux), c’est plus une diminution de puissance des tests (de différenciation, en particulier) qu’une augmentation qu’il faut attendre (augmentation des variances d’estimation), comme cela peut être illustré par les immenses intervalles de confiance de F_ST obtenus par bootstrap sur les loci chez la levure opportuniste Candida albicans (voir la figure 1 dans Nébavi et al., 2006).

Le problème des tests répétés

53Comme nous l’avons déjà vu, le but d’un test statistique est d’évaluer la probabilité avec laquelle le hasard permet d’expliquer nos données si celles-ci proviennent d’une population respectant l’hypothèse nulle. Si cette probabilité est inférieure à un seuil choisi α, on décide que les données dévient significativement de ce que l’on attend sous H₀. Par conséquent, et par définition, pour un seuil choisi de α = 0,05 (le plus classique), on s’attend à ce que sous H₀ 5 % des tests soient significatifs par hasard. Autrement dit, si j’échantillonne 100 fois dans une population panmictique et que j’effectue un test du F_IS pour chaque échantillon, je m’attends à trouver en moyenne cinq tests significatifs au seuil α = 0,05 (si la taille des échantillons et si le polymorphisme des loci sont suffisants).

54Par conséquent, la répétition de tests pose un problème. Plusieurs méthodes existent pour résoudre le problème des tests répétés et dont l’application dépend de la question posée et du type de tests répétés.

Les tests répétés sont indépendants

55Ces tests répétés peuvent correspondre à différents cas de figure dont voici une liste non exhaustive :

• je voudrais combiner différents tests (de la même H₀) trouvés dans la littérature pour lesquels je n’ai pas les données brutes ;
• je cherche à savoir si le F_IS de chaque locus dévie significativement de 0 dans un sens ou dans l’autre ;
• je dispose de données de structuration de plusieurs sites comparables, sur plusieurs années et je cherche à combiner les P-values obtenues lorsque j’ai testé la significativité du F_ST dans chacun de ces jeux de données d’années différentes ;
• je compare la différenciation entre deux catégories d’individus (mâles versus femelles ; parasites d’hôtes d’espèces différentes ou de sexes différents, etc.) dans plusieurs sites (je souhaite combiner l’information de tous les sites).

56Dans tous les cas, je peux chercher à savoir si un signal global existe ou je peux désirer identifier quels tests sont significatifs.

Tester si un signal global existe

57On peut alors combiner les k tests de quatre façons différentes : le test binomial et sa version généralisée (Teriokhin et al., 2007 ; De Meeûs et al., 2009), la procédure de Fisher (Fisher, 1970), le test SGM (Goudet, 1999) et la transformation Z de Stouffer (Whitlock, 2005).

58On peut procéder à un test binomial pour un nombre d’essais correspondant au nombre de tests et un attendu correspondant au seuil α. Pour α = 0,05, la structure du test est la suivante :

• H₀ : la proportion de tests significatifs observés n’est pas différente de 0,05 ;
• H₁ : la proportion observée de tests significatifs est supérieure à l’attendu 0,05 (test unilatéral).

59La plupart des logiciels de statistiques font le test binomial et son application est assez simple. La loi binomiale concerne les cas où on ne peut avoir que deux possibilités : vrai ou faux, présence ou absence, noir ou blanc ou, comme ici, significatif ou non. Elle est définie par le nombre d’essais (ou taille de l’échantillon) k, les probabilités complémentaires d’état de l’événement p et q = 1 – p pour significatif et non significatif respectivement et k’ le nombre de fois où l’événement « significatif » a effectivement été observé parmi les N essais. Dans notre cas, k correspond donc au nombre de tests que l’on souhaite combiner, et k’ au nombre de tests significatifs au seuil de 5 % parmi ces k tests. On souhaite avoir la probabilité d’obtenir par hasard un nombre de tests significatifs aussi grand ou plus grand que k’. Cette probabilité est :

60où k! = k(k – 1)(k – 2)…(k – k + 2)

61Donc si on a dix tests dont cinq sont significatifs, on a P = 0,00006 (valeur hautement significative donc). Pour un seul test significatif observé sur 10, cette P-value devient 0,4. Il existe depuis peu une version généralisée (Binomial généralisé) de ce test (Teriokhin et al., 2007) implémentée par le logiciel MultiLocus V2.2 (De Meeûs et al., 2009). La philosophie de ce test est décrite en détail dans l’aide qui accompagne le logiciel et je ne reviendrai donc pas dessus.

62La procédure de Fisher (Fisher, 1970), qu’il ne faut pas confondre avec le test exact du même auteur car cela n’a pas de rapport, propose la formule suivante :

63où P_i correspond à la P-value obtenue au ième test.

64Cette expression suit normalement une loi du χ² (Chi-2) avec 2k degrés de liberté (ddl), dont on peut donc extraire la P-value associée à partir d’une table du χ², d’un logiciel ou en tapant la formule LOI.CHIDEUX(χ²_obs;2*k) sous Excel.

65Le test de randomisation SGM de symétrie autour de 0,5 de la moyenne géométrique (la moyenne géométrique correspond à la racine kième du produit des k P-values entre elles) (Goudet, 1999) est implémenté par le logiciel SGM distribué sur demande par l’auteur lui-même.

66Le test de transformation Z de Stouffer (Whitlock, 2005) consiste en la transformation des P-values en leur équivalent Z, avec par exemple la commande Excel LOI.NORMALE.INVERSE(P_i;0;1) ou LOI.NORMALE.STANDARD.INVERSE(P_i) (mettre 0,9999 pour les P_i = 1) qui donne un Z_i pour chaque P_i que l’on combine en la statistique Z_s :

67La P-value globale est obtenue en comparant cette statistique à la loi normale, avec par exemple la commande Excel LOI.NORMALE.STANDARD(Z_s).

68De mon espérience, la procédure binomiale généralisée est plus performante dans plus de situations. Il faut juste prendre garde à choisir un seuil de k’ = k quand k < 4 ou k’ = k/2 dans les autres cas (De Meeûs, 2014). Par ailleurs, d’une façon qui ne concerne pas vraiment notre propos ici, le test binomial n’exige pas de connaître la P-value exacte des tests à combiner (même si cela est préférable), ce qui peut représenter un avantage certain lorsque l’on combine des données de la littérature.

69On pourra trouver une discussion plus théorique de ce type de problèmes dans la littérature (Goudet, 1999 ; Whitlock, 2005 ; De Meeûs et al., 2009). La procédure de Fisher sera inadéquate dans certaines configurations de distribution des P-values (en U, en cloche, en L, ou en J) autres que la distribution uniforme. Il est en effet important de faire attention à cela et ne pas être esclave de ses données. La présence d’une P-value = 0 devrait en toute rigueur interdire l’utilisation de la procédure de Fisher.

70La procédure de Fisher répond davantage à la question : y a-t-il au moins un test significatif ?

71Le SGM est quant à lui très (trop) conservateur, une propriété qui pourrait s’avérer utile dans le cadre des méta-analyses (sur jeux de données publiées) où le biais de publication en faveur des résultats significatifs pourrait être ainsi partiellement corrigé.

Déterminer quels sont les tests significatifs, procédures de type Bonferroni

72Une autre configuration pourrait nous amener à rechercher lesquels, parmi ces k tests, sont réellement significatifs. Ce peut être le cas si on recherche un marqueur de sous-dominance au milieu de plusieurs marqueurs (quels loci sont déficitaires en hétérozygotes ?). Dans ce cas, il n’y a pas d’autre solution que de procéder à une correction de Bonferroni (très conservatrice) (Holm, 1979 ; Rice, 1989) ou, de façon moins conservatrice (sauf pour la plus basse P-value), le Bonferroni séquentiel.

73Il faut ordonner les N P-values de chaque test de la plus petite à la plus grande. La plus petite des probabilités est multipliée par N, la deuxième plus petite par N – 1, la troisième par N – 2, etc. jusqu’à tomber sur un test non significatif après correction. Les tests significatifs sont ceux dont la P-value ainsi corrigée reste inférieure au seuil choisi α (= 0,05). On peut aussi plus classiquement corriger le seuil lui-même par le nombre de tests correspondants : α’₁ = α/N, α’₂ = α/(N – 1), α’₃ = α/(N – 2), etc., ce qui revient au même, mais personnellement je préfère disposer d’une P-value « exacte ». Le tableau 4 donne un exemple d’une série de 10 tests triés par ordre croissant de leur P-value, leur P-value corrigée par la procédure du Bonferroni séquentiel, le seuil corrigé correspondant à α = 0,05, ainsi que la décision statistique concernant chaque test, compte tenu du nombre de tests effectués.

74Dans cet exemple on décide que trois des dix tests ont donné une valeur déviant significativement de ce qui est attendu sous H₀. On remarquera que cette procédure est très conservatrice. Il faut le savoir au moment d’échantillonner et ne pas lésiner sur le nombre d’individus génotypés. De faibles échantillons ne permettront jamais d’obtenir des P-values suffisamment basses pour supporter une procédure de Bonferroni. C’est ce que les statisticiens appellent le risque de seconde espèce (β) ou risque de se tromper en acceptant H₀. Alors que α correspond au risque de première espèce, ou risque de se tromper en rejetant H₀. Nous avons déjà évoqué ces concepts en p. 71. Il existe une procédure moins conservatrice et simple à calculer avec R (R-Core-Team, 2020), c’est la procédure de Benjamini et Hochberg (1995) (BH), que je recommande maintenant d’utiliser en lieu et place de la beaucoup trop conservatrice procédure de Bonferroni.

Les tests répétés ne sont pas indépendants

75C’est typiquement le cas des tests de déséquilibre de liaison par paire de loci. C’est aussi le cas de tests de différenciation par paire de sous-échantillons. Ici encore, il est nécessaire de distinguer le cas où une réponse globale est souhaitée du cas où on recherche quels tests sont significatifs.

Tester si un signal global existe

76Je ne pense pas qu’un test global puisse être appliqué dans le cas de séries non indépendantes. À défaut, j’utilise la correction de Benjamini et Yekutieli (2001) (voir le paragraphe suivant) et regarde la proportion de tests qui restent significatifs, ainsi que le seuil de la plus grande des p-values significatives obtenue.

Déterminer quels sont les tests significatifs, procédure de Benjamini et Yekutieli

77Ici, s’il est souhaité de déterminer quelles paires de loci sont statistiquement associées ou quelles paires de sous-échantillons sont différenciées, dans ce cas il convient d’utiliser la procédure de Benjamini et Yekutieli (BY) (Benjamini et Yekutieli, 2001) (très facile sous R avec la commande « p.adjust »).

Le cas des déséquilibres de liaison

78Les déséquilibres de liaison représentent le cas de figure le plus fréquent. Ici, en fonction de la taille de chaque sous-échantillon et du degré de polymorphisme des loci, le seuil de BY peut s’avérer impossible à atteindre (trop conservateur). Dans ce cas, il est plus raisonnable de ne prendre en compte que les loci les plus polymorphes, c’est-à-dire ceux pour lesquels les chances de détecter quelque chose sont les plus grandes. Par exemple, on peut écarter les loci dont un allèle atteint ou dépasse la fréquence de 90 %, on peut même être plus sévère en fonction des circonstances, car de tels loci ne présenteront qu’exceptionnellement des P-values suffisamment faibles alors qu’ils contribuent à l’augmentation de la sévérité du BY. Comme suggéré précédemment, ce qui est la plupart du temps recherché c’est si les loci ne sont pas trop liés. Il n’y a de toutes les façons pas d’agrément général sur la meilleure procédure et il est donc laissé libre choix aux empiristes de décider si une correction plus ou moins sévère ou pas de correction doit être appliquée. Les tests Multilocus (Agapow et Burt, 2001) ont été spécifiquement conçus pour tester un effet global, tel que celui attendu sous régime clonal de reproduction. Des études de simulations (De Meeûs et Balloux, 2004) ont suggéré que la mesure la plus précise (ou plutôt la moins mauvaise) serait le coefficient de corrélation (Agapow et Burt, 2001). Ce coefficient se base sur l’indice d’association I_A (Brown et al., 1980 ; Maynard-Smith et al., 1993 ; Haubold et al., 1998), mais contrairement à ce dernier est indépendant du nombre de loci étudiés dans l’analyse. Cette mesure est également utilisée comme statistique dans les tests de randomisation implémentés dans Multilocus. Par rapport au test bi-locus cette procédure permet l’obtention directe d’un test global sur l’ensemble des loci, mais ne peut être réalisée que sous-échantillon par sous-échantillon. Il peut cependant n’être significatif qu’à cause d’une seule paire de loci. Le test bi-loci de Fstat (basé sur le G) permet d’obtenir un test sur l’ensemble des sous-échantillons, mais pour chaque paire de loci prise une à une. Là encore, le G global se montre plus performant que les alternatives (De Meeûs et al., 2009 ; De Meeûs, 2014).

Tester la corrélation entre distances

79Il s’agit ici de regarder si la différenciation génétique (distance génétique) que l’on observe entre les sous-populations de notre échantillon peut être due aux distances géographiques qui les séparent, ou à d’autres distances. Ces autres distances peuvent être écologiques (différences de températures moyennes de janvier, pluviométriques, etc.), ou même génétiques si on souhaite comparer les distances génétiques hôtes et parasites, par exemple.

80Dans tous les cas, on cherche à corréler deux matrices de distances entre elles. Comme pour les déséquilibres de liaison, les mesures de ces matrices ne sont pas indépendantes, ce qui empêche de procéder à des tests classiques de corrélation ou de régression. Nous allons donc effectuer un test de Mantel (1967).

81La valeur d’un paramètre d’association, ou d’un coefficient de corrélation, entre les deux matrices est calculée à partir des données réelles, puis comparée à la série de pseudo-valeurs obtenues par permutation aléatoire de l’ordre des populations dans l’une des deux matrices de distances. À partir de là, la suite devient similaire à tout autre test par permutation. On pourra consulter la réponse 8 pour plus de détails sur le test de Mantel.

Distances génétiques et géographiques

82Ce cas de figure a été étudié en profondeur par Rousset (1997). Ici, les sous-échantillons peuvent être distribués de deux façons différentes qui requièrent chacune une analyse qui lui est particulière. Cependant, dans les deux cas, la matrice des distances génétiques doit contenir une mesure corrigée de la différenciation entre paire de sous-populations, à savoir θ/(1 – θ) (voir Rousset, 1997), θ étant l’estimateur du F_ST (voir p. 55). On sait en effet par l’équation (21) que :

83Sachant que, dans le cas d’un isolement par la distance, c’est-à-dire quand la différenciation augmente avec l’éloignement géographique des individus, et si chaque dème connaît un fonctionnement raisonnablement similaire (à peu près même taille et même système de reproduction) on voit bien que Q_S, la probabilité d’identité de gènes entre deux individus de la même sous-population, sera à peu près la même d’une sous-population à l’autre alors que Q_T, probabilité d’identité entre dèmes, sera une fonction décroissante de la distance entre dèmes. On voit donc bien que puisque Q_T se trouve au numérateur, mais aussi au dénominateur du F_ST, la relation entre F_ST et la distance géographique ne peut pas être linéaire. Par contre, comme on peut le voir, le rapport F_ST/(1 – F_ST) ne subit pas ce problème, en effet :

84À partir de là, deux cas de figures sont à distinguer en fonction du schéma de dispersion des organismes étudiés, indépendamment du plan d’échantillonnage, qui peut suivre un tracé en une ou deux dimensions.

Les sous-échantillons sont alignés en une seule dimension

85Certains organismes sont susceptibles de se distribuer en une seule dimension et seront donc échantillonnés comme tels. C’est typiquement ce qui se passe pour des organismes côtiers comme des bivalves (moules), des cirripèdes (balanes et anatifs), des patelles, des algues (fucus vésiculeux, laminaires), ou les organismes vivant le long des cours d’eau (mouches tsé-tsé, lymnées). La matrice des distances par paire de sous-échantillons est alors comparée à celle des θ/(1 – θ) (car c’est bien l’estimateur non biaisé de Weir et Cockerham qu’il faut utiliser). Si le test de Mantel est significatif, on peut alors utiliser la pente b de la droite de régression du F_ST/(1 – F_ST) = a + bD_G, où D_G est la distance géographique, afin d’estimer le produit Dσ² de la densité d’adultes reproducteurs dans un site (D) et la moyenne des carrés de la distance axiale (moitié de la distance entre deux points) entre parents et descendants adultes (σ²). En effet, Rousset (1997) montre qu’alors :

86Cette méthodologie fut appliquée pour l’escargot intertidal (côtier) Bendicium vitatum (voir Rousset, 1997 pour une réanalyse) ou la tique d’oiseaux marins Ixodes uriae (McCoy et al., 2003). Une revue sur les glossines est parue en 2019 (De Meeûs et al., 2019b).

Les sous-échantillons sont distribués sur deux dimensions

87Dans ce cas, le test de Mantel doit être effectué entre la matrice des Log népériens des distances géographiques par paire de populations et celle des θ/(1 – θ) (voir Rousset, 1997). Si le test est significatif, la pente de la régression F_ST/(1 – F_ST) ≈ a + bLn(D_G) va permettre d’estimer Dσ² avec l’équation (Rousset, 1997) :

88Si l’un des deux paramètres D ou σ peut être estimé, même approximativement, de façon indépendante, on obtient un pouvoir d’inférence relativement puissant ici (voir Koffi et al., 2006a ; Bouyer et al., 2009 ; De Garine-Wichatitsky et al., 2009 pour illustration). Le lecteur trouvera également des commentaires utiles dans deux articles plus récents (Séré et al., 2017 ; De Meeûs et al., 2019b).

89Les mêmes procédures peuvent être appliquées entre individus entre lesquels un équivalent du F_ST/(1 – F_ST) (a) (calculé dans Genepop) et développé par Rousset (2000), Leblois et al. (2003) et Leblois et al. (2004) peut être régressé contre les distances entre individus (directe pour une dimension, en Log pour deux dimensions), ce qui conduit aux mêmes possibilités d’inférences que celles décrites ci-dessus. Watts et al. (2007) proposent une statistique e en principe plus puissante lorsque le voisinage (4Dσ² ou 4πDσ²) est grand. Nous verrons cela plus en détail dans la partie pratique de ce manuel.

90Dans le cas particulier de deux dimensions, Rousset (1997) montre que le nombre d’immigrants présents dans un sous-échantillon peut directement être tiré de la pente de la régression F_ST/(1 – F_ST) ≈ a + bLn(D_G), Nm = 1/2πb.

91Il faut enfin ajouter qu’une procédure de bootstrap sur les loci peut permettre de calculer un intervalle de confiance de la pente b qui, s’il ne contient pas le 0, permet de décider si la pente est significativement positive (voir à ce titre les ajouts que j’ai apporté aux analyses des données réelles pour la réédition de ce manuel).

Autres distances

92On peut souhaiter vérifier si la différenciation entre sites est corrélée à une différence écologique entre sites ou tester s’il existe une corrélation entre différenciation génétique des sous-échantillons des hôtes et des parasites qui les infestent. Comme nous l’avons vu précédemment, le F_ST a été défini dans le cadre d’un modèle en îles. De fait, il ne se comporte pas idéalement par paire de populations (fortes variances, voir Balloux et Goudet, 2002) et on lui préférera d’autres mesures pour les tests de Mantel telles que la distance de corde (chord distance) de Cavalli-Sforza et Edwards (Cavalli-Sforza et Edwards, 1967) ou la distance d’allèles partagés (shared allelic distance) (Bowcock et al., 1994) (déjà discuté en p. 62-63). Pour la construction d’arbres (dendrogrammes), il semble aussi que les distances de cordes donnent de meilleurs résultats (Takezaki et Nei, 1996).

93En fait, la performance de différentes mesures et leur choix vont dépendre des situations rencontrées, même si en principe toutes les distances devraient aboutir en théorie à des résultats concordants. Ceci peut être illustré par la corrélation que Prugnolle et al. (2005) ont montrée entre les distances génétiques entre infra-populations5 de schistosomes et celles mesurées entre les rats qui les portaient (ou leur apparentement si on préfère) en Guadeloupe. Dans l’article, c’était la distance de Cavalli-Sforza et Edwards (1967) qui avait été utilisée entre infra-populations de schistosomes et la « shared allele distance » (Bowcock et al., 1994) entre les individus rats. Le logiciel MSA (Dieringer et Schlötterer, 2003, téléchargeable à http://i122server.vu-wien.ac.at/) calcule cette distance. La corrélation obtenue était très significative (P-value = 0,0005), mais De Meeûs et al. (2007a) ont montré que si le F_ST est utilisé pour les deux matrices, la corrélation n’est plus significative (P-value = 0,15) et elle l’est beaucoup moins (P-value = 0,0113) quand c’est Cavalli-Sforza et Edwards qui est utilisé pour les deux matrices. Le choix d’une statistique n’est donc pas entièrement neutre. Ajoutons enfin que d’autres mesures d’apparentement entre individus existent, telles que l’estimateur de Queller et Goodnight (1989) ou de Wang (2002) dont nous avons déjà parlé à propos des tests de pangamie (p. 79-80) ou, plus récemment, de Kalinowski et al. (2006) qui pourrait être encore plus puissant. Notons que FreeNA (Chapuis et Estoup, 2007) calcule D_CSE et les F_ST par paires ainsi que l’intervalle de confiance de bootstrap de ces derniers.

Tester les biais de dispersion de certaines catégories d’individus

94Dans les populations naturelles, il se peut qu’un sexe disperse davantage que l’autre sexe ou que les individus parasités dispersent plus ou moins bien que les individus sains. Dans ce cas, il existe plusieurs statistiques (mesures) qui peuvent être comparées (Goudet et al., 2002). Je ne parlerai ici que de trois d’entre elles et dans le cas d’un biais de dispersion sexe-spécifique.

95L’indice d’assignement, dont nous avons déjà parlé en p. 66-67, consiste à calculer la probabilité qu’un individu a d’appartenir à la sous-population où il a été échantillonné, compte tenu de son génotype à tous les loci génotypés et de celui de l’ensemble des individus de son sous-échantillon. Afin de tenir compte du degré de polymorphisme qui peut beaucoup varier d’un site à l’autre, il faut corriger cette probabilité. Ceci est fait en soustrayant à cette valeur la valeur moyenne obtenue sur l’ensemble des individus du sous-échantillon concerné, après une transformation Log afin de minimiser les risques d’erreurs dus aux petites valeurs (voir Favre et al,. 1997 pour plus de détails). Cet indice se note AI_c. Sa distribution sur l’ensemble des populations est nécessairement centrée sur 0. Et une valeur négative signifie que l’individu est moins bien assigné que la moyenne à son propre sous-échantillon. On calcule ensuite la moyenne de cet indice sur les mâles et la moyenne sur les femelles de l’ensemble de l’échantillon. La statistique suivante est calculée :

96où les signes – et + désignent la catégorie qui disperse le moins et le plus respectivement, la barre désignant la moyenne, s² la variance et Nb le nombre total d’individus de la catégorie considérée, observés dans l’ensemble des sous-échantillons.

97La moyenne du sexe le moins dispersant (donc mieux assigné) doit être supérieure à celle du sexe le plus dispersant.

98La deuxième statistique qui nous intéresse correspond à :

99La variance de l’indice d’assignement du sexe le plus dispersant doit être supérieure à celle du sexe le moins dispersant.

100La troisième statistique dépend de la différence des F_ST estimés pour chaque catégorie :

101Δθ = θ(-) – θ(+) (56)

102La différenciation mesurée sur la catégorie d’individus les moins dispersants doit être plus élevée que celle mesurée pour la catégorie la plus vagile.

103Ensuite, l’appartenance à une catégorie (mâle ou femelle) est re-distribuée au hasard pour chaque individu de chaque sous-échantillon, en gardant les individus dans leur sous-échantillon, et en conservant la même proportion de chaque catégorie (même sexe-ratio) et la statistique est mesurée. Cette randomisation est répétée un grand nombre de fois afin d’obtenir une distribution des valeurs possibles sous H₀ (pas de différence de dispersion) à laquelle la valeur observée est comparée. Les tests peuvent être unilatéraux ou bilatéraux. Dans ce dernier cas, ce sont les valeurs absolues des différences [dans (53) et (55)] ou le ratio de la plus grande sur la plus petite valeur de chaque randomisation qui sont utilisés. Ces mesures et randomisations sont toutes implémentées dans Fstat (menu “biased dispersal”). Ces procédures ont été utilisées avec succès pour mettre en évidence, dans les populations suisses de la tique Ixodes ricinus, un biais de dispersion sexe-spécifique, les femelles représentant le sexe peu ou pas dispersant (De Meeûs et al., 2002a), et un biais de dispersion pathogène spécifique, les tiques infectées par le spirochète Borrelia afzelii dispersant très peu ou pas du tout (De Meeûs et al., 2004b). De même, Prugnolle et al. (2002) ont pu mettre en évidence une structure génétique spécifique du sexe chez le trématode Schistosoma mansoni infectant des rats en Guadeloupe.

104Dans certains cas, l’échantillonnage ne permet pas de tester une différence entre sexes ou entre catégories d’individus, par randomisation, notamment pour tester une différence de F_ST. Dans ce cas, une alternative moins puissante existe et permet de comparer H_s, F_IS ou le déséquilibre de liaison entre catégories d’individus dans un seul échantillon. Il suffit d’utiliser les loci (ou les paires de loci pour les déséquilibres de liaison) comme des répliquats (plus ou moins indépendants d’ailleurs) et de faire un test de comparaison pour données appariées, le critère d’appariement correspondant donc au locus (ou la paire de loci). Comme la distribution de telles données a toutes les chances de ne pas suivre une loi normale, il est conseillé ici de procéder à un test de rang de Wilcoxon pour données appariées (Wilcoxon signed ranks test for paired data) (Siegel et Castellan, 1 988).

Tester la différence entre groupes

105Ce cas de figure se présente lorsque différents types de sites doivent être comparés. C’est typiquement le cas si on souhaite comparer différents paramètres génétiques, tels que H_s, F_IS, F_ST ou d’autres, entre infra-populations trouvées dans des hôtes mâles et celles trouvées dans des hôtes femelles. Ce peut être aussi le cas entre des sites de différentes natures tels que des prés et des bois dans un paysage de bocages (par exemple, H₀ : les populations de bois sont-elles plus structurées que celles de pré, ou plus pamictiques, etc.). Ce peut également être le cas pour comparer des parasites trouvés sur des espèces hôtes différentes. Les procédures suivent toujours la même philosophie. Ici, le paramètre d’intérêt est moyenné sur l’ensemble des sous-échantillons de chaque catégorie. Soit x_i cette valeur moyenne pour les sous-échantillons du groupe i. Pour un test unilatéral avec deux groupes, on calcule juste la différence (x₁ – x₂) (x₁ étant la plus grande). Pour les autres cas de figure, la statistique utilisée sera :

106où ng représente le nombre de groupes à comparer.

107Ensuite, les échantillons de chaque groupe sont randomisés (permutations aléatoires des échantillons dans les différents groupes en gardant le nombre d’échantillon par groupe constant) un grand nombre de fois (10 000) et la statistique (x₁ – x₂) ou celle définie en (57) est recalculée pour chaque randomisation. La valeur observée est ensuite comparée à la distribution des valeurs obtenues par randomisation, la P-value du test correspondant (encore une fois) à la proportion de fois qu’une valeur aussi grande ou plus grande a été observée au cours des randomisations. Cette procédure est implémentée dans Fstat (menu “Comparison among groups”).

108Comme précédemment, l’échantillonnage peut ne pas permettre de procéder à ce test sans qu’il soit pour autant impossible de tester des différences de F_IS, de H_s ou de déséquilibres de liaison. Ici aussi, les loci (ou paires de loci) peuvent être utilisés comme répliquats pour un test de rangs pour données appariées (voir par exemple Nébavi et al., 2006).

Analyses multivariées

109Les analyses multivariées permettent souvent une représentation didactique de l’organisation générale de la variabilité génétique globale des échantillons génotypés. Dans certains cas, ils permettent également des analyses statistiques et des inférences. Il en existe plusieurs types, de même nature, mais offrant des possibilités différentes.

Analyse factorielle des correspondances (AFC)

110Cette analyse, introduite par Benzécri (1973), a été adaptée aux données génétiques diploïdes par She et al. (1987). L’AFC place chaque individu dans un hyper-espace à K dimensions (K étant le nombre total d’allèles présents sur l’ensemble des loci) et les projette sur les plans définis par les axes orthogonaux (donc indépendants) expliquant le mieux la dispersion des points (même principe que celui d’une régression). Une mesure de la pertinence des axes ainsi définis est représentée par le pourcentage d’inertie de chaque axe. Comme il y a K axes, un axe représentant 100/K % d’inertie ne veut rien dire. L’inertie est donc proportionnelle non seulement à la quantité d’information que l’axe correspondant représente, mais est aussi fonction du nombre total d’axes (plus il y a d’axes et moins chaque axe peut avoir une très forte inertie). L’AFC est une procédure qui peut s’avérer utile pour classer les individus en fonction de leur proximité génétique.

Exemples

111L’utilisation de l’AFC s’est avérée payante pour analyser la présence de trématodes parasites dans une zone d’hybridation de leur hôte (moule de bouchot, Mytilus edulis) avec une autre espèce (moule d’Espagne, M. galloprovincialis) incompatible pour le parasite (Coustau et al., 1991) ou, de façon plus spectaculaire, dans le cas du monogène Diplozoon gracile, spécifique du poisson Barbus meridionalis, en zone d’hybridation avec B. barbus, un hôte moins favorable au parasite, comme présenté dans la figure 11.

112Cette technique peut également être utilisée pour détecter une structure cachée dans un échantillon comme celle qui proviendrait d’un effet Wahlund (déficits en hétérozygotes à tous les loci non expliqués par le système de reproduction), comme cela a été réalisé dans Solano et al. (2000) (voir le paragraphe suivant). GENETIX 4.05.4 (développé par Belkhir et al. et téléchargeable gratuitement à http://www.univ-montp2.fr/~genetix/genetix/genetix.htm) offre une interface extrêmement conviviale, en français qui plus est (assez rare pour être souligné), pour produire des AFC en deux ou même trois dimensions (pas nécessairement les plus faciles à lire en ce qui me concerne).

Figure 11

Recommandations et astuces pour les utilisateurs de l’AFC

113Quand on procède à une AFC (FCA ou FA en anglais), le programme génère différents fichiers tels que celui contenant les coordonnées des individus sur les différents axes. On peut être tenté d’utiliser ces coordonnées, qui sont donc des données ordinales continues issues de données qualitatives disjointes, pour procéder à des analyses de type analyse de variance (Anova) ou régression. Vérifier si les coordonnées des individus sur le premier axe de l’AFC sont expliquées plus ou moins bien par telle ou telle autre variable écologique peut en effet représenter une perspective séduisante. Je sais que beaucoup de personnes considèrent qu’il n’est pas valide de procéder à ce genre d’analyses à partir d’une AFC (alors qu’ils considèrent que cela est possible à partir d’une ACP, traitée plus loin) qui transforme des données discrètes bornées (0, 1 ou 2) en données continues de distribution incertaine. Mon opinion à ce sujet est que si on peut s’en passer on évite les ennuis, mais sinon je ne vois pas vraiment où est le problème à partir du moment où certaines précautions sont prises, comme de vérifier la distribution des données avant de procéder à une Anova.

114Il existe aussi une astuce à connaître par rapport au fait que ce type d’analyse est très sensible à la présence d’individus porteurs d’un allèle rare (outliers en anglais). En effet, les individus porteurs d’un allèle rare vont tirer le nuage vers eux. Le résultat est néfaste, car les autres individus se retrouveront compactés dans un nuage trop dense pour qu’on puisse y détecter quoi que ce soit. Cela va aussi remettre sérieusement en cause toute utilisation des coordonnées, car les coordonnées de chaque individu seront alors conditionnées majoritairement par la position de quelques individus exceptionnels. Il est souvent nécessaire de retirer plusieurs individus de l’analyse et parfois même un grand nombre. Dans l’échantillon de Nyafaro (Burkina Faso) dans Solano et al. (2000), près de 42 % des individus ont dû être ainsi écartés de l’analyse afin de pouvoir déceler une sous-structure dans les individus restants (60 sur les 97).

Analyse en composantes principales (ACP)

115Une ACP (PCA en anglais) suit le même principe que l’AFC sauf que ce sont des données ordinales continues qui sont utilisées au lieu de données disjonctives. Ici, ce sont des groupes d’individus (sous-échantillons) qui seront positionnés dans un hyperespace de K dimensions. Les coordonnées de chaque groupe sur chacun des axes principaux peuvent être utilisées pour des analyses statistiques supplémentaires telles que des analyses de variance ou autres régressions comme dans Nébavi et al. (2006). C’est une procédure fort utile pour positionner des sous-échantillons les uns par rapport aux autres en fonction de leur appartenance à un groupe écologique particulier comme des sous-échantillons de tiques d’oiseaux marins sur différentes espèces hôtes, comme on peut le voir dans la figure 12 (voir aussi McCoy et al., 2003, 2005).

Figure 12

116Le logiciel PCA-GEN ver. 1.2 (développé par J. Goudet librement téléchargeable à http://www2.unil.ch/popgen/softwares/pcagen.htm) permet cette analyse à partir de données au format Fstat (mais avec un format limité à deux caractères par allèle). Ce logiciel, en plus de fournir les graphiques en deux dimensions de la projection des points selon les axes demandés et leur pourcentage d’inertie, fournit également des tests de significativité de ces axes selon la méthode du bâton brisé (broken stick), une technique empirique appliquée à l’ACP (Frontier, 1976 ; Legendre et Legendre, 1998 ; King et Jackson, 1999) qui correspond davantage à un critère qu’à un test réel. Une explication plus détaillée de cette technique peut être consultée en réponse 9 à la fin de ce manuel. PCA-GEN propose aussi une procédure de permutations des génotypes complets entre sous-échantillons afin de tester la significativité de chaque axe (basé sur le pourcentage d’inertie). Comme seuls les génotypes complets sont permutés, il est donc important de disposer de jeux de données suffisamment complets si on souhaite que cette procédure ait un minimum de sens.

Analyse canonique des correspondances (ACC)

117L’ACC (CCA en anglais) est malheureusement implémentée actuellement par un logiciel commercial qui s’appelle CANOCO (Ter Braak, 1986, 1987 ; Ter Braak et Šmilauer, 2002). Il s’agit d’une méthode complexe d’ordination des données visant à directement corréler des tableaux de données multivariées. L’ordination des données couplées aux techniques de régression suivies de tests par permutation des données offre une méthode sophistiquée pour corréler les données génétiques à des variables environnementales. Elle offre également l’opportunité d’obtenir une projection en deux dimensions des centroïdes (barycentres) des données génétiques de chaque sous-échantillon défini, autour de laquelle une ellipse correspondant à l’intervalle de confiance à 95 % de cette projection peut également être dessinée. L’utilisation de l’ACC est rare, mais peut se montrer efficace ou au moins illustrative (Škalamera et al., 1999 ; Angers et al., 1999). Le logiciel ADE-4 permet aussi ce genre d’approches (Chessel et al., 2004). Je discuterai d’autres algorithmes disponibles dans un paragraphe spécifique (voir plus bas « Commentaires sur les algorithmes bayésiens de clusterisation »).

Construction d’arbres

118Construire des dendrogrammes censés relier les différents individus ou sous-échantillons en fonction de leur proximité génétique procure un moyen assez élégant et relativement simple de représenter les données génétiques suivant un schéma hiérarchique. Ce mode de représentation est d’ailleurs si populaire qu’innombrables sont les études qui l’utilisent. Un des champs d’application privilégié de la construction d’arbres peut être trouvé dans les études d’épidémiologie moléculaire d’organismes clonaux (voir Taylor et al., 1999 pour revue).

119Plusieurs méthodes existent. Pour des données type microsatellites ou allozymes, à cause de l’homoplasie, il n’est pas raisonnable d’espérer obtenir quelque chose ayant valeur phylogénétique. Selon Takezaki et Nei (1996), la méthode par NJTREE (neighbor-joining tree) basée sur une matrice de distances de corde (type Cavalli-Sforza et Edwards, 1967) paraît l’option la plus pertinente. Le logiciel MEGA 5 (Tamura et al., 2011a) (Kumar et al., 2004 ; Tamura et al., 2011b), librement téléchargeable de http://www.megasoftware.net/, offre une interface conviviale pour construire simplement un tel arbre à partir d’une demi-matrice de distances par paire. La méthode UPGMA, qui fait l’hypothèse d’une horloge moléculaire (les branches ont nécessairement la même longueur) est également très souvent utilisée. Je n’ai pas d’opinion bien tranchée sur la question et je doute que l’UPGMA donne souvent des résultats forts différents du NJTREE. Mais comme certains auteurs ont fait des analyses comparatives théoriques ayant conduit à mettre en avant NJTREE et les distances de corde, je préfère d’instinct suivre leur recommandation. L’avantage d’utiliser MEGA est que les arbres générés peuvent être sauvés dans le presse-papier (clipboard) et collés dans un logiciel de graphique où, après dissociation on peut travailler tranquillement la figure obtenue.

Trouver une sous-structure cachée

120Dans certaines situations, il n’existe aucun indice visible qui permettrait de subdiviser un échantillon en plusieurs sous-unités objectives. Dans de telles situations, les stratégies d’échantillonnage peuvent se montrer inefficaces à représenter une réalité biologique ou écologique existante. En effet, si un facteur biologique et/ou écologique contribue fortement à l’élaboration de l’architecture génotypique des individus étudiés, on s’attend à ce qu’un tel phénomène laisse une signature génétique visible dans l’échantillon sous la forme d’un déficit en hétérozygotes (effet Wahlund). Le site d’échantillonnage peut, par exemple, correspondre à une aire de nourrissage d’individus provenant d’unités de reproductions très séparées. Il existe différentes méthodes permettant de regrouper les individus d’un échantillon par affinité génétique en différents groupes (sous-populations putatives) en utilisant leur génotype multilocus. Par exemple, de très importants déficits en hérérozygotes sont souvent trouvés pour les loci microsatellites des mouches tsé-tsé. En particulier, Solano et al. (2000) sur Glossina palpalis gambiensis avaient mis en évidence de très importants F_IS. Ces forts déficits ne pouvant être expliqués facilement, et en tous les cas pas en totalité, par la présence d’allèles nuls ou par la dominance d’allèles courts (voir p. 110-116), une structure cachée a été recherchée. À l’époque où ces données ont été analysées, une AFC avait été réalisée et avait permis d’identifier des sous-groupes de tsé-tsé où le déficit en hétérozygotes se retrouvait inférieur au déficit initial (individus regroupés), ce qui confirmait que ces déficits provenaient bien d’un effet Wahlund (dont l’origine exacte reste elle-même à identifier) (fig. 13).

Figure 13

121D’autres méthodes, basées sur les statistiques pseudo-bayésiennes et des simulations de Monte-Carlo par chaîne de Markhov, sont maintenant disponibles. Elles permettent d’inférer la vraisemblance avec laquelle certains individus peuvent être regroupés et donc considérés comme appartenant à la même sous-population (cluster), ce qui peut conduire à la détection d’une structure cachée. Différentes méthodes avec différents logiciels sont disponibles. On peut citer ici les deux principaux que sont STRUCTURE de Pritchard et al. (2002) (Pritchard et al., 2000 ; Falush et al., 2003), librement téléchargeable à http://pritch.bsd.uchicago.edu/software/structure2_1.html, et BAPS 4 de Corander et al. (2006) (Corander et al., 2003, 2004 ; téléchargeable à http://www.rni.helsinki.fi/~jic/bapspage.html).

122BAPS a notamment permis de détecter une structure cachée dans deux échantillons (séparés dans le temps) de Glossina palpali palpalis de Bonon (Côte d’Ivoire) (Ravel et al., 2007).

123Dans les deux algorithmes (STRUCTURE et BAPS), l’hypothèse de panmixie est supposée dans chaque cluster que le logiciel cherche à construire. Cependant, la notion de panmixie telle qu’affirmée par les auteurs n’est pas claire et nous verrons que les clusters obtenus ne sont pas nécessairement conformes à Hardy-Weinberg. D’une manière générale, il est très difficile de savoir ce qui est fait et les différences entre ces méthodes mériteraient d’être mieux explorées dans différentes situations de populations structurées (y compris l’absence de structure) et pour différents systèmes de reproduction. Cela signifie qu’il ne faut en aucun cas être esclave du résultat fourni par ces méthodes et garder la tête froide en toute circonstance. D’une manière générale, BAPS est beaucoup plus facile d’utilisation et beaucoup plus rapide, mais produit davantage de clusters que ce qu’il y a en réalité (Latch et al., 2006). INSTRUCT (Gao et al., 2007) fonctionne comme STRUCTURE, mais avec prise en compte de l’autofécondation. Il ne marche qu’en ligne à http://cbsuapps.tc.cornell.edu/InStruct.aspx, ce qui peut poser un problème pour les pays mal connectés. Enfin, il s’agit d’un domaine en pleine évolution et lorsque ce manuel paraîtra, d’autres logiciels avec d’autres options seront parus. En particulier, dans la seconde partie de ce manuel, nous utiliserons aussi un autre de ces logiciels plus récent, FLOCK Duchesne et Turgeon (2009), qui bien que différent de BAPS, donne des résultats très comparables et est quasiment aussi facile d’utilisation. Avant analyse, j’engage les lecteurs souhaitant aller plus loin de vérifier sous Google, en recherchant “admixture AND population structure”, par exemple.

Commentaires sur les algorithmes bayésiens de clusterisation

124Dans cette réédition, il m’a semblé utile de rajouter un paragraphe dédié aux problèmes liés à l’utilisation des logiciels de clusterisation bayésienne, eu égard à leur extraordinaire popularité, et leur utilisation quasi systématique, et le plus souvent sans but précis. Il en existe aujourd’hui un nombre suffisant pour dérouter à peu près tout le monde. Il y a manifestement une niche fructueuse à exploiter en termes de publications et de citations. Je ne vais me focaliser que sur trois des plus populaires ou des plus aisés d’utilisation. À titre d’exemple, l’article original de l’algorithme et du logiciel STRUCTURE (Pritchard et al., 2000) a été cité 28 801 fois à ce jour (30/07/2020), soit 1 140 fois par an, et son complément (Evanno et al., 2005) 16 223 fois (1 081 fois par an). Ceux de BAPS (Corander et al., 2003) 899 fois, (Corander et al., 2004) 489 fois, et (Corander et al., 2006) 221 fois pour les principales. Enfin, DAPC (Jombart et al., 2010) a été cité 2 681 fois (268 fois par an).

125Un premier point à préciser est que les algorithmes utilisés sont très difficiles (voire impossible) à comprendre. Je soupçonne les auteurs eux-mêmes de ne pas bien maîtriser tous les résultats que peuvent générer leurs logiciels, en particulier chez des organismes non-modèles présentant des problèmes d’amplification divers (stuttering, dominance d’allèles courts, allèles nuls), des régimes de reproduction qui ne sont pas strictement panmictiques (dioécie, croisements fermés, clonalité), des modèles de structuration spatiale différents du modèle en îles discontinues, ou des effets temporels (présence d’individus appartenant à différentes cohortes). Pour STRUCTURE et BAPS, par exemple, l’hypothèse de base est que les vraies sous-populations sont à l’équilibre de Hardy-Weinberg et indemnes de tout déséquilibre de liaison, mais ils n’évaluent jamais la validité des partitions trouvées quand ces hypothèses ne sont pas respectées.

126Un second point qui me semble primordial est la question à quoi ces procédures sont censées répondre. C’est très loin d’être clair. Personnellement, si l’échantillonnage spatiotemporel dont je dispose est représentatif de ce qui se passe dans la nature, je ne vois absolument aucun intérêt à utiliser ce genre de méthodes. Avec les statistiques et estimateurs classiques, si je ne rencontre pas de problèmes (voir plus bas), je vais être en mesure de faire des inférences nettement plus fiables et parlantes que les approches bayésiennes en termes de système de reproduction, tailles des sous-populations, dispersion, biais de dispersion de certaines catégories d’individus, biais de structuration de certains groupes de sous-populations, ou description des différents niveaux de structuration hiérarchiques et/ou croisés. Il n’est pas rare de trouver dans des publications des calculs de F_ST et des tests de significativité entre clusters, ce qui n’a strictement aucun sens biologique. Je ne suis pas le seul à le constater (Meirmans, 2015). Toutes ces procédures visent en effet à ranger ensemble les individus les plus semblables génétiquement et à maximiser la différence entre groupes ainsi formés (clusters). En quoi sera-t-il étonnant de trouver une différenciation significative entre ces groupes ? Si je classe des boules rouges ensemble et des boules noires ailleurs, est-il vraiment pertinent de mesurer et tester la différence de couleur entre ces groupes ?

127De mon expérience, ces méthodes convergent rarement vers la même partition, ce qui à mon avis représente un vrai problème. Plus grave encore, pour un même logiciel, la partition peut varier en fonction des loci utilisés (Manangwa et al., 2019). L’existence de déséquilibres de liaisons, qui comme nous l’avons vu sont forcément présents partout dans la nature, sont connus pour affecter les résultats de STRUCTURE (Kaeuffer et al., 2007) et probablement de BAPS, et un schéma d’isolement par la distance conduit à des partitions erronées (Frantz et al., 2009).

128Vous allez me dire que je suis un Tartuffe, et vous aurez parfaitement raison, car j’ai moi-même utilisé ces procédures à plusieurs occasions. En fait, je n’utilise les méthodes de clusterisation bayésiennes que dans deux situations très précises. Pour résumer, comme je le dis dans mes cours, « je n’utilise les méthodes bayésiennes que quand je suis désespéré ». La principale situation, celle qui motive le plus souvent chez moi l’utilisation d’une de ces méthodes, est la présence de F_IS positifs inexpliqués dans mes sous-échantillons, éventuellement accompagnés d’une proportion prohibitive de paires de loci en déséquilibre de liaison significatif (par exemple 15-20 %). Dans ce cas j’utilise BAPS, non pas parce que les résultats attendus sont plus vrais que les autres, mais parce que c’est de loin le plus facile d’utilisation, tout simplement. Si je trouve une partition qui me permet de comprendre ce qui se passe, comme la coexistence avérée d’espèces cryptiques (Manangwa et al., 2019), je suis content. Si le résultat est plus ambigu, mais permet d’observer une baisse significative des F_IS et des déséquilibres de liaisons dans les clusters, et en fonction des configurations, je peux tenter de trouver des interprétations qui resteront spéculatives dans une plus (Ravel et al., 2007) ou moins (Chevillon et al., 2007) large mesure. Si la procédure ne permet pas d’interprétation biologique, il reste les problèmes techniques d’amplification corrigibles par les méthodes disponibles ou en rejetant les loci responsables (De Meeûs et al., 2019). Enfin, dans le cas où rien de ce qui précède ne fonctionne, il convient d’être extrêmement prudent dans les interprétations, voire d’abandonner le jeu de données. La seconde motivation vient des referees d’articles qui très souvent exigent (sans jamais donner de raison convaincante) que l’on procède à une analyse par STRUCTURE et, souvent également, avec DAPC, même lorsqu’une analyse BAPS est déjà présente. La plupart des referees en question considèrent en effet BAPS comme beaucoup moins performant, sans donner d’arguments convaincants. Dans ce cas, malgré mon caractère têtu qui n’est plus à démontrer, je préfère dorénavant présenter les trois méthodes pour avoir la paix, même si je trouve cette redondance parfaitement ridicule.

129Plus généralement, BAPS et DAPC seront plus efficaces pour trouver des partitions de groupes peu différentiés, alors que STRUCTURE va rechercher plutôt les entités fortement disjointes. STRUCTURE et ses dérivés sont particulièrement bien adaptés pour déterminer l’appartenance d’individus à de telles entités et caractériser leurs intermédiaires éventuels. Ils se montrent par exemple très efficaces pour classer deux espèces d’une zone hybride et calculer le degré d’introgression des différents hybrides, ou pour déterminer l’introgression ou l’appartenance à différentes races domestiques (Flori et al., 2019 ; Arteaga et al., 2020).

130Enfin, STRUCTURE n’est pas très simple à utiliser et doit être suivi d’une analyse par STRUCTURE Harvester (Earl et Vonholdt, 2012) avec à la clé environ 1 000 fichiers de sorties dont on ne sait que faire. DAPC est une procédure du package AdeGenet (Jombart, 2008) de R dont l’utilisation est pour moi un véritable cauchemar (documentation totalement obscure, succession de commandes nombreuses et parfaitement non intuitives). BAPS est très simple à utiliser et ne génère qu’un fichier de sortie, très simple à comprendre.

Estimer des effectifs efficaces

131Nous avons déjà vu, à propos des tests d’isolement par la distance en p. 91-92, que certains paramètres démographiques sont extrapolables à partir des données génétiques. Il existe d’autres situations où certains paramètres, tels que l’effectif efficace ou le taux de migration, peuvent être inférés à partir de données séparées dans le temps et/ou dans l’espace (sans isolement par la distance).

Définition de l’effectif efficace d’une population

132L’effectif efficace d’une population, aussi appelé effectif génétique et généralement noté N_e, est censé représenter avec quelle vitesse une population perd sa diversité génétique par dérive génétique. En effet, la fraction 1/N_e (1/(2N_e) pour des diploïdes) donne la probabilité que deux gènes pris au hasard dans la population descendent d’un seul et même gène ancêtre des parents. 1/N_e représente aussi la probabilité pour deux gamètes qui s’unissent de provenir du même parent. Ce phénomène est appelé coalescence. Cette coalescence récurrente de certains gènes implique donc que d’autres gènes ne contribuent pas au pool des générations suivantes. Certains gènes sont donc perdus, ce qui signifie également que la diversité génétique s’érode. Le ratio entre la taille réelle de la population N_c (aussi appelée taille de recensement ou census size en anglais) et cet effectif efficace représente donc une mesure de la dynamique de la diversité génétique d’une population focale par rapport à une population dite idéale. Une population idéale perdrait sa diversité génétique à la vitesse 1/N_c par génération de telle sorte que son effectif efficace serait égal à son effectif de recensement. Une telle situation correspond donc à des populations monoïques à générations non chevauchantes se reproduisant de façon panmictique dans un environnement constant sans sélection, migration ni mutation. À titre d’exemple, une population composée de N_c = 100 individus dioïques avec un seul mâle (N_m = 1) et 99 femelles (N_f = 99) aurait un effectif efficace de (Hartl et Clark, 1989 : 86) :

133Ce qui donne un N_e ≈ 4, soit 25 fois plus petit que la taille de recensement (N_c = 100). Ceux qui le souhaitent trouveront une démonstration de ceci en réponse 10. On comprend aisément qu’une telle population perd sa diversité à une vitesse très grande. D’autres facteurs influencent l’effectif efficace. En particulier, la subdivision des populations est susceptible d’augmenter l’effectif efficace d’une population, car une population subdivisée conservera en général mieux sa diversité génétique qu’une population homogène. Pour caricaturer, une population totalement subdivisée en sous-unités stables n’échangeant aucun migrant aura un effectif efficace infini, car la diversité génétique se trouvera rapidement fixée à l’échelle globale quand chaque sous-population se retrouve fixée pour l’un ou l’autre des allèles présents (et donc quand la diversité est totalement perdue à une échelle locale). Les lecteurs qui n’auront pas encore jeté ce manuel au feu trouveront une excellente revue de Criscione et Blouin (2005) sur le calcul des effectifs efficaces dans le cas des organismes parasites.

134Enfin, il faut également signaler qu’il existe plusieurs définitions de l’effectif efficace avec, dans la plupart des situations, des conséquences négligeables sur les valeurs estimées. Citons l’effectif efficace de consanguinité qui, comme son nom l’indique, mesure la probabilité d’identité par descendance des gènes, l’effectif efficace de variance qui s’attache à analyser la variation des fréquences alléliques (leur amplitude plus exactement) d’une génération à l’autre, l’effectif efficace de valeur propre qui s’intéresse à l’évolution de l’hétérozygotie au cours du temps à l’aide de certaines propriétés des matrices et de leur algèbre, et enfin l’effectif efficace de coalescence qui s’intéresse au temps qu’il faut mettre pour retrouver l’ancêtre commun de deux représentants du même gène dans une population.

Méthodes de calcul de l’effectif efficace des populations naturelles

135Deux familles de méthodes permettent d’inférer l’effectif efficace de populations étudiées, à l’aide de marqueurs moléculaires. Les études basées sur les fluctuations temporelles des fréquences alléliques, au cours des générations, permettent d’estimer ce que les spécialistes appellent l’effectif efficace de variance de populations échantillonnées de façon répétée au cours de leur cycle de vie (Waples, 1989). Le logiciel MACLEEPS 1.1 (Anderson et al., 2000) (téléchargeable à http://www.stat.washington.edu/thompson/Genepi/Mcleeps.shtml) réalise une estimation de N_e par maximum de vraisemblance en utilisant la variation des fréquences des allèles entre générations. Il est donc nécessaire de connaître le temps de génération de l’espèce étudiée. L’algorithme utilisé fait l’hypothèse que la sélection, la migration et la mutation ont un impact négligeable comparé à la dérive. Un intervalle de confiance à 95 % est également calculé (Anderson et al., 2000). Cela est également réalisé par le logiciel NeEstimator (Peel et al., 2004) (voir plus bas).

136Les études des variations spatiales des fréquences des allèles permet d’estimer ce que les spécialistes (encore eux) appellent l’effectif efficace de consanguinité. Les estimations produites par ESTIM 1.2 appartiennent à cette catégorie (Vitalis et Couvet, 2001a) (téléchargeable gratuitement : http://www.t-de-meeus.fr/ProgMeeusGB.html). Ce logiciel utilise l’information monolocus fournie par le F_ST et celle offerte par le déséquilibre d’identité par paire de loci afin de pouvoir en tirer N_e et m (effectif efficace et taux de migration) sans avoir besoin de connaître le modèle ni le taux de mutation ou de migration (Vitalis et Couvet, 2001b, c). Cependant, les taux d’autofécondation et de recombinaison entre les loci utilisés doivent eux être connus.

137La façon d’utiliser la méthode temporelle (effectif efficace de variance) et la méthode spatiale de Vitalis et Couvet (2001a) (effectif efficace de consanguinité), une comparaison ainsi que les problèmes possibles que l’on peut rencontrer en les utilisant, sont discutés dans Meunier et al. (2004b). ESTIM donne souvent des effectifs peu informatifs (0 ou infini). Il fournit également, quand le calcul est possible, les intervalles de confiance à 95 % des valeurs estimées.

138Il existe une autre méthode utilisant l’information des déséquilibres de liaison entre loci, dans une seule population (Bartley et al., 1992). Cette méthode, ainsi que celle de Waples (échantillons temporels), est implémentée par le logiciel NeEstimator Version 1.3 (logiciel non publié de Peel D., Ovenden J. R., Peel S. L., 2004, NeEstimator Version 1.3 : software for estimating effective population size. Queensland Government, Department of Primary Industries and Fisheries) téléchargeable gratuitement à http://www.dpi.qld.gov.au/fishweb/11629.html. Ce logiciel propose aussi une méthode basée sur les niveaux d’hétérozygotie observés (Luikart et Cornuet, 1999), mais à priori moins précise que la méthode proposée par Balloux (2004) (voir plus bas). Dans tous les cas, et quand le calcul est possible, le logiciel donne les intervalles de confiance à 95 %. La méthode de Bartley étant biaisée quand les échantillons sont de taille inférieure au N_e (England et al., 2006 ; Waples, 2006), Waples et Do proposent LDNe (Waples et Do, 2008) qui donne rarement des résultats utilisables, mais corrige le biais dû aux faibles échantillons et est donc utile pour vérifier qu’on ne sous-estime pas les N_e. Notons que la dernière version de NeEstimator (Do et al., 2014) propose plusieurs méthodes temporelles et non temporelles, avec corrections dans certains cas de l’effet des données manquantes (Peel et al., 2013). Balloux (2004) propose un estimateur corrigé par rapport à la méthode de Luikart et Cornuet et facile à calculer en utilisant l’estimateur de Weir et Cockerham : N_e = 1/(– 2F_IS) – F_IS/(1 + F_IS), dans le cas de populations dioiques, ou auto-incompatibles.

139Enfin, il est important de signaler qu’une approche synthétique, alliant les deux types d’informations (spatiale et temporelle), est également disponible (Wang et Whitlock, 2003). Un logiciel appelé MLNE estimant N_e et m en même temps peut être téléchargé gratuitement à partir de http://www.zoo.cam.ac.uk/ioz/software.htm.

Détection de goulots d’étranglement

140Ce paragraphe figure ici car les notions d’effectifs efficaces de goulot d’étranglement (bottleneck en anglais) et de biologie de la conservation sont étroitement liées. Une population qui subit une forte réduction d’effectif (goulot d’étranglement) va avoir tendance à présenter une réduction simultanée du nombre d’allèles par locus et de leur diversité génétique (que nous avons plus haut appelée H_s). Durant un goulot d’étranglement, le nombre d’allèles est réduit plus fortement que la diversité génétique. Il en résulte qu’une population ayant subi un goulot d’étranglement récent présentera une diversité génétique supérieure à celle attendue à l’équilibre mutation/dérive compte tenu du nombre d’allèles observés, sous l’hypothèse d’une taille constante de la population. Plusieurs modèles de mutation peuvent être utilisés selon les situations. Selon Cornuet et Luikart (1996), dans le cas de microsatellites il vaut mieux utiliser le modèle de mutation SMM ou à deux phases, même si c’est avec un IAM que la détection semble la plus aisée. Il s’agit donc de faire un choix entre ce qui, de la détection ou de la non-détection d’un goulot d’étranglement, est plus ou moins grave, ce qui dépend évidemment du contexte. La détection et les tests de significativité de cet excès d’hétérozygotie (signature d’un goulot d’étranglement) sont mis en œuvre dans le logiciel Bottleneck (Piry et al., 1997) (voir Cornuet et Luikart, 1996).

141Dans une population à l’équilibre mutation/dérive dont la taille n’a pas varié depuis un temps raisonnable, il y a autant de chance d’observer un excès qu’un déficit de diversité génétique, par rapport à l’attendu, aux différents loci. Afin de détecter si le nombre d’excès observé dépasse significativement ce qui est attendu sous cette hypothèse nulle, on peut utiliser trois tests (décrits par les auteurs dans l’aide du logiciel), mais le plus commode et le plus puissant est le test de Wilcoxon.

142Dans leur article, Cornuet et Luikart (1996) montrent (voir leur figure 3) que la détection d’une telle signature ne semble possible que dans certaines conditions, qui dépendent du degré de polymorphisme observé, du nombre de générations écoulées depuis le dernier goulot d’étranglement (qu’on cherche à détecter) et de l’effectif efficace de la population (celui qu’elle acquiert après l’événement de goulot d’étranglement). Par exemple, avec des loci raisonnablement polymorphes (microsatellites), des échantillons inférieurs à 40 individus et moins de 10 loci, la détection d’un goulot d’étranglement n’est possible que si ce dernier a eu lieu dans une fourchette de temps définie par les limites 0,025 x 2 x N_e et 2,5 x 2 x N_e générations et où N_e représente l’effectif efficace qui s’est mis en place après le goulot d’étranglement. Par conséquent, la connaissance de ce temps de générations τ depuis le dernier goulot d’étranglement probable peut offrir une manière détournée d’estimer une fenêtre probable pour N_e. Ici, cette fenêtre serait de [τ/5N_e, τ/0,5N_e]. C’est ce principe qui a permis d’estimer grossièrement les effectifs efficaces probables de la tique du bétail Rhipicephalus (Boophilus) microplus récemment introduite en Nouvelle-Calédonie comme très grands malgré des traitements acaricides soutenus dans les élevages bovins de l’île (Koffi et al., 2006a).

143Enfin, il n’est pas inutile de signaler ici que le test de Bottleneck aura tendance à donner des résultats légèrement significatifs quand les populations étudiées sont de petites tailles. Dans ce cas il est utile, voire indispensable, d’obtenir des informations sur la taille des populations étudiées, par l’utilisation de méthodes d’estimation d’effectifs efficaces, par exemple. Par ailleurs, il ne faudra pencher en faveur d’un goulot d’étranglement que si les P-values sont très significatives et/ou si au moins deux, et encore mieux les trois, modèles de mutation convergent vers cette conclusion.

144L’utilisation de plusieurs méthodes pour estimer N_e pourra être d’un grand secours pour convaincre les referees toujours réticents (si les valeurs obtenues par différentes méthodes convergent, bien entendu).

Le cas spécial des allèles nuls

Présentation générale

145Les allèles nuls correspondent à des allèles qu’on ne peut pas déceler avec la méthode de détection biochimique utilisée. Ils sont invisibles à l’état hétérozygote, car récessifs par rapport aux autres allèles, et mal détectés à l’état homozygotes (blancs), car il est souvent difficile de séparer les cas où la manipulation a échoué (mauvaise amplification, matériel dégradé, etc.) des cas où on a réellement à faire à un homozygote nul. Les allèles nuls sont fréquemment rencontrés dans les études de génétique des populations naturelles, bien que fréquemment ignorés. Il est même probable que bon nombre de déficits en hétérozygotes documentés dans de nombreux articles soient en fait dus à ce phénomène, alors que d’autres causes sont privilégiées dans les articles en question. Les allèles nuls peuvent être fréquents même dans le cas des allozymes (Gaffney, 1994 ; Nébavi et al., 2006), où on ne les attend pourtant guère, car ils correspondent dans cette circonstance à des enzymes non fonctionnels, bien qu’indispensables à la vie (pour la plupart). On augurerait donc ici une moindre valeur sélective des allèles nuls, à moins qu’un mécanisme permette de les garder à l’état hétérozygote le plus fréquemment possible, comme cela peut être le cas chez les organismes clonaux (Nébavi et al., 2006). C’est un problème rencontré typiquement chez les marqueurs microsatellites (Paetkau et Strobeck, 1995 ; Pemberton et al., 1995 ; Brookfield, 1996). Une mutation dans la séquence flanquante, au niveau des séquences correspondant à un des primers, empêche la bonne amplification de cet allèle. Il apparaîtra « blanc » (aucun signal) à l’état homozygote et sera dominé par les allèles avec lesquels il sera hétérozygote. Les hétérozygotes pour ce type d’allèles apparaissent donc homozygotes pour l’autre allèle.

146Une conséquence importante de la présence d’allèles nuls, en plus d’une surestimation du déficit local en hétérozygotes, sera une surestimation des mesures de subdivisions (Chapuis et Estoup, 2007). Nous verrons que ce biais peut être facilement et efficacement corrigé par le logiciel FeeNA (Chapuis et Estoup, 2007).

Détecter la présence d’allèles nuls

147Nous savons maintenant que la présence d’allèles nuls à un locus va provoquer des déficits en hétérozygotes inexplicables biologiquement. Normalement, dans ce cas, on s’attend à ce que les différents loci donnent une mesure différente (variance entre loci) (De Meeûs et al., 2002a ; Hurtrez-Boussès et al., 2004). Par ailleurs, s’il y a structuration entre sous-échantillons, la fréquence de ces allèles nuls, aux loci concernés, devrait changer d’un sous-échantillon à l’autre et provoquer une variance des déficits (F_IS) entre sous-échantillons, mais seulement pour des niveaux de différenciation génétique élevés. Ensuite, il existe des procédures plus ou moins complexes pour estimer, à chaque locus et dans chaque sous-échantillon, la fréquence d’allèles nuls nécessaires pour expliquer les déficits observés (Brookfield, 1996). Le logiciel Micro-checker V 2.2.3. (Van Oosterhout et al., 2004), téléchargeable librement de http://www.microchecker.hull.ac.uk/), permet de faire ces estimations pour chaque locus et chaque sous-échantillon. Ces fréquences estimées d’allèles nuls peuvent ensuite permettre d’évaluer la proportion attendue d’individus blancs, sous l’hypothèse de panmixie et si ces allèles nuls expliquent la totalité du déficit.

Trucs et astuces pour tester la présence des allèles nuls

148Il est important d’insister encore sur le fait que, si tous les loci convergent vers le même déficit en hétérozygotes (tous présentent un F_IS comparable à celui des autres), il n’est alors pas nécessaire d’invoquer les allèles nuls, mais plus parcimonieusement une cause biologique (autofécondation, effet Wahlund). Dans ce qui suit, nous partons donc du principe qu’une forte variance entre loci a été observée.

149Il faut tout d’abord savoir que Micro-checker est conçu spécifiquement pour les microsatellites. Avant de procéder à l’analyse, il est demandé le type de motif pour chaque marqueur. Si vous n’êtes pas sûr de vos données (quelques mutants atypiques d’un pas différent du motif de base), il vaut en général mieux adopter l’option mononucléotidique pour tous les loci. Cependant, si vous avez défini le motif de certains de vos loci comme mononucléotides, il faudra s’en souvenir au moment d’interpréter les tests de stuttering, en observant les déficits à une et deux répétitions sur les graphiques de sorties de MicroChecker. Ensuite, parmi les résultats que propose Micro-Checker, il faut garder, pour chaque locus et chaque sous-échantillon, la fréquence des allèles nuls, ainsi que la présence ou non de stuttering. Il vaut mieux utiliser la méthode 2 de Brookfield (1996) qui tient compte des données manquantes (doubles nuls). La fréquence attendue d’homozygotes blancs sera, sous l’hypothèse de croisement au hasard, égale à p_nul². On peut comparer par un test binomial cette fréquence attendue aux nombres de blancs effectivement observés à ce locus dans le sous-échantillon concerné. Ce test permet de vérifier si les allèles nuls expliquent raisonnablement les déficits observés aux loci concernés. On peut utiliser le logiciel R qui est gratuit avec la commande binom.test.

150MicroChecker ne peut pas travailler avec des sous-échantillons trop petits (il renvoie un message d’erreur dans ce cas). Si la plupart de vos sous-échantillons renvoient un message d’erreur dû à la faiblesse des effectifs, il existe encore une alternative. Vous pouvez effectuer la régression du F_IS observé par locus et sous-échantillon en fonction du nombre de blancs observés par locus et sous-échantillon (ou plus simplement par locus). Si le test de corrélation de Spearman correspondant est significatif, c’est que les allèles nuls expliquent une partie du F_IS. Cette méthode, en calculant le R² de la régression (proportion de la variance expliquée par la régression), permet aussi d’appréhender à quel point les allèles nuls expliquent les données. En cas d’allèles nuls, on attends également une corrélation positive entre F_IS et F_ST et une erreur standard de jackknife du F_IS au moins deux fois plus grande que celle du F_ST. Tout ces trucs et astuces sont détaillés dans plusieurs de mes articles (Séré et al., 2017 ; De Meeûs, 2018 ; De Meeûs et al., 2019 ; Manangwa et al., 2019), y compris dans le cas spécial des organismes clonaux (Séré et al., 2014).

151Hormis chez les clones, toutes ces méthodes font l’hypothèse qu’il y a grosso modo panmixie. Si la fréquence de nuls ne permet pas d’expliquer correctement tous vos déficits en hétérozygotes, en particulier si aucun locus n’est conforme à l’attendu panmictique, il se peut que d’autres phénomènes soient en cause. Si par exemple, il y a de l’autofécondation ou effet Wahlund, en plus des allèles nuls, nous ne pourrons pas expliquer les données à l’aide des seuls allèles nuls. Dans le cas de l’autofécondation, il existe un logiciel permettant d’estimer le taux d’autofécondation en tenant compte des allèles nuls ou autre problème (dominance partielle, dominance des allèles courts). Il s’agit de RMES (David et al, 2007), qui fait l’hypothèse d’équilibre de liaison entre loci et utilise les déséquilibres d’hétérozygotie par paire de loci, ce qui peut poser un problème dans les petites populations très autofécondantes. Le logiciel est librement téléchargeable à http://www.cefe.cnrs.fr/genetique-et-ecologie-evolutive/patrice-david.

152Cependant, précisons que les méthodes de régressions décrites ci-dessus doivent fonctionner quel que soit le système de reproduction.

153Finalement, la seule méthode pour avoir une estimation fiable du déficit en hétérozygotes local de base réel des sous-échantillons (F_IS) est soit d’éliminer les loci touchés et recalculer ce déficit, soit de prendre l’intercept de la régression du F_IS par le nombre de données manquantes par locus ou leur fréquence. Pour le F_ST ou la distance de corde de Cavalli-Sforza et Edwards (1967), les algorithmes ENA et INA respectivement, implémentés dans le logiciel FreeNA (Chapuis et Estoup, 2007), rétablissent assez fidèlement la valeur réelle de ces estimateurs, tant que la fréquence de ces allèles nuls ne dépasse pas 40 % (Séré et al., 2017).

Le cas très spécial de la dominance des allèles courts

Point de vue théorique

154La dominance des allèles courts, ou « short allele dominance » ou encore « large allele dropout », (Wattier et al., 1998 ; De Meeûs et al., 2004a), peut survenir plus ou moins fréquemment en fonction du modèle biologique. La logique qui se cache derrière ce terme est la suivante. Si, par un mécanisme qui reste à mettre en évidence, une compétition existe, au cours de la PCR, entre les deux portions d’ADN correspondant aux deux allèles d’un même locus devant être amplifiés, alors il semble logique que ce soit l’allèle le plus court (s’il y a une différence de taille entre les deux, bien entendu) qui sera le mieux amplifié.

155Une tentative de modélisation du phénomène peut être trouvée dans De Meeûs et al. (2004a). Dans ce modèle, on suppose une population panmictique et un locus pour lequel la PCR favorisera l’amplification de l’allèle le plus court de façon proportionnelle à la différence de taille qui l’oppose à l’autre allèle, ainsi qu’à un paramètre α variant entre 0 (pas de dominance) et 1 (dominance totale des allèles courts). Si les allèles existant à ce locus se rangent dans un ordre de tailles croissantes de s₁ à s_n et que la fréquence d’un allèle quelconque i est notée p_i, on peut poser que la fréquence observée d’hétérozygotes pour cet allèle avec un autre allèle j sera de :

156et

157Dans les équations (59) et (60), on voit bien que le biais sera maximal pour l’écart maximal de taille, c’est-à-dire pour un hétérozygote pour les allèles 1 et n, et minimal entre deux allèles les plus proches. En utilisant ces deux équations, on en déduit que la proportion observée d’hétérozygotes pour l’allèle i avec tous les autres allèles sera de :

158soit :

159L’équation (62) nous donne donc l’hétérozygotie attendue sous panmixie moins la proportion des génotypes hétérozygotes erronément interprétés comme homozygotes pour le plus court des allèles. La proportion d’homozygotes observés pour l’allèle i sera donc celle attendue sous panmixie plus la proportion d’individus portant l’allèle i et un allèle plus long que j et interprétés comme homozygotes pour cet allèle. Cette homozygotie observée sera donc de :

160L’utilisation des équations (1), (62) et (63) nous permet alors d’estimer la fréquence erronément observée de l’allèle i dans l’échantillon comme :

161En utilisation l’équation (6), nous pouvons alors estimer le déficit artificiel d’hétérozygotes observé par rapport aux attendus panmictiques :

162Dans la figure 14, il n’est pas inintéressant de constater que la relation entre taille des allèles et leur déficit en hétérozygotes n’est ni linéaire ni monotone et dépend de la distribution des fréquences des allèles (le tableau 5 décrit les différentes distributions utilisées), mais globalement on s’attend quand même à observer une décroissance de F_IS en fonction de la taille des allèles.

163On constate aussi que ce phénomène modifie également l’estimation des fréquences des allèles.

Figure 14

Du point de vue pratique : détection de la dominance des allèles courts

164Il existe une procédure de détection de la dominance des allèles courts dans le logiciel Micro-Checker, appelée ici « large allele drop-out ». Mais cette procédure ne teste le phénomène que dans chaque sous-échantillon pris séparément. Il en résulte un manque de puissance. On peut tester l’existence d’une dominance d’allèles courts sur l’ensemble des sous-échantillons en utilisant une approche de corrélation et/ou de régression.

Méthode fastidieuse de régression multiple

165C’est la méthode ancienne que j’ai développée dans mon article de 2004 (De Meeûs et al., 2004), que je préfère éviter aujourd’hui car plus longue, fastidieuse et possiblement moins robuste (test paramétrique).

166Tout d’abord, pour le locus étudié, il faut récupérer le F_IS de chaque allèle dans chaque sous-échantillon. Fstat ne le fait malheureusement pas automatiquement. Il faut créer autant de fichiers Fstats qu’il y a de sous-échantillons et, dans chacun de ces fichiers, il faut créer une population fictive fixée (un seul allèle présent) pour les locus dont on veut les F_IS par allèle. Fstat n’aime en effet pas travailler sur une seule population. Une fois qu’on a fait calculer ces F_IS par le logiciel, on a tout ce qui est nécessaire pour effectuer une régression linéaire généralisée ou GLiM. GLiM (Generalised Linear Model) est une forme de régression qui permet d’analyser des données de n’importe quelle forme (gaussiennes, poissoniennes, logistiques pour les plus utilisées) en fonction de n’importe quel type de variable (facteur catégoriel, logique, ordinal discontinu ou continu). Cette régression doit donc être de la forme F_IS = S + T + Cte, avec S pour le sous-échantillon, T la taille de l’allèle et Cte une constante. La régression est idéalement pondérée par le degré de polymorphisme de chaque allèle : p_i(1 – p_i), où p_i est la fréquence de l’allèle. Une analyse de variance sur le modèle permet ensuite de tester si l’effet de la taille des allèles, corrigé de l’effet des sous-échantillons, est significatif ou non. Attention, il faut que la relation entre taille des allèles et F_IS soit négative. Les relations positives doivent donc être ignorées.

167Nous verrons tout ceci en détail dans la mise en pratique de toutes ces connaissances dans la deuxième partie.

Méthode rapide de corrélation

168C’est la méthode du F_IT que j’ai adoptée dernièrement et qui est décrite dans un article de 2019 (Manangwa et al., 2019). Cette méthode est très rapide, elle fait moins d’hypothèses sur la distribution des données et procède à un test unilatéral (corrélation attendue négative), donc plus puissant que la régression (pente différente de 0 donc bilatérale). Il s’agit de prendre le F_IT (plus puissant que le F_IS) par allèle calculé sur l’ensemble des sous-échantillons et de procéder à une corrélation non-paramétrique de Spearman. Il faut quand même vérifier que le résultat (significatif ou non) ne dépend pas de quelques allèles rares. En cas de doute, pour le vérifier, on peut procéder à une régression du F_IT (voire du F_IS pour être sûr) en fonction de la taille des allèles pondérée par leur polymorphisme p_i(1 – p_i). Comme on réalise autant de tests qu’il y a de loci testés, il peut être avisé de procéder à une correction de BH. Mais c’est à vous de décider si vous souhaitez être certains de corriger ou d’exclure les loci avec possible dominance des allèles courts.

169Il existe deux méthodes pour régler le problème. La première consiste à retourner vers les profils homozygotes des loci incriminés et à essayer de détecter des micro-pics d’allèles plus longs que celui déjà vu. Cette méthode peut entièrement corriger le problème. Comme il n’existe pas de méthode analytique pour le moment, la seconde méthode consiste à simplement éliminer les loci non corrigibles. Il ne faut pas hésiter à le faire car la dominance d’allèles courts biaise potentiellement tous les estimateurs de génétique des population (Manangwa et al., 2019 ; De Meeûs et al., 2019).

Le cas du « stuttering »

170Le mot anglais « stuttering » se traduit par bégaiement. Si les amorces de la PCR ne s’accrochent pas très bien, il peut y avoir glissement plus ou moins fréquent de la Taq et ajout ou retrait d’un motif dans le produit amplifié. La récurrence du phénomène tend à produire des produits de différentes tailles, dont la fréquence diminue avec l’écart par rapport à l’allèle de départ (voir la figure 4 de l’article De Meeûs et al., 2019a), avec un avantage pour les produits plus courts (Shinde et al., 2003). Il doit donc exister une relation entre déficit observé et différence de taille entre allèles, le déficit devant être plus important pour les individus hétérozygotes pour des allèles de tailles proches. Ceci est détecté automatiquement par Micro-Checker. Ici, il sera utile de porter le nombre de randomisations de MicroChecker à 10 000, d’ignorer les sorties textes et de ne regarder que les sorties graphiques pour déterminer s’il y a stuttering. Il y a stuttering significatif quand la fréquence observée d’hétérozygotes entre allèles d’une seule répétition de différence est significativement en dessous de l’intervalle de confiance de cette valeur. Dans le cas où le locus microsatellite est imparfait, vous serez obligés de le signaler comme mononucléotide. Dans ce cas, il faudra regarder les hétérozygotes entre allèles différant d’une répétition et deux répétitions (pour un dinucléotide) pour tester le stuttering. Un peu plus de détails sont consultables dans mon article de 2019 (De Meeûs et al., 2019a).

171Comme décrit en détail dans De Meeûs et al. (2019a), une méthode qui peut s’avérer très efficace pour corriger l’effet du stuttering est le regroupement d’allèles de taille proche pour fabriquer des allèles synthétiques (AS), en prenant garde de ne pas regrouper dans un même AS uniquement des allèles rares (de fréquence inférieure à 0,05) (afin d’éviter des excès artificiels d’hétérozygotes). Si cette procédure est efficace (chute significative du F_IS), il faut garder ce re-codage pour la suite des analyses, sinon il faut garder le codage initial.