Chapitre 2 Adaptation locale

Cette première partie fera dans un premier temps un état de l’art des méthodes destinées à identifier des locus impliqués dans des processus d’adaptation locale. Nous présentons différentes méthodes classiques de scan génomique pour la sélection. Ensuite, nous présentons l’utilisation de l’Analyse en Composantes Principales en génétique des populations et nous montrons comment l’utiliser pour faire des scans à sélection. Par souci de clarté, nous ne considérons ici que des espèces diploïdes, bien qu’une grande partie des résultats présentés ici puisse être adaptée au cas d’espèces haploïdes. Les locus seront par ailleurs supposés bi-alléliques, c’est-à-dire que pour un locus donné, au plus deux allèles sont observés sur ce locus à l’échelle de la population étudiée.

2.1 L’état de l’art pour les scans génomiques

2.1.1 Modèles démographiques

Afin de mieux comprendre l’heuristique des méthodes de scan génomique présentées ici, nous donnons dans ce paragraphe une brève description des modèles démographiques fréquemment utilisés en génétique des populations. En effet l’idée de sélection dans une population est généralement relative à (au moins) une autre population, et l’histoire démographique de ces populations joue un rôle important sur la distribution théorique des fréquences alléliques.

2.1.1.1 Modèle en îles

Dans un modèle en îles, les différentes populations échangent entre elles des individus au cours du temps (Figure 2.1). La proportion d’individus échangés est appelée taux de migration. De forts taux de migration vont avoir tendance à homogénéiser les variations génétiques entre les populations. Des faibles taux de migration vont en revanche conduire à une différenciation plus forte. Les différences de taux de migration peuvent par exemple être expliquées par l’existence de barrières naturelles (Landguth, Cushman, Murphy, & Luikart, 2010).

2.1.1.2 Modèle star-like

Le modèle star-like suppose l’existence d’une population ancestrale de laquelle sont issues différentes populations (Figure 2.1, panels B et C). Contrairement au modèle en îles, les populations évoluent de façon indépendante sans s’échanger d’individus, et se différencient éventuellement sous l’effet de la dérive génétique et de mutations aléatoires. Il existe un modèle de divergence instantanée dans lequel la quantité de dérive est la même pour toutes les branches de l’arbre de populations (Figure 2.1, panel C).

Modèles démographiques. A. Représentation schématique d’un modèle en îles à trois populations. B. Représentation schématique d’un modèle star-like à trois populations. Les longueurs de branches correspondent à la dérive génétique depuis la divergence initiale et sont ici différentes les unes des autres. Dans un modèle de divergence, la quantité de dérive est proportionnelle au temps de divergence divisé par la taille efficace de la sous-population. C. Représentation schématique d’un modèle star-like à trois populations où les trois branches ont subi la même quantité de dérive génétique. D. Représentation schématique d’un modèle de divergence présentant une structure hiérarchique.

Figure 2.1: Modèles démographiques. A. Représentation schématique d’un modèle en îles à trois populations. B. Représentation schématique d’un modèle star-like à trois populations. Les longueurs de branches correspondent à la dérive génétique depuis la divergence initiale et sont ici différentes les unes des autres. Dans un modèle de divergence, la quantité de dérive est proportionnelle au temps de divergence divisé par la taille efficace de la sous-population. C. Représentation schématique d’un modèle star-like à trois populations où les trois branches ont subi la même quantité de dérive génétique. D. Représentation schématique d’un modèle de divergence présentant une structure hiérarchique.

2.1.2 L’indice de fixation

Indépendamment de l’histoire démographique, un allèle sélectionné voit généralement sa prévalence augmenter au sein d’une population. Si bien que l’observation d’une fréquence allélique anormalement élevée dans une population relativement aux autres donne à suggérer que cet allèle a été favorablement sélectionné. Dans l’optique de détecter des signaux de sélection, il semble alors naturel de proposer une statistique testant si au moins deux populations présentent des fréquences alléliques significativement différentes l’une de l’autre (Holsinger & Weir, 2009). L’indice de fixation, ou encore $F_{ST}$ en abrégé, est une statistique basée sur cette heuristique.

Définition 2.1 (Indice de fixation) Pour un locus donné, dénotant $N$ le nombre de populations considérées et ($p_1, p_2, \dots, p_N$) les fréquences d’un des deux allèles existant, la $F_{ST}$ est définie par la relation suivante (Wright, 1943) :

\[\begin{equation} F_{ST} = \frac{\frac{1}{N-1}\sum_{i = 1}^N (p_i - \bar{p})^2}{\bar{p}(1 - \bar{p})} \end{equation}\] où $\bar{p} = \frac{1}{N-1}\sum_{i = 1}^N p_i$.

Dans l’expression de la $F_{ST}$ donnée en 2.1, le numérateur correspond à la variance génétique interpopulationnelle tandis que le dénominateur correspond à la variance génétique mesurée à l’échelle de la métapopulation⁶. La $F_{ST}$ peut être vue comme la réduction de variance génétique due à la structure de populations. Autrement dit, on regarde si le fait de grouper les individus dans des populations a une incidence ou non sur la variance génétique.

2.1.3 Test de Lewontin-Krakauer

Comme expliqué plus haut, la détection de locus sous sélection passe généralement par la caractérisation d’un modèle décrivant l’évolution de locus sous l’effet de processus neutres tels que la dérive génétique. Dans le cas des scans génomiques, il s’agit d’estimer la distribution neutre de la statistique de test (calculée en chaque locus), afin d’identifier les locus qui s’en écartent le plus. Suivant ce principe, Lewontin et Krakauer proposent pour la $F_{ST}$ un test d’adéquation du $\chi^2$ (Lewontin & Krakauer, 1973).

Définition 2.2 (Test de Lewontin-Krakauer) Notant $N$ le nombre de populations, la statistique de test introduite par Lewontin & Krakauer (1973), dénotée $T_{LK}$, a pour expression :

\[\begin{equation} T_{LK} = \frac{N - 1}{\bar{F}_{ST}} F_{ST} \end{equation}\]

Dans un scénario de divergence instantanée (Figure 2.1), sous l’hypothèse que les fréquences alléliques sont distribuées selon une loi normale ou binomiale, $T_{LK}$ suit une loi du $\chi^2$ à $N-1$ degrés de libertés (Lewontin & Krakauer, 1973). En effet, notant $p = (p_1, p_2, \dots, p_N)$, en utilisant la définition 2.1 de la $F_{ST}$ et en remarquant que :

\[\begin{equation} (N-1)F_{ST} = \frac{1}{\bar{p}(1 - \bar{p})}\sum_{i = 1}^N (p_i - \bar{p})^2 = \left(\frac{p - \bar{p}}{\sqrt{\bar{p}(1-\bar{p})}}\right) \left(\frac{p - \bar{p}}{\sqrt{\bar{p}(1-\bar{p})}}\right)^T, \tag{2.1} \end{equation}\]

il est possible de réécrire $T_{LK}$ sous la forme d’une somme quadratique de lois normales. Cependant, les contraintes sur le modèle démographique sous-jacent sont extrêmement fortes et dans certaines situations elles ne seront pas vérifiées. Par exemple, dans les exemples de la figure 2.1, l’approximation $\chi^2$ est correcte pour le modèle en îles (Figure 2.1, panel A) et pour le modèle star-like (Figure 2.1, panel C). En revanche, pour le panel C où la longueur des branches est différente, l’approximation $\chi^2$ n’est plus correcte parce que les variances des quantités $p_i - \bar{p}$ ne sont pas identiques pour différentes valeurs de $i$. Des variantes de ce test ont donc été proposées pour s’adapter à des modèles de structure de populations plus flexibles (Bonhomme et al., 2010; Excoffier, Hofer, & Foll, 2009; Whitlock & Lotterhos, 2015).

2.1.3.1 Estimation du nombre de degrés de liberté effectif

Une manière de s’adapter à des modèles plus flexibles est de garder la statistique de test de l’équation (2.1) et d’améliorer l’approximation $\chi^2$. Pour améliorer l’approximation $\chi^2$, Whitlock & Lotterhos (2015) proposent d’approcher la statistique $T_{LK}$ avec un $\chi^2$ à $\text{df}$ degrés de libertés au lieu de $N-1$ degrés de libertés où $\text{df}$ est un paramètre à estimer. Pour estimer $\text{df}$, Whitlock & Lotterhos (2015) dérivent un modèle de vraisemblance basé sur la distribution de la $F_{ST}$. Partant de la densité d’une variable aléatoire suivant un $\chi^2$ à $\text{df}$ degrés de libertés, la densité de la $F_{ST}$ s’écrit :

\[\begin{equation} f(F_{ST}) = \frac{\text{df}}{\bar{F}_{ST}} \times \frac{1}{2^{\frac{\text{df}}{2}}\Gamma\left(\frac{\text{df}}{2}\right)} \times \left(\frac{\text{df}}{\bar{F}_{ST}}F_{ST}\right)^{-1+\frac{\text{df}}{2}} \tag{2.2} \end{equation}\]

L’estimation de $\text{df}$ se fait en maximisant la fonction de log-vraisemblance $\sum_{i=1}^p \log(f(F_{ST}^i))$ (où $F_{ST}^i$ désigne la $F_{ST}$ observée pour l’allèle $i$). La correction du test de Lewontin-Krakauer consiste alors à tester l’adéquation de $\frac{\text{df}}{\bar{F}_{ST}}F_{ST}$ à un $\chi^2$ à $\text{df}$ degrés de liberté.

2.1.3.2 Dérivation du test de Lewontin-Krakauer dans le cas de populations structurées

Dans le cas de structure hiérarchique (Figure 2.1, panel D), les quantités $p_i - \bar{p}$ ne sont plus indépendantes entre elles, ce qui remet en cause l’approximation $\chi^2$. Pour tenir compte de la structure hiérarchique, Bonhomme et al. (2010) proposent de corriger la statistique $T_{LK}$ pour l’apparentement génétique des populations, modélisé par une matrice $\mathcal{F} = (f_{ij})_{1 \leq i,j \leq N} \in \mathcal{M}_N(\mathbb{R})$, où $f_{ij}$ mesure la corrélation entre $p_i$ et $p_j$ et peut être interprétée comme la probabilité qu’un individu de la population $i$ et un individu de la population $j$ aient hérité de cet allèle d’un même ancêtre commun.

La statistique de test $T_{F-LK}$, correspondant au test de Lewontin-Krakauer corrigé pour l’apparentement génétique $\mathcal{F}$, garde une forme analogue à celle développée en (2.1) :

\[\begin{equation} T_{F-LK} = \left(\frac{p - \hat{p}_0}{\sqrt{\hat{p}_0(1-\hat{p}_0)}}\right) \mathcal{F}^{-1} \left(\frac{p - \hat{p}_0}{\sqrt{\hat{p}_0(1-\hat{p}_0)}}\right)^T \tag{2.3} \end{equation}\]

où $\hat{p}_0$ désigne un estimateur de la fréquence $p_0$ de l’allèle dans la population ancestrale. Dans le cas où les locus sont soumis uniquement à la dérive génétique, $T_{F-LK} \sim \chi_{N - 1}^2$ (Bonhomme et al., 2010). En pratique, la méthode implémentée dans le logiciel hapflk cherche à estimer les paramètres $p_0$ et $\mathcal{F}$.

2.1.4 Le modèle $F$

Une autre idée consiste à affirmer qu’en l’absence de sélection, dans un modèle en îles (Figure 2.1, panel A), la proportion d’allèles immigrants⁷ doit être la même pour tous les locus (Beaumont & Balding, 2004). Cette proportion d’allèles immigrants mesure la dérive génétique subie par la population qui intègre ces allèles (Villemereuil & Gaggiotti, 2015). En effet, une population échangeant moins d’allèles avec les autres populations se retrouverait alors plus différenciée. Les locus susceptibles d’être sous sélection sont ceux présentant une proportion d’allèles migrants anormalement basse (Bazin, Dawson, & Beaumont, 2010; Petry, 1983). Beaumont & Balding (2004) proposent alors un modèle de régression logistique pour la $F_{ST}$ (Balding, 2003; Balding, Bishop, & Cannings, 2008), en la modélisant par des effets $\alpha$ spécifiques au locus (taux de mutation, sélection) et des effets $\beta$ spécifiques à la population (taille de population efficace, taux de migration), c’est le modèle $F$ :

