ChIP-chip versus ChIP-seq: Lessen voor experimenteel ontwerp en data-analyse

Analyse van input DNA-profielen

Om de technologische verschillen tussen sequencing-gebaseerde en microarray-gebaseerde ChIP-gegevens te begrijpen, analyseerden we eerst de profielen van cross-linked en gesoniceerde DNA-fragmenten (input DNA) gegenereerd door microarray (INPUT-chip) en high-throughput sequencing (INPUT-seq). Aangezien het input DNA profiel onafhankelijk moet zijn van het antilichaam gebruikt voor ChIP, kan deze vergelijking inzicht geven in de specifieke verschillen tussen deze twee profileringstechnologieën. Wij verkregen INPUT-chipgegevens uit het achtergrondkanaal van onze tweekanaals microarraygegevens. Terwijl dit microarray-platform gebruik maakt van competitieve hybridisatie, is aangetoond dat de twee kanalen in onze Agilent-microarray relatief onafhankelijk zijn, aangezien verzadiging in een van beide kanalen zeer zeldzaam is. Van alle INPUT-chip profielen die we geëxtraheerd, presenteren we alleen de analyse van acht representatieve profielen hier (twee van elk van de vier ontwikkelingstijdstippen), omdat de meeste van de INPUT-chip profielen zijn zeer vergelijkbaar (Additional file 2: Figuur S1). De acht INPUT-chip profielen werden vervolgens vergeleken met de negen INPUT-seq profielen verzameld in deze studie (Additional file 1: Tabel S3).

Een van de meest opvallende observaties is dat INPUT-chip en INPUT-seq profielen aanzienlijk lijken te verschillen, ook al werd hetzelfde input-DNA materiaal gebruikt voor microarray hybridisatie en sequencing (figuur 1). De relatieve omvang en locatie van de pieken lijken consistent te zijn over de INPUT-chip profielen van meerdere experimenten. De patronen in de negen INPUT-seq profielen lijken echter meer variabel te zijn. We kunnen visueel identificeren veel regio’s die inconsistent signaal verrijking hebben over meerdere INPUT-seq profielen (gemarkeerd in figuur 1a). Een clustering analyse werd uitgevoerd om deze waarneming te kwantificeren. We vonden dat alle acht INPUT-chip profielen geclusterd dicht bij elkaar (figuur 1b). Dit resultaat toont aan dat de achtergrond DNA distributie gemeten uit microarray en high-throughput sequencing verschillend is. Alle INPUT-chip en zeven van de negen INPUT-seq profielen positief gecorreleerd met genomische GC-gehalte op het genoom-brede niveau (figuur 1b), alsmede rond de transcriptie startplaatsen (TSS) en transcriptie eindplaatsen (TES) (figuur 1c). De sterkte van de correlatie met GC is zeer consistent onder INPUT-chip profielen, maar zeer variabel onder de INPUT-seq profielen (figuur 1b-c en Additional file 2: figuur S2). Met name de INPUT-seq profielen verkregen op E-16-20 h (E16) en E-20-24 h (E20) niet correleren met GC-gehalte.

