Distribution des étiquettes autour des positions de liaison aux protéines

En général, l’immunoprécipitation sélectionne un ensemble de fragments d’ADN qui se chevauchent autour des positions liées. Le séquençage à haut débit identifie des étiquettes courtes (∼35 pb pour les plates-formes Solexa ou solides) aux extrémités 5 ‘ des fragments de l’un ou l’autre des brins d’ADN. Les positions des étiquettes sont ensuite déterminées en les alignant sur l’ensemble du génome, les alignements ambigus étant généralement écartés. La distribution spatiale résultante des occurrences d’étiquettes autour d’une position de liaison stable montrera donc des pics distincts de densité d’étiquettes sur les brins positifs et négatifs (Fig. 1b, c). La distance entre les pics doit refléter la taille de la région protégée, bien qu’elle puisse également être influencée par la distribution de taille des fragments d’ADN. Cette distance ne présente pas de forte dépendance au nombre de balises dans les pics (Tableau supplémentaire 1 en ligne).

Une signature à l’échelle du génome de ce motif d’étiquettes peut être évaluée en calculant la corrélation croisée des densités d’étiquettes de brins positifs et négatifs, en décalant les brins les uns par rapport aux autres en augmentant la distance. Tous les ensembles de données examinés présentent un pic clair dans le profil de corrélation croisée des brins, correspondant à la taille prédominante de la région protégée (Fig. 1d et Fig. supplémentaires. 1 en ligne). L’amplitude du pic reflète la fraction d’étiquettes dans l’ensemble de données qui apparaît conformément au modèle d’étiquettes de liaison attendu. Dans un cas idéal, lorsque toutes les étiquettes séquencées participent à de tels motifs de liaison, l’amplitude de corrélation atteint une valeur maximale. Inversement, l’amplitude diminue à mesure que les positions des étiquettes sont randomisées (Fig. 2 en ligne).

Utilisation d’alignements d’étiquettes de qualité variable

Bien que certaines étiquettes s’alignent parfaitement avec le génome de référence, d’autres ne s’alignent que partiellement, avec des lacunes ou des non-concordances. Les étiquettes mal alignées peuvent résulter de problèmes expérimentaux tels que la contamination des échantillons, correspondre à des régions polymorphes ou non assemblées du génome ou refléter des erreurs de séquençage. Pour la plate-forme Solexa, les erreurs de séquençage sont plus abondantes vers les extrémités 3′ des fragments séquencés, ce qui entraîne fréquemment des alignements partiels qui n’incluent que les portions des étiquettes proches des extrémités 5′. Nous estimons que cette augmentation des fréquences de discordance vers les terminaisons 3′ représente 41 à 75 % de toutes les discordances observées dans les ensembles de données examinés (fig. 3 en ligne). Comme il n’est pas rare que > 50% du total des balises ne donnent lieu qu’à un alignement partiel, l’inclusion de balises partiellement alignées mais toujours informatives est importante pour optimiser l’utilisation de tout ensemble de données11,12. Nous avons donc choisi d’utiliser la longueur de l’appariement et le nombre de nucléotides couverts par des non-appariements et des lacunes pour classer la qualité de l’alignement des étiquettes (Tableau 1 et Tableau supplémentaire 2 en ligne).

Compte tenu d’une classification des étiquettes par qualité d’alignement, nous proposons d’utiliser le profil de corrélation croisée de brins pour déterminer si une classe particulière d’étiquettes doit être inclus dans une analyse plus approfondie. Un ensemble d’étiquettes informatives sur les positions de liaison devrait augmenter l’amplitude de la corrélation croisée, tandis qu’un ensemble d’étiquettes cartographié aléatoirement devrait la diminuer (Fig. 2). Utiliser cette approche pour l’ensemble de données NRSF (Fig. 2), nous avons constaté que les alignements avec des correspondances couvrant au moins 18 pb et zéro décalage amélioraient le profil de corrélation croisée. Cependant, seules les correspondances complètes (25 pb) doivent être prises en compte pour les étiquettes avec deux incompatibilités. L’utilisation de ce critère pour accepter les étiquettes a augmenté leur nombre sur l’ensemble des étiquettes parfaitement alignées de 27% pour l’ensemble de données NRSF, de 30% pour l’ensemble de données CTCF et de 36% pour l’ensemble de données STAT1 (Fig. supplémentaire. 4 en ligne). L’incorporation de ces étiquettes a amélioré la sensibilité et la précision des positions de fixation identifiées (fig. 5 en ligne).

Contrôle de la distribution des étiquettes d’arrière-plan