\[\begin{equation} \log \left( \frac{F_{ST}}{1 - F_{ST}} \right) = \alpha + \beta \tag{2.4} \end{equation}\]

La décomposition de la $F_{ST}$, en une somme d’effets locus-spécifique et population-spécifique, permet d’identifier les locus à forte $F_{ST}$ qui présentent un effet qui n’est pas partagé par les autres locus. L’un des défauts de ce modèle est qu’il ne prend pas en compte la structure hiérarchique (Foll, Gaggiotti, Daub, Vatsiou, & Excoffier, 2014). Le modèle $F$ est implémenté dans le logiciel Bayescan.

La liste des méthodes présentées ici n’est pas exhaustive, mais elle met en avant plusieurs défauts qui sont communs à la plupart des méthodes de scan génomique.

Dans le cas de la $F_{ST}$, la nécessité de travailler avec des fréquences alléliques populationnelles impose de travailler à l’échelle des populations plutôt qu’à l’échelle des individus. Les conséquences directes d’une telle nécessité sont :

l’assignation arbitraire de chaque individu à une population. La présence d’individus métissés peut s’avérer problématique.
la supposition que la structure de populations n’est pas continue.

Par ailleurs, les méthodes de scan génomique basées sur des méthodes bayésiennes telles que Bayescan ou la première version de PCAdapt (Duforet-Frebourg, Bazin, & Blum, 2014; Foll & Gaggiotti, 2008) sont connues pour être computationnellement lourdes, et peuvent nécessiter plusieurs jours de calcul même dans le cas de jeux de données comportant seulement quelques milliers de SNPs.

2.2 L’Analyse en Composantes Principales en génétique des populations

L’Analyse en Composantes Principales est un outil incontournable pour l’analyse de données génétiques. Elle est notamment connue pour sa capacité à retrouver la structure génétique, et donne la possibilité de ne garder qu’un nombre réduit de variables tout en résumant l’essentiel de la variation génétique. Nous présentons ici l’Analyse en Composantes Principales ainsi que ses applications en génétique des populations.

2.2.1 Principe de l’Analyse en Composantes Principales

L’Analyse en Composantes Principales est une méthode statistique consistant à transformer un ensemble de variables possiblement corrélées en un nouvel ensemble de variables orthogonales et donc non corrélées appelées composantes principales. De plus, elle est définie de telle sorte que la première composante principale maximise la variance, ce qui signifie que parmi toutes les droites affines de l’espace de départ, la première composante correspond à celle où la projection orthogonale des observations présente la dispersion (ou la variance) la plus grande (Figure 2.2). De la même manière, la deuxième composante principale correspond à la droite affine maximisant la variance, sous la contraite d’être orthogonale à la composante principale précédente (Figure 2.2). Les composantes principales suivantes se déduisent donc des précédentes en suivant ce schéma itératif. La contrainte d’orthogonalité impose que le nombre de composantes principales soit inférieur au nombre de variables et au nombre d’observations. En pratique, le calcul des composantes principales repose sur la diagonalisation de la matrice de covariance ou de la matrice de corrélation.

Analyse en Composantes Principales de données distribuées selon une Gaussienne multivariée. La droite rouge correspond à l’axe de projection maximisant la variance, et donc par définition, à la première composante principale. La droite verte correspond à la deuxième composante principale et se déduit de la première grâce à la contraite d’orthogonalité et au fait qu’il n’y a ici que deux variables.

Figure 2.2: Analyse en Composantes Principales de données distribuées selon une Gaussienne multivariée. La droite rouge correspond à l’axe de projection maximisant la variance, et donc par définition, à la première composante principale. La droite verte correspond à la deuxième composante principale et se déduit de la première grâce à la contraite d’orthogonalité et au fait qu’il n’y a ici que deux variables.

2.2.2 Apparentement génétique interindividuel

L’utilisation de l’Analyse en Composantes Principales en génétique des populations a été popularisée par Cavalli-Sforza (P. Menozzi, Piazza, & Cavalli-Sforza, 1978). En génétique des populations, l’Analyse en Composantes Principales est réalisée à partir de la diagonalisation d’une matrice de covariance particulière, appelée matrice d’apparentement génétique (G. McVean, 2009). Nous avons vu un peu plus haut la notion d’apparentement génétique interpopulationnel ainsi que l’intérêt de corriger la $F_{ST}$ pour celui-ci. Depuis l’apparition des données génomiques, il est possible de définir des mesures de similarité génétique à l’échelle de l’individu à partir des données de génotype, contrairement à l’apparentement génétique interpopulationnel qui est défini à partir des fréquences alléliques. La figure 2.3 compare les deux matrices d’apparentement pour une même simulation d’un modèle démographique à trois populations, et met en évidence le fait que l’apparentement génétique peut s’apprécier à une échelle plus fine. Il existe cependant différentes définitions de l’apparentement génétique interindividuel (D. Speed & Balding, 2015). Nous utiliserons la définition basée sur la corrélation allélique, retenue par Galinsky et al. (2016) et G. Chen, Lee, Zhu, Benyamin, & Robinson (2016), auquel cas l’apparentement génétique entre l’individu $i$ et l’individu $j$ est donnée par la quantité suivante :

\[\begin{equation} G_{RM, ij} = \frac{1}{p} \sum_{k = 1}^p \frac{(G_{ki} - 2p_k) \times (G_{kj} - 2p_k)}{2p_k(1-p_k)} \tag{2.5} \end{equation}\]

À gauche une matrice d’apparentement génétique interpopulationnelle. À droite une matrice d’apparentement génétique interindividuelle. Ces matrices d’apparentement ont été estimées à partir d’une simulation d’un modèle en îles à trois populations.

Figure 2.3: À gauche une matrice d’apparentement génétique interpopulationnelle. À droite une matrice d’apparentement génétique interindividuelle. Ces matrices d’apparentement ont été estimées à partir d’une simulation d’un modèle en îles à trois populations.

2.2.3 Applications en génétique des populations

Visualisation

En génétique des populations, l’ACP est devenu un outil de visualisation extrêmement utilisé. Cela s’explique notamment par sa capacité à rendre compte de la structure de populations à l’aide d’un faible nombre d’axes principaux (Figure 2.4), appelés aussi axes de variation génétique (Price et al., 2006).

ACP réalisée sur la phase 3 du jeu de données HapMap à l’aide de la librarie R pcadapt (Gibbs et al., 2003). Le jeu de données HapMap est un jeu de données humaines incluant une grande diversité de populations. Les deux premières composantes principales distinguent trois groupes génétiques correspondant aux populations africaines, asiatiques et européennes.

Figure 2.4: ACP réalisée sur la phase 3 du jeu de données HapMap à l’aide de la librarie R pcadapt (Gibbs et al., 2003). Le jeu de données HapMap est un jeu de données humaines incluant une grande diversité de populations. Les deux premières composantes principales distinguent trois groupes génétiques correspondant aux populations africaines, asiatiques et européennes.

Correction pour la structure de populations

En génétique associative, où l’on cherche à détecter les gènes associés à un phénotype en comparant des individus porteurs et non-porteurs du phénotype, les composantes principales servent par exemple à corriger pour la structure de populations pour éviter les associations dues à de la différenciation génétique entre les individus porteurs et non-porteurs du phénotype (Price et al., 2006).

Structure géographique

Novembre et al. (2008) ont montré que ces axes de variation génétique pouvaient également être interprétés en terme d’axes géographiques. L’Analyse en Composantes Principales a été réalisée sur un échantillon d’individus européens⁸ issus du jeu de données POPRES (Nelson et al., 2008). En figure 2.5, nous observons en effet que la projection des individus sur les deux premiers axes de l’ACP reflète de façon particulièrement frappante la disposition géographique des différentes populations.

Figure 2.5: ACP réalisée sur le jeu de données POPRES à l’aide de la librarie R pcadapt (Novembre et al., 2008).

Ascendance génétique

Une autre particularité de l’ACP réside dans la possibilité d’inférer les coefficients de métissage ou coefficients d’ascendance à partir des composantes principales (Ma & Amos, 2012; G. McVean, 2009). Un coefficient de métissage quantifie pour un individu donné la proportion de son génôme provenant d’un groupe génétique spécifique (appelé aussi population source ou population ancestrale). L’un des premiers articles à établir un lien entre l’ACP et les coefficients de métissage global fut sur l’interprétation généalogique de l’ACP par G. McVean (2009) (Figure 2.6). Nous reviendrons sur cette notion dans le chapitre consacré à l’introgression.

$Coefficients de métissage et ACP (G. McVean, 2009). A. Chaque ligne de la figure correspond à un chromosome. Chacune de ces lignes représente la projection des individus issus des populations CEU, ASW et YRI sur la première composante principale. B. Chaque point représente un individu correspondant à la moyenne des projections précédentes, ramenées à l’intervalle $[0, 1]$ à l’aide d’une transformation affine.$

Figure 2.6: Coefficients de métissage et ACP (G. McVean, 2009). A. Chaque ligne de la figure correspond à un chromosome. Chacune de ces lignes représente la projection des individus issus des populations CEU, ASW et YRI sur la première composante principale. B. Chaque point représente un individu correspondant à la moyenne des projections précédentes, ramenées à l’intervalle $[0, 1]$ à l’aide d’une transformation affine.

Malgré l’importance et la popularité de l’Analyse Composantes Principales en génétique des populations, ce n’est que récemment que son utilisation a été étendue aux scans génomiques (Duforet-Frebourg et al., 2015; Galinsky et al., 2016). Dans la partie qui suit, nous proposons de nouvelles statistiques basées sur l’Analyse en Composantes Principales et démontrons en quoi elles généralisent les tests classiques de différenciation.

2.3 Statistiques de test basées sur l’Analyse en Composantes Principales

Nous présentons dans cette partie des statistiques basées sur l’Analyse en Composantes Principales dans le cadre des scans génomiques. Les raisons pour lesquelles nous nous sommes intéressés à l’ACP sont multiples. Tout d’abord, l’ACP est particulièrement adaptée au traitement de données en grande dimension, justement parce qu’elle permet de les résumer à l’aide d’un nombre réduit de variables. De plus, comme expliqué précédemment, l’ACP permet de retrouver la structure génétique de façon non paramétrique et sans prior sur l’appartenance individuelle. Grâce à cette propriété, nous montrons la possibilité d’étendre la $F_{ST}$ au cas de populations structurées sans information populationnelle a priori. L’idée de notre démarche repose sur l’hypothèse que les axes principaux reflètent la structure génétique et que les marqueurs génétiques les plus corrélés à ces axes sont des candidats crédibles pour l’adaptation locale. Cette partie sera dédiée à la présentation des méthodes statistiques développées à partir de cette hypothèse de travail. Leur présentation sera accompagnée de validations numériques conduites sur des simulations ainsi que de justifications théoriques quant aux similarités qu’elles présentent avec les méthodes classiques de scan génomique.

Dans l’article 1 (Duforet-Frebourg et al., 2015), nous introduisons la communalité en tant que statistique de test pour la détection de signaux d’adaptation locale. La communalité est une notion empruntée à l’analyse factorielle et s’interprète comme la proportion de variance expliquée par le modèle à facteurs. Nous justifierons également d’un point de vue théorique les observations établissant la correspondance entre l’indice de fixation et la communalité. Pour ce faire, nous montrerons que la $F_{ST}$ peut se réécrire sous la forme d’une statistique de communalité pour un modèle à facteurs discrets, nous invitant de ce fait à considérer la communalité comme une extension de la $F_{ST}$ au cas continu. Cette généralisation est particulièrement intéressante lorsqu’il est difficile de définir des populations de façon claire comme cela peut être le cas en présence d’individus métissés. Cependant, de la même manière que la $F_{ST}$, la communalité n’est pas adaptée au cas de populations structurées.