We merken ook op dat INPUT-seq met een hogere sequencing diepte (>4 miljoen gemapte leest) de neiging om nauwer te clusteren dan die met een lagere sequencing diepte, wat suggereert dat er een relatie kan zijn tussen sequencing diepte en input DNA variabiliteit. Om deze hypothese te testen, genereerden we 11 extra INPUT-seq profielen door subsampling sequencing leest van de meest diep gesequenced input DNA monster (AdultMale; AM) op verschillende bemonsterings proportie (Figuur 1d en Additional file 2: Figuur S3). Zoals verwacht, hebben profielen met een hogere sequencing diepte de neiging om sterker samen te clusteren, en hun correlatie met GC-inhoud variatie is meer consistent. Echter, de GC-inhoud correlatie wordt alleen veel zwakker bij een zeer lage sequencing diepte (<2 miljoen leest; Figuur 1d). Dit geeft aan dat lage sequencing diepte is niet de enige factor die invloed INPUT-seq kwaliteit. Bovendien kunnen sommige INPUT-seq met relatief lage sequencing diepte (E0 en AF, <4 miljoen leest) geven consistente input DNA-profielen. Dit impliceert dat INPUT-seq variabiliteit kan ook worden toegeschreven aan andere experimentele factoren. Hoewel verdere studies nodig zijn om het volledige scala van experimentele factoren die de variabiliteit van de input DNA-bibliotheken te ontleden, kan worden beïnvloed door variaties in de monstervoorbereiding (bijvoorbeeld verschillende chromatine voorbereiding en sonicatie), run-to-run variatie van de sequencer, sequencer-to-sequencer variatie voor hetzelfde model, en een groot aantal andere variabelen in experimenten. De hoge variabiliteit tussen INPUT-seq profielen is inderdaad een kritische kwestie, omdat grote variabiliteit bijdraagt tot instabiliteit van de dichtheid schatting in een ChIP-seq profiel, die downstream data-analyse zal beïnvloeden. Zoals zal worden aangetoond in de volgende delen van dit document, kan een INPUT-seq met ongewoon zwakke correlatie met GC-gehalte invloed hebben op de bouw van gemiddelde profielen op belangrijke genomische locaties. Het is dus noodzakelijk om de sequentie van de input DNA tot voldoende diepte en om na te gaan dat het verkregen profiel consistent is met die van soortgelijke experimenten.

Genomic dekking is een andere belangrijke overweging bij de keuze tussen ChIP-chip en ChIP-seq. De genomische dekking van ChIP-chip wordt beperkt door de microarray probe ontwerp, en de dekking van ChIP-seq is afhankelijk van sequencing diepte. De genomische dekking die door onze Agilent-microarray wordt bereikt, bedraagt ongeveer 70%. Met behulp van de sub-bemonsterde INPUT-seq gegevens, tonen we aan dat INPUT-seq over het algemeen een hogere genomische dekking geeft bij sequencing diepte zo laag als een miljoen leest. Deze trend geconstrueerd uit de willekeurig gesubsamplede gegevens bevestigt de waargenomen genomische dekking van de andere acht echte INPUT-seq datasets (figuur 1e).

Vergelijking van profielkenmerken

We vergeleken vervolgens de kenmerken van ChIP-chip en ChIP-seq profielen. Om de profielen gegenereerd door de twee technologieën te vergelijken, hebben we het genoom verdeeld in 1 kb niet-overlappende bins en gedefinieerd het verrijkingsniveau op elke bin als het gemiddelde van de log ratio van het IP-kanaal over het invoerkanaal (zie de Methoden sectie voor details). We verwijzen naar een signaal verdeling van een ChIP profiel als de verdeling van de verrijking waarden van alle bins. Ten eerste hebben wij getracht de gemiddelde signaal-ruisverhouding voor profielen gegenereerd door de twee technologieën te karakteriseren. Wij gebruikten de (afgeknotte) scheefheid van het signaaldichtheidsprofiel na het verwijderen van signalen uit de hoogste en laagste 5% van de verdeling als een maat voor de signaal-ruisverhouding van een profiel. Scheefheid is een maat voor de asymmetrie van een verdeling en een positieve scheefheid geeft aan dat de staart aan de rechterkant langer is, wat een goede signaal-ruisverhouding impliceert. In bijna alle gevallen heeft een ChIP-seq profiel een hogere skewness dan het overeenkomstige ChIP-chip profiel van dezelfde biologische conditie (figuur 2 en Additional file 1: tabel S4). We merken op dat het verschil van skewness is afhankelijk van de IP-factor die kan te wijten zijn aan verschillende antilichaam kwaliteit en prevalentie van histon-modificatie of binding gebeurtenissen. Dezelfde conclusie kan worden getrokken, zelfs als een andere bin grootte werd gebruikt (Additional file 2: figuur S4). Onze resultaten bevestigden de algemene waarneming dat ChIP-seq meestal een meer onderscheidend signaal profiel dan ChIP-chip produceert.

