ChIP-Chip versus ChIP-seq: Lektionen für Versuchsplanung und Datenanalyse

Analyse von Input-DNA-Profilen

Um die technologischen Unterschiede zwischen Sequenzierungs- und Microarray-basierten ChIP-Daten zu verstehen, analysierten wir zunächst die Profile von vernetzten und beschallten DNA-Fragmenten (Input-DNA), die durch Microarray (INPUT-chip) und Hochdurchsatzsequenzierung (INPUT-seq) erzeugt wurden. Da das Input-DNA-Profil unabhängig von dem für ChIP verwendeten Antikörper sein sollte, kann dieser Vergleich Aufschluss über die spezifischen Unterschiede zwischen diesen beiden Profilierungstechnologien geben. Wir haben INPUT-Chip-Daten aus dem Hintergrundkanal unserer Zweikanal-Mikroarray-Daten gewonnen. Während diese Microarray-Plattform eine kompetitive Hybridisierung verwendet, haben sich die beiden Kanäle in unserem Agilent-Microarray als relativ unabhängig erwiesen, da eine Sättigung in beiden Kanälen sehr selten ist. Von allen INPUT-Chip-Profilen, die wir extrahiert haben, stellen wir hier nur die Analyse von acht repräsentativen Profilen vor (zwei aus jedem der vier Entwicklungszeitpunkte), da die meisten INPUT-Chip-Profile sehr ähnlich sind (Zusatzdatei 2: Abbildung S1). Die acht INPUT-Chip-Profile wurden dann mit den neun INPUT-seq-Profilen verglichen, die in dieser Studie gesammelt wurden (Zusätzliche Datei 1: Tabelle S3).

Eine der auffälligsten Beobachtungen ist, dass INPUT-Chip- und INPUT-seq-Profile sich deutlich zu unterscheiden scheinen, obwohl dasselbe DNA-Inputmaterial für die Microarray-Hybridisierung und Sequenzierung verwendet wurde (Abbildung 1). Die relative Größe und Lage der Peaks scheinen bei den INPUT-Chip-Profilen aus mehreren Experimenten konsistent zu sein. Allerdings scheinen die Muster in den neun INPUT-seq-Profilen variabler zu sein. Wir können visuell viele Regionen identifizieren, die über mehrere INPUT-seq-Profile hinweg eine inkonsistente Signalanreicherung aufweisen (in Abbildung 1a hervorgehoben). Um diese Beobachtung zu quantifizieren, wurde eine Clustering-Analyse durchgeführt. Wir fanden heraus, dass alle acht INPUT-Chip-Profile eng beieinander lagen (Abbildung 1b). Dieses Ergebnis zeigt, dass die von Microarray und Hochdurchsatzsequenzierung gemessene Hintergrund-DNA-Verteilung unterschiedlich ist. Alle INPUT-Chip- und sieben von neun INPUT-seq-Profilen korrelierten positiv mit dem genomischen GC-Gehalt auf genomweiter Ebene (Abbildung 1b) sowie um die Transkriptionsstartstellen (TSS) und Transkriptionsendstellen (TES) (Abbildung 1c). Die Stärke der Korrelation mit dem GC-Gehalt ist bei den INPUT-Chip-Profilen sehr einheitlich, bei den INPUT-seq-Profilen jedoch sehr unterschiedlich (Abbildung 1b-c und Zusatzdatei 2: Abbildung S2). Insbesondere korrelieren die INPUT-seq-Profile, die bei E-16-20 h (E16) und E-20-24 h (E20) erhalten wurden, nicht mit dem GC-Gehalt.