L’article 2 présente une nouvelle statistique de test basée sur la distance robuste de Mahalanobis pour pallier au problème de la structure de populations (Luu et al., 2017). Enfin nous établirons le lien entre cette nouvelle statistique et le test de Lewontin-Krakauer corrigé pour l’apparentement génétique (2.3), ce qui permettra de conclure quant à la généralisation de la statistique de test $T_{F-LK}$ par cette nouvelle statistique.

Ce travail a notamment abouti sur le développement d’une librairie R implémentant ces statistiques, appelée pcadapt (Duforet-Frebourg et al., 2015; Luu et al., 2017). L’aspect computationnel de ces méthodes sera cependant traité dans le chapitre correspondant, et permettra notamment de discuter des problématiques liées à la présence de données manquantes ainsi que de la complexité algorithmique.

2.3.1 La communalité

Nous présentons dans ce paragraphe un résumé des travaux relatifs à l’article 1 (Duforet-Frebourg et al., 2015). Dans cet article, nous y abordons la possibilité de réaliser des scans génomiques pour la sélection en utilisant l’Analyse en Composantes Principales (ACP). Nous expliquons comment l’indice de différenciation génétique, communément appelé $F_{ST}$, peut être vu comme la proportion de variance expliquée par les composantes principales. La corrélation entre les variants génétiques et les composantes principales donne un cadre conceptuel permettant la détection de variants génétiques impliqués dans les processus d’adaptation locale sans qu’il n’y ait besoin de définir de populations a priori.

Définitions

La première approche présentée repose sur la corrélation des locus avec chaque axe principal :

Définition 2.3 (Corrélation à un axe principal) Soit $G \in \mathcal{M}_{np}(\{0, 1, 2\})$, où $n$ désigne le nombre d’individus et $p$ le nombre de locus. Soit $U \Sigma V^T$ la décomposition en valeurs singulières de rang $K$ de $\tilde{G}$. Notant $\sqrt{\lambda_1}, \dots, \sqrt{\lambda_K}$ les élements diagonaux de $\Sigma$, la corrélation $\rho_{jk}$ du locus $j \in [|1, p|]$ avec le $k$-ième axe principal est donnée par la formule ci-dessous (Cadima & Jolliffe, 1995) :

\[\begin{equation} \rho_{jk} = \frac{\sqrt{\lambda_{k}}V_{jk}}{\sqrt{n-1}}. \end{equation}\]

La seconde approche statistique présentée est une approche multivariée et propose de considérer la somme quadratique de ces corrélations ainsi que la somme quadratique des loadings. L’intérêt de cette approche par rapport à la précédente est de considérer les locus comme des variables multidimensionnelles, généralisant ainsi l’approche composante par composante.

Définition 2.4 (Communalité) Reprenant les notations de la définition 2.3, la communalité $h^2$ est définie pour chaque locus $j$ par la formule ci-dessous :

\[\begin{equation} \begin{split} h_j^2 &= \sum_{k=1}^K \rho_{jk}^2 \\ &= \frac{1}{n-1}\sum_{k=1}^K\sqrt{\lambda_{k}}V_{jk} \\ &\simeq ||\tilde{G}_{.,j}^TU||_2^2. \end{split} \end{equation}\]

Le choix de cette statistique est fondé sur l’interprétation usuelle de la communalité en tant que proportion de variance expliquée par les $K$ premières composantes principales.

Résultats

La $F_{ST}$ peut être interprétée comme une proportion de variance expliquée. Pour chacun des scénarios démographiques étudiés (modèle en île et modèle de divergence à 3 populations), nous avons montré numériquement que la $F_{ST}$ et la communalité sont corrélées à plus de $90\%$ (Duforet-Frebourg et al., 2015). Dans la suite de ce manuscript, nous explicitons un modèle à facteurs où la $F_{ST}$ correspond à la proportion de variance expliquée, ce qui permet de réécrire la $F_{ST}$ sous une forme analytique analogue à celle de la communalité $h^2$.

Proposition 2.1 Soient $N$ le nombre de populations, $n_j$ le nombre d’individus appartenant à la $j$-ème population, $n$ le nombre total d’individus et $\delta_{ji}$ le symbole de Kronecker tel que $\delta_{ji} = 1$ si l’individu $i$ appartient à la $j$-ème population et $0$ sinon. Notant $U_{\delta} = \left(\frac{\delta_{ji}}{\sqrt{2(N-1)}n_j}\right)_{1 \leq i \leq n, \; 1\leq j \leq N} \in \mathcal{M}_{nN}(\mathbb{R})$, il existe une matrice $L \in \mathcal{M}_{Np}(\mathbb{R})$ telle que :

\[\tilde{G} = U_{\delta}L.\]

Preuve. En constatant que $U_{\delta}^TU_{\delta} = \text{Diag}\left(\frac{1}{2(N-1)n_1}, \dots, \frac{1}{2(N-1)n_N}\right)$, nous pouvons définir $L = (U_{\delta}^TU_{\delta})^{-1}U_{\delta}^T\tilde{G}$, si bien que $\tilde{G} = U_{\delta}L$.

En considérant cette factorisation matricielle $\tilde{G} = U_{\delta}L$ comme un modèle à facteurs et en reprenant la définition 2.4, la communalité pour ce modèle s’écrit $||\tilde{G}_{.,j}^T U_{\delta}||_2^2$.

Proposition 2.2 \[F_{ST} = ||\tilde{G}_{.,j}^T U_{\delta}||_2^2.\]

Preuve. En reprenant les notations de la proposition 2.1 et en observant que $p_k = \frac{\sum_{i=1}^n \delta_{ki} G_{ij}}{\sum_{i=1}^n 2\delta_{ki}}$, \[\begin{equation} \begin{split} F_{ST} & = \frac{1}{N-1}\frac{\sum_{k=1}^N (p_k - \bar{p})^2}{\bar{p}(1-\bar{p})} \\ & = \frac{1}{(N-1)\bar{p}(1-\bar{p})}\sum_{k=1}^N \left(\frac{\sum_{i=1}^n \delta_{ki} G_{ij}}{\sum_{i=1}^n 2\delta_{ki}} - \bar{p}\right)^2\\ & = \frac{1}{N-1}\sum_{k=1}^N \left(\frac{\sum_{i=1}^n \delta_{ki} (G_{ij} - 2\bar{p})}{2n_k\sqrt{\bar{p}(1-\bar{p})}}\right)^2\\ \end{split} \end{equation}\]

Or $\tilde{G}_{ij} = \frac{G_{ij} - 2\bar{p}}{\sqrt{2\bar{p}(1-\bar{p})}}$ :

\[\begin{equation} \begin{split} F_{ST} &= \sum_{k=1}^N \left(\sum_{i=1}^n \frac{\delta_{ki}}{\sqrt{2(N-1)}n_k} \tilde{G_{ij}}\right)^2 \\ &= ||\tilde{G}_{.,j}^T U_{\delta}||_2^2. \end{split} \end{equation}\]

La proposition 2.2 permet d’identifier la $F_{ST}$ comme une statistique de communalité dans le cas particulier où les facteurs $U_{\delta}$ sont discrets. La statistique $h^2$ utilise les scores continus de l’ACP en guise de facteurs, ce qui suggère l’utilisation de la communalité en tant que généralisation de la $F_{ST}$. En réalité, ceci ne nous permet pas tout à fait de conclure puisque les matrices $U$ et $U_{\delta}$ peuvent en principe être bien différentes. Dans le cas de deux populations, nous montrons en annexe que la matrice $U$ peut être assimilée à $U_{\delta}$, ce qui permet d’établir le lien entre la $F_{ST}$ et la communalité.

L’approche basée sur l’ACP permet de retrouver des signaux d’adaptation connus. Les statistiques de corrélation et de communalité ont été calculées sur deux jeux de données humaines présentant une structure de populations discrète (1000 Genomes) et une structure de populations continue (POPRES). Les approches basées sur l’ACP permettent de s’affranchir des considérations liées au caractère continu ou discret de la structure de populations, et par la même occasion de la difficulté de définir des populations. L’analyse de ces jeux de données a permis la détection de différentes régions du génôme connues pour être impliquées dans des processus d’adaptation locale (Figure 2.7), achevant de valider l’utilisation de l’ACP en tant qu’outil statistique pour les scans à sélection.

Analyse du jeu données POPRES réalisée avec pcadapt en utilisant la statistique de la communalité. Les points rouges correspondent aux 100 locus ayant les valeurs de communalité les plus élevées.

Figure 2.7: Analyse du jeu données POPRES réalisée avec pcadapt en utilisant la statistique de la communalité. Les points rouges correspondent aux 100 locus ayant les valeurs de communalité les plus élevées.

Limites

De par sa définition, la communalité est dépendante de la proportion de variance expliquée par chaque composante principale retenue et donc de l’amplitude relative des valeurs singulières $\sqrt{\lambda_1}, \dots, \sqrt{\lambda_K}$ de $\tilde{G}$. Dans des cas de figures où $\lambda_1 \gg \lambda_2$, $h^2$ est équivalente à $\rho_1^2$, excluant la possibilité de détecter d’éventuels signaux d’adaptation sur les composantes suivantes. Or ceci s’explique très bien en exploitant le caractère géométrique des statistiques telles que la communalité qui sont distribuées selon un $\chi^2$. Puisqu’il s’agit de sommes quadratiques de lois normales, les courbes de niveau décrites par ces statistiques sont des ellipsoïdes. Les courbes de niveau peuvent être vues comme le pendant géométrique des seuils de $p$-valeur. Si l’on s’intéresse aux courbes de niveau de la communalité, nous nous apercevons en effet qu’elle aura tendance à favoriser la détection de SNPs corrélés avec la première composante, malgré la présence éventuelle de SNPs fortement corrélés avec les composantes suivantes. Pour l’illustrer, plaçons-nous dans le cas $K=2$ (supposant que les deux premières composantes principales aient été retenues). Notant $X$ (resp. $Y$) la variable aléatoire associée aux loadings de la première (resp. seconde) composante principale. La communalité s’écrit $h^2 = \lambda_1 X^2 + \lambda_2 Y^2 = \left(\frac{X}{\sqrt{\lambda_1}}\right)^2 + \left(\frac{Y}{\sqrt{\lambda_2}}\right)^2$ avec $\lambda_1 > \lambda_2$. En dimension $2$, les lignes de niveaux de $h^2$ décrivent des ellipses dont le petit axe est orienté selon $X$ (et paramétré par $\frac{1}{\sqrt{\lambda_1}}$) (Figure 2.8), justifiant ainsi qu’il est plus facile de se trouver en-dehors de l’ellipse en ayant une forte valeur de $X$ qu’une forte valeur de $Y$.

$À gauche les lignes de niveau de la communalité avec $\lambda_1 = 2\lambda_2$ où l’axe des abscisse et l’axe des ordonnées correspondent aux loadings de chaque composante principale. À droite les lignes de niveau de la distance de Mahalanobis.$

Figure 2.8: À gauche les lignes de niveau de la communalité avec $\lambda_1 = 2\lambda_2$ où l’axe des abscisse et l’axe des ordonnées correspondent aux loadings de chaque composante principale. À droite les lignes de niveau de la distance de Mahalanobis.

L’interprétation de la $F_{ST}$ en termes de proportion de variance expliquée fait de la communalité une statistique de test intéressante d’un point de vue conceptuel, permettant de considérer la communalité comme une extension de la $F_{ST}$. Cependant, elle garde de la $F_{ST}$ les défauts liés à la non-prise en compte de l’apparentement génétique interpopulationnel (cf. équation (2.1) (Bonhomme et al., 2010)). Afin de remédier à ce problème, nous proposons une statistique de test basée sur la distance robuste de Mahalanobis.