Naar aanleiding daarvan hebben we de verrijkingsregio’s binnen elk ChIP-profiel gekarakteriseerd. Om een eerlijke vergelijking uit te voeren, zouden we graag een algoritme gebruiken dat peak calling uitvoert op ChIP-seq en ChIP-chip data met dezelfde criteria. Momenteel zijn veel veelgebruikte peak calling algoritmen specifiek ontworpen voor het analyseren van ChIP-chip of ChIP-seq data, maar niet voor beide. Om deze beperking te overwinnen, hebben we pieken geïdentificeerd van zowel ChIP-chip en ChIP-seq profielen met behulp van dezelfde genoom-scanning heuristiek (zie de Methoden sectie). Onze resultaten geven aan dat we bijna altijd ontdekken een groter aantal pieken en smallere pieken met behulp van gegevens die zijn gegenereerd uit ChIP-seq in vergelijking met ChIP-chip bij het analyseren van hetzelfde biologische monster, en deze conclusie is consistent, ongeacht de strengheid van de identificatiecriteria gebruikt (figuur 2 en Additional file 2: figuur S5). In de praktijk kunnen we waarschijnlijk een nog groter aantal smalle pieken in ChIP-seq gegevens identificeren als we expliciet gebruik maken van streng-specifieke informatie binnen de piek-oproep procedure (naast alleen het verschuiven van elke lees in de richting van zijn 5 ‘einde met een constant aantal baseparen), zodat de huidige analyse biedt een ondergrens op de effectiviteit van ChIP-seq in vergelijking met ChIP-chip. Alles bij elkaar genomen, onze resultaten tonen aan dat ChIP-seq een hogere ruimtelijke resolutie en signaal-ruisverhouding biedt.

Genoom-brede signaal reproduceerbaarheid binnen en tussen technologieën

Daarnaast hebben we de reproduceerbaarheid tussen ChIP-chip en / of ChIP-seq profielen geschat op het genoom-brede niveau (1 kb bins). Om vertekeningen als gevolg van verschillen in genomische dekking en sequentie mapping (figuur 1e) te voorkomen, hebben we uitsluiten genomische regio’s die geen microarray probes en regio’s met ongewoon hoge variabiliteit over meerdere INPUT-seq profielen. De Pearson correlatiecoëfficiënt, r, werd gebruikt als een maat voor correlatie, omdat het gevoeliger dan de Spearman correlatiecoëfficiënt voor het vergelijken van de staart van twee signaal distributies, die bijzonder belangrijk is bij het analyseren van ChIP verrijking signaal profielen. De correlatie tussen ChIP-chip replicaatparen en tussen ChIP-seq replicaatparen is over het algemeen hoog (mediaan r = 0,85 en 0,82, respectievelijk), wat erop wijst dat beide technologieën reproduceerbare resultaten kunnen produceren. Zoals verwacht, is de cross-platform correlatie tussen replicaatparen van ChIP-chip en ChIP-seq profielen bescheidener (mediaan r = 0,41; Additional file 1: Tabel S5). Vergelijkbare conclusies kunnen worden bereikt, zelfs als we gebruik maken van verschillende bin maten voor de berekening van inter-profiel correlatie (Additional file 2: Figuur S6). Een representatieve scatter plot vergelijking van elk paar van technologieën wordt getoond in figuur 3b-d. We zien ook een positieve correlatie tussen de skewness en inter-profiel reproduceerbaarheid (Additional file 2: figuur S7), wat suggereert dat meer gevoelige antilichamen kunnen produceren meer consistente profielen tussen de twee technologieën.

Bouw van het gemiddelde signaal profiel bij TSS en TES