Wir stellen auch fest, dass INPUT-seq mit höherer Sequenzierungstiefe (>4 Millionen gemappte Reads) dazu neigen, sich enger zu gruppieren als solche mit geringerer Sequenzierungstiefe, was darauf hindeutet, dass es eine Beziehung zwischen Sequenzierungstiefe und Variabilität der Input-DNA geben könnte. Um diese Hypothese zu testen, generierten wir 11 zusätzliche INPUT-seq-Profile durch Subsampling von Sequenzierungs-Reads aus der am tiefsten sequenzierten Input-DNA-Probe (AdultMale; AM) mit unterschiedlichem Sampling-Anteil (Abbildung 1d und Additional file 2: Abbildung S3). Wie erwartet, neigen Profile mit höherer Sequenzierungstiefe dazu, sich stärker zu gruppieren, und ihre Korrelation mit der Variation des GC-Gehalts ist konsistenter. Allerdings wird die Korrelation mit dem GC-Gehalt erst bei einer sehr geringen Sequenzierungstiefe (<2 Millionen Reads; Abbildung 1d) deutlich schwächer. Dies deutet darauf hin, dass eine geringe Sequenzierungstiefe nicht der einzige Faktor ist, der die INPUT-seq-Qualität beeinflusst. Darüber hinaus können einige INPUT-seq mit relativ geringer Sequenzierungstiefe (E0 und AF, <4 Millionen Reads) konsistente Input-DNA-Profile liefern. Dies bedeutet, dass die INPUT-seq-Variabilität auch auf andere experimentelle Faktoren zurückgeführt werden kann. Obwohl weitere Studien erforderlich sind, um das gesamte Spektrum der experimentellen Faktoren zu untersuchen, die die Variabilität der Input-DNA-Bibliotheken beeinflussen, könnte sie durch Variationen bei der Probenvorbereitung (z. B. unterschiedliche Chromatinvorbereitung und Sonikation), Variationen von Lauf zu Lauf des Sequenzers, Variationen von Sequenzer zu Sequenzer für dasselbe Modell und eine Vielzahl anderer Variablen in den Experimenten beeinflusst werden. Die hohe Variabilität zwischen den INPUT-seq-Profilen ist in der Tat ein kritischer Punkt, da eine große Variabilität zur Instabilität der Dichteschätzung in einem ChIP-seq-Profil beiträgt, was sich auf die nachfolgende Datenanalyse auswirkt. Wie in den folgenden Abschnitten dieses Papiers gezeigt wird, kann eine INPUT-seq mit ungewöhnlich schwacher Korrelation mit dem GC-Gehalt die Konstruktion von Durchschnittsprofilen an wichtigen genomischen Stellen beeinträchtigen. Es ist daher unbedingt erforderlich, die Input-DNA in ausreichender Tiefe zu sequenzieren und sicherzustellen, dass das erhaltene Profil mit denen aus ähnlichen Experimenten übereinstimmt.

Die genomische Abdeckung ist ein weiterer wichtiger Aspekt bei der Wahl zwischen ChIP-Chip und ChIP-seq. Die genomische Abdeckung von ChIP-Chip ist durch das Microarray-Sondendesign begrenzt, und die Abdeckung von ChIP-seq hängt von der Sequenziertiefe ab. Die von unserem Agilent-Mikroarray erreichte genomische Abdeckung liegt bei etwa 70 %. Anhand der sub-sampled INPUT-seq-Daten zeigen wir, dass INPUT-seq im Allgemeinen eine höhere genomische Abdeckung bei einer Sequenzierungstiefe von nur einer Million Reads bietet. Dieser Trend, der sich aus den zufällig unterabgetasteten Daten ergibt, bestätigt die beobachtete genomische Abdeckung der anderen acht echten INPUT-seq-Datensätze (Abbildung 1e).

Vergleich der Profileigenschaften