La signification statistique du regroupement des étiquettes observé pour une position de liaison présumée aux protéines dépend du motif d’arrière-plan attendu. Le modèle le plus simple suppose que la densité de l’étiquette d’arrière-plan est répartie uniformément le long du génome et indépendamment entre les brins 11. En plus de l’échantillon de puce NRSF, Johnson et al.2 ont séquencé un échantillon d’entrée de contrôle, fournissant une évaluation expérimentale de la distribution des étiquettes de fond. Nous avons constaté que la distribution des étiquettes d’arrière-plan présente un degré de clustering significativement plus élevé que prévu à partir d’un processus de Poisson homogène suggéré par le modèle simple susmentionné (P< 10-6, supplémentaire Fig. 6 en ligne).

Notre examen de la densité des balises d’entrée identifie trois principaux types d’anomalies de fond. Le premier type se traduit par des pics singuliers de densité d’étiquette à une seule position chromosomique de plusieurs ordres de grandeur supérieurs à la densité environnante (Fig. 3 bis). De tels pics se produisent généralement à la même position sur les deux brins chromosomiques. Le deuxième type d’anomalie se traduit par des amas non uniformes et larges (> 1 000 pb) de densité de tag accrue apparaissant sur l’un ou l’autre des brins ou les deux (Fig. 3b). Le troisième type présente de petits groupes de densité d’étiquettes spécifiques aux brins ressemblant au motif attendu d’une position stable de liaison aux protéines, bien qu’il montre généralement une plus petite séparation entre les pics de brins (Fig. 3c). Un ensemble similaire d’anomalies peut être observé dans le séquençage d’entrée d’autres organismes (données non présentées).

Le premier type d’anomalie peut être facilement détecté et éliminé en raison de son écart extrême par rapport à la densité de balise environnante. Cependant, les autres types d’anomalies, en particulier la troisième, sont difficiles à distinguer au sein des données de la puce. Cela indique que le séquençage du matériel d’entrée est essentiel pour tenir compte correctement de la distribution des balises d’arrière-plan. Le séquençage d’une expérience témoin simulée (anticorps non spécifique ou pas d’anticorps) peut également être nécessaire.

Pour contrôler la distribution inégale de l’arrière-plan, les méthodes de liaison proposées ci-dessous soustraient la densité d’étiquettes de fond redimensionnée avant de déterminer les positions de liaison, si de telles données sont disponibles. De plus, seules les positions de liaison dans des régions présentant des rapports PuCE/étiquette d’entrée significatifs sont acceptées2. L’effet de ces corrections d’arrière-plan sera caractérisé dans les sections qui suivent.

Méthodes de détection de liaison et couverture relative des sites de liaison

Nous avons examiné cinq méthodes différentes d’appel de positions de liaison, dont deux algorithmes précédemment publiés (CSP, XSET) et trois méthodes qui nous sont propres. En bref, le procédé de localisateur de pics ChIPSeq (CSP) identifie les régions d’enrichissement significatif par rapport au profil d’entrée et détermine les positions de liaison comme celles avec le plus grand nombre de balises dans ces régions2. La méthode de l’ensemble étendu (XSET) étend les étiquettes de brins positifs et négatifs selon la longueur attendue du fragment d’ADN et détermine les positions de liaison comme celles avec le plus grand nombre de fragments qui se chevauchent.11.

Nos méthodes tirent parti du motif d’étiquette spécifique aux brins observé aux positions de fixation (Fig. 1c). La première méthode de ce type, la densité d’étiquettes de fenêtre (WTD), est similaire à XSET, mais note les positions en fonction du nombre d’étiquettes spécifiques au brin en amont et en aval de la position examinée (Fig. 4 bis). La deuxième méthode, l’appariement des pics de brin (MSP), détermine les pics locaux de densité d’étiquette spécifique au brin et identifie les positions entourées de pics de brin positifs et négatifs d’une magnitude comparable à la distance attendue (Fig. 4b). La troisième méthode, la corrélation d’étiquettes en miroir (MTC), analyse le génome pour identifier des positions présentant des motifs d’étiquettes de brins positifs et négatifs prononcés qui se reflètent les uns les autres (Fig. 4c). Le code source est disponible en ligne (Code Source supplémentaire) et un package R à jour peut être téléchargé à l’adresse http://compbio.med.harvard.edu/Supplements/ChIP-seq.