Bouw van het gemiddelde ChIP-signaal profielen rond belangrijke genomische kenmerken zoals TSS en TES is een veel voorkomende manier om signaalverrijking te visualiseren rond deze kenmerken. Daarom onderzochten we de reproduceerbaarheid van de gemiddelde TSS en TES profielen (2 kb omhoog en 2 kb stroomafwaarts) voor elk paar replicaat ChIP-profielen (Additional file 2: figuur S8). De gemiddelde profielen van de meeste repliceren paren zijn zeer consistent. Er zijn echter een paar paren die significant verschillend zijn, met name de profielen van H3K27Me3 en H3K9Me3 in zowel stadium E-16-20 h en E-20-24 h (Additional file 2: figuren S8c en S8g). Zonder externe validatie is het onmogelijk om te bepalen of de gemiddelde signaalprofielen gegenereerd door ChIP-chip of ChIP-seq nauwkeuriger zijn. Niettemin leidden twee lijnen van bewijsmateriaal ons om te geloven dat de gemiddelde signaalprofielen van ChIP-chip waarschijnlijk nauwkeuriger waren. Ten eerste hadden alle drie ChIP-chip replicaten op deze tijdstippen zeer consistente gemiddelde profielen. Ten tweede, de ChIP-seq gemiddelde signaal profielen op deze biologische omstandigheden leek op de trend van GC-gehalte variatie op TSS en TES (figuur 1c). De ongewoon lage correlaties tussen GC-inhoud en de INPUT-seq profielen van E-16-20 h en E-20-24 h (figuur 1b en Additional file 2: figuur S2b) bracht ons tot de hypothese dat de waargenomen discrepantie was te wijten aan een verkeerde voorstelling van GC-inhoud variatie door de respectieve INPUT-seq profielen. Zowel H3K27Me3 en H3K9Me3 zijn repressieve merktekens die meestal uitgeput zijn bij TSSs en TESs, zodat elke variatie in achtergrondaftrekking waarschijnlijk veel meer uitgesproken is dan andere histonmerken die een sterke signaalverrijking hebben bij deze genomische kenmerken. Om onze hypothese te testen, hebben we de overeenkomstige INPUT-seq achtergrond vervangen door de INPUT-seq van het AdultFemale monster, aangezien dit de hoogste correlatie heeft met GC-gehalte variatie. Na de vervanging, de gemiddelde signaal profielen gegenereerd door ChIP-seq en ChIP-chip op deze twee ontwikkelingsstadia eens (figuur 4 en Additional file 2: figuur S9). Dit resultaat is opvallend omdat het laat zien dat het gebruik van verschillende INPUT-seq als negatieve controle van dezelfde ChIP-seq profiel kan leiden tot wezenlijk verschillende interpretatie van de gegevens.

Effect van het gebruik van verschillende invoerprofielen bij de normalisatie van ChIP-seq-gegevens

Na de impact van INPUT-seq bij de constructie van gemiddelde TSS- en TES-profielen te hebben waargenomen, vroegen we ons af of het gebruik van verschillende INPUT-seq-profielen voor achtergrondnormalisatie de resultaten van ChIP-seq-piekoproeping aanzienlijk beïnvloedt. We gebruikten SPP om pieken op te roepen voor 10 van onze ChIP-seq monsters (CBP, H3K9Ac, H3K9Me3, H3K27Ac, H3K27Me3 op E16-20 h en E20-24 h) waarbij elk ChIP-profiel genormaliseerd werd tegen vier verschillende INPUT-seq als achtergrond (de input van het overeenstemmende tijdstip, AdultFemale, AdultMale, en E-4-8 h). Deze INPUT-seq profielen werden gekozen omdat ze verschillende sequencing diepte en correlatie met GC-gehalte (figuur 1b). Een vergelijking van het aantal pieken en de mediane piekbreedte wordt getoond in figuur 5. We zagen een groot verschil in het aantal pieken dat werd opgeroepen voor elk ChIP-seq monster wanneer verschillende INPUT-seq werd gebruikt als achtergrond. In het extreme geval (E-16-24 h, H3K9Me3 ChIP), kan het aantal pieken verandert van nul tot bijna 40.000 bij een FDR van 5% (figuur 5a). In het algemeen werden meer statistisch significante pieken (FDR < 0,05) gedetecteerd bij normalisatie tegen een diep gesequenced input DNA monster (AdultMale en E-4-8 h in dit experiment), hoewel de absolute omvang van het verschil varieert tussen ChIP datasets. Het verschil in piekaantal wijst waarschijnlijk op een verschil in detectievermogen. Voor elk ChIP monster berekenden we het aandeel van de overlap tussen elk paar van de piek sets gegenereerd door vier verschillende input DNA achtergrond (dat wil zeggen, zes vergelijkingen per ChIP monster). Wij vonden dat het gemiddelde aandeel van overlapping met betrekking tot de kleinere piekreeks ongeveer 95% bedraagt, wat erop wijst dat de verschillen in aantal gedetecteerde pieken waarschijnlijk te wijten zijn aan een verschillend vermogen om zwakkere pieken op te roepen. Wij stelden vast dat de sterke pieken (d.w.z. die met een lage detectie-FDR) vaker in verschillende piekreeksen werden gedetecteerd (zie Additional file 2: figuur S10 voor een voorbeeld). De mediane breedte van de gedetecteerde pieken wordt ook beïnvloed door het gebruik van verschillende INPUT-seq als achtergrond (figuur 5b). Deze analyse toonde aan dat de normalisatie met behulp van verschillende INPUT-seq kan een belangrijke, en ondergewaardeerde, invloed hebben op piek calling.