Anschließend verglichen wir die Eigenschaften von ChIP-Chip- und ChIP-seq-Profilen. Um die von den beiden Technologien erzeugten Profile zu vergleichen, teilten wir das Genom in nicht überlappende 1-kb-Bins ein und definierten das Anreicherungsniveau in jedem Bin als den Durchschnitt des logarithmischen Verhältnisses des IP-Kanals zum Eingangskanal (siehe Abschnitt „Methoden“ für Einzelheiten). Wir bezeichnen die Signalverteilung eines ChIP-Profils als seine Verteilung der Anreicherungswerte aller Bins. Zunächst wollten wir das durchschnittliche Signal-Rausch-Verhältnis der mit den beiden Technologien erzeugten Profile charakterisieren. Wir verwendeten die (abgeschnittene) Schiefe des Signaldichteprofils nach Entfernung der Signale aus den höchsten und niedrigsten 5 % der Verteilung als Maß für das Signal-Rausch-Verhältnis eines Profils. Die Schiefe ist ein Maß für die Asymmetrie einer Verteilung, und eine positive Schiefe zeigt an, dass der Schwanz auf der rechten Seite länger ist, was auf ein gutes Signal-Rausch-Verhältnis schließen lässt. In fast allen Fällen hat ein ChIP-seq-Profil eine höhere Skewness als das entsprechende ChIP-Chip-Profil der gleichen biologischen Bedingung (Abbildung 2 und Additional file 1: Tabelle S4). Wir stellen fest, dass der Unterschied in der Schiefe vom IP-Faktor abhängt, der auf die unterschiedliche Qualität der Antikörper und die Prävalenz von Histonmodifikationen oder Bindungsereignissen zurückzuführen sein könnte. Dieselbe Schlussfolgerung lässt sich auch dann ziehen, wenn eine andere Bin-Größe verwendet wurde (Zusätzliche Datei 2: Abbildung S4). Unsere Ergebnisse bestätigen die allgemeine Beobachtung, dass ChIP-seq in der Regel ein ausgeprägteres Signalprofil als ChIP-Chip erzeugt.

Als nächstes haben wir die Anreicherungsregionen innerhalb jedes ChIP-Profils charakterisiert. Um einen fairen Vergleich durchzuführen, möchten wir einen Algorithmus verwenden, der Peak-Calling für ChIP-seq- und ChIP-Chip-Daten nach denselben Kriterien durchführt. Derzeit sind viele gängige Peak-Calling-Algorithmen speziell für die Analyse von ChIP-Chip- oder ChIP-seq-Daten konzipiert, aber nicht für beide. Um diese Einschränkung zu überwinden, haben wir Peaks sowohl aus ChIP-Chip- als auch aus ChIP-seq-Profilen mithilfe derselben Genom-Scan-Heuristik identifiziert (siehe Abschnitt Methoden). Unsere Ergebnisse zeigen, dass wir fast immer eine größere Anzahl von Peaks und schmalere Peaks entdecken können, wenn wir Daten aus ChIP-seq im Vergleich zu ChIP-chip verwenden, wenn wir dieselbe biologische Probe analysieren, und diese Schlussfolgerung ist unabhängig von der Strenge der verwendeten Identifizierungskriterien konsistent (Abbildung 2 und Zusatzdatei 2: Abbildung S5). In der Praxis können wir wahrscheinlich eine noch größere Anzahl von schmalen Peaks in ChIP-seq-Daten identifizieren, wenn wir im Rahmen des Peak-Calling-Verfahrens explizit strangspezifische Informationen nutzen (neben der Verschiebung jedes Reads in Richtung seines 5′-Endes um eine konstante Anzahl von Basenpaaren), so dass die aktuelle Analyse eine untere Grenze für die Effektivität von ChIP-seq im Vergleich zu ChIP-Chip darstellt. Insgesamt zeigen unsere Ergebnisse, dass ChIP-seq eine höhere räumliche Auflösung und ein besseres Signal-Rausch-Verhältnis bietet.

Genomweite Signalreproduzierbarkeit innerhalb und zwischen Technologien

