ChIP-chip versus ChIP-seq : Leçons pour la conception expérimentale et l’analyse des données

Analyse des profils d’ADN d’entrée

Pour comprendre les différences technologiques entre les données ChIP basées sur le séquençage et celles basées sur les microréseaux, nous avons d’abord analysé les profils des fragments d’ADN réticulés et soniqués (ADN d’entrée) générés par les microréseaux (INPUT-chip) et le séquençage à haut débit (INPUT-seq). Puisque le profil de l’ADN d’entrée devrait être indépendant de l’anticorps utilisé pour le ChIP, cette comparaison peut donner un aperçu des différences spécifiques entre ces deux technologies de profilage. Nous avons obtenu des données INPUT-chip à partir du canal de fond de nos données de microarray à deux canaux. Alors que cette plateforme de biopuces utilise l’hybridation compétitive, les deux canaux de notre biopuce Agilent se sont avérés relativement indépendants, la saturation de l’un ou l’autre canal étant très rare. Parmi tous les profils INPUT-chip que nous avons extraits, nous ne présentons ici que l’analyse de huit profils représentatifs (deux de chacun des quatre points de temps de développement), car la plupart des profils INPUT-chip sont très similaires (fichier supplémentaire 2 : Figure S1). Les huit profils INPUT-chip ont ensuite été comparés aux neuf profils INPUT-seq recueillis dans cette étude (Additional file 1 : Table S3).

L’une des observations les plus frappantes est que les profils INPUT-chip et INPUT-seq semblent être sensiblement différents, même si le même matériel d’ADN d’entrée a été utilisé pour l’hybridation des puces et le séquençage (Figure 1). L’ampleur relative et l’emplacement des pics semblent être cohérents dans les profils INPUT-chip provenant de plusieurs expériences. En revanche, les schémas des neuf profils INPUT-seq semblent être plus variables. Nous pouvons identifier visuellement de nombreuses régions qui présentent un enrichissement du signal incohérent dans plusieurs profils INPUT-seq (mis en évidence sur la figure 1a). Une analyse de regroupement a été effectuée pour quantifier cette observation. Nous avons constaté que les huit profils INPUT-chip se regroupent étroitement les uns avec les autres (figure 1b). Ce résultat montre que la distribution de l’ADN de fond mesurée à partir de la biopuce et du séquençage à haut débit est différente. Tous les profils INPUT-chip et sept des neuf profils INPUT-seq ont été corrélés positivement avec la teneur génomique en GC au niveau du génome entier (figure 1b), ainsi qu’autour des sites de début de transcription (TSS) et des sites de fin de transcription (TES) (figure 1c). La force de la corrélation avec le GC est très cohérente entre les profils INPUT-chip, mais très variable entre les profils INPUT-seq (figure 1b-c et fichier supplémentaire 2 : figure S2). Notamment, les profils INPUT-seq obtenus à E-16-20 h (E16) et E-20-24 h (E20) ne présentent pas de corrélation avec la teneur en GC.