2.3.2 La distance robuste de Mahalanobis

Ce paragraphe résume de façon détaillée les résultats obtenus dans l’article 2, assortis de justifications qui ne figurent pas nécessairement dans l’article.

Définition

Nous choisissons d’utiliser les coefficients de régression standardisés plutôt que les loadings pour calculer la distance de Mahalanobis. Ce choix sera discuté dans la suite de ce paragraphe.

Définition 2.5 (Coefficients de régression standardisés) Soit $U \Sigma V^T$ la décomposition en valeurs singulières de rang $K$ de $\tilde{G}$. Notant $\epsilon = \tilde{G} - U \Sigma V^T \in \mathcal{M}_{np}(\mathbb{R})$, les coefficients de régression standardisés $z$, appelés encore $z$-scores, sont définis pour chaque locus $j$ par la relation ci-dessous (Saporta, 2006) :

\[\begin{equation} z_j = \frac{U^TG_{.,j}}{\sqrt{\frac{||\epsilon_{.,j}||_2^2}{n - K - 1}}}. \end{equation}\]

La distance de Mahalanobis est une statistique classique de détection de valeurs aberrantes pour des données multivariées, ce qui signifie qu’elle permet de déterminer si une observation $x$ provient effectivement d’un ensemble d’observations $X$. Dans le cas de la communalité par exemple, nous cherchions à détecter les locus significativement isolés de l’ensemble des locus neutres, à partir de la distribution des corrélations. La distance de Mahalanobis apparaît donc comme une alternative naturelle à la communalité pour la détection de locus sous sélection.

Définition 2.6 (Distance de Mahalanobis) Soit $z \in \mathbb{R}^K$, et $Z = (z_1, \dots, z_p) \in \mathcal{M}_{pK}(\mathbb{R})$ une matrice représentant un ensemble de $p$ $z$-scores $K$-dimensionnels. La distance de Mahalanobis de $z$ relativement à cet ensemble est donnée par la formule ci-dessous :

\[\begin{equation} D^2(z) = (z - \bar{z})^T \Sigma^{-1} (z - \bar{z}), \end{equation}\]

où $\bar{z}$ (resp. $\Sigma$) désigne la moyenne des $z$-scores (resp. la matrice de covariance).

Si les $z$-scores $z_i$ sont distribués selon une loi normale multidimensionnelle de moyenne $\bar{z}$ et de matrice de covariance $\Sigma$ définie positive, alors $D^2 \sim \chi_K^2$. En pratique, nous utilisons le facteur d’inflation génomique $\lambda = \text{median}(D^2) / \text{median}(\chi^2_K)$ et approchons $D^2 / \lambda$ par un $\chi^2$ à $K$ degrés de liberté (François, Martins, Caye, & Schoville, 2016).

Contrairement aux loadings, les coefficients de régression standardisés prennent en compte la variance résiduelle. Une comparaison de leur utilisation avec la distance de Mahalanobis est présentée dans la suite de ce paragraphe.

Estimation robuste de la matrice de covariance

La statistique de test $T_{F-LK}$ de Bonhomme et al. (2010) et la $F_{ST}$ généralisée proposée par Ochoa & Storey (2016) sont basées sur une estimation de la matrice d’apparentement génétique $\mathcal{F}$ par une méthode des moments. Pour des données gaussiennes multivariées, cela consiste à estimer la moyenne et la matrice de covariance. Les estimateurs classiques de moyenne sont connus pour être sensibles à la présence de données aberrantes (Figure 2.9). Un estimateur est dit robuste s’il n’est pas ou s’il est peu affecté par la présence de données aberrantes. La distance de Mahalanobis peut être rendue robuste si l’estimation de la moyenne $\bar{x}$ et de la matrice de covariance $\Sigma$ se fait à l’aide d’estimateurs robustes. Nous montrons ici la nécessité d’utiliser un estimateur robuste pour la covariance et justifions notre choix d’estimateur. Nous proposons une comparaison géométrique (Figure 2.9) de deux estimateurs robustes de la matrice de covariance : MCD (Covariance de Déterminant Minimal (Rousseeuw, 1985)) et OGK (Gnanadesikan-Kettenring Orthogonalisé (Maronna & Zamar, 2002)). La procédure de comparaison est la suivante :

les matrices de covariance sont estimées pour les méthodes OGK et MCD sur deux jeux de données simulées ($96,7\%$ de locus neutres pour le modèle de divergence contre $92,8\%$ de locus neutres pour le modèle en îles). L’estimateur de covariance classique est également inclus dans la comparaison.
les matrices de covariance ainsi estimées sont ensuites représentées par des ellipses relativement à un niveau de confiance fixé ici à $95\%$, signifiant que les ellipses de covariance (Figure 2.9) sont censées recouvrir $95\%$ des observations neutres.
nous regardons ensuite quelle ellipse contient le moins d’observations aberrantes, tout en couvrant au moins $95\%$ des locus neutres.

La figure 2.9 montre que l’utilisation d’une méthode d’estimation robuste est nécessaire, même pour des jeux de données présentant moins de $10\%$ de données aberrantes. Quant à la comparaison entre l’estimateur OGK et l’estimateur MCD, les deux méthodes semblent effectivement tenir compte de la présence de données aberrantes pour ajuster le calcul de la matrice de covariance. Nous notons néanmoins que l’estimateur OGK retient moins de données aberrantes tout en recouvrant une proportion de locus neutres visiblement supérieure à $95\%$ au sein de son ellipse de confiance à $95\%$. Notre choix d’estimateur s’est donc porté sur l’estimateur OGK et a été réimplémenté dans la librairie pcadapt.

$Comparaison des estimations de la matrice de covariance. Les matrices de covariance peuvent être interprétées géométriquement en termes d’ellipses. Une ellipse de confiance au seuil $\alpha$ est censée contenir une proportion $\alpha$ de l’ensemble des observations effectivement issues de la distribution (correspondant aux observations neutres). Nous constatons que l’estimateur classique de covariance surestime les valeurs propres de la matrice de covariance même lorsque la proportion de données aberrantes est faible.$

Figure 2.9: Comparaison des estimations de la matrice de covariance. Les matrices de covariance peuvent être interprétées géométriquement en termes d’ellipses. Une ellipse de confiance au seuil $\alpha$ est censée contenir une proportion $\alpha$ de l’ensemble des observations effectivement issues de la distribution (correspondant aux observations neutres). Nous constatons que l’estimateur classique de covariance surestime les valeurs propres de la matrice de covariance même lorsque la proportion de données aberrantes est faible.

Loadings et coefficients de régression standardisés

Nous justifions ici de façon empirique l’utilisation des coefficients de régression standardisés ($z$-scores) au détriment des loadings pour calculer la distance de Mahalanobis. La procédure de comparaison est la suivante :

la distance de Mahalanobis et les $p$-valeurs associées sont calculées à partir des loadings et des coefficients de régression standardisés.
pour chacun des seuils de significativité (ici $5\%$, $10\%$ et $20\%$), nous déterminons l’ensemble des SNPs présentant une $p$-valeur ajustée⁹ inférieure à ce seuil et calculons le taux de fausses de découvertes et la puissance relativement à cet ensemble.

Cette procédure de comparaison nous permet par exemple d’évaluer si le taux de fausses découvertes est bien contrôlé ou encore si une méthode est trop conservative ou non. En figure 2.10, nous constatons que les deux statistiques sont bien contrôlées, c’est-à-dire que pour chacun des seuils de significativité, aucune des deux statistiques ne présentent un taux de fausses découvertes supérieur à ce seuil. La figure 2.10 met en évidence le côté relativement conservatif de la distance de Mahalanobis calculée à partir des loadings, ce qui a pour effet de diminuer la puissance du test.

$Comparaison des distances de Mahalanobis calculées à partir des loadings et des $z$-scores. Les taux de fausses de découvertes et les puissances sont calculées puis moyennées sur l’ensemble des simulations de modèles en îles utilisées dans l’article 2. La ligne en pointillé représente le taux de fausses découvertes attendu pour chaque seuil de significativité.$

Figure 2.10: Comparaison des distances de Mahalanobis calculées à partir des loadings et des $z$-scores. Les taux de fausses de découvertes et les puissances sont calculées puis moyennées sur l’ensemble des simulations de modèles en îles utilisées dans l’article 2. La ligne en pointillé représente le taux de fausses découvertes attendu pour chaque seuil de significativité.

Rappel des résultats principaux de l’article 2

Nous avons développé une méthode de scan à sélection valable pour différentes structures de populations, aussi bien adaptée au cas de populations discrètes qu’au cas de populations continues. À l’aide de simulations reproduisant différents scénarios démographiques, exhibant des structures de populations discrètes et continues, nous montrons que l’utilisation de la distance robuste de Mahalanobis permet de pallier aux défauts de la communalité. En termes de sensibilité statistique et de contrôle du taux de fausses découvertes, notre méthode affiche de meilleurs résultats en comparaison à d’autres méthodes de scan à sélection (OutFLANK, Bayescan, FLK), quelque soit le modèle démographique sous-jacent, à l’exception du modèle en îles où les performances sont équivalentes. Nous étudions également l’impact du caractère discret ou continu de la structure de populations sur les méthodes de scan à sélection et montrons que notre méthode est moins sensible à la présence d’individus métissés, contrairement aux méthodes basées sur la $F_{ST}$.

Une généralisation de $T_{F-LK}$. La similarité des résultats numériques produits par pcadapt et FLK suggèrent l’existence d’un lien entre les deux méthodes et nous justifions ce résultat en annexe. La distance robuste de Mahalanobis calculée à partir de l’ACP généralise la statistique de test $T_{F-LK}$ car elle peut être utilisée dans le cas de populations continues.

2.4 Article 1

Abstract

To characterize natural selection, various analytical methods for detecting candidate genomic regions have been developed. We propose to perform genome-wide scans of natural selection using principal component analysis (PCA). We show that the common $F_{ST}$ index of genetic differentiation between populations can be viewed as the proportion of variance explained by the principal components. Considering the correlations between genetic variants and each principal component provides a conceptual framework to detect genetic variants involved in local adaptation without any prior definition of populations. To validate the PCA-based approach, we consider the 1000 Genomes data (phase 1) considering 850 individuals coming from Africa, Asia, and Europe. The number of genetic variants is of the order of 36 millions obtained with a low-coverage sequencing depth (3x). The correlations between genetic variation and each principal component provide well-known targets for positive selection (EDAR, SLC24A5, SLC45A2, DARC), and also new candidate genes (APPBPP2, TP1A1, RTTN, KCNMA, MYO5C) and noncoding RNAs. In addition to identifying genes involved in biological adaptation, we identify two biological pathways involved in polygenic adaptation that are related to the innate immune system (beta defensins) and to lipid metabolism (fatty acid omega oxidation). An additional analysis of European data shows that a genome scan based on PCA retrieves classical examples of local adaptation even when there are no well-defined populations. PCA-based statistics, implemented in the PCAdapt R package and the PCAdapt fast open-source software, retrieve well-known signals of human adaptation, which is encouraging for future whole-genome sequencing project, especially when defining populations is difficult.

Significance Statement

Positive natural selection or local adaptation is the driving force behind the adaption of individuals to their environment. To identify genomic regions responsible for local adaptation, we propose to consider the genetic markers that are the most related with population structure. To uncover genetic structure, we consider principal component analysis that identifies the primary axes of variation in the data. Our approach generalizes common approaches for genome scan based on measures of population differentiation. To validate our approach, we consider the human 1000 Genomes data and find well-known targets for positive selection as well as new candidate regions. We also find evidence of polygenic adaptation for two biological pathways related to the innate immune system and to lipid metabolism.