Weiterhin haben wir die Reproduzierbarkeit zwischen ChIP-Chip- und/oder ChIP-seq-Profilen auf genomweiter Ebene (1 kb-Bins) geschätzt. Um Verzerrungen aufgrund von Unterschieden in der Genomabdeckung und der Sequenzzuordnung (Abbildung 1e) zu vermeiden, haben wir Genomregionen, die keine Microarray-Sonden enthalten, und Regionen mit ungewöhnlich hoher Variabilität über mehrere INPUT-seq-Profile hinweg ausgeschlossen. Der Pearson-Korrelationskoeffizient r wurde als Maß für die Korrelation verwendet, da er für den Vergleich des Schwanzes von zwei Signalverteilungen empfindlicher ist als der Spearman-Korrelationskoeffizient, was bei der Analyse von ChIP-Anreicherungssignalprofilen besonders wichtig ist. Die Korrelation zwischen ChIP-Chip-Replikatpaaren und zwischen ChIP-Seq-Replikatpaaren ist im Allgemeinen hoch (Median r = 0,85 bzw. 0,82), was darauf hindeutet, dass beide Technologien reproduzierbare Ergebnisse liefern können. Wie erwartet ist die plattformübergreifende Korrelation zwischen ChIP-Chip- und ChIP-Seq-Profilen bescheidener (Median r = 0,41; Zusatzdatei 1: Tabelle S5). Ähnliche Schlussfolgerungen lassen sich auch dann ziehen, wenn wir unterschiedliche Bin-Größen für die Berechnung der Korrelation zwischen den Profilen verwenden (Zusätzliche Datei 2: Abbildung S6). Ein repräsentatives Streudiagramm zum Vergleich der einzelnen Technologiepaare ist in Abbildung 3b-d dargestellt. Wir beobachten auch eine positive Korrelation zwischen der Schiefe und der Reproduzierbarkeit zwischen den Profilen (Zusätzliche Datei 2: Abbildung S7), was darauf hindeutet, dass empfindlichere Antikörper konsistentere Profile zwischen den beiden Technologien erzeugen können.

Konstruktion des durchschnittlichen Signalprofils an TSS und TES

Die Konstruktion von durchschnittlichen ChIP-Signalprofilen um wichtige genomische Merkmale wie TSS und TES ist eine gängige Methode zur Visualisierung der Signalanreicherung um diese Merkmale. Daher untersuchten wir die Reproduzierbarkeit der durchschnittlichen TSS- und TES-Profile (2 kb aufwärts und 2 kb abwärts) für jedes Paar von ChIP-Replikatprofilen (Zusatzdatei 2: Abbildung S8). Die durchschnittlichen Profile der meisten Replikatpaare sind sehr konsistent. Es gibt jedoch einige wenige Paare, die sich signifikant unterscheiden, insbesondere die Profile von H3K27Me3 und H3K9Me3 in den Stadien E-16-20 h und E-20-24 h (Zusätzliche Datei 2: Abbildungen S8c und S8g). Ohne externe Validierung ist es unmöglich festzustellen, ob die durch ChIP-Chip oder ChIP-seq erzeugten durchschnittlichen Signalprofile genauer sind. Zwei Indizien lassen jedoch vermuten, dass die durchschnittlichen Signalprofile von ChIP-Chip mit größerer Wahrscheinlichkeit genau sind. Erstens wiesen alle drei ChIP-Chip-Replikate zu diesen Zeitpunkten sehr konsistente Durchschnittsprofile auf. Zweitens ähnelten die durchschnittlichen ChIP-seq-Signalprofile unter diesen biologischen Bedingungen dem Trend der Variation des GC-Gehalts am TSS und TES (Abbildung 1c). Die ungewöhnlich niedrigen Korrelationen zwischen den GC-Gehalten und den INPUT-seq-Profilen von E-16-20 h und E-20-24 h (Abbildung 1b und Zusatzdatei 2: Abbildung S2b) veranlassten uns zu der Hypothese, dass die beobachtete Diskrepanz auf eine falsche Darstellung der GC-Gehaltsvariation durch die jeweiligen INPUT-seq-Profile zurückzuführen ist. Sowohl H3K27Me3 als auch H3K9Me3 sind repressive Markierungen, die in der Regel an TSSs und TESs verarmt sind, so dass jede Variation bei der Hintergrundsubtraktion wahrscheinlich viel ausgeprägter ist als andere Histonmarkierungen, die eine starke Signalanreicherung an diesen genomischen Merkmalen aufweisen. Um unsere Hypothese zu testen, haben wir den entsprechenden INPUT-seq-Hintergrund durch den INPUT-seq aus der AdultFemale-Probe ersetzt, da dieser die höchste Korrelation mit der Variation des GC-Gehalts aufweist. Nach der Ersetzung stimmen die durchschnittlichen Signalprofile, die durch ChIP-seq und ChIP-Chip in diesen beiden Entwicklungsstadien erzeugt wurden, überein (Abbildung 4 und Zusatzdatei 2: Abbildung S9). Dieses Ergebnis ist bemerkenswert, da es zeigt, dass die Verwendung unterschiedlicher INPUT-seq als Negativkontrolle desselben ChIP-seq-Profils zu einer erheblich unterschiedlichen Interpretation der Daten führen kann.

Auswirkungen der Verwendung unterschiedlicher Input-Profile bei der Normalisierung von ChIP-seq-Daten

Nachdem wir die Auswirkungen von INPUT-seq bei der Erstellung von durchschnittlichen TSS- und TES-Profilen beobachtet hatten, fragten wir, ob die Verwendung unterschiedlicher INPUT-seq-Profile für die Hintergrundnormalisierung die ChIP-seq-Peak-Calling-Ergebnisse signifikant beeinflusst. Wir verwendeten SPP, um Peaks für 10 unserer ChIP-seq-Proben (CBP, H3K9Ac, H3K9Me3, H3K27Ac, H3K27Me3 bei E16-20 h und E20-24 h) zu rufen, wobei jedes ChIP-Profil gegen vier verschiedene INPUT-seq als Hintergrund normalisiert wurde (der Input vom passenden Zeitpunkt, AdultFemale, AdultMale und E-4-8 h). Diese INPUT-seq-Profile wurden ausgewählt, weil sie eine unterschiedliche Sequenzierungstiefe und Korrelation mit dem GC-Gehalt aufweisen (Abbildung 1b). Ein Vergleich der Anzahl der Peaks und der medianen Peakbreite ist in Abbildung 5 dargestellt. Wir beobachteten einen großen Unterschied in der Anzahl der Peaks, die für jede ChIP-seq-Probe aufgerufen wurden, wenn verschiedene INPUT-seq als Hintergrund verwendet wurden. Im Extremfall (E-16-24 h, H3K9Me3 ChIP) kann sich die Anzahl der Peaks bei einer FDR von 5 % von Null auf fast 40.000 ändern (Abbildung 5a). Im Allgemeinen wurden mehr statistisch signifikante Peaks (FDR < 0,05) entdeckt, wenn gegen eine tief sequenzierte Eingangs-DNA-Probe normalisiert wurde (AdultMale und E-4-8 h in diesem Experiment), obwohl die absolute Größe des Unterschieds zwischen den ChIP-Datensätzen variiert. Der Unterschied in der Anzahl der Peaks deutet wahrscheinlich auf einen Unterschied in der Nachweisleistung hin. Für jede ChIP-Probe berechneten wir den Anteil der Überschneidungen zwischen jedem Paar von Peak-Sets, die durch vier verschiedene Eingangs-DNA-Hintergründe erzeugt wurden (d. h. sechs Vergleiche pro ChIP-Probe). Wir fanden heraus, dass der mittlere Anteil der Überlappung in Bezug auf den kleineren Peaksatz etwa 95 % beträgt, was darauf hindeutet, dass die Unterschiede in der Anzahl der erkannten Peaks wahrscheinlich auf die unterschiedliche Leistung bei der Erkennung schwächerer Peaks zurückzuführen sind. Wir beobachteten, dass die starken Peaks (d. h. diejenigen mit niedriger FDR) mit größerer Wahrscheinlichkeit in verschiedenen Peak-Sets erkannt wurden (siehe Zusatzdatei 2: Abbildung S10 für ein Beispiel). Die mediane Breite der erkannten Peaks wird auch durch die Verwendung verschiedener INPUT-seq als Hintergrund beeinflusst (Abbildung 5b). Diese Analyse zeigte, dass die Normalisierung mit unterschiedlichen INPUT-seq einen signifikanten und unterschätzten Einfluss auf das Peak-Calling haben kann.

Die Bewertung von Abweichungen aufgrund der Verwendung verschiedener Peak-Caller

Eine weitere wichtige Quelle für Abweichungen bei der Analyse von ChIP-Chip- und ChIP-seq-Profilen ist auf die Verwendung unterschiedlicher Analysealgorithmen zurückzuführen. Bis heute wurden zahlreiche öffentlich zugängliche ChIP-Chip- und ChIP-seq-Analysetools entwickelt, die alle unterschiedliche Methoden für Tag-Shifting, Profilnormalisierung, Glättung, Peak-Identifizierung und Berechnung der Falschentdeckungsrate verwenden. Es ist daher nicht allzu überraschend, dass verschiedene Peak-Caller recht unterschiedliche Ergebnisse bei der Identifizierung von Bindungsstellen liefern können, insbesondere wenn es um Peaks mit schwachen Signalen geht. Anhand unseres Kompendiums von ChIP-Chip- und ChIP-seq-Datensätzen konnten wir beurteilen, wie viel Abweichung bei der Peak-Identifizierung auf den Einsatz unterschiedlicher Profilierungstechnologien und die Verwendung unterschiedlicher Peak-Caller zurückzuführen ist. In dieser Studie haben wir unsere ChIP-Chip-Profile mit zwei Peak-Callern analysiert: MA2C und Splitter, und analysierten unsere ChIP-seq-Profile mit zwei weiteren Peak-Callern: MACS und SPP (siehe Zusatzdatei 1: Tabelle S8). Diese Peak-Caller wurden ausgewählt, weil sie weit verbreitet und öffentlich zugänglich sind und in früheren Vergleichsstudien im Allgemeinen eine gute Leistung gezeigt haben. Wir berechneten die Überlappung der obersten 1.000 Peaks von vier der Faktoren (CBP, H3K4Me1, H3K4Me3 und H3K27Me3) über mehrere Entwicklungsstadien hinweg. Die vier IP-Faktoren wurden ausgewählt, da sie repräsentative Profile mit breiten Peaks (CBP und H3K27Me3) und schmalen Peaks (H3K4Me1 und H3K4Me3) darstellen. Wir stellen hier nur die Ergebnisse des Vergleichs der obersten 1.000 Peaks vor, da dies eine biologisch sinnvolle Anzahl von Anreicherungsstellen mit hoher Vertrauenswürdigkeit in diesen Profilen ist. Die allgemeine Schlussfolgerung dieser Analyse ist robust gegenüber einer Vielzahl von Peak-Calling-Schwellenwerten (Zusatzdatei 2: Abbildung S11). Die Übereinstimmung zwischen zwei Peak-Sets wurde anhand des durchschnittlichen Anteils der sich überschneidenden Peaks gemessen. Wie in Abbildung 6 dargestellt, ergaben die Vergleiche auf der Grundlage von H3K4Me1- und H3K4Me3-Profilen die erwarteten Ergebnisse, bei denen die plattforminterne Konkordanz höher ist als die plattformübergreifende Konkordanz (d. h., Peak-Sets, die von zwei Peak-Callern auf demselben Profil erzeugt wurden, sind konkordanter als Peak-Sets, die von zwei Peak-Callern auf zwei Profilen erzeugt wurden). Bei der Analyse der Profile von H3K27Me3 und CBP kann die Intra-Plattform-Konkordanz jedoch genauso niedrig sein wie die Inter-Plattform-Konkordanz, was bedeutet, dass die Unterschiede bei den Peak-Calling-Algorithmen genauso groß sein können wie die Verwendung unterschiedlicher Profilierungstechnologien für einige IP-Faktoren. Die Beobachtung, dass aktuelle Peak-Calling-Algorithmen für ChIP-Profile mit breiten Domänen (CBP und H3K27Me3) weniger übereinstimmende Ergebnisse liefern als für solche mit scharfen Spitzen (H3K4Me1 und H3K4Me3), könnte darauf hindeuten, dass sie bei der Identifizierung breiter Anreicherungsregionen weniger konsistent sind, was ein interessantes Thema für weitere Untersuchungen sein könnte.