Nous remarquons également que les INPUT-seq avec une profondeur de séquençage plus élevée (>4 millions de lectures mappées) ont tendance à se regrouper plus étroitement que ceux avec une profondeur de séquençage plus faible, ce qui suggère qu’il pourrait y avoir une relation entre la profondeur de séquençage et la variabilité de l’ADN d’entrée. Pour tester cette hypothèse, nous avons généré 11 profils INPUT-seq supplémentaires en sous-échantillonnant les lectures de séquençage de l’échantillon d’ADN d’entrée le plus profondément séquencé (AdultMale ; AM) à différentes proportions d’échantillonnage (figure 1d et fichier supplémentaire 2 : figure S3). Comme prévu, les profils ayant une profondeur de séquençage plus élevée ont tendance à se regrouper plus fortement, et leur corrélation avec la variation du contenu en GC est plus cohérente. Toutefois, la corrélation avec la teneur en GC ne devient beaucoup plus faible qu’à une profondeur de séquençage très faible (<2 millions de lectures ; figure 1d). Cela indique qu’une faible profondeur de séquençage n’est pas le seul facteur affectant la qualité de INPUT-seq. En outre, certaines séquences INPUT-seq avec une profondeur de séquençage relativement faible (E0 et AF, <4 millions de lectures) peuvent donner des profils d’ADN d’entrée cohérents. Cela implique que la variabilité de INPUT-seq peut également être attribuée à d’autres facteurs expérimentaux. Bien que des études supplémentaires soient nécessaires pour disséquer l’ensemble des facteurs expérimentaux affectant la variabilité des bibliothèques d’ADN d’entrée, celle-ci pourrait être affectée par des variations dans la préparation des échantillons (par exemple, différentes préparations de chromatine et sonication), des variations d’une exécution à l’autre du séquenceur, des variations d’un séquenceur à l’autre pour le même modèle, et une foule d’autres variables dans les expériences. La grande variabilité entre les profils INPUT-seq est en effet un problème critique, car une grande variabilité contribue à l’instabilité de l’estimation de la densité dans un profil ChIP-seq, ce qui affectera l’analyse des données en aval. Comme nous le montrerons dans les sections suivantes de cet article, un INPUT-seq présentant une corrélation anormalement faible avec le contenu GC peut avoir un impact sur la construction de profils moyens à des emplacements génomiques importants. Il est donc impératif de séquencer l’ADN d’entrée à une profondeur suffisante et de s’assurer que le profil obtenu est cohérent avec ceux d’expériences similaires.

La couverture génomique est une autre considération clé lors du choix entre ChIP-chip et ChIP-seq. La couverture génomique de ChIP-chip est limitée par la conception de la sonde du microarray, et la couverture de ChIP-seq dépend de la profondeur de séquençage. La couverture génomique obtenue par notre micropuce Agilent est d’environ 70 %. En utilisant les données INPUT-seq sous-échantillonnées, nous montrons que INPUT-seq fournit généralement une couverture génomique supérieure à une profondeur de séquençage aussi faible qu’un million de lectures. Cette tendance construite à partir des données sous-échantillonnées de manière aléatoire corrobore la couverture génomique observée des huit autres ensembles de données INPUT-seq réels (figure 1e).

Comparaison des caractéristiques des profils

Nous avons ensuite comparé les caractéristiques des profils ChIP-chip et ChIP-seq. Pour comparer les profils générés par les deux technologies, nous avons divisé le génome en bacs non chevauchants de 1 kb et défini le niveau d’enrichissement à chaque bac comme la moyenne du rapport logarithmique du canal IP sur le canal d’entrée (voir la section Méthodes pour plus de détails). Nous faisons référence à la distribution du signal d’un profil ChIP comme étant sa distribution des valeurs d’enrichissement de tous les bacs. Nous avons d’abord cherché à caractériser le rapport signal/bruit moyen des profils générés par les deux technologies. Nous avons utilisé l’asymétrie (tronquée) du profil de densité du signal après avoir retiré les signaux des 5% les plus élevés et les plus bas de la distribution comme mesure du rapport signal/bruit d’un profil. L’asymétrie est une mesure de l’asymétrie d’une distribution et une asymétrie positive indique que la queue du côté droit est plus longue, ce qui implique un bon rapport signal/bruit. Dans presque tous les cas, un profil ChIP-seq présente une asymétrie plus élevée que son profil ChIP-chip correspondant de la même condition biologique (figure 2 et fichier supplémentaire 1 : tableau S4). Nous notons que la différence d’asymétrie dépend du facteur IP, qui pourrait être dû à la qualité différente des anticorps et à la prévalence de la modification des histones ou des événements de liaison. La même conclusion peut être tirée même si une taille de bin différente a été utilisée (fichier supplémentaire 2 : figure S4). Nos résultats ont confirmé l’observation générale selon laquelle ChIP-seq produit généralement un profil de signal plus distinctif que ChIP-chip.