Beoordeling van variatie als gevolg van het gebruik van verschillende peak callers

Een andere belangrijke bron van variatie in de analyse van ChIP-chip- en ChIP-seq-profielen is afkomstig van het gebruik van verschillende analysealgoritmen. Een groot aantal openbaar beschikbare ChIP-chip en ChIP-seq analyse-instrumenten zijn ontwikkeld tot op heden, en alle van hen maken gebruik van verschillende methoden voor tag shifting, profiel normalisatie, afvlakken, piek identificatie, en berekening van false discovery rate. Het is daarom niet al te verwonderlijk dat verschillende peak callers heel verschillende resultaten kunnen genereren in termen van identificatie van bindingsplaatsen, vooral wanneer het gaat om pieken met zwakke signalen. Met behulp van ons compendium van ChIP-chip en ChIP-seq datasets, konden we nagaan hoeveel variatie in piekidentificatie kan worden toegeschreven aan het gebruik van verschillende profiling technologie en het gebruik van verschillende peak callers. In deze studie analyseerden wij onze ChIP-chip profielen met behulp van twee peak callers: MA2C en Splitter en onze ChIP-seq-profielen geanalyseerd met twee andere piek-callers: MACS en SPP (zie aanvullend bestand 1: tabel S8). Deze peak callers werden gekozen omdat ze veel gebruikt worden, publiek beschikbaar zijn, en over het algemeen goede prestaties laten zien in eerdere vergelijkende studies. We berekenden de overlap van de top 1.000 pieken van vier van de factoren (CBP, H3K4Me1, H3K4Me3, en H3K27Me3) over meerdere ontwikkelingsstadia. De vier IP-factoren werden gekozen omdat het representatieve profielen waren met brede pieken (CBP en H3K27Me3) en smalle pieken (H3K4Me1 en H3K4Me3). Wij presenteren hier alleen de resultaten van de vergelijking van de top 1000 pieken, aangezien dit een biologisch redelijk aantal aanverrijkingsplaatsen met hoge betrouwbaarheid in deze profielen is. De algemene conclusie van deze analyse is robuust tegen een verscheidenheid van piek calling drempels (Additional file 2: figuur S11). Overeenstemming tussen twee piek sets werd gemeten door het gemiddelde aandeel van overlappende pieken. Zoals blijkt uit figuur 6, de vergelijkingen op basis van profielen van H3K4Me1 en H3K4Me3 leverde verwachte resultaten, waarin de intra-platform concordantie is hoger dan cross-platform concordantie (dat wil zeggen, piek sets gegenereerd door twee piek bellers op hetzelfde profiel zijn meer concordant dan piek sets gegenereerd door twee piek bellers op twee profielen). De intra-platform concordantie kan echter even laag zijn als de inter-platform concordantie wanneer de profielen van H3K27Me3 en CBP worden geanalyseerd, hetgeen impliceert dat de variatie in algoritmen voor piekafroepingen even groot kan zijn als het gebruik van verschillende profilingtechnologieën voor sommige IP-factoren. De waarneming dat de huidige “peak calling”-algoritmen minder eensluidende resultaten opleveren voor ChIP-profielen met brede domeinen (CBP en H3K27Me3) dan die met scherpe pieken (H3K4Me1 en H3K4Me3) kan erop wijzen dat zij minder consistent zijn bij het identificeren van brede verrijkingsgebieden, wat een interessant onderwerp voor verder onderzoek kan zijn.