Introduction

Because of the flood of genomic data, the ability to understand the genetic architecture of natural selection has dramatically increased. Of particular interest is the study of local positive selection which explains why individuals are adapted to their local environment. In humans, the availability of genomic data fostered the identification of loci involved in positive selection (Barreiro, Laval, Quach, Patin, & Quintana-Murci, 2008; Grossman et al., 2013; Pickrell et al., 2009; Sabeti et al., 2007). Local positive selection tends to increase genetic differentiation, which can be measured by difference of allele frequencies between populations (Colonna et al., 2014; Nielsen, 2005; Sabeti et al., 2006). For instance, a mutation in the DARC gene that confers resistance to malaria is fixed in sub-Saharan African populations whereas it is absent elsewhere (Hamblin, Thompson, & Di Rienzo, 2002). In addition to the variants that confer resistance to pathogens, genome scans also identify other genetic variants, and many of these are involved in human metabolic phenotypes and morphological traits (Barreiro et al., 2008; Hancock et al., 2010).

In order to provide a list of variants potentially involved in natural selection, genome scans compute measures of genetic differentiation between populations and consider that extreme values correspond to candidate regions (G. Luikart, England, Tallmon, Jordan, & Taberlet, 2003). The most widely used index of genetic differentiation is the $F_{ST}$ index which measures the amount of genetic variation that is explained by variation between populations (Excoffier, Smouse, & Quattro, 1992). However the $F_{ST}$ statistic requires to group individuals into populations which can be problematic when ascertainment of population structure does not show well-separated clusters of individuals (e.g., Novembre et al. (2008)). Other statistics related to $F_{ST}$ have been derived to reduce the false discovery rate (FDR) obtained with $F_{ST}$ but they also work at the scale of populations (Bonhomme et al., 2010; Fariello, Boitard, Naya, SanCristobal, & Servin, 2013; Günther & Coop, 2013). Grouping individuals into populations can be subjective, and important signals of selection may be missed with an inadequate choice of populations (W.-Y. Yang, Novembre, Eskin, & Halperin, 2012). We have previously developed an individual-based approach for selection scan based on a Bayesian factor model but the Markov chain Monte Carlo (MCMC) algorithm required for model fitting does not scale well to large data sets containing a million of variants or more (Duforet-Frebourg et al., 2014).

We propose to detect candidates for natural selection using principal component analysis (PCA). PCA is a technique of multivariate analysis used to ascertain population structure (N. Patterson, Price, & Reich, 2006). PCA decomposes the total genetic variation into $K$ axes of genetic variation called principal components. In population genomics, the principal components can correspond to evolutionary processes such as evolutionary divergence between populations (G. McVean, 2009). Using simulations of an island model and of a model of population fission followed by isolation, we show that the common $F_{ST}$ statistic corresponds to the proportion of variation explained by the first K principal components when K has been properly chosen. With this point of view, the $F_{ST}$ of a given variant is obtained by summing the squared correlations of the first $K$ principal components opening the door to new statistics for genome scans. At a genome-wide level, it is known that there is a relationship between $F_{ST}$ and PCA (G. McVean, 2009), and our simulations show that the relationship also applies at the level of a single variant.

The advantages of performing a genome scan based on PCA are multiple: it does not require to group individuals into populations, the computational burden is considerably reduced compared with genome scan approaches based on MCMC algorithms (Foll & Gaggiotti, 2008); Riebler, Held, & Stephan (2008); Günther & Coop (2013); Duforet-Frebourg et al. (2014)], and candidate single nucleotide polymorphisms (SNPs) can be related to different evolutionary events that correspond to the different principal components. Using simulations and the 1000 Genomes data, we show that PCA can provide useful insights for genome scans. Looking at the correlations between SNPs and principal components provides a novel conceptual framework to detect genomic regions that are candidates for local adaptation.

New Method

New Statistics for Genome Scan

We denote by $Y$ the $(n \times p)$ centered and scaled genotype matrix where $n$ is the number of individuals and $p$ is the number of loci. The new statistics for genome scan are based on PCA. The objective of PCA is to find a new set of orthogonal variables called the principal components, which are linear combinations of (centered and standardized) allele counts, such that the projections of the data onto these axes lead to an optimal summary of the data. To present the method, we introduce the truncated singular value decomposition (SVD) that approximates the data matrix $Y$ by a matrix of smaller rank

\[\begin{equation} Y \approx U \Sigma V^T, \tag{2.6} \end{equation}\]

where $U$ is a $(n \times K)$ orthonormal matrix, $V$ is a $(p \times K)$ orthonormal matrix, $\Sigma$ is a diagonal $(K \times K)$ matrix and $K$ corresponds to the rank of the approximation. The solution of PCA with $K$ components can be obtained using the truncated SVD: the $K$ columns of $V$ contain the coefficients of the new orthogonal variables, the $K$ columns of $U$ contain the projections (called “scores”) of the original variables onto the principal components and capture population structure (supplementary fig. B.1), and the squares of the elements of $\Sigma$ are proportional to the proportion of variance explained by each principal component (Jolliffe, 1986). We denote the diagonal elements of $\Sigma$ by $\sqrt{\lambda_k}, \; k = 1, \dots, K$ where the $\lambda_k$’s are the ranked eigenvalues of the matrix $YY^T$. Denoting by $V_{jk}$, the entry of $V$ at the $j^{th}$ line and $k^{th}$ column, then the correlation $\rho_{jk}$ between the $j^{th}$ SNP and the $k^{th}$ principal component is given by $\rho_{jk} = \sqrt{\lambda_{k}}V_{jk}/\sqrt{n-1}$ (Cadima & Jolliffe, 1995). In the following, the statistics $\rho_{jk}$ are referred to as “loadings” and will be used for detecting selection. The second statistic we consider for genome scan corresponds to the proportion of variance of a SNP that is explained by the first $K$ PCs. It is called the communality in exploratory factor analysis because it is the variance of observed variables accounted for by the common factors, which correspond to the first $K$ PCs. Because the principal components are orthogonal to each other, the proportion of variance explained by the first $K$ principal components is equal to the sum of the squared correlations with the first $K$ principal components. Denoting by $h_j^2$ the communality of the $j^{th}$ SNP, we have

\[\begin{equation} h_j^2 = \sum_{k=1}^K \rho_{jk}^2. \tag{2.7} \end{equation}\]

The last statistic we consider for genome scans sums the squared of normalized loadings. It is defined as ${h^{\prime}_j}^2 = \sum_{k=1}^K V_{jk}^2$. Compared to the communality $h^2$, the statistic ${h^{\prime}_j}^2$ should theoretically give the same importance to each PC because the normalized loadings are on the same scale as we have $\sum_{j=1}^K V_{jk}^2 = 1$, for $k = 1, \dots, K$.

Numerical Computations

The method of selection scan should be able to handle a large number $p$ of genetic variants. In order to compute truncated SVD with large values of $p$, we compute the $n \times n$ covariance matrix $\Omega = YY^T/(p-1)$. The covariance matrix $\Omega$ is typically of much smaller dimension than the $p \times p$ covariance matrix. Considering the $n \times n$ covariance matrix $\Omega$ speeds up matrix operations. Computation of the covariance matrix is the most costly operation and it requires a number of arithmetic operations proportional to $p n^2$. After computing the covariance matrix $\Omega$, we compute its first $K$ eigenvalues and eigenvectors to find $\Sigma^2/(p-1)$ and $U$. Eigenanalysis is performed with the dsyevr routine of the linear algebra package LAPACK (Anderson et al., 1999). The matrix $V$, which captures the relationship between each SNPs and population structure, is obtained by the matrix operation $V^T = \Sigma^{-1} U^T Y$. The software PCAdapt fast, process data as a stream and never store in order to have a very low memory access whatever the size of the data.

Results

Island Model

To investigate the relationship between communality $h^2$ and $F_{ST}$, we consider an island model with three islands. We use $K = 2$ when performing PCA because there are three islands. We choose a value of the migration rate that generates a mean $F_{ST}$ value (across the 1,400 neutral SNPs) of 4%. We consider five different simulations with varying strengths of selection for the 100 adaptive SNPs. In all simulations, the $R^2$ correlation coefficient between $h^2$ and $F_{ST}$ is larger than 98%. Considering as candidate SNPs the 1% of the SNPs with largest values of $F_{ST}$ or of $h^2$, we find that the overlap coefficient between the two sets of SNPs is comprised between 88% and 99%. When varying the strength of selection for adaptive SNPs, we find that the relative difference of FDRs obtained with $F_{ST}$ (top 1%) and with $h^2$ (top 1%) is smaller than 5%. The similar values of FDR obtained with $h^2$ and with $F_{ST}$ decrease for increasing strength of selection (supplementary fig. B.2).

Divergence Model

To compare the performance of different PCA-based summary statistics, we simulate genetic variation in models of population divergence. The divergence models assume that there are three populations, $A$, $B_1$ and $B_2$ with $B_1$ and $B_2$ being the most related populations (figs. 2.11 and 2.12). The first simulation scheme assumes that local adaptation took place in the lineages corresponding to the environments of populations $A$ and $B_1$ (fig. 2.11). The SNPs, which are assumed to be independent, are divided into three groups: 9,500 SNPs evolve neutrally, 250 SNPs confer a selective advantage in the environment of $A$, and 250 other SNPs confer a selective advantage in the environment of $B_1$. Genetic differentiation, measured by pairwise $F_{ST}$, is equal to 14% when comparing population $A$ to the other ones and is equal to 5% when comparing populations $B_1$ and $B_2$. Performing PCA with $K = 2$ shows that the first component separates population $A$ from $B_1$ and $B_2$ whereas the second component separates $B_1$ from $B_2$ (supplementary fig. B.1). The choice of $K = 2$ is evident when looking at the scree plot because the eigenvalues, which are proportional to the proportion of variance explained by each PC, drop beyond $K = 2$ and stay almost constant as $K$ further increases (supplementary fig. B.3).

$Repartition of the 1% top-ranked SNPs for each PCA-based statistic under a divergence model with two types of adaptive constraints. Thicker and colored lineages correspond to lineages where adaptation took place. The squared loadings with PC1 $\rho_{j1}^2$ pick a large proportion of SNPs involved in selection in population $A$ whereas the squared loadings with PC2 $\rho_{j2}^2$ pick SNPs involved in selection in population $B_1$. This difference is reflected in the different repartition of the top-ranked SNPs for the communality $h^2$ and the statistic ${h^{\prime}}^2$.$

Figure 2.11: Repartition of the 1% top-ranked SNPs for each PCA-based statistic under a divergence model with two types of adaptive constraints. Thicker and colored lineages correspond to lineages where adaptation took place. The squared loadings with PC1 $\rho_{j1}^2$ pick a large proportion of SNPs involved in selection in population $A$ whereas the squared loadings with PC2 $\rho_{j2}^2$ pick SNPs involved in selection in population $B_1$. This difference is reflected in the different repartition of the top-ranked SNPs for the communality $h^2$ and the statistic ${h^{\prime}}^2$.

$Repartition of the 1% top-ranked SNPs of each PCA-based statistic under a divergence model with four types of adaptive constraints. Thicker and colored lineages correspond to lineages where adaptation occurred. The different types of SNPs picked by the squared loadings $\rho_{j1}^2$ and $\rho_{j2}^2$ are also found when comparing the communality $h^2$ and the statistic ${h^{\prime}}^2$.$

Figure 2.12: Repartition of the 1% top-ranked SNPs of each PCA-based statistic under a divergence model with four types of adaptive constraints. Thicker and colored lineages correspond to lineages where adaptation occurred. The different types of SNPs picked by the squared loadings $\rho_{j1}^2$ and $\rho_{j2}^2$ are also found when comparing the communality $h^2$ and the statistic ${h^{\prime}}^2$.