Puis, nous avons caractérisé les régions d’enrichissement dans chaque profil ChIP. Pour effectuer une comparaison équitable, nous aimerions utiliser un algorithme qui effectue le peak calling sur les données ChIP-seq et ChIP-chip en utilisant les mêmes critères. Actuellement, de nombreux algorithmes d’appel de pics couramment utilisés sont spécifiquement conçus pour analyser les données ChIP-chip ou ChIP-seq, mais pas les deux. Pour surmonter cette limitation, nous avons identifié les pics à partir des profils ChIP-chip et ChIP-seq en utilisant la même heuristique de balayage du génome (voir la section Méthodes). Nos résultats indiquent que nous pouvons presque toujours découvrir un plus grand nombre de pics et des pics plus étroits en utilisant les données générées par ChIP-seq par rapport à ChIP-chip lors de l’analyse du même échantillon biologique, et cette conclusion est cohérente quelle que soit la rigueur des critères d’identification utilisés (figure 2 et fichier additionnel 2 : figure S5). Dans la pratique, nous pouvons probablement identifier un nombre encore plus important de pics étroits dans les données ChIP-seq si nous utilisons explicitement les informations spécifiques au brin dans la procédure d’appel de pic (en plus de décaler chaque lecture vers son extrémité 5′ d’un nombre constant de paires de bases), de sorte que l’analyse actuelle fournit une limite inférieure de l’efficacité de ChIP-seq par rapport à ChIP-chip. Dans l’ensemble, nos résultats démontrent que ChIP-seq offre une résolution spatiale et un rapport signal/bruit plus élevés.

Reproductibilité du signal à l’échelle du génome au sein et entre les technologies

En outre, nous avons estimé la reproductibilité entre les profils ChIP-chip et/ou ChIP-seq à l’échelle du génome (bins de 1 kb). Pour éviter les biais dus aux différences de couverture génomique et de cartographie des séquences (figure 1e), nous avons exclu les régions génomiques qui ne contiennent aucune sonde de microarray et les régions présentant une variabilité inhabituellement élevée entre plusieurs profils INPUT-seq. Le coefficient de corrélation de Pearson, r, a été utilisé comme mesure de la corrélation, car il est plus sensible que le coefficient de corrélation de Spearman pour comparer la queue de deux distributions de signaux, ce qui est particulièrement important dans l’analyse des profils de signaux d’enrichissement ChIP. La corrélation entre les paires de réplicats ChIP-chip et entre les paires de réplicats ChIP-seq est généralement élevée (r médian = 0,85 et 0,82, respectivement), ce qui indique que les deux technologies peuvent produire des résultats reproductibles. Comme prévu, la corrélation interplateforme entre les paires de répliques de profils ChIP-chip et ChIP-seq est plus modeste (médiane r = 0,41 ; fichier supplémentaire 1 : tableau S5). Des conclusions similaires peuvent être obtenues même si nous utilisons des tailles de bin différentes pour calculer la corrélation inter-profils (fichier supplémentaire 2 : figure S6). La figure 3b-d présente un diagramme de dispersion représentatif comparant chaque paire de technologies. Nous observons également une corrélation positive entre l’asymétrie et la reproductibilité inter-profils (fichier supplémentaire 2 : figure S7), ce qui suggère que des anticorps plus sensibles peuvent produire des profils plus cohérents entre les deux technologies.

Construction du profil de signal moyen au SCT et au TES