Bien qu’une liste complète des véritables sites de liaison ne soit connue pour aucun des ensembles de données examinés, les trois protéines présentent des spécificités de séquence de liaison connues. Bien que les méthodes de détection de liaison décrites dans ce travail ne reposent pas sur des informations de séquence, nous avons utilisé des instances de motif de séquence à score élevé pour évaluer les performances relatives de différentes méthodes de détection de liaison. Ce faisant, nous supposons uniquement que les instances de motif à score élevé contiennent un sous-ensemble représentatif de positions de reliure réelles et n’exigent pas que tous les motifs à score élevé soient liés, ou que tous les sites de reliure réels présentent une signature de motif. Nous avons évalué la performance en utilisant des motifs de séquence canonique pour la liaison par NRSF et CTCF14,15, et le motif de site activé par le gamma (GAZ) comme prédicteur de la liaison de STAT15,11. Les méthodes de détection de liaison fournissent des scores de magnitude de crête associés aux positions de liaison identifiées, permettant ainsi la hiérarchisation des positions de liaison déterminées par chaque méthode.

Pour comparer la sensibilité des différentes méthodes, nous avons sélectionné un nombre croissant de positions de liaison supérieures renvoyées par chaque méthode et examiné la fraction des occurrences de motifs pour lesquelles une position de liaison a été identifiée (Fig. 4d). Nous avons constaté que 89% des correspondances de motifs NRSF sélectionnées ayant obtenu le meilleur score coïncidaient avec les positions de liaison détectées. Le taux de couverture du motif dépasse nettement celui attendu de la prédiction aléatoire, permettant de comparer les performances relatives des différents procédés de détection de liaison. À l’exception du MSP et du CSP, toutes les méthodes permettent d’obtenir une couverture de motif également élevée. La méthode CSP fonctionne moins bien pour les positions de liaison les plus proéminentes (top 500), tandis que l’approche MSP fonctionne mal sur toute la plage. Les analyses de liaison STAT1 et CTCF montrent des résultats analogues en termes de performances relatives des différentes méthodes (Fig. 7 en ligne). Ces résultats sont également confirmés par l’analyse de loci de liaison validés par PCR à partir de la littérature2, 11, 15 (Figures supplémentaires. 8 et 9 en ligne). Nous notons que les ensembles de test validés par motif et PCR ne représentent qu’une fraction des véritables sites de liaison. Comme cette fraction est plus petite pour CTCF et STAT1, des ensembles plus importants de positions de liaison supérieures sont utilisés pour illustrer la couverture des ensembles de tests par différentes méthodes.

Les méthodes de soustraction d’arrière-plan décrites dans la section précédente améliorent la couverture du motif NRSF, atteignant le même niveau de couverture jusqu’à 11 % moins de positions de fixation supérieures (fig. 10 en ligne). Les corrections ont peu d’effet sur les 1 500 positions de liaison supérieures, qui sont associées à un nombre d’étiquettes plus élevé que tous les pics faussement positifs résultant d’un arrière-plan irrégulier. Les positions faussement positives basées sur le contexte sont généralement de plus petite ampleur et commencent à influencer les prédictions à mesure que des positions plus contraignantes sont prises en compte.

Précision des positions de liaison

Pour évaluer la précision spatiale avec laquelle les positions de liaison aux protéines sont identifiées par différentes méthodes, nous avons analysé les distances entre les positions prédites et les emplacements des frappes de motifs à score élevé (Fig. 5 bis). Pour l’ensemble de données NRSF, la méthode WTD prédit les positions de liaison avec la plus grande précision, avec > 60% des pics prédits situés à moins de 10 pb du centre du motif (Fig. 5b et Fig. supplémentaires. 11a en ligne). Il est suivi par les méthodes XSET, MTC et MSP, avec CSP appelant440% des pics à moins de 10 pb des motifs. Les corrections de fond ont un effet limité sur la précision des positions prédites, seule la méthode WTD montrant une amélioration de 3% pour les positions de liaison fortes (données non montrées).

Pour les prédictions CTCF et STAT1, cependant, la méthode MTC atteint une meilleure précision que WTD (Fig. 5c, d et Fig. supplémentaires. 11 ter, c). La différence peut s’expliquer par les propriétés de la distribution de balises immédiatement proche du centre de la région protégée. Contrairement à WTD et XSET, la méthode MTC ne prend pas en compte les balises dans la région centrale (30 pb) lors de la notation des positions de liaison. La modification de la méthode MTC pour prendre en compte de telles positions réduit la précision des positions de liaison déterminées à un niveau similaire aux prédictions WTD. En examinant la distribution globale des positions des étiquettes par rapport aux coups de motif à score élevé, nous avons constaté que CTCF et STAT1 présentaient des pics inattendus de densité d’étiquettes immédiatement adjacents (entre 10 et 15 pb) à la position du motif (Fig. supplémentaire. 12 en ligne). Ce motif, dans lequel de petits ensembles d’étiquettes de brin négatives apparaissent immédiatement en amont de la région protégée et sont reflétées par les étiquettes de brin positives immédiatement en aval, peut résulter d’interactions de réticulation se produisant au-delà de la région protégée centrale (Fig. 1b, ligne brisée). De ce fait, les méthodes de détection des pics qui prennent en compte les balises près de la région centrale ont tendance à appeler des positions de 15 à 20 pb en amont ou en aval du site de liaison véritable.