We investigate the relationship between the communality statistic $h^2$, which measures the proportion of variance explained by the first two PCs, and the $F_{ST}$ statistic. We find a squared Pearson correlation coefficient between the two statistics larger than 98.8% in the simulations corresponding to figures 2.11 and 2.12 (supplementary fig. B.4). For these two simulations, we look at the SNPs in the top 1% (respectively, 5%) of the ranked lists based on $h^2$ and $F_{ST}$, and we find an overlap coefficient always larger than 93% for the lists provided by the two different statistics (respectively, 95%). Providing a ranking of the SNPs almost similar to the ranking provided by $F_{ST}$ is therefore possible without considering that individuals originate from predefined populations.

We then compare the performance of the different statistics based on PCA by investigating if the top-ranked SNPs (top 1%) manage to pick SNPs involved in local adaptation (fig. 2.11). The squared loadings $\rho_{j1}^2$ with the first PC pick SNPs involved in selection in population $A$ (39% of the top 1%), a few SNPs involved in selection in $B_1$ (9%), and many false positive SNPs (FDR of 53%). The squared loadings with the second PC $\rho_{j2}^2$ pick less false positives (FDR of 12%) and most SNPs are involved in selection in $B_1$ (88%) with just a few involved in selection in $A$ (1%). When adaptation took place in two different evolutionary lineages of a divergence tree between populations, a genome scan based on PCA has the nice property that outlier loci correlated with PC1 or with PC2 correspond to adaptive constraints that occurred in different parts of the tree.

Because the communality $h^2$ gives more importance to the first PC, it picks preferentially the SNPs that are the most correlated with PC1. There is a large overlap of 72% between the 1% top-ranked lists provided by $h^2$ and $\rho_{j1}^2$. Therefore, the communality statistic $h^2$ is more sensitive to ancient adaptation events that occurred in the environment of population $A$. In contrast, the alternative statistic ${h^{\prime}}^2$ is more sensitive to recent adaptation events that occurred in the environment of population $B_1$. When considering the top-ranked 1% of the SNPs, ${h^{\prime}}^2$ captures only one SNP involved in selection in $A$ (1% of the top 1%) and 88 SNPs related to adaptation in $B_1$ (88% of the top 1%). The overlap between the 1% top-ranked lists provided by ${h^{\prime}}^2$ and by $\rho_{j2}^2$ is of 86%.

The ${h^{\prime}}^2$ statistic is mostly influenced by the second principal component because the distribution of squared loadings corresponding to the second PC has a heavier tail, and this result holds for the two divergence models and for the 1000 Genomes data (supplementary fig. B.5). To summarize, the $h^2$ and ${h^{\prime}}^2$ statistics give too much importance to PC1 and PC2, respectively, and they fail to capture in an equal manner both types of adaptive events occurring in the environment of populations $A$ and $B_1$.

We also investigate a more complex simulation in which adaptation occurs in the four branches of the divergence tree (fig. 2.12). Among the 10,000 simulated SNPs, we assume that there are four sets of 125 adaptive SNPs with each set being related to adaptation in one of the four branches of the divergence tree. Compared with the simulation of figure 2.11, we find the same pattern of population structure (supplementary fig. B.1). The squared loadings $\rho_{j1}^2$ with the first PC mostly pick SNPs involved in selection in the branch that predates the split between $B_1$ and $B_2$ (51% of the top 1%), SNPs involved in selection in the environment of population $A$ (9%), and false positive SNPs (FDR of 38%). Except for false positives (FDR of 14%), the squared loadings $\rho_{j2}^2$ with the second PC rather pick SNPs involved in selection in $B_1$ and $B_2$ (42% for $B_1$ and 44% for $B_2$). Once again, there is a large overlap between the SNPs picked by the communality $h^2$ and by $\rho_1^2$ (92% of overlap) and between the SNPs picked by ${h^{\prime}}^2$ and $\rho_2^2$ (93% of overlap). Because the first PC discriminates population $A$ from $B_1$ and $B_2$ (supplementary fig. B.1), the SNPs most correlated with PC1 correspond to SNPs related to adaptation in the (red and green) branches that separate $A$ from populations $B_1$ and $B_2$. In contrast, the SNPs that are most correlated to PC2 correspond to SNPs related to adaptation in the two (blue and yellow) branches that separate population $B_1$ from $B_2$ (fig. 2.12).

We additionally evaluate to what extent the results are robust with respect to some parameter settings. When considering the 5% of the SNPs with most extreme values of the statistics instead of the top 1%, we also find that the summary statistics pick SNPs related to different evolutionary events (supplementary fig. B.6). The main difference being that the FDR increases considerably when considering the top 5% instead of the top 1% (supplementary fig. B.6). We also consider variation of the selection coefficient ranging from $s = 1.01$ to $s = 1.1$ ($s = 1.025$ corresponds to the simulations of figs. 2.11 and 2.12). As expected, the FDR of the different statistics based on PCA is considerably reduced when the selection coefficient increases (supplementary fig. B.7).

In the divergence model of figure 2.11, we also compare the FDRs obtained with the statistics $h^2$, ${h^{\prime}}^2$, and with a Bayesian factor model implemented in the software PCAdapt (Duforet-Frebourg et al., 2014). For the optimal choice of $K = 2$, the statistic ${h^{\prime}}^2$ and the Bayesian factor model provide the smallest FDR (supplementary fig. B.8). However, when varying the value of $K$ from $K = 1$ to $K = 6$, we find that the communality $h^2$ and the Bayesian approach are robust to overspecification of $K$ ($K > 3$) whereas the FDR obtained with ${h^{\prime}}^2$ increases importantly as $K$ increases beyond $K = 2$ (supplementary fig. B.8).

We also consider a more general isolation-with-migration model. In the divergence model where adaptation occurs in two different lineages of the population tree (fig. 2.11), we add constant migration between all pairs of populations. We assume that migration occurred after the split between $B_1$ and $B_2$. We consider different values of migration rates generating a mean $F_{ST}$ of 7.5% for the smallest migration rate to a mean $F_{ST}$ of 0% for the largest migration rate. We find that the $R^2$ correlation between $F_{ST}$ and $h^2$ decreases as a function of the migration rate (supplementary fig. B.9). For $F_{ST}$ values larger than 0.5%, $R^2$ is larger than 97%. The squared correlation $R^2$ decreases to 47% for the largest migration rate. Beyond a certain level of migration rate, population structure, as ascertained by principal components, is no more described by well-separated clusters of individuals (supplementary fig. B.10) but by a more clinal or continuous pattern (supplementary fig. B.10) explaining the difference between $F_{ST}$ and $h^2$. However, the FDRs obtained with the different statistics based on PCA and with $F_{ST}$ evolve similarly as a function of the migration rate. For both types of approaches, the FDR increases for larger migration with almost no true discovery (only one true discovery in the top 1% lists) when considering the largest migration rate.

The main results obtained under the divergence models can be described as follows. The principal components correspond to different evolutionary lineages of the divergence tree. The communality statistic $h^2$ provides similar list of candidate SNPs than $F_{ST}$ and it is mostly influenced by the first principal component which can be problematic if other PCs also convey adaptive events. To counteract this limitation, which can potentially lead to the loss of important signals of selection, we show that looking at the squared loadings with each of the principal components provide adaptive SNPs that are related to different evolutionary events. When adding migration rates between lineages, we find that the main results are unchanged up to a certain level of migration rate. Above this level of migration rate, the relationship between $F_{ST}$ and $h^2$ does not hold anymore and genome scans based on either PCA or $F_{ST}$ produce a majority of false positives.

1000 Genomes Data

Since we are interested in selective pressures that occurred during the human diaspora out of Africa, we decide to exclude individuals whose genetic makeup is the result of recent admixture events (African Americans, Columbians, Puerto Ricans, and Mexicans). The first three principal components capture population structure whereas the following components separate individuals within populations (fig. 2.13 and supplementary fig. B.11). The first and second PCs ascertain population structure between Africa, Asia, and Europe (fig. 2.13) and the third principal component separates the Yoruba from the Luhya population (supplementary fig. B.11). The decay of eigenvalues suggests to use $K = 2$ because the eigenvalues drop between $K = 2$ and $K = 3$ where a plateau of eigenvalues is reached (supplementary fig. B.3).

$PCA with $K = 2$ applied to the 1000 Genomes data. The sampled populations are the following: British in England and Scotland (GBR), Utah residents with Northern and Western European ancestry (CEU), Finnish in Finland (FIN), Iberian populations in Spain (IBS), Toscani in Italy (TSI), Han Chinese in Bejing (CHB), Southern Han Chinese (CHS), Japanese in Tokyo (JPT), Luhya in Kenya (LWK), Yoruba in Nigeria (YRI).$

Figure 2.13: PCA with $K = 2$ applied to the 1000 Genomes data. The sampled populations are the following: British in England and Scotland (GBR), Utah residents with Northern and Western European ancestry (CEU), Finnish in Finland (FIN), Iberian populations in Spain (IBS), Toscani in Italy (TSI), Han Chinese in Bejing (CHB), Southern Han Chinese (CHS), Japanese in Tokyo (JPT), Luhya in Kenya (LWK), Yoruba in Nigeria (YRI).

When performing a genome scan with PCA, there are different choices of statistics. The first choice is the $h^2$ communality statistic. Using the three continents as labels, there is a squared correlation between $h^2$ and $F_{ST}$ of $R^2 = 0.989$. To investigate if $h^2$ is mostly influenced by the first PC, we determine if the outliers for the $h^2$ statistics are related with PC1 or with PC2. Among the top 0.1% of SNPs with the largest values of $h^2$, we find that 74% are in the top 0.1% of the squared loadings $\rho_{j1}^2$ corresponding to PC1 and 20% are in the top 0.1% of the squared loadings $\rho_{j2}^2$ corresponding to PC2. The second possible choice of summary statistics is the ${h^{\prime}}^2$ statistic. Investigating the repartition of the 0.1% outliers for ${h^{\prime}}$, we find that 0.005% are in the top 0.1% of the squared loadings $\rho_{j1}^2$ corresponding to PC1 and 85% are in the top 0.1% of the squared loadings $\rho_{j2}^2$ corresponding to PC2. The ${h^{\prime}}^2$ statistic is mostly influenced by the second PC because the distribution of the $V_{j2}^2$ (normalized squared loadings) has a longer tail than the corresponding distribution for PC1 (supplementary fig. B.5). Because the $h^2$ statistic is mostly influenced by PC1 and ${h^{\prime}}^2$ is mostly influenced by PC2, confirming the results obtained under the divergence models, we rather decide to perform two separate genome scans based on the squared loadings $\rho_{j1}^2$ and $\rho_{j2}^2$.

The two Manhattan plots based on the squared loadings for PC1 and PC2 are displayed in figures 2.14 and 2.15. Because of linkage disequilibrium (LD), Manhattan plots generally produce clustered outliers. To investigate if the top 0.1% outliers are clustered in the genome, we count–for various window sizes–the proportion of contiguous windows containing at least one outlier. We find that outlier SNPs correlated with PC1 or with PC2 are more clustered than expected if they would have been uniformly distributed among the 36,536,154 variants (supplementary fig. B.12). Additionally, the clustering is larger for the outliers related to the second PC as they cluster in fewer windows (supplementary fig. B.12). As the genome scan for PC2 captures more recent adaptive events, it reveals larger genomic windows that experienced fewer recombination events.

$Manhattan plot for the 1000 Genomes data of the squared loadings $\rho_{j1}^2$ with the first principal component. For sake of presentation, only the top-ranked SNPs (top 0.1%) are displayed and the 100 top-ranked SNPs are colored in red.$

Figure 2.14: Manhattan plot for the 1000 Genomes data of the squared loadings $\rho_{j1}^2$ with the first principal component. For sake of presentation, only the top-ranked SNPs (top 0.1%) are displayed and the 100 top-ranked SNPs are colored in red.

$Manhattan plot for the 1000 Genomes data of the squared loadings $\rho_{j2}^2$ with the second principal component. For sake of presentation, only the top-ranked SNPs (top 0.1%) are displayed and the 100 top-ranked SNPs are colored in red.$

Figure 2.15: Manhattan plot for the 1000 Genomes data of the squared loadings $\rho_{j2}^2$ with the second principal component. For sake of presentation, only the top-ranked SNPs (top 0.1%) are displayed and the 100 top-ranked SNPs are colored in red.

The 1000 Genome data contain many low-frequency SNPs; 82% of the SNPs have a minor allele frequency smaller than 5%. However, these low-frequency variants are not found among outlier SNPs. There are no SNP with a minor allele frequency smaller than 5% among the 0.1% of the SNPs most correlated with PC1 or with PC2.

The 100 SNPs that are the most correlated with the first PC are located in 24 genomic regions. Most of the regions contain just one or a few SNPs except a peak in the gene APPBP2 that contains 33 out of the 100 top SNPs, a peak encompassing the RTTN and CD226 genes containing 17 SNPS and a peak in the ATP1A1 gene containing seven SNPs (fig. 2.14). Confirming a larger clustering for PC2 outliers, the 100 SNPs that are the most correlated with PC2 cluster in fewer genomic regions. They are located in 14 genomic regions including a region overlapping with EDAR contains 44 top hits, two regions containing eight SNPs and located in the pigmentation genes SLC24A5 and SLC45A2, and two regions with seven top hit SNPs, one in the gene KCNMA1 and another one encompassing the RGLA/MYO5C genes (fig. 2.15).

We perform Gene Ontology (GO) enrichment analyses using Gowinda for the SNPs that are the most correlated with PC1 and PC2. For PC1, we find, among others, enrichment ($\text{FDR} \leq 5 \%$) for ontologies related to the regulation of arterial blood pressure, the endocrine system and the immunity response (interleukin production, response to viruses). For PC2, we find enrichment ($\text{FDR} \leq 5 \%$) related to olfactory receptors, keratinocyte and epidermal cell differentiation, and ethanol metabolism. We also search for polygenic adaptation by looking for biological pathways enriched with outlier genes (Daub et al., 2013). For PC1, we find one enriched ($\text{FDR} \leq 5 \%$) pathway consisting of the beta defensin pathway. The beta defensin pathway contains mainly genes involved in the innate immune system consisting of 36 defensin genes and of two Toll-Like receptors (TLR1 and TLR2). There are additionally two chemokine receptors (CCR2 and CCR6) involved in the beta defensin pathway. For PC2, we also find one enriched pathway consisting of fatty acid omega oxidation ($\text{FDR} \leq 5 \%$). This pathway consists of genes involved in alcohol oxidation (CYP, ALD, and ALDH genes). Performing a less stringent enrichment analysis which can find pathways containing overlapping genes, we find more enriched pathways: the beta defensin and the defensin pathways for PC1 and ethanol oxidation, glycolysis/gluconeogenesis and fatty acid omega oxidation for PC2.

To further validate the proposed list of candidate SNPs involved in local adaptation, we test for an enrichment of genic or nonsynonymous SNP among the SNPs that are the most correlated with the PC. We measure the enrichment among outliers by computing odds ratio (Fagny et al., 2014; Kudaravalli, Veyrieras, Stranger, Dermitzakis, & Pritchard, 2008). For PC1, we do not find significant enrichments (table 2.1) except when measuring the enrichment of genic regions compared with nongenic regions (OR = 10.18 for the 100 most correlated SNPs, $P \leq 5 \%$ using a permutation procedure). For PC2, we find an enrichment of genic regions among outliers as well as an enrichment of nonsynonymous SNPs (table 2.1). By contrast with the enrichment of genic regions for SNPs extremely correlated with the first PC, the enrichment for the variants extremely correlated with PC2 outliers is significant when using different thresholds to define outliers (table 2.1).

Table 2.1: Enrichment Measured with Odds Ratio (OR) of the Variants Most Correlated with the Principal Components Obtained from the 1000 Genomes Data. Enrichment significant at the 1% (resp. 5%) level are indicated with ** (resp. *).
	Top 0.1%	Top 0.01%	Top 0.005%	Top 100 SNPs
pc1–genic/nogenic	1.60*	1.24	1.09	1.93
pc1–nonsyn/all	1.70	1.18	2.42	10.07*
pc1–UTR/all	1.37	0.80	1.65	3.44
pc2–genic/nogenic	1.51*	2.27	4.73**	4.44*
pc2–nonsyn/all	1.72	4.66*	7.40	12.18*
pc2–UTR/all	1.68	4.01*	3.36	2.73

Discussion

The promise of a fine characterization of natural selection in humans fostered the development of new analytical methods for detecting candidate genomic regions (Vitti, Grossman, & Sabeti, 2013). Population-differentiation based methods such as genome scans based on $F_{ST}$ look for marked differences in allele frequencies between population (Holsinger & Weir, 2009). Here, we show that the communality statistic $h^2$, which measures the proportion of variance of a SNP that is explained by the first $K$ principal components, provides a similar list of outliers than the $F_{ST}$ statistic when there are $K + 1$ populations. In addition, the communality statistic $h^2$ based on PCA can be viewed as an extension of $F_{ST}$ because it does not require to define populations in advance and can even be applied in the absence of well-defined populations.

To provide an example of genome scans based on PCA when there are no clusters of populations, we additionally consider the POPRES data consisting of 447,245 SNPSs typed for 1,385 European individuals (Nelson et al., 2008). The scree plot indicates that there are $K = 2$ relevant clusters (supplementary fig. B.3). The first principal component corresponds to a Southeast–Northwest gradient and the second one discriminates individuals from Southern Europe along a East–West gradient (Jay, Sjödin, Jakobsson, & Blum, 2012; Novembre et al., 2008) (fig. 2.16). Considering the 100 SNPs most correlated with the first PC, we find that 75 SNPs are in the lactase region, 18 SNPs are in the HLA region, 5 SNPs are in the ADH1C gene, 1 SNP is in HERC2, and another is close to the LOC283177 gene (fig. 2.17). When considering the 100 SNPs most correlated with the second PC, we find less clustering than for PC1 with more peaks (supplementary fig. B.13). The regions that contain the largest number of SNPs in the top 100 SNPs are the HLA region (41 SNPs) and a region close to the NEK10 gene (10 SNPs), which is a gene potentially involved in breast cancer (Ahmed et al., 2009). The genome scan retrieves well-known signals of adaption in humans that are related to lactase persistence (LCT) (Bersaglieri et al., 2004), immunity (HLA), alcohol metabolism (ADH1C) (Han et al., 2007), and pigmentation (HERC2) (Wilde et al., 2014). The analysis of the POPRES data shows that genome scan based on PCA can be applied when there is a clinal or continuous pattern of population structure without well-defined clusters of individuals.

$PCA with $K = 2$ applied to the POPRES data.$

Figure 2.16: PCA with $K = 2$ applied to the POPRES data.

$Manhattan plot for the POPRES data of the squared loadings $\rho_{j1}^2$ with the first principal component. For sake of presentation, only the top-ranked SNPs (top 5%) are displayed and the 100 top-ranked SNPs are colored in red.$

Figure 2.17: Manhattan plot for the POPRES data of the squared loadings $\rho_{j1}^2$ with the first principal component. For sake of presentation, only the top-ranked SNPs (top 5%) are displayed and the 100 top-ranked SNPs are colored in red.

When there are clusters of populations, we have shown with simulations that genome scans based on $F_{ST}$ can be reproduced with PCA. Genome scans based on PCA have the additional advantage that a particular axis of genetic variation, which is related to adaptation, can be pinpointed. Bearing some similarities with PCA, performing a spectral decomposition of the kinship matrix has been proposed to pinpoint populations where adaptation took place (Fariello et al., 2013). However, despite of some advantages, the statistical problems related to genome scans with $F_{ST}$ remain. The drawbacks of $F_{ST}$ arise when there is hierarchical population structure or range expansion because $F_{ST}$ does not account for correlations of allele frequencies among subpopulations (Bierne, Roze, & Welch, 2013; Lotterhos & Whitlock, 2014). An alternative presentation of the issues arising with $F_{ST}$ is that it implicitly assumes either a model of instantaneous divergence between populations or an island-model (Bonhomme et al., 2010). Deviations from these models severely impact FDRs (Duforet-Frebourg et al., 2014). Viewing $F_{ST}$ from the point of view of PCA provides a new explanation about why $F_{ST}$ does not provide an optimal ranking of SNPs for detecting selection. The statistic $F_{ST}$ or the proposed $h^2$ communality statistic are mostly influenced by the first principal component and the relative importance of the first PC increases with the difference between the first and second eigenvalues of the covariance matrix of the data. Because the first PC can represent ancient adaptive events, especially under population divergence models (G. McVean, 2009), it explains why $F_{ST}$ and the communality $h^2$ are biased toward ancient evolutionary events. Following recent developments of $F_{ST}$-related statistics that account for hierarchical population structure (Bonhomme et al., 2010; Foll et al., 2014; Günther & Coop, 2013), we proposed an alternative statistic ${h^{\prime}_j}^2$, which should give equal weights to the different PCs. However, analyzing simulations and the 1000 Genomes data show that ${h^{\prime}_j}^2$ do not properly account for hierarchical population structure because outliers identified by ${h^{\prime}_j}^2$ are almost always related to the last PC kept in the analysis. To avoid to bias data analysis in favor of one principal component, it is possible to perform a genome scan for each principal component.

In addition to ranking the SNPs when performing a genome scan, a threshold should be chosen to extract a list of outlier SNPs. We do not have addressed the question of how to choose the threshold and rather used empirical threshold such as the 99% quantile of the distribution of the test statistic (top 1%). If interested in controlling the FDR, we can assume that the loadings $\rho_{kj}$ are Gaussian with zero mean (Galinsky et al., 2016). Because of the constraints imposed on the loadings when performing PCA, the variance of the $\rho_{kj}$’s is equal to the proportion of variance explained by the $k^{th}$ PC, which is given by $\lambda_k / (p \times (n - 1))$ where $\lambda_k$ is the $k^{th}$ eigenvalue of the matrix $YY^T$. Assuming a Gaussian distribution for the loadings, the communality can then be approximated by a weighted sum of chi-square distribution. Approximating a weighted sum of chi-square distribution with a chi-square distribution, we have (Yuan & Bentler, 2010)

\[\begin{equation} h^2 \times K / c \leadsto \chi^2_K, \tag{2.8} \end{equation}\]

where $c = \sum_{k=1}^K \lambda_k / (p \times (n - 1))$ is the proportion of variance explained by the first $K$ PCs. The chi-square approximation of equation (2.8) bears similarity with the approximation of Lewontin & Krakauer (1973) that states that $F_{ST} \times (n_{\text{pops}}-1)/\bar{F}_{ST}$ follows a chi-square approximation with $(n_{\text{pops}}-1)$ degrees of freedom where $\bar{F}_{ST}$ is the mean $F_{ST}$ over loci and $(n_{\text{pops}}-1)$ is the number of populations. In the simulations of an island model and of a divergence model, quantile-to-quantile plots indicate a good fit to the theoretical chi-square distribution of expression (2.8) (supplementary fig. B.14). When using the chi-square approximation to compute $P$-values, we evaluate if FDR can be controlled using Benjamini–Hochberg correction (Benjamini and Hochberg 1995). We find that the actual proportion of false discoveries corresponds to the target FDR for the island model but the procedure is too conservative for the divergence model (supplementary fig. B.15). For instance, when controlling FDR at a level of 25%, the actual proportion of false discoveries is of 15%. A recent test based on $F_{ST}$ and a chi-square approximation was also found to be conservative (Whitlock & Lotterhos, 2015). Analysing the phase 1 release of the 1000 Genomes data demonstrates the suitability of a genome scan based on PCA to detect signals of positive selection. We search for variants extremely correlated with the first PC, which corresponds to differentiation between Africa and Eurasia and with the second PC, which corresponds to differentiation between Europe and Asia. For variants most correlated with the second PC, there is a significant enrichment of genic and nonsynonymous SNPs whereas the enrichment is less detectable for variants related to the first PC. The enrichment analysis confirms that positive selection may favor local adaptation of human population by increasing differentiation in genic regions especially in nonsynonymous variants (Barreiro et al., 2008). Consistent with LD, we find that candidate variants are clustered along the genome with a larger clustering for variants correlated with the Europe–Asia axis of differentiation (PC2). The difference of clustering illustrates that statistical methods based on LD for detecting selection will perform differently depending on the time frame under which adaptation had the opportunity to occur (Sabeti et al., 2006). The fact that population divergence, and its concomitant adaptive events, between Europe and Asia is more recent that the out-of-Africa event is a putative explanation of the difference of clustering between PC1 and PC2 outliers. Explaining the difference of enrichment between PC1 and PC2 outliers is more difficult. The weaker enrichment for PC1 outliers can be attributed either to a larger number of false discoveries or to a larger importance of other forms of natural selection such as background selection (Hernandez et al., 2011).

When looking at the 100 SNPs most correlated with PC1 or PC2, we find genes for which selection in humans was already documented (9/24 for PC1 and 5/14 for PC2). Known targets for selection include genes involved in pigmentation (MATP, OCA2 for PC1 and SLC45A2, SLC24A5, and MYO5C for PC2), in the regulation of sweating (EDAR for PC2), and in adaptation to pathogens (DARC, SLC39A4, and VAV2 for PC1). A 100 kb region in the vicinity of the APPBPP2 gene contains one-third of the 100 SNPs most correlated with PC1. This APPBPP2 region is a known candidate for selection and has been identified by looking for miRNA binding sites with extreme population differentiation (J. Li et al., 2012). APPBPP2 is a nervous system gene that has been associated with Alzheimer disease, and it may have experienced a selective sweep (Williamson et al., 2007). For some SNPs in APPBPP2, the differences of allele frequencies between Eurasiatic population and sub-Saharan populations from Africa are of the order of 90% (http://popgen.uchicago.edu/ggv/, last accessed December 2015) calling for a further functional analysis. Moreover, looking at the 100 SNPs most correlated with PC1 and PC2 confirms the importance of noncoding RNA (FAM230B, D21S2088E, LOC100133461, LINC00290, LINC01347, LINC00681), such as miRNA (MIR429), as a substrate for human adaptation (Grossman et al., 2013; J. Li et al., 2012). Among the other regions with a large number of candidate SNPs, we also found the RTTN/CD226 regions, which contain many SNPs correlated with PC1. In different selection scans, the RTTN genes has been detected (Barreiro et al., 2008; Carlson et al., 2005), and it is involved in the development of the human skeletal system (Wu & Zhang, 2010). An other region with many SNPs correlated with PC1 contains the ATP1A1 gene involved in osmoregulation and associated with hypertension (Gurdasani et al., 2015). The regions containing the largest number of SNPs correlated with PC2 are well-documented instances of adaptation in humans and includes the EDAR, SLC24A5, and SLC45A2 genes. The KCNMA1 gene contains seven SNPs correlated with PC2 and is involved in breast cancer and obesity (Jiao et al., 2011; Oeggerli et al., 2012). As for KCNMA1, the MYO5C has already been reported in selection scans although no mechanism of biological adaption has been proposed yet (H. Chen, Patterson, & Reich, 2010; Fumagalli et al., 2010). To summarize, the list of most correlated SNPs with the PCs identifies well-known genes related to biological adaptation in humans (EDAR, SLC24A5, SLC45A2, DARC), but also provides candidate genes that deserve further studies such as the APPBPP2, TP1A1, RTTN, KCNMA1, and MYO5C genes, as well as the ncRNAs listed above.

We also show that a scan based on PCA can also be used to detect more subtle footprints of positive selection. We conduct an enrichment analysis that detects polygenic adaptation at the level of biological pathways (Daub et al., 2013). We find that genes in the beta-defensin pathway are enriched in SNPs correlated with PC1. The beta-defensin genes are key components of the innate immune system and have evolved through positive selection in the catarrhine primate lineages (Hollox & Armour, 2008). As for the HLA complex, some beta-defensin genes (DEFB1, DEFB127) show evidence of long-term balancing selection with major haplotypic clades coexisting since millions of years (Cagliani et al., 2008; Hollox & Armour, 2008). We also find that genes in the omega fatty acid oxidation pathways are enriched in SNPs correlated with PC2. This pathway was also found when investigating polygenic adaptation to altitude in humans (Foll et al., 2014). The proposed explanation was that omega oxidation becomes a more important metabolic pathway when beta oxidation is defective, which can occur in case of hypoxia (Foll et al., 2014). However, this explanation is not valid in the context of the 1000 Genomes data when there are no populations living in hypoxic environments. Proposing phenotypes on which selection operates is complicated by the fact that the omega fatty acid oxidation pathway strongly overlaps with two other pathways: ethanol oxidation and glycolysis. Evidence of selection on the alcohol dehydrogenase locus have already been provided (Han et al., 2007) with some authors proposing that a lower risk for alcoholism might have been beneficial after rice domestication in Asia (Y. Peng et al., 2010). This hypothesis is speculative and we lack a confirmed biological mechanism explaining the enrichment of the fatty acid oxidation pathway. More generally, the enrichment of the beta-defensin and of the omega fatty acid oxidation pathways confirms the importance of pathogenic pressure and of metabolism in human adaptation to different environments (Barreiro & Quintana-Murci, 2010; Daub et al., 2013; Fumagalli et al., 2011; Hancock et al., 2008).

In conclusion, we propose a new approach to scan genomes for local adaptation that works with individual genotype data. Because the method is efficiently implemented in the software PCAdapt fast, analyzing 36,536,154 SNPs took only 502 min using a single core of an Intel(R) Xeon(R) (E5-2650, 2.00GHz, 64 bits). Even with low-coverage sequence data (3x), PCA-based statistics retrieve well-known examples of biological adaptation which is encouraging for future whole-genome sequencing project, especially for nonmodel species, aiming at sampling many individuals with limited cost.

Materials and Methods

Simulations of an Island Model

Simulations were performed with ms (Hudson, 2002). We assume that there are three islands with 100 sampled individuals in each of them. There is a total of 1,400 neutral SNPs, and 100 adaptive SNPs. SNPs are assumed to be unlinked. To mimic adaptation, we consider that adaptive SNP have a migration rate smaller than the migration rate of neutral SNPs ($4N_0m = 4$ for neutral SNPs) (Bazin et al. 2010). The strength of selection is equal to the ratio of the migration rates of neutral and adaptive SNPs. Adaptation is assumed to occur in one population only. The ms command lines for neutral and adaptive SNPs are given below (assuming an effective migration rate of $4N_0m = 0.1$ for adaptive SNPs).

The values of migrations rates we consider for adaptive SNPs are $4N_0m = 0.04, 0.1, 0.4, 1, 2$.

Simulations of Divergence Models

We assume that each population has a constant effective population size of $N_0 = 1000$ diploid individuals, with 50 individuals sampled in each population. The genotypes consist of 10,000 independent SNPs. The simulations were performed in two steps. In the first step, we used the software ms to simulate genetic diversity (Hudson, 2002) in the ancestral population. We kept only variants with a minor allele frequency larger than 5% at the end of the first step. The second step was performed with simuPOP (B. Peng & Kimmel, 2005) and simulations were started using the allele frequencies generated with ms in the ancestral population. Looking forward in time, we consider that there are 100 generations between the initial split and the following split between the two $B$ subpopulations, and 200 generations following the split between the two $B$ subpopulations. We assume no migration between populations. In the simulation of figure 2.11, we assume that 250 SNPs confer a selective advantage in the branch leading to population $A$ and 250 other SNPs confer a selective advantage in the branch leading to population $B_1$. We consider an additive model for selection with a selection coefficient of $s = 1.025$ for heterozygotes. For the simulation of figure 2.12, we assume that there are four nonoverlapping sets of 125 adaptive SNPs with each set being related to adaptation in one of the four branches of the divergence tree. A SNP can confer a selective advantage in a single branch only.

When including migration, we consider that there are 200 generations between the initial split and the following split between the two $B$ subpopulations, and 100 generations following the split between the two $B$ subpopulations. We consider migration rates ranging from 0.2% to 5% per generation. Migration is assumed to occur only after the split between $B_1$ and $B_2$. The migration rate is the same for the three pairs of populations. To estimate the $F_{ST}$ statistic, we consider the estimator of Weir & Cockerham (1984).

1000 Genomes Data

We downloaded the 1000 Genomes data (phase 1 v3) (Consortium & others, 2012). We kept low-coverage genome data and excluded exomes and triome data to minimize variation in read depth. Filtering the data resulted in a total of 36,536,154 SNPs that have been typed on 1,092 individuals. Because the analysis focuses on biological adaptation that took place during the human diaspora out of Africa, we removed recently admixed populations (Mexican, Columbian, PortoRican, and AfroAmerican individuals from the Southwest of the United States). The resulting data set contains 850 individuals coming from Asia (two Han Chinese and one Japanese populations), Africa (Yoruba and Luhya), and Europe (Finish, British in England and Scotland, Iberian, Toscan, and Utah residents with Northern and Western European ancestry).

Enrichment Analyses

We used Gowinda (Kofler & Schlötterer, 2012) to test for enrichment of GO. A gene is considered as a candidate if there is at least one of the most correlated SNPs (top 1%) that is mapped to the gene (within an interval of 50 kb upstream and downstream of the gene). Enrichment was computed as the proportion of genes containing at least one outlier SNPs among the genes of the given GO category that are present in the data set. In order to sample a null distribution for enrichment, Gowinda performs resampling without replacement of the SNPs. We used the -gene option of Gowinda that assumes complete linkage within genes.

We performed a second enrichment analysis to determine if outlier SNPs are enriched for genic regions. We computed odds ratio (Kudaravalli et al., 2008)

\[\begin{equation} \text{OR} = \frac{\text{Pr(genic|outlier)}}{\text{Pr(not genic|outlier)}} \frac{\text{Pr(not genic|not outlier)}}{\text{Pr(genic|not outlier)}} \tag{2.9} \end{equation}\]

We implemented a permutation procedure to test if an odds ratio is significantly larger than 1 (Fagny et al., 2014). The same procedure was applied when testing for enrichment of UTR regions (untranslated regions) and of nonsynonymous SNPs.

Polygenic Adaptation

To test for polygenic adaptation, we determined whether genes in a given biological pathway show a shift in the distribution of the loadings (Daub et al., 2013). We computed the SUMSTAT statistic for testing if there is an excess of selection signal in each pathway (Daub et al., 2013). We applied the same pruning method to take into account redundancy of genes within pathways. The test statistic is the squared loading standardized into a $z$-score (Daub et al., 2013). SUMSTAT is computed for each gene as the sum of test statistic of each SNP belonging to the gene. Intergenic SNPs are assigned to a gene provided they are situated 50kb up- or downstream. We downloaded 63,693 known genes from the UCSC website and we mapped SNPs to a gene if a SNP is located within a gene transcript or within 50kb of a gene. A total of 18,267 genes were mapped with this approach. We downloaded 2,681 gene sets from the NCBI Biosystems database. After discarding genes that were not part of the aforementioned gene list, removing gene sets with less than 10 genes and pooling nearly identical gene sets, we kept 1,532 sets for which we test if there was a shift of the distribution of loadings.

Acknowledgments

This work has been supported by the LabEx PERSYVAL-Lab (ANR-11-LABX-0025-01) and the ANR AGRHUM project (ANR-14-CE02-0003-01). POPRES data were obtained from dbGaP (accession number phs000145.v1.p1).