La construction de profils de signal ChIP moyens autour de caractéristiques génomiques importantes telles que le SCT et le TES est un moyen courant de visualiser l’enrichissement du signal autour de ces caractéristiques. Par conséquent, nous avons étudié la reproductibilité des profils moyens de SCT et de TES (2 kb en amont et 2 kb en aval) pour chaque paire de profils ChIP répliqués (fichier supplémentaire 2 : figure S8). Les profils moyens de la plupart des paires répliquées sont très cohérents. Cependant, quelques paires sont significativement différentes, notamment les profils de H3K27Me3 et H3K9Me3 aux deux stades E-16-20 h et E-20-24 h (fichier supplémentaire 2 : figures S8c et S8g). Sans validation externe, il est impossible de déterminer si les profils de signaux moyens générés par ChIP-chip ou ChIP-seq sont plus précis. Néanmoins, deux séries de preuves nous ont amenés à penser que les profils de signaux moyens générés par ChIP-chip étaient plus susceptibles d’être précis. Premièrement, les trois répliques de ChIP-chip à ces points temporels présentaient des profils moyens très cohérents. Deuxièmement, les profils de signaux moyens de ChIP-seq dans ces conditions biologiques ressemblaient à la tendance de la variation de la teneur en GC à TSS et TES (figure 1c). Les corrélations inhabituellement faibles entre les teneurs en GC et les profils INPUT-seq de E-16-20 h et E-20-24 h (figure 1b et fichier complémentaire 2 : figure S2b) nous ont incités à émettre l’hypothèse que l’écart observé était dû à une mauvaise représentation de la variation de la teneur en GC par les profils INPUT-seq respectifs. H3K27Me3 et H3K9Me3 sont des marques répressives qui sont généralement appauvries au niveau des SCT et des TES, de sorte que toute variation dans la soustraction du fond est probablement beaucoup plus prononcée que d’autres marques d’histones qui ont un fort enrichissement du signal au niveau de ces caractéristiques génomiques. Pour vérifier notre hypothèse, nous avons remplacé le fond INPUT-seq correspondant par l’INPUT-seq de l’échantillon AdultFemale, car il présente la corrélation la plus élevée avec la variation du contenu en GC. Après le remplacement, les profils de signaux moyens générés par ChIP-seq et ChIP-chip à ces deux stades de développement concordent (figure 4 et fichier supplémentaire 2 : figure S9). Ce résultat est frappant car il montre que l’utilisation de différents INPUT-seq comme contrôle négatif du même profil ChIP-seq peut conduire à une interprétation substantiellement différente des données.

Effet de l’utilisation de différents profils d’entrée dans la normalisation des données ChIP-seq

Ayant observé l’impact de l’INPUT-seq dans la construction des profils TSS et TES moyens, nous avons demandé si l’utilisation de différents profils INPUT-seq pour la normalisation du fond affecte significativement les résultats d’appel de pics ChIP-seq. Nous avons utilisé SPP pour appeler des pics pour 10 de nos échantillons ChIP-seq (CBP, H3K9Ac, H3K9Me3, H3K27Ac, H3K27Me3 à E16-20 h et E20-24 h) où chaque profil ChIP a été normalisé par rapport à quatre INPUT-seq différents comme fond (l’entrée du point de temps correspondant, AdultFemale, AdultMale, et E-4-8 h). Ces profils INPUT-seq ont été choisis parce qu’ils ont une profondeur de séquençage différente et une corrélation avec la teneur en GC (figure 1b). Une comparaison du nombre de pics et de la largeur médiane des pics est présentée à la figure 5. Nous avons observé une grande différence dans le nombre de pics appelés pour n’importe quel échantillon ChIP-seq lorsque différents profils INPUT-seq étaient utilisés comme arrière-plan. Dans le cas extrême (E-16-24 h, ChIP H3K9Me3), le nombre de pics peut varier de zéro à près de 40 000 pour un FDR de 5 % (figure 5a). En général, davantage de pics statistiquement significatifs (FDR < 0,05) ont été détectés lors de la normalisation par rapport à un échantillon d’ADN d’entrée profondément séquencé (AdultMale et E-4-8 h dans cette expérience), bien que l’ampleur absolue de la différence varie entre les ensembles de données ChIP. La différence dans le nombre de pics indique probablement une différence dans le pouvoir de détection. Pour chaque échantillon ChIP, nous avons calculé la proportion de chevauchement entre chaque paire d’ensembles de pics générés par quatre fonds d’ADN d’entrée différents (c’est-à-dire six comparaisons par échantillon ChIP). Nous avons constaté que la proportion moyenne de chevauchement par rapport à l’ensemble de pics le plus petit est d’environ 95%, ce qui indique que les différences dans le nombre de pics détectés sont probablement dues à une puissance différente pour appeler les pics plus faibles. Nous avons observé que les pics forts (c’est-à-dire ceux dont la FDR de détection est faible) étaient plus souvent détectés dans les différents ensembles de pics (voir le fichier supplémentaire 2 : figure S10 pour un exemple). La largeur médiane des pics détectés est également affectée par l’utilisation de différentes séquences INPUT comme fond (figure 5b). Cette analyse a montré que la normalisation utilisant différents INPUT-seq peut avoir un impact significatif, et sous-estimé, sur l’appel de pics.

Évaluer la variation due à l’utilisation de différents appelants de pics

Une autre source importante de variation dans l’analyse des profils ChIP-chip et ChIP-seq provient de l’utilisation de différents algorithmes d’analyse. Un grand nombre d’outils d’analyse ChIP-chip et ChIP-seq accessibles au public ont été développés à ce jour, et tous utilisent des méthodes différentes pour le décalage des marqueurs, la normalisation des profils, le lissage, l’identification des pics et le calcul du taux de fausse découverte. Il n’est donc pas trop surprenant de constater que différents appelants de pics peuvent générer des résultats très différents en termes d’identification du site de liaison, en particulier lorsqu’il s’agit de pics à faible signal. En utilisant notre recueil de données ChIP-chip et ChIP-seq, nous avons pu évaluer dans quelle mesure la variation dans l’identification des pics peut être attribuée à l’utilisation de différentes technologies de profilage et à l’utilisation de différents appelants de pics. Dans cette étude, nous avons analysé nos profils ChIP-chip en utilisant deux appelants de pics : MA2C et Splitter et analysé nos profils ChIP-seq en utilisant deux autres appelants de pics : MACS et SPP (voir le fichier supplémentaire 1 : tableau S8). Ces appelants de pics ont été choisis parce qu’ils sont largement utilisés, disponibles publiquement et ont généralement montré de bonnes performances dans des études comparatives précédentes. Nous avons calculé le chevauchement des 1 000 premiers pics de quatre des facteurs (CBP, H3K4Me1, H3K4Me3 et H3K27Me3) à travers plusieurs stades de développement. Les quatre facteurs IP ont été choisis car il s’agissait de profils représentatifs contenant des pics larges (CBP et H3K27Me3) et des pics étroits (H3K4Me1 et H3K4Me3). Nous ne présentons ici que les résultats de la comparaison des 1 000 pics les plus élevés, car il s’agit d’un nombre biologiquement raisonnable de sites d’enrichissement de haute confiance dans ces profils. La conclusion générale de cette analyse est robuste face à une variété de seuils d’appel de pics (fichier supplémentaire 2 : figure S11). La concordance entre deux ensembles de pics a été mesurée par la proportion moyenne de pics qui se chevauchent. Comme le montre la figure 6, les comparaisons basées sur les profils de H3K4Me1 et H3K4Me3 ont donné les résultats attendus, à savoir que la concordance intra-plateforme est supérieure à la concordance inter-plateforme (c’est-à-dire que les ensembles de pics générés par deux appelants de pics sur le même profil sont plus concordants que les ensembles de pics générés par deux appelants de pics sur deux profils). Cependant, la concordance intra-plateforme peut être aussi faible que la concordance inter-plateforme lors de l’analyse des profils de H3K27Me3 et de CBP, ce qui implique que la variation des algorithmes d’appel de pic peut être aussi importante que l’utilisation de différentes technologies de profilage pour certains facteurs IP. L’observation selon laquelle les algorithmes actuels d’appel de pic produisent des résultats moins concordants pour les profils de ChIP avec des domaines larges (CBP et H3K27Me3) que ceux avec des pics aigus (H3K4Me1 et H3K4Me3) peut suggérer qu’ils sont moins cohérents pour identifier de larges régions d’enrichissement, ce qui peut être un sujet intéressant pour une étude plus approfondie.

.