Estimation de la fréquence des allèles et de la cartographie d’association à l’aide de données de séquençage de nouvelle génération

L’allèle mineur est l’allèle le moins fréquent dans la population à un site variable. Nous décrivons d’abord deux approches principales pour estimer la fréquence des allèles mineurs (CRG) à un site particulier du génome. La première approche consiste à déduire des génotypes individuels et à traiter ces génotypes inférés comme étant complètement exacts lors de l’estimation du CRG. Nous examinons ensuite la performance d’un cadre de vraisemblance qui prend directement en compte l’incertitude dans l’attribution des génotypes. Tout au long de notre travail, nous supposons que tous les sites de ségrégation sont bialléliques.

Estimation du CRG à partir de génotypes appelés

Une façon d’estimer le CRG à partir de données de séquençage de nouvelle génération consiste à appeler d’abord un génotype pour chaque individu à l’aide de données de séquençage, puis à utiliser ces génotypes comme s’ils étaient les vrais. C’était l’approche traditionnellement utilisée pour les données de génotype et les données de séquençage de Sanger. Il n’est pas clair dans quelle mesure il fonctionnera lorsqu’il sera appliqué aux données de séquençage de nouvelle génération.

Une approche du maximum de vraisemblance peut être utilisée pour déduire le génotype de chaque individu à partir des données de séquençage de la prochaine génération. À chaque site j, pour chaque individu i, la probabilité pour chacun des trois génotypes possibles (en supposant que nous connaissions l’allèle mineur) est donnée comme suit:

(1)

où D i, j est les données de séquençage observées chez l’individu i au site j, g i, j ∈ {0, 1, 2} est le nombre d’allèles mineurs contenus dans le génotype de chaque individu, et et contrôlent respectivement les erreurs de séquençage et les qualités de base de lecture. Les données de séquençage observées pour chaque individu peuvent être considérées comme l’alignement des lectures au site j en tenant compte des scores de qualité de lecture. Ceci est représenté comme la probabilité de génotype et se trouve dans le fichier de probabilité de génotype (GLF) qui est produit dans de nombreux programmes qui analysent les données de séquençage de nouvelle génération, tels que SOAPsnp et MAQ.

Pour attribuer un génotype à un individu particulier, la probabilité de chacun des trois génotypes possibles peut être calculée pour l’individu. Le génotype avec la probabilité la plus élevée peut alors être attribué. Cependant, les chercheurs préfèrent souvent un critère d’appel plus rigoureux et n’attribueront pas de génotype à un individu à moins que le génotype le plus probable ne soit beaucoup plus probable que le second. Ici, les trois génotypes possibles sont triés par leurs probabilités : , où g(k) correspond au génotype avec la k th plus grande probabilité. Avec un seuil f donné, on peut appeler le génotype g(1) if . Sinon, un génotype n’est pas appelé et le génotype de l’individu est considéré comme manquant. Une valeur seuil commune de f est de 1, ce qui indique que le génotype le plus probable est au moins 10 fois plus probable que le deuxième plus probable. Notez que ce type de filtrage peut entraîner une plus grande confiance pour le génotype « appelé », mais il en résulte également plus de données manquantes.

Estimateur du maximum de vraisemblance de la fréquence des allèles

Au lieu d’estimer le CRG à partir des génotypes appelés, une méthode du maximum de vraisemblance (ML) introduite par Kim et al. (voir aussi Lynch pour une approche similaire) estime directement les MAF et prend en compte l’incertitude du génotype. Plus précisément, étant donné un allèle mineur, la probabilité d’observer les données de séquence à chaque individu i est obtenue en additionnant les probabilités correspondant aux trois génotypes possibles.

Supposons que les trois probabilités de génotype définies dans l’équation 1 sont disponibles. En utilisant la même notation que ci-dessus, soient respectivement D j et p j les données de séquençage observées au site j et le CRG correspondent. La probabilité de génotype étant donné que la fréquence des allèles mineurs peut être calculée en supposant l’équilibre de Hardy-Weinberg (HWE). Ensuite, en supposant l’indépendance entre les individus, la probabilité du CRG à ce locus est un produit de toutes les probabilités calculées sur tous les N individus :

(2)

L’estimation ML de p j peut être calculée soit en maximisant directement la probabilité pour un espace de paramètres restreint à l’aide du Broyden-Fletcher – Méthode Goldfarb-Shanno (BFGS) ou en utilisant l’algorithme de maximisation des attentes (EM). Lors de l’utilisation de l’algorithme EM, l’attente postérieure d’un génotype est calculée pour chaque individu et la moyenne de ces postérieurs est mise à jour à plusieurs reprises. Notre implémentation de BFGS était plus rapide que l’algorithme EM. Par exemple, pour obtenir des estimations à partir de 100 000 sites, les BFG ont pris environ 16 secondes, mais les EM ont pris environ 100 secondes. Cependant, la différence de vitesse peut être spécifique à la mise en œuvre. Dans notre cas, pour les deux méthodes, nous avons cessé de mettre à jour les paramètres lorsque l’augmentation de la probabilité était inférieure à 0,001.

Estimateur du maximum de vraisemblance avec un allèle mineur incertain

En pratique, le deuxième nucléotide le plus courant chez les individus peut souvent être utilisé comme allèle mineur. Cependant, pour les SNP rares (par exemple, CRG < 1%), il est difficile de déterminer quel allèle est l’allèle mineur, car les quatre nucléotides peuvent apparaître dans certaines lectures en raison d’erreurs de séquençage. Pour faire face à cette situation, nous décrivons maintenant un cadre de vraisemblance qui tient compte de l’incertitude dans la détermination de l’allèle mineur.

Supposons que pour le site j on connaisse l’allèle majeur M. Notez que décider lequel des deux allèles communs est susceptible d’être le plus important n’est pas important car nous nous intéressons principalement à l’estimation des fréquences à des SNP rares. De plus, pour les allèles de fréquences intermédiaires (environ 50%), la distinction entre allèle majeur et allèle mineur est moins importante. Assignez les trois autres nucléotides non majeurs m1, m2 et m3. La probabilité introduite dans l’équation 2 suppose un allèle majeur fixe M et un allèle mineur fixe m. Par conséquent, pour tenir compte de l’incertitude dans la désignation de l’allèle mineur, la fonction de vraisemblance peut être modifiée comme suit:

(3)

De plus, en supposant que l’un des trois allèles mineurs possibles est également probable, nous obtenons :

(4)

div>

. Étant donné que peut être très petit avec de gros ensembles de données (par exemple, avec de nombreuses personnes), il est utile de calculer la probabilité dans l’échelle log. Ordonnez les trois probabilités log conditionnelles comme (l(1), l(2), l(3)), où l(1) est la plus grande. Ensuite,

Test G utilisant des génotypes appelés pour la cartographie d’association

Dans les études d’association, les SNP montrant des différences significatives dans la fréquence des allèles entre les cas et les témoins sont censés être associés au phénotype d’intérêt. La cartographie d’association peut être réalisée à l’aide des données d’études de séquençage de nouvelle génération. Nous discutons d’abord des approches qui nécessitent d’appeler des génotypes individuels, puis effectuons un test d’association à l’aide des génotypes appelés. Dans cette approche, un génotype est d’abord appelé pour chaque individu. Les génotypes peuvent être filtrés ou non filtrés. En supposant l’indépendance entre les individus et les HWE, un tableau de contingence 2 × 2 peut être construit en comptant le nombre d’allèles majeurs et mineurs dans les cas et les témoins. Ceci conduit au test d’indépendance du rapport de vraisemblance bien connu, le test G :

(5)

où O k, h est la fréquence observée dans une cellule, et E k, h est la fréquence attendue dans l’hypothèse nulle dans laquelle la fréquence de l’allèle est la même entre les cas et les contrôles. Le test du chi carré bien connu de Pearson est asymptotiquement équivalent au test G. Si la table est générée à partir de génotypes vrais, alors la statistique G suit asymptotiquement une distribution du chi carré avec 1 degré de liberté (χ2(1)). Cependant, dans nos études, nous construisons la statistique G en utilisant des génotypes « appelés », donc HWE peut ne pas tenir en raison de la sur- et de la sous-appel des hétérozygotes. De plus, la construction de la statistique de test en comptant les génotypes « appelés » au lieu des génotypes « observés » introduit probablement une variabilité supplémentaire. Par conséquent, la théorie statistique peut ne plus être valable. Notez que lorsqu’un génotype n’est pas appelé pour un certain individu, les données sont considérées comme manquantes et ne sont pas incluses dans le tableau 2 × 2.

Test de rapport de vraisemblance tenant compte de l’incertitude dans les génotypes observés pour la cartographie d’association

Au lieu d’appeler des génotypes, le cadre de vraisemblance permet une incertitude dans les génotypes et les tests à chaque site j si la fréquence des allèles est la même entre les cas et les témoins. Formellement, on calcule la vraisemblance des hypothèses H O : p j, 1 = p j, 2 (= p j, 0) et H A : p j,1 ≠ p j, 2, où p j,1 et p j, 2 sont les CRG dans les cas et les contrôles, respectivement.

En supposant que les allèles mineurs (m) et majeurs (M) sont connus, la probabilité de la fréquence des allèles mineurs peut être calculée comme décrit dans l’équation 2, et la statistique du test du rapport de vraisemblance est calculée comme suit:

(6)

et sont les données observées pour les cas et les contrôles, respectivement, et et sont les MLE des MAF dans les cases et les contrôles, respectivement.

Si l’allèle mineur est inconnu, la probabilité sous l’hypothèse nulle est calculée comme dans l’équation 3, et la statistique LRT est modifiée comme suit:

(7)

où D j est les données observées pour les cas et les témoins, et est la fréquence de l’allèle sous l’hypothèse nulle. Les autres notations sont les mêmes que dans l’équation 6.

Estimation du CRG dans des données simulées

Nous comparons les estimations de la fréquence des allèles sur des données simulées en utilisant des génotypes vrais (True), appelés génotypes sans filtrage (Call NF), appelés génotypes avec filtrage (f = 1; Call F), et la méthode du maximum de vraisemblance (ML). Pour les SNP rares, le type d’allèle mineur n’est souvent pas apparent. Lors de l’appel de génotypes, le deuxième nucléotide le plus courant est supposé être l’allèle mineur. La méthode ML intègre directement l’incertitude dans la détermination de l’allèle mineur et, sauf indication contraire, les résultats utilisant la méthode des allèles mineurs inconnus (équation 3) sont présentés. Notez que la méthode ML de l’allèle mineur inconnu fonctionne de manière similaire à la méthode ML de l’allèle mineur connue mais la première est meilleure pour les SNP très rares (fichier supplémentaire 1).

Nous avons d’abord évalué dans quelle mesure les différentes approches ont pu estimer le CRG chez 200 individus à travers une gamme de profondeurs de séquençage pour 1 000 SNP avec un CRG réel de 5%. La figure 1 montre les diagrammes des distributions des CRG estimées en utilisant les quatre approches différentes. Comme prévu, pour des données de couverture plus élevées, telles qu’une profondeur individuelle de 12×, toutes les méthodes fonctionnent aussi bien que lorsque les génotypes sont connus avec certitude (Vrai). Cependant, lorsque la profondeur diminue, les estimations du CRG obtenues en appelant d’abord les génotypes deviennent biaisées. Par exemple, le CRG médian estimé à l’aide de la méthode Call F est de 5,3 % pour une couverture de 6 × et de 12,5 % pour une couverture de 2×. La raison du biais à la hausse est qu’il devient plus difficile d’appeler hétérozygotes car les vrais hétérozygotes ressemblent souvent à des erreurs de séquençage. Par conséquent, plus d’hétérozygotes que d’homozygotes mineurs ont tendance à avoir des génotypes manquants. Cependant, le biais global dans les estimations du CRG à partir des génotypes appelés n’est pas toujours dans une direction (données non présentées). Fait intéressant, le biais semble être pire pour la méthode Call F que pour la méthode Call NF. Ce schéma peut sembler contre-intuitif car le filtrage des appels de génotype semble diminuer la probabilité d’appeler une erreur de séquençage un hétérozygote. Cependant, la méthode Call F entraîne également une plus grande quantité de données manquantes car de nombreux homozygotes pour l’allèle majeur ne seront pas appelés en raison d’erreurs de séquençage. Ainsi, dans ce cas, appeler des génotypes sans filtrage semble être la meilleure stratégie que de filtrer les génotypes lorsqu’on essaie d’estimer le CRG.

Figure 1
figure1

Estimations de la fréquence des allèles aux sites avec un CRG vrai de 5% pour différentes profondeurs de couverture. À chaque profondeur, 1 000 sites ont été simulés en utilisant 200 individus, et à chaque site, une estimation de la fréquence des allèles est calculée en utilisant: (1) les génotypes vrais (True); (2) les génotypes appelés sans filtrage (Appel NF); (3) les génotypes appelés avec filtrage (appel F); et (4) la méthode du maximum de vraisemblance (ML). Pour plus de détails sur les méthodes d’estimation, voir Méthodes.

Les résultats sont radicalement différents pour la nouvelle méthode ML. Cette méthode fournit des estimations impartiales du CRG (médiane d’environ 4,9 %) sur une gamme de profondeurs. Même à 2×, les estimations ne montrent qu’une variance légèrement plus grande que celles basées sur les génotypes réels.

Nous avons également comparé l’erreur quadratique moyenne estimée (MSE; Espérance () des différentes estimations du CRG sur une plage de profondeurs de séquençage (Figure 2). La méthode ML a un MSE inférieur aux méthodes d’appel avec 50 ou 200 individus. En particulier, le MSE calculé sur la base de la méthode Call F est beaucoup plus élevé que ceux des autres méthodes surtout lorsque la profondeur diminue. L’ESM des estimations du CRG basées sur les génotypes réels reflète la limite inférieure de l’ESM et n’est pas constante d’une profondeur à l’autre en raison de la variance d’échantillonnage et d’une taille d’échantillon finie. En utilisant 50 individus, le MSE s’approche de 0,0005 avec une profondeur croissante et lorsqu’on utilise une taille d’échantillon de 200 individus, il s’approche de 0,0013 avec une profondeur croissante.

Figure 2
figure2

Erreur quadratique moyenne (MSE; Attendu ) de quatre erreurs différentes types d’estimateurs de fréquence des allèles pour différentes tailles d’échantillon (panneau gauche et droit) et profondeurs de couverture (axe des abscisses). À chaque profondeur, MSE a été calculé à partir des estimations de fréquence des allèles effectuées en utilisant quatre méthodes différentes: True, Call NF, Call F et ML (pour plus de détails sur les méthodes, voir la légende de la figure 1).

Dans l’ensemble, la nouvelle méthode ML surpasse les méthodes d’appel de génotype.

Estimation d’une distribution de MAF à partir de données simulées

Nous examinons ensuite comment les différentes approches d’estimation ont été utilisées pour estimer la proportion de SNP à différentes fréquences dans la population (similaire au spectre de fréquences du site, mais basé sur la fréquence de l’allèle de la population au lieu de la fréquence de l’échantillon). Ici, nous avons simulé 20 000 SNP où la distribution des véritables CRG suivait la distribution stationnaire standard pour une taille de population effective de 10 000 (voir Méthodes). Notez qu’en pratique, cependant, il est très difficile de distinguer un SNP très rare d’une erreur de séquençage. Par conséquent, à des fins de comparaison avec des données réelles, nous avons éliminé les SNP dont le CRG estimé était inférieur à 2%. La figure 3 montre la proportion de SNP tombant dans chaque bac de fréquences différentes après avoir exclu ces SNP avec un CRG estimé < 2%.

Figure 3
figure3

Distribution des fréquences des allèles des SNP simulés en supposant la distribution stationnaire standard des fréquences des allèles. À chaque profondeur (chaque panneau), 20 000 SNP ont été simulés, et pour chaque SNP, des estimations du CRG ont été obtenues en utilisant quatre méthodes différentes (voir la légende de la figure 1). Ensuite, pour chaque méthode (chaque couleur), seuls les sites avec des fréquences d’allèles estimées > 2% sont utilisés pour générer chaque histogramme (axe des abscisses).

Comme prévu, avec une profondeur de couverture élevée, telle que 10 × par individu, toutes les méthodes fournissent des distributions CRG estimées similaires à la distribution attendue basée sur les génotypes réels (figure 3). Avec une profondeur de couverture moins profonde, par exemple moins de 4 × par individu, les distributions des CRG obtenues par les méthodes d’appel de génotypes s’écartent considérablement de la distribution prévue des CRG basée sur les génotypes réels (figure 3). En particulier, ces méthodes surestiment la proportion de SNP basse fréquence. Par exemple, la proportion attendue de SNP dans le deuxième bac (CRG estimé entre 2 et 4 %) est de 18 %. La proportion correspondante basée sur la méthode d’appel NF à une profondeur de 4 × est de 26%, ce qui est 1,4 fois plus élevé que prévu. La surestimation de la proportion de SNP à basse fréquence est due à la confusion des erreurs de séquençage avec les vrais hétérozygotes, ce qui entraîne une surestimation des génotypes hétérozygotes. L’ampleur de cette inflation diffère selon les seuils de filtrage, mais un seuil plus élevé n’augmente ou ne diminue pas nécessairement l’inflation.

L’image est entièrement différente pour la méthode ML. La distribution estimée du CRG obtenue à partir de la nouvelle méthode ML suit de près la distribution réelle même avec de faibles profondeurs de couverture. Ici, il n’y a presque pas d’excès de SNPS basse fréquence. A une profondeur de 4×, la proportion de SNP dans le deuxième bac de l’histogramme est de 18,4%, ce qui est très proche de la proportion attendue (18%). Ainsi, des estimations plus fiables du spectre de fréquences peuvent être faites à partir de données à faible couverture en utilisant notre approche de vraisemblance plutôt qu’en utilisant les approches d’appel de génotype.

Cartographie d’association dans les données simulées

Nous comparons les performances des méthodes qui traitent les génotypes inférés comme des génotypes vrais dans les tests d’association (à l’aide d’un test G) à notre test du rapport de vraisemblance (LRT) qui tient compte de l’incertitude dans les génotypes. Nous examinons la distribution de la statistique de test sous l’hypothèse nulle de l’absence de différence de fréquence d’allèle entre les cas et les témoins. Nous comparons également la puissance des différentes approches.

Avec des tailles d’échantillon raisonnablement grandes, la théorie asymptotique standard suggère que, dans l’hypothèse nulle, la statistique G et la statistique LRT suivent une distribution du chi carré avec un degré de liberté (χ2(1)). Par conséquent, nous avons comparé la distribution nulle de la statistique G calculée sur la base des méthodes d’appel ainsi que la statistique LRT à la distribution χ2(1) en utilisant des graphiques QQ (Figure 4). Nous avons simulé 5 000 SNP à travers une variété de profondeurs de séquençage dans 500 cas et témoins où le CRG utilisé pour simuler les génotypes était de 5% dans les deux cas et les témoins. La distribution de la statistique G calculée en utilisant les génotypes vrais montre une très bonne correspondance avec une distribution χ2(1). Cependant, la distribution de la statistique G calculée à partir des génotypes appelés s’écarte sensiblement d’une distribution χ2(1). Appeler des génotypes puis traiter ces génotypes comme étant précis produit un vaste excès de signaux faussement positifs si les valeurs de p sont calculées à l’aide d’une distribution du χ2(1). Par exemple, à une profondeur de 2×, 11% des SNP avaient une valeur p inférieure à 5%, par rapport aux 5% attendus. L’effet est causé par une augmentation de la variance, due à la sur-concentration des homozygotes en hétérozygotes, dans le test allélique utilisé ici pour détecter l’association. Les tests génotypiques tels que le test de tendance d’Armitage, qui sont robustes aux écarts par rapport à l’équilibre de Hardy-Weinberg, ne montrent pas une augmentation similaire du taux de faux positifs (fichier supplémentaire 2). Conformément à cette observation, le filtrage des génotypes appelés entraîne une diminution de la fraction des tests significatifs lors de l’utilisation du test G, bien que le filtrage ne résout pas complètement le problème. D’autre part, la statistique du TLR ne montre qu’un très léger écart par rapport à une distribution du χ2(1) pour des profondeurs de couverture de 2× ou 5×.

Figure 4
figure4

QQ – tracés comparant la distribution nulle de la statistique de test d’intérêt avec une distribution χ2(1). Chaque colonne correspond à une statistique de test différente: (1) G-statistique calculée en utilisant les génotypes vrais (True); (2) G-statistique calculée en utilisant des génotypes appelés sans filtrage (Appel NF); (3) Statistique G calculée à l’aide de génotypes appelés avec filtrage (Appel F); et (4) la statistique de test du rapport de vraisemblance avec allèle mineur inconnu (LRT). En supposant 500 cas et 500 témoins, dans l’hypothèse nulle, un ensemble de 5 000 sites a été simulé avec un CRG de 5% avec une profondeur de séquençage de 2 × (panneaux supérieurs) et 5 × (panneaux inférieurs). Le facteur « Inflation » est indiqué dans le coin supérieur gauche de chaque figure.

Nous avons également généré des courbes de caractéristiques de fonctionnement du récepteur (ROC) pour chacun des différents tests d’association. Ces courbes montrent la puissance du test à différents taux de faux positifs. Étant donné que les distributions de certaines des statistiques de test ne suivent pas la distribution du χ2(1) sous l’hypothèse nulle, pour faire une comparaison équitable, nous avons obtenu la valeur critique pour chaque taux de faux positifs sur la base de la distribution empirique nulle. La puissance est calculée comme la fraction des loci de maladie simulés dont la statistique dépasse la valeur critique. Dans l’ensemble, nous constatons que le TLR fonctionne mieux que le test G basé sur l’une ou l’autre méthode d’appel du génotype (figure 5). Par example, à un taux de faux positifs de 5% et avec une profondeur de séquençage de 5×, la puissance de détection d’un locus de maladie avec un CRG de 1% et un risque relatif (RR) de 2 est de 51% avec le TLR, mais la puissance chute à 33% en utilisant la méthode d’appel sans filtrage et à 34% en utilisant la méthode d’appel avec filtrage. En particulier, à faible profondeur, le test G appliqué aux génotypes appelés avec filtrage fonctionne très mal (colonne la plus à gauche de la figure 5). Si l’on compare la puissance du LRT au test de tendance d’Armitage en utilisant des génotypes appelés, on constate que le LRT a également une puissance supérieure au test de tendance d’Armitage (fichier supplémentaire 3). Cela suggère que si l’on souhaite utiliser des génotypes appelés, les filtrer en fonction de la confiance des appels peut entraîner une perte de puissance.

Figure 5
figure5

Courbes de caractéristiques de fonctionnement du récepteur (ROC) de quatre tests d’association. Pour la définition des quatre statistiques, voir la légende de la figure 4. En supposant 500 cas et 500 témoins, un ensemble de 20 000 sites ont été simulés sous la valeur nulle et sous l’alternative à des profondeurs de séquençage individuelles de 2 ×, 5× et 10× (trois colonnes). À chaque taux de faux positifs (axe des abscisses), la valeur critique correspondante a été calculée à l’aide de la distribution empirique nulle. Le taux positif réel (puissance; axe des ordonnées) a été obtenu en calculant la fraction des sites responsables avec des statistiques de test qui dépassent la valeur critique.

Application à des données réelles

Nous avons analysé 200 exomes de témoins pour une étude d’association de maladies qui ont été séquencés à l’aide de la technologie Illumina à une profondeur par individu de 8×. Nous avons utilisé les probabilités de génotype générées par le programme « SOAPsnp » pour notre inférence. Pour plus de détails, voir Méthodes.

Tout d’abord, nous avons exploré la précision des estimations du CRG à partir de données de séquençage de nouvelle génération pour 50 SNP en les comparant aux CRG estimées à partir de données de génotype de séquenom. Les estimations utilisant la méthode ML et la méthode d’appel du génotype sans filtrage sont fortement corrélées aux estimations faites à partir des données du génotype du séquenome (c.-à-d. une petite différence normalisée entre les deux estimations de la figure 6). Cependant, les estimations basées sur l’appel de génotype avec filtrage montrent une mauvaise correspondance avec les fréquences estimées à partir des données de génotype de séquenome, en particulier lorsque la profondeur de séquençage est faible. Fait intéressant, il existe un SNP où le CRG estimé à partir des données de reséquençage est très différent de l’estimation obtenue à partir des données du génotype du séquenome, même si la profondeur de séquençage est très élevée (14 ×). Plus précisément, le CRG estimé à partir des données du génotype du séquenome est de 22,5%, mais est de 17,2% lorsqu’il est estimé à l’aide de l’approche ML. L’examen individuel montre que chez de nombreuses personnes, le génotype hautement soutenu basé sur les données de séquençage diffère des génotypes de séquenome. Étant donné que ce SNP est couvert par de nombreuses lectures chez ces individus et que les bases de lecture observées ont des scores de haute qualité (> Q20), il est probable que la différence soit due à des erreurs de génotypage des séquénomes. Notez qu’il existe quelques SNP dans lesquels les MAF estimées de l’approche d’appel du génotype sans filtrage semblent mieux correspondre aux MAF estimées du génotypage des séquénomes que les estimations de l’approche ML. Par exemple, pour un SNP, le CRG estimé est de 25,7 % à partir des données du génotype du séquenome, 25.9% de la méthode d’appel du génotype sans filtrage et 27,2% de la méthode ML. Cependant, l’inspection individuelle révèle qu’il existe quelques individus pour lesquels le génotype appelé à partir des données de séquençage diffère du génotype du séquenome. Dans ces cas, les erreurs des génotypes appelés se sont annulées, donnant l’apparence d’une meilleure correspondance avec les données du génotype du séquenome. Par conséquent, pour ces SNP, il est difficile de dire quelle méthode fonctionne le mieux.

Figure 6
figure6

Estimations de la fréquence des allèles calculées à partir de 200 individus à l’aide de données de séquençage de nouvelle génération par rapport aux données de génotype du séquenome. Pour chaque site, seuls les individus possédant à la fois des données de génotype de séquenome et des données de séquençage ont été utilisés pour estimer la fréquence des allèles. Pour les données de séquençage, les estimations du CRG ont été obtenues en utilisant trois méthodes différentes (Appel NF; Appel F; et ML). La différence normalisée pour chaque estimation a été calculée comme suit : , où et sont les MAF estimées à partir des données de séquençage et des données de génotype de séquenome, respectivement, et n est le nombre d’individus utilisés pour l’estimation. Chaque site est classé dans l’un des quatre bacs en fonction de la profondeur de couverture individuelle moyenne (couleur): moins de 4×, plus de 4× mais moins de 8×, plus de 8× mais moins de 16× et plus de 16×.

Nous avons ensuite examiné la distribution des MAF calculées à l’aide de plusieurs approches sur une gamme de profondeurs de séquençage à partir de nos données de séquençage d’exomes de nouvelle génération (figure 7). Nous avons éliminé les SNP avec un CRG estimé < 2% car il est difficile de distinguer ces SNP à très basse fréquence des erreurs de séquençage dans cet ensemble de données. Nous avons également supprimé les sites dans lesquels il y avait une différence significative (valeur p inférieure à 10-5 en utilisant un test de somme de rang) dans le score de qualité des bases de lecture entre les allèles mineurs et majeurs. Ces sites sont susceptibles d’être des SNP artificiels qui peuvent se produire en raison d’une cartographie incorrecte ou de biais inconnus introduits au cours de la procédure expérimentale. Ensuite, nous avons classé chaque site dans des bacs en fonction de la profondeur de couverture. Le nombre de SNP dans chaque bac est indiqué dans le tableau 1. Lorsque la profondeur moyenne est inférieure à 9 ×, les distributions des MAF estimées basées sur les méthodes d’appel du génotype sont très différentes de celles basées sur la méthode ML. Plus précisément, les approches d’appel du génotype donnent lieu à un grand excès de SNP basse fréquence (CRG compris entre 2% et 4%). Ce schéma reflète ce qui a été observé dans nos études de simulation (figure 3). De plus, pour les méthodes d’appel du génotype, la distribution de fréquence des allèles change radicalement à mesure que la profondeur du séquençage change. Par conséquent, comme discuté précédemment, lorsque la profondeur n’est pas très élevée, les méthodes d’appel de génotypage sont susceptibles d’inclure beaucoup de faux SNP qui sont des erreurs de séquençage. Ces erreurs apparaissent comme un excès de SNP basse fréquence dans la distribution de fréquences. La distribution basée sur la méthode ML est plus stable d’une profondeur à l’autre, mais il existe toujours un excès de SNP à faible fréquence d’allèle avec une profondeur inférieure à 9 × par rapport à la proportion de SNP à basse fréquence à des profondeurs plus importantes.

Figure 7
figure7

Distribution de la fréquence des allèles mineurs estimée à partir des exomes de 200 individus séquencés. Pour chaque site, la fréquence des allèles mineurs a été estimée à l’aide de quatre méthodes différentes: (1) la méthode ML avec un allèle mineur inconnu, (2) la méthode ML avec un allèle mineur connu ou fixe, (3) appeler des génotypes sans filtrage (Appel NF) et (4) appeler des génotypes avec filtrage (Appel F). Chaque site est classé en bacs en fonction de la profondeur de couverture. De plus, dans chaque histogramme, les sites dont le CRG estimé est inférieur à 2 % ne sont pas pris en compte. Pour le nombre de SNP utilisés pour cette analyse, voir le tableau 1.

Tableau 1 Nombre de SNP avec un CRG estimé supérieur à 2% en utilisant une méthode particulière (ligne) dans chaque bac (colonne) définie par la profondeur de séquence moyenne entre les individus.

Enfin, nous avons utilisé ces données de reséquençage d’exome pour simuler une étude d’association cas-témoins. Pour examiner la distribution des statistiques du test d’association dans l’hypothèse nulle, nous avons attribué au hasard 100 individus à un groupe de cas et les 100 autres au groupe témoin. Pour tous les SNP sur le chromosome 2 avec des estimations du CRG > 2% (sur la base de la méthode ML des allèles mineurs inconnus), nous avons testé les différences de fréquence des allèles entre les cas et les témoins en calculant la statistique G en utilisant des génotypes appelés avec et sans filtrage ainsi que la statistique LRT. La figure 8 montre les graphiques QQ comparant les distributions des statistiques de test à la distribution standard χ2(1). Comme on le voit dans les études de simulation, la distribution nulle de la statistique G calculée en appelant des génotypes sans filtrage s’écarte sensiblement de la distribution du χ2(1). Cependant, la distribution nulle de la statistique LRT suit de près la distribution χ2(1). Le facteur d’inflation est de 1,01, ce qui implique que la statistique du TLR fonctionne bien lorsqu’elle est appliquée aux données réelles.

Figure 8
figure8

QQ – plots comparant les statistiques de test d’association pour les différences de fréquence des allèles entre 100 cas et 100 contrôles à une distribution χ2(1). Les phénotypes ont été assignés au hasard à des individus dans l’ensemble de données de rééquençage de l’exome, de sorte qu’il y ait 100 cas et 100 témoins. Pour chaque site, trois statistiques ont été calculées: la statistique G utilisant des génotypes appelés sans filtrage (Appel NF), la statique G utilisant des génotypes appelés avec filtrage (Appel F) et la statistique LRT. Pour minimiser l’inclusion de faux SNP, les sites dont les estimations du CRG de ML sont inférieures à 2 % sont éliminés. À des fins d’affichage, les résultats des sites sur le chromosome 2 sont affichés. Notez que le facteur d’inflation est affiché dans le coin supérieur gauche de chaque graphique QQ.

Related Posts

Laisser un commentaire

Votre adresse e-mail ne sera pas publiée. Les champs obligatoires sont indiqués avec *