Positions statistiquement significatives

Les méthodes de détection de liaison devraient limiter les positions de liaison résultantes à celles qui ne sont pas susceptibles d’avoir eu lieu par hasard. Le niveau de signification statistique souhaité est généralement donné en termes de taux de fausses découvertes (FDR) ou de nombre de positions faussement positives attendues (valeur E).

Les méthodes de détection peuvent ensuite utiliser la distribution de balises d’arrière-plan pour déterminer le score de position de liaison minimal satisfaisant le niveau de signification spécifié. De nombreux appels faussement positifs proviennent des grandes régions anormales décrites précédemment. Ces erreurs systématiques peuvent être filtrées avant la détermination des seuils de signification. Sur la base des données d’échantillon d’entrée pour le NRSF, nous avons trouvé un total de 2 755 positions de liaison pour le seuil FDR de 0,01 en utilisant la méthode WTD. Cela correspond étroitement au nombre de pics supérieurs requis pour obtenir une couverture maximale des positions de motifs à haut pointage utilisées dans les sections précédentes (Fig. 4d).

En l’absence d’une estimation empirique de la distribution des étiquettes de fond, il peut être possible de s’appuyer sur un modèle analytique. Le modèle le plus simple de ce type est un processus de Poisson spatial où les étiquettes sont réparties uniformément dans les régions accessibles du génome11. Cependant, étant donné que les distributions d’étiquettes de fond réelles présentent un degré important de regroupement d’étiquettes, ce seuil basé sur Poisson est nettement inférieur à celui obtenu à partir de la mesure empirique de l’arrière-plan, ce qui entraîne une surestimation du nombre de positions de liaison significatives (9 206 contre 2 755 pour un FDR de 0,01). La comparaison avec les calculs du FDR basés sur les entrées révèle que le modèle basé sur Poisson sous-estime les FDR entre 8 et 20 fois, selon le FDR cible (Tableau supplémentaire 3 en ligne).

Une estimation plus précise des seuils statistiques peut être obtenue en tenant compte du degré de regroupement présent dans la distribution des étiquettes d’arrière-plan. Une approche simple consiste à utiliser une randomisation qui maintient ensemble les balises se produisant à la même position ou à proximité, au lieu de leur attribuer des positions indépendantes, comme cela est fait en utilisant le modèle de Poisson. Le nombre de positions significatives déterminées à l’aide de tels modèles de randomisation avec différentes tailles de bac est indiqué dans le tableau supplémentaire 3. Pour le FDR de 0.01, un modèle de randomisation qui maintient ensemble des étiquettes se produisant exactement à la même position dans le génome donne un nombre comparable de positions de liaison au NRSF (2 985). Nous avons utilisé cette randomisation pour déterminer le nombre de positions de liaison statistiquement significatives pour les ensembles de données CTCF (2 3981 positions pour un FDR de 0,01) et STAT1 (44 921 positions pour un FDR de 0,01). L’appariement du nombre de positions de liaison pour des valeurs de FDR plus strictes nécessite des blocs de randomisation d’étiquettes plus grands (tableau supplémentaire 3), ce qui indique que les stratégies de randomisation simples ne peuvent pas tenir compte correctement des propriétés de regroupement d’arrière-plan.

Test de profondeur de séquençage suffisante

Pour évaluer si la profondeur de séquençage a atteint un point de saturation au-delà duquel aucun site de liaison supplémentaire n’est détecté, nous avons analysé comment l’ensemble des sites de liaison prédits a changé lorsque seul un sous-ensemble de données de balise a été utilisé pour la prédiction. En échantillonnant des fractions croissantes des données de l’étiquette, nous avons déterminé les positions de liaison et comparé ces prédictions avec l’ensemble des sites de liaison de référence identifiés à partir des données complètes (Fig. 6a et Fig. supplémentaires. 13 en ligne).

Si la profondeur de séquençage a dépassé le point de saturation, il serait possible d’arriver à l’ensemble de référence en utilisant uniquement un sous-ensemble des données de balise. Cependant, nous avons constaté qu’aucun des trois ensembles de données n’atteignait un tel point de saturation (asymptote horizontale) et que la fraction des positions de liaison concordantes diminuait lorsque même une petite fraction des données d’étiquette était omise. Cela indique que des sites de liaison supplémentaires sont continuellement identifiés avec une profondeur de séquençage croissante. La tendance observée se maintient pour une gamme de seuils FDR (fig. 13): bien que la pente de la courbe de saturation puisse être réduite en fixant un seuil de FDR considérablement plus strict, il en résulte un nombre nettement plus réduit de sites de liaison.

Pour comprendre les propriétés de la couverture du site de liaison, nous avons examiné le nombre de balises associées à des motifs de séquence à score élevé (Fig. 6b et Fig. supplémentaires. 14 en ligne). Dans les trois ensembles de données, la distribution des nombres d’étiquettes a montré une plage dynamique très large. Alors que certaines positions comportaient des centaines de balises, d’autres dépassaient à peine les chiffres de fond attendus. De plus, ces distributions semblaient continues en ce sens qu’elles ne montraient pas de sous-populations distinctes de positions de liaison. Cela suggère que l’augmentation de la profondeur de séquençage peut permettre de distinguer un plus grand nombre de positions de liaison faibles sans seuil qualitatif qui définirait un ensemble complet de sites de liaison.

Comme des positions de liaison plus prononcées sont identifiées en utilisant une profondeur de séquençage plus petite, une expérience d’une profondeur donnée peut saturer la détection des positions de liaison qui dépassent un certain taux d’enrichissement de l’étiquette par rapport au fond. Nous appelons ce rapport d’enrichissement le rapport d’enrichissement saturé minimal (MSER). Les critères de saturation qui définissent la pente maximale acceptable de la courbe de saturation (Fig. 6a) peut être formulé comme une exigence de stabilité de l’ensemble des sites de liaison prédits. Par exemple, nous exigeons un accord de 99% dans l’ensemble des positions de liaison lorsque l’ensemble de données est réduit de 105 balises. En utilisant les données d’étiquette d’entrée NRSF pour déterminer les intervalles de confiance pour le rapport d’enrichissement de chaque position de liaison, nous avons constaté que la profondeur de séquençage obtenue était suffisante pour saturer la détection de positions de liaison avec des rapports d’enrichissement d’étiquette significativement supérieurs à 7,5 (valeur P < 0,05; Fig. 6a et Fig. supplémentaires. 15 en ligne). Des 2 755 positions de liaison NRSF détectées à un FDR de 0.01, 1,879 (68%) avait des taux d’enrichissement nettement supérieurs à la valeur MSER de 7,5 (Fig. supplémentaire. 13). Nous notons qu’une valeur de MSER particulière n’implique pas que toutes les positions de liaison réelles de cet enrichissement par pli ont été découvertes; elle indique plutôt que de nouvelles positions de liaison avec un enrichissement significativement supérieur à la valeur de MSER sont détectées à un rythme suffisamment lent. Une plage potentielle de rapports d’enrichissement réels peut être évaluée à partir des intervalles de confiance d’enrichissement calculés pour chaque position de liaison (fig. 16 en ligne). Comme l’estimation des intervalles de confiance du rapport d’enrichissement dépend également de la quantité d’informations disponibles sur la distribution des étiquettes de fond, des ensembles de données d’entrée de couverture génomique similaire doivent être utilisés lors de la comparaison de différentes valeurs MSER.

Pour des raisons pratiques, il est important de pouvoir prédire le nombre de balises nécessaires pour saturer la détection de pics au-dessus d’un taux d’enrichissement cible donné. La relation entre le nombre de balises et le MSER s’installe dans une dépendance qui peut être extrapolée à l’aide d’un modèle log-log (Fig. 6c). Nous prédisons, par exemple, que 1.2 × 106 balises supplémentaires seraient nécessaires pour atteindre la saturation dans la détection des positions de liaison NRSF avec un enrichissement sur le fond significativement plus élevé que le double (valeur P < 0,05). Les valeurs et extrapolations MSER dépendent des critères de saturation et des méthodes utilisées pour calculer les intervalles de confiance d’enrichissement (Fig. 17 en ligne).

L’augmentation de la profondeur de séquençage est également susceptible d’entraîner une précision accrue des positions de liaison déterminées. À l’aide de l’ensemble de données NRSF, nous avons analysé comment la distance moyenne entre les positions de liaison détectées et les motifs de séquence dépend du nombre de balises utilisées pour les prédictions. Nos résultats montrent que la précision s’est en effet améliorée avec l’augmentation du nombre de balises (fig. 18 en ligne). L’amélioration, cependant, était mineure: la précision ne diminuait que de plusieurs paires de bases même lorsque le nombre de balises était réduit de moitié.