Heading

Estimation et partitionnement basé sur l’âge-.based partitioning of individual genomic inbreeding levels in the Belgian Blue Beef cattle population

Nous avons commencé par utiliser un modèle Mix14R (avec R k allant de 2 à 8192) pour estimer la proportion du génome appartenant à différentes classes HBD pour les 634 taureaux BBB (Fig. 1a), ce qui permet d’estimer le coefficient de consanguinité par rapport aux différentes populations de base comme expliqué dans Méthodes (Fig. 1b). En considérant toutes les classes de DBC, la fraction du génome qui est DBC (correspondant au coefficient de consanguinité estimé avec la population de base la plus éloignée) était égale à 0,303 en moyenne (allant de 0,258 à 0,375), avec une contribution majeure des classes de DBC avec des taux de R k élevés (R k > 256) qui représentent 71,4% de la proportion totale de DBC en moyenne. Ces petits ROH reflètent l’histoire de la population (consanguinité de fond et déséquilibre de liaison associés à la taille effective passée de la population (Ne)) mieux que la variation individuelle. Les classes associées à des taux R k plus faibles (c’est-à-dire avec des segments DHB plus longs) représentaient une proportion plus faible de la proportion totale de DHB (le coefficient de consanguinité moyen était égal à 0,054 et 0,087 en incluant les classes DHB avec R k ≤ 32 et R k ≤ 256, respectivement, et en fixant la population de base il y a environ 16 ou 128 générations) mais présentaient plus de variation entre les individus. Par exemple, le coefficient de consanguinité associé aux ancêtres communs remontant jusqu’à environ quatre générations (correspondant aux classes HBD avec R k ≤ 8) variait de 0,000 à 0,137. Pour les taureaux nés de 1980 à 2010, le pourcentage du génome dans les segments HBD a augmenté de 3,3 % (+ 0,11 % par an), c’est-à-dire d’environ 28 à 31 % (voir le fichier supplémentaire 1 : Fig. S1a). Cependant, la tendance pour les classes DBC plus récentes (R k ≤ 32) était plus prononcée (voir fichier additionnel 1 : Fig. S1b), c’est-à-dire de presque 0 à 6% (+ 0,20% par an) et correspondait plus étroitement à la tendance observée avec les coefficients de consanguinité basés sur le pedigree (voir fichier additionnel 1 : Fig. S1c). Les taureaux nés avant 1980 présentaient peu de signes d’autozygotie récente par rapport aux taureaux modernes.

Pour évaluer la contribution de chaque classe HBD au pourcentage du génome dans les segments HBD et à sa variation dans le bétail BBB, nous avons divisé la fraction totale du génome en classes de DBC en quatre classes principales (classes de DBC très récentes avec R k = 2 à 8, classes de DBC récentes avec R k = 16 à 64, classes de DBC anciennes avec R k = 128 à 512, et classes de DBC très anciennes avec R k = 1024 à 8192), chaque groupe ayant trois classes de DBC sauf le dernier avec quatre classes de DBC. La fraction moyenne du génome associée à chacune de ces classes principales (ordonnée de récente à ancienne) était égale à 0,027 (SD = 0,029), 0,041 (SD = 0,019), 0,054 (SD = 0,013) et 0,180 (SD = 0,011). Notons que les proportions élevées de segments HBD très récents sont mécaniquement associées à des proportions plus faibles de segments HBD très anciens (r = – 0,407) car les segments HBD récents masquent les segments HBD plus anciens. Bien que le pourcentage du génome dans les classes HBD associées à des ancêtres communs récents ne représente que 22,6% de l’autozygotie totale, il présente plus de variation individuelle que celui des classes plus anciennes (plus de 50% de la variance totale est associée aux classes HBD très récentes). En ajustant un modèle linéaire, nous avons estimé que les classes HBD très récentes représentent 59% de la variation totale de l’autozygotie et que l’ajout de classes HBD récentes au modèle augmente cette valeur à 83%. De même, les corrélations entre les coefficients de consanguinité mesurés par rapport à différentes populations de base (par exemple, en incluant différentes classes HBD dans le calcul) avec les coefficients de consanguinité estimés en utilisant toutes les classes HBD ont augmenté brusquement de 0,16 pour les estimations basées sur la première classe (R k = 2) à 0,77 pour les coefficients de consanguinité estimés en incluant les classes HBD avec un R k ≤ 8 et à 0,90 avec un R k ≤ 16, puis ne se sont améliorées que marginalement en ajoutant plus de classes HBD (Fig. 2). La diminution de la corrélation observée à R k = 1024 résulte du fait que l’autozygotie ancienne est concentrée à R k = 1024 pour certains individus et à R k = 2048 pour d’autres.

Comparaison des résultats pour les bovins BBB avec ceux d’autres races

Pour déterminer si des niveaux et des modèles d’autozygotie comparables sont également observés dans d’autres races d’origine européenne, nous avons appliqué le même modèle à 10 races génotypées avec la même matrice (Fig. 3). Dans la plupart de ces races, les coefficients de consanguinité estimés par rapport aux différentes populations de base ont augmenté modérément jusqu’à FG-256 (par exemple, la classe HBD avec R k ≤ 256 incluse dans l’estimation) et plus fortement avec les populations de base plus anciennes (FG-512 à FG-2048), qui comprennent beaucoup plus de générations d’ancêtres. De grandes différences dans les coefficients de consanguinité ont été observées avec des populations de base relativement récentes (FG-64, il y a environ 32 générations), allant de 0,013 et 0,042 chez les Piémontais et les Limousins à 0,164 et 0,200 chez les bovins Jersey et Hereford. Certains individus Hereford présentaient des coefficients de consanguinité extrêmes estimés avec des populations de base récentes (voir fichier additionnel 2), c’est-à-dire jusqu’à 40 % pour FG-8 (il y a environ quatre générations). Une partie des individus Hereford de cet ensemble de données provient de la lignée Hereford 1, une lignée consanguine, ce qui indique que notre modèle capture correctement les événements extrêmes mais aussi que les individus génotypés inclus dans cette étude ne sont pas nécessairement représentatifs de la race.

Estimation des coefficients de consanguinité et des probabilités de DHB avec différentes densités de SNP

Nous avons ajusté un modèle Mix14R en utilisant différentes densités de SNP, à savoir.e., de LD (6844 SNP) à HD (601 226 SNP) sur l’ensemble de données 634 BBB et même à WGS (5 653 911 SNP) pour les 50 individus séquencés sur le génome entier. Les coefficients de consanguinité estimés moyens mesurés par rapport aux différentes populations de base (Fig. 4) et Additional file 3 : Fig. S2 étaient similaires pour tous les panels de SNP pour les populations de base les plus récentes (FG-32). Pour les populations de base plus anciennes, moins d’autozygotie a été capturée avec le panel LD avec des différences marquées pour les classes HBD anciennes qui ont été capturées uniquement avec les panels HD ou WGS. Une tendance similaire a été observée avec le panel 50 K mais les coefficients de consanguinité moyens étaient similaires à ceux du panel HD jusqu’à FG-256 (environ 128 générations en arrière). Les coefficients de consanguinité moyens estimés en utilisant la population de base la plus éloignée et les panels LD, 50 K et HD étaient égaux à 0,060, 0,093 et 0,303, respectivement (lorsqu’ils étaient estimés sur les 50 individus séquencés seulement, ces valeurs étaient égales à 0,047, 0,101 et 0,309, respectivement, et à 0,359 avec le panel WGS). La population de base est donc fonction des plus petits segments HBD qui peuvent être capturés par le panel utilisé. Les corrélations entre ces coefficients de consanguinité estimés avec différents panels étaient élevées, c’est-à-dire 0,934 (LD-HD), 0,944 (LD-50 K) et 0,975 (50 K-HD). Malgré les coefficients de consanguinité beaucoup plus faibles obtenus avec le panel 50 K, il capture essentiellement toute la variation individuelle obtenue avec un panel HD, en accord avec l’observation précédente que la plupart de la variation était associée aux classes récentes de HBD.

Nous avons ensuite utilisé l’algorithme de Viterbi pour identifier les segments HBD avec différents panels de SNP (tableau 1). L’algorithme de Viterbi classe chaque position de SNP comme HBD ou non-HBD tandis que l’algorithme forward-backward fournit la probabilité locale de HBD. Comme prévu, des segments plus nombreux et plus courts de HBD sont capturés avec des panels de densité plus élevée. Avec le panel HD, une proportion limitée de segments extrêmement petits (quelques kb) a été capturée. La longueur de la majorité des segments était comprise entre 10 et 500 kb, plus de la moitié étant inférieure à 100 kb, mais ces segments n’ont pas nécessairement la plus grande contribution au pourcentage total du génome dans les classes HBD puisque les classes avec des segments moins nombreux mais plus longs peuvent représenter une grande proportion de l’autozygotie. Nous avons également observé des segments HBD extrêmement longs (> 50 Mb), ce qui a confirmé la présence d’une autozygotie récente (le segment HBD le plus long faisait plus de 90 Mb). En moyenne, chacun des 634 taureaux présentait 4,25 segments HBD d’une longueur supérieure à 10 Mb et associés à un ancêtre commun présent environ cinq générations en arrière. Le nombre de ces segments HBD variait de 0 à 14 par individu. Soixante et un taureaux avaient même un ou plusieurs (jusqu’à trois) segments HBD de plus de 50 Mb. Avec les panels 50 K et LD, plus de 99 % des segments identifiés avaient une longueur supérieure à 100 et 500 kb, respectivement (avec un pic dans les classes de 1 à 5 Mb et de 5 à 10 Mb, respectivement), et seule une fraction des segments a été capturée par rapport à l’utilisation du panel HD. En particulier, la grande majorité des segments HBD plus courts que 1 Mb n’ont pas été identifiés. À des densités de SNP plus faibles, les plus petits segments ne sont tout simplement pas capturés parce qu’ils ne contiennent aucun SNP ou trop peu. Les segments de taille intermédiaire peuvent ne pas atteindre des probabilités élevées de DHB en raison d’un plus petit nombre de SNP dans le segment. Inversement, la longueur de certains segments DHB peut être surestimée lors de l’utilisation du panel LD, par exemple lorsqu’il n’y a pas assez de SNP pour identifier les petits segments non DHB qui flanquent les segments DHB. La figure 5a illustre l’identification des segments HBD pour un chromosome. Elle montre que (1) plus de segments ont été identifiés à une densité plus élevée, (2) les probabilités HBD étaient plus élevées avec des cartes plus denses, (3) l’algorithme de Viterbi a déclaré certaines positions de SNP comme HBD bien qu’elles n’aient eu que des probabilités HBD modérées, et (4) les limites des segments HBD ont varié avec la densité du panel. De même, la figure 5b représente les segments HBD qui ont été identifiés sur le chromosome 5 de Bos taurus (BTA) pour 50 individus avec l’algorithme de Viterbi avec différentes densités de SNP. Les résultats sont en accord avec ceux rapportés dans le tableau 1. De plus grandes proportions du génome ont été déclarées HBD avec le panel HD et les petits segments HBD ont représenté la plus grande partie de la différence avec les résultats des panels de densité inférieure. Néanmoins, nous avons observé que certains segments HBD de quelques Mb de long n’ont pas été identifiés avec une densité de SNP plus faible (et encore plus avec le panel LD). Comme pour la Fig. 5a, la longueur de certains segments HBD est surestimée lorsque le panel LD est utilisé. Nous avons également comparé les probabilités de DBC locales estimées en utilisant le panel LD ou le panel 50 K avec les classes de DBC locales déduites en utilisant le panel HD et l’algorithme de Viterbi (Fig. 6). Les probabilités HBD étaient élevées pour les classes HBD récentes et diminuaient pour les ancêtres communs plus éloignés. Comme prévu, le panel LD n’était efficace que pour les ancêtres communs les plus récents (la probabilité de DBC était de 0,90 ou plus lorsque R k < 16 et ~ 0,50 pour R k = 32), tandis que le panel 50 K permettait la capture d’une autozygotie plus ancienne (la probabilité de DBC était de 0,90 ou plus lorsque R k < 64 et ~ 0,50 pour R k = 128). D’autres résultats concernant l’âge (ou la longueur) des segments HBD qui peuvent être capturés avec différentes densités de SNP sont décrits dans Druet et Gautier .

Comparaison des modèles

Modèles qui estiment les taux R k des classes HBD (modèles KR)

Pour les différentes densités de SNP testées et pour chaque individu, nous avons utilisé le BIC (voir ) pour sélectionner le modèle KR avec le meilleur support statistique (i.e., avec le nombre optimal de classes K, avec K – 1 classes HBD et une classe non-HBD) après avoir estimé le(s) taux R k pour chaque individu avec chaque modèle testé. Pour chaque panel de SNP, le tableau 2 montre le nombre de fois où un modèle a été sélectionné comme le meilleur pour l’individu analysé. Plus la densité de SNP augmente, plus le nombre de générations passées peut être exploré et le K optimal augmente en conséquence. Dans la plupart des cas, les modèles avec une classe HBD sont préférés pour le panel LD, les modèles avec deux classes HBD pour le panel 50 K, les modèles avec trois classes HBD pour les panels HD et WGS (bien que le modèle avec quatre classes HBD soit aussi souvent sélectionné pour ce dernier, c’est-à-dire pour 23 des 50 individus). Avec ces modèles optimaux, la première classe HBD capture l’autozygotie la plus récente (R k de 15 à 20), la deuxième classe HBD capture l’autozygotie qui est associée à des ancêtres communs de quelques centaines de générations en arrière et les classes ultérieures sont associées à des R k plus élevés (> 1000) (tableau 2). Les corrélations des coefficients de consanguinité estimés avec ces modèles KR sélectionnés avec ceux obtenus avec le modèle Mix14R complet (allant de 0,981 à 1,000) et la comparaison des coefficients de consanguinité estimés moyens indiquent qu’avec ces modèles KR réduits, nous pouvons capturer efficacement l’autozygotie à l’échelle du génome. Avec les modèles 1R et des densités de SNP faibles ou modérées, nous avons observé une légère sous-estimation des coefficients de consanguinité par rapport au modèle Mix14R et des corrélations légèrement inférieures (toujours supérieures à 0,98). Les taux de R k estimés pour chaque individu avec ces panels ont une valeur médiane plus faible (respectivement 15 et 41 avec les panels LD et 50 K) que les taux de R k estimés avec des panels de densité plus élevée (médiane R k > 1000) pour lesquels la contribution des petits ROH est beaucoup plus importante. Par conséquent, certains petits fragments n’ont pas été capturés par le modèle à faible densité alors qu’à densité plus élevée, les coefficients de consanguinité sont presque identiques aux estimations obtenues avec le modèle Mix14R. Les modèles contenant deux classes HBD ou plus ont capturé la même quantité d’autozygotie que le modèle Mix14R, quelle que soit la densité de SNP. Bien que le coefficient de consanguinité soit correctement estimé avec un modèle 1R (une classe HBD et une classe non-HBD avec le même taux) avec les données WGS, les segments HBD identifiés ont tendance à être plus petits puisque les taux R k estimés sont plus élevés (c’est-à-dire que les longueurs attendues des fragments sont plus petites) comme le montre le fichier additionnel 4 : Fig. S3. En effet, le modèle 1R produit plus de segments de 10 à 100 kb que le modèle Mix14R, mais moins de segments de plus de 100 kb. Ainsi, avec un modèle 1R, les longs segments HBD pourraient être coupés en fragments plus petits en présence de SNP hétérozygotes (peut-être des erreurs de séquençage), alors qu’avec des modèles incluant une ou plusieurs classes HBD associées à des ancêtres communs récents (avec de petits taux R k), ces segments HBD sont identifiés comme un seul fragment long et récent (parce que la pénalité pour terminer et commencer un nouveau segment est plus élevée). La figure 7 illustre ce phénomène à l’aide d’un exemple. En effet, nous avons observé un long segment avec des probabilités HBD élevées bien qu’il y ait de multiples positions où la probabilité du génotype hétérozygote est non nulle (mais ceci est limité par rapport aux régions flanquantes). Avec le modèle Mix14R, on considère qu’il s’agit d’un long segment et la probabilité locale de DBC reste supérieure à 0,99 pour toute la région (sauf pour une région avec cinq SNP hétérozygotes consécutifs). Avec le modèle 1R, les probabilités HBD chutent de manière répétée en raison de ces SNP éventuellement hétérozygotes et le segment HBD le plus long est coupé en plusieurs fragments plus petits (sur la base des résultats de l’algorithme de Viterbi). Notez qu’avec le panel HD, cet individu est homozygote pour les 13 009 SNP qui sont inclus dans ce segment de 56,1 Mb de long. Comme dans la figure 5, nous remarquons que l’algorithme de Viterbi classe certaines positions avec une faible probabilité estimée de DHB comme DHB.

Modèles avec des taux R k prédéfinis de classes HBD (modèles MixKR)

Par rapport aux modèles KR, les modèles MixKR présentent l’avantage d’utiliser les mêmes classes HBD pour tous les individus (les taux R k des classes HBD ne sont pas estimés individuellement mais prédéfinis par l’utilisateur) et facilitent les comparaisons entre individus (par exemple, comparer deux individus avec une seule classe HBD mais avec R k = 8 pour le premier et R k = 96 pour le second ne serait pas facile – les R k estimés vont de 4 à 1000). Plusieurs de ces modèles MixKR (avec K = 2, 3 et 4) ont été testés avec le panel LD (tableau 3) pour évaluer si les modèles réduits avec des taux prédéfinis de classes HBD sont efficaces. Pour sélectionner ces taux prédéfinis, soit nous avons utilisé les médianes des taux estimés obtenus à partir de modèles avec le même nombre de classes (voir section précédente), soit nous avons sélectionné quelques classes du modèle MixKR afin de couvrir la gamme des valeurs estimées (par exemple, une classe pour les segments HBD récents et une pour les segments HBD anciens). En accord avec les observations précédentes sur les modèles KR, les comparaisons des coefficients de consanguinité estimés avec ceux obtenus avec le modèle Mix14R indiquent que les modèles avec une seule classe HBD sous-estiment légèrement les coefficients de consanguinité et entraînent des corrélations plus faibles (> 0,96) que les modèles avec deux classes HBD ou plus (> 0,99). La présence de plusieurs classes HBD (> 2) permet de mieux évaluer les contributions des différentes générations passées (ex, R k = 8 vs 64) mais ne fournit pas de meilleures estimations du coefficient de consanguinité à l’échelle du génome.

Comparaison avec d’autres estimateurs de coefficients de consanguinité

Les moyennes et les plages des coefficients de consanguinité estimés avec différentes méthodes et le panel HD sont dans le tableau 4 et leurs corrélations sont dans le tableau 5, et dans le fichier additionnel 5 : Tableaux S1 et S2 pour les autres panels. Comme notre modèle, les modèles basés sur l’homozygotie observée et le ROH ont donné lieu à des coefficients de consanguinité élevés (respectivement, 0,644 et 0,151 en moyenne) alors que les autres estimateurs génomiques ont donné lieu à des coefficients de consanguinité centrés sur 0 et incluant des valeurs négatives. Il faut noter que des valeurs plus élevées sont obtenues en moyenne (0,268) lorsqu’on utilise des règles moins strictes pour identifier les ROH (par exemple, des fenêtres de 20 SNP et au moins 10 SNP par ROH). Nous avons observé des corrélations très élevées entre les estimations basées sur les HMM et les mesures basées sur l’homozygotie (r = 0,95 avec FHOM et FExHOM, ces deux mesures présentant une corrélation de 1 et étant essentiellement les mêmes) ou sur les ROH (r = 0,95 avec FROH), ce qui suggère qu’avec un grand nombre de SNP, les heuristiques simples (ignorant les fréquences alléliques, l’espacement des SNP, etc.) sont efficaces (FHOM et FROH étant fortement corrélées, r = 0,97). La corrélation entre les FHOM estimés avec les panels LD et HD est égale à 0,890, ce qui est légèrement inférieur à la corrélation entre les estimations obtenues avec le HMM pour ces deux panels (r = 0,934), ce qui indique que les estimateurs globaux fonctionnent toujours correctement avec 6844 SNP dans cette population. Les méthodes ROH basées sur des règles sont moins efficaces à des densités de SNP plus faibles car elles ne capturent que les fragments les plus longs (5 Mb ou plus et 20 Mb en moyenne) avec les paramètres utilisés dans l’étude actuelle (la taille des fenêtres par défaut dans plink). En fait, les estimateurs basés sur les ROH sont rarement utilisés avec le panel LD chez les bovins, bien que davantage de segments HBD puissent être identifiés avec des règles moins strictes, au prix d’un taux accru de faux positifs. À faible densité de SNP, le cadre HMM fournit encore des probabilités globales et locales correctes de HBD bien que les segments HBD ne soient pas identifiés sans ambiguïté .

Les corrélations des estimations du GRM traditionnel avec nos estimations sont modérément élevées (r = 0.73) et plus faibles avec les estimateurs d’homozygotie (r = 0,63) et les estimateurs basés sur le ROH (0,61). Le FGRM a été calculé avec la formule proposée par , qui divise toutes les contributions des SNP par le même poids. Lorsqu’il est estimé avec la formule alternative, qui divise chaque contribution SNP par son propre poids 2f i (1 – f i ) (f i étant la fréquence du SNP i) comme dans Amin et al. , les corrélations étaient plus faibles (c’est-à-dire 0,48 avec FG, 0,34 avec FHOM et 0,33 avec FROH). L’estimateur basé sur les corrélations unifiées entre les gamètes proposées par Yang et al. a présenté des corrélations relativement élevées avec FG et FGRM (respectivement, 0,90 et 0,92) et des corrélations légèrement plus faibles avec les autres estimateurs (r = 0,87 et 0,85 avec FHOM et FROH, respectivement).

Les corrélations de ces estimations avec les coefficients de consanguinité du pedigree (en ne considérant que les individus nés après 1999 pour augmenter la profondeur du pedigree) sont également dans le tableau 5. Les corrélations globales étaient modérées, avec les valeurs les plus élevées pour les corrélations avec les mesures basées sur l’homozygotie et le ROH (0,55 pour les deux mesures) et des valeurs légèrement inférieures pour celles avec l’estimateur basé sur le HMM (0,46), tandis que nous avons observé une faible relation avec le FGRM (0,29) et une corrélation modérée avec le FUNI (0,45). Nous avons également comparé la FPED et les coefficients de consanguinité estimés avec notre modèle par rapport à différentes populations de base (Fig. 8) et avons constaté que les corrélations augmentaient jusqu’à FG-32 (capturant la consanguinité à partir d’ancêtres remontant à environ 16 générations) et présentaient ensuite un plateau de FG-32 à FG-256 atteignant un maximum à r = 0,56 (c’est-à-dire légèrement meilleur que les estimateurs basés sur l’homozygotie). Cette tendance était attendue puisque la FPED est estimée pour un nombre limité de générations en arrière dans le temps. Le nombre équivalent moyen de générations connues estimé avec PEDIG était de 6,3 pour les taureaux nés après 1999 (il est passé de 5,5 pour les taureaux nés en 2000 à 7,5 pour ceux nés en 2010) correspondant en moyenne à FG-16. L’ajout de la classe HBD R k = 32 permet de capturer les contributions de certaines branches plus anciennes du pedigree et les plus petits segments HBD hérités d’ancêtres communs dans le pedigree.

.