ChIP-chip versus ChIP-seq: Lessons for experimental design and data analysis

Analiza profilelor ADN de intrare

Pentru a înțelege diferențele tehnologice dintre datele ChIP bazate pe secvențiere și cele bazate pe microarray, am analizat mai întâi profilele fragmentelor de ADN reticulate și sonicate (ADN de intrare) generate de microarray (INPUT-chip) și de secvențiere de mare randament (INPUT-seq). Deoarece profilul ADN de intrare ar trebui să fie independent de anticorpul utilizat pentru ChIP, această comparație poate oferi o perspectivă asupra diferențelor specifice dintre aceste două tehnologii de profilare. Am obținut date INPUT-chip din canalul de fond al datelor noastre de microarray cu două canale. În timp ce această platformă de microarray utilizează hibridizarea competitivă, s-a demonstrat că cele două canale din microarray-ul nostru Agilent sunt relativ independente, deoarece saturația în oricare dintre canale este foarte rară . Dintre toate profilurile INPUT-chip pe care le-am extras, prezentăm aici doar analiza a opt profiluri reprezentative (două din fiecare dintre cele patru momente de dezvoltare), deoarece majoritatea profilurilor INPUT-chip sunt foarte asemănătoare (Fișier suplimentar 2: Figura S1). Cele opt profiluri INPUT-chip au fost apoi comparate cu cele nouă profiluri INPUT-seq colectate în acest studiu (Fișier suplimentar 1: Tabelul S3).

Una dintre cele mai izbitoare observații este că profilurile INPUT-chip și INPUT-seq par a fi substanțial diferite, chiar dacă același material ADN de intrare a fost utilizat pentru hibridizarea microarray și secvențiere (Figura 1). Magnitudinea relativă și localizarea vârfurilor par a fi consecvente între profilurile INPUT-chip din mai multe experimente. Cu toate acestea, modelele din cele nouă profiluri INPUT-seq par a fi mai variabile. Putem identifica vizual multe regiuni care prezintă o îmbogățire inconsecventă a semnalului în mai multe profiluri INPUT-seq (evidențiate în figura 1a). S-a efectuat o analiză de grupare pentru a cuantifica această observație. Am constatat că toate cele opt profiluri INPUT-chip s-au grupat strâns unul lângă celălalt (Figura 1b). Acest rezultat arată că distribuția ADN de fond măsurată din microarray și secvențierea de mare randament este diferită. Toate profilurile INPUT-chip și șapte din cele nouă profiluri INPUT-seq s-au corelat pozitiv cu conținutul GC genomic la nivelul întregului genom (figura 1b), precum și în jurul situsurilor de început de transcriere (TSS) și al situsurilor de sfârșit de transcriere (TES) (figura 1c). Forța corelației cu GC este foarte consistentă între profilurile INPUT-chip, dar foarte variabilă între profilurile INPUT-seq (Figura 1b-c și Fișierul suplimentar 2: Figura S2). În special, profilurile INPUT-seq obținute la E-16-20 h (E16) și E-20-24 h (E20) nu se corelează cu conținutul de GC.

De asemenea, observăm că INPUT-seq cu o adâncime de secvențiere mai mare (>4 milioane de citiri cartografiate) tind să se grupeze mai strâns decât cele cu o adâncime de secvențiere mai mică, sugerând că poate exista o relație între adâncimea de secvențiere și variabilitatea ADN-ului de intrare. Pentru a testa această ipoteză, am generat 11 profile INPUT-seq suplimentare prin subeșantionarea citirilor de secvențiere din eșantionul de ADN de intrare cel mai profund secvențiat (AdultMale; AM) la diferite proporții de eșantionare (Figura 1d și Fișierul suplimentar 2: Figura S3). După cum era de așteptat, profilurile cu o adâncime de secvențiere mai mare tind să se grupeze mai puternic, iar corelația lor cu variația conținutului GC este mai consistentă. Cu toate acestea, corelația conținutului GC devine mult mai slabă doar la o adâncime de secvențiere foarte mică (<2 milioane de citiri; Figura 1d). Acest lucru indică faptul că adâncimea scăzută de secvențiere nu este singurul factor care afectează calitatea INPUT-seq. Mai mult decât atât, unele INPUT-seq cu o adâncime de secvențiere relativ scăzută (E0 și AF, <4 milioane de citiri) pot oferi profile de ADN de intrare consistente. Acest lucru implică faptul că variabilitatea INPUT-seq poate fi atribuită și altor factori experimentali. Deși sunt necesare studii suplimentare pentru a diseca întreaga gamă de factori experimentali care afectează variabilitatea bibliotecilor de ADN de intrare, aceasta ar putea fi afectată de variațiile în pregătirea probelor (de exemplu, pregătirea diferită a cromatinei și sonicare), variația de la o rulare la alta a secvențiatorului, variația de la secvențiator la secvențiator pentru același model și o multitudine de alte variabile în cadrul experimentelor. Variabilitatea ridicată între profilurile INPUT-seq este într-adevăr o problemă critică, deoarece variabilitatea mare contribuie la instabilitatea estimării densității într-un profil ChIP-seq, ceea ce va afecta analiza datelor din aval. După cum se va arăta în secțiunile ulterioare ale acestei lucrări, un INPUT-seq cu o corelație neobișnuit de slabă cu conținutul de GC poate avea un impact asupra construirii profilelor medii în locații genomice importante. Prin urmare, este imperativ să se secvențieze ADN-ul de intrare la o profunzime suficientă și să se verifice dacă profilul obținut este în concordanță cu cele din experimente similare.

Cuprinderea genomică este un alt considerent cheie atunci când se alege între ChIP-chip și ChIP-seq. Acoperirea genomică a ChIP-chip este limitată de designul sondei microarray, iar acoperirea ChIP-seq depinde de adâncimea de secvențiere. Acoperirea genomică obținută de microarray-ul nostru Agilent este de aproximativ 70%. Utilizând datele INPUT-seq subeșantionate, arătăm că INPUT-seq oferă, în general, o acoperire genomică mai mare la o adâncime de secvențiere de până la un milion de citiri. Această tendință construită din datele subeșantionate aleatoriu coroborează acoperirea genomică observată a celorlalte opt seturi de date INPUT-seq reale (Figura 1e).

Compararea caracteristicilor profilurilor

Apoi am comparat caracteristicile profilurilor ChIP-chip și ChIP-seq. Pentru a compara profilurile generate de cele două tehnologii, am împărțit genomul în intervale de 1 kb care nu se suprapun și am definit nivelul de îmbogățire la fiecare interval ca fiind media raportului logaritmic al canalului IP față de canalul de intrare (a se vedea secțiunea Metode pentru detalii). Ne referim la distribuția semnalului unui profil ChIP ca fiind distribuția sa de valori de îmbogățire a tuturor bins-urilor. În primul rând, am urmărit să caracterizăm raportul mediu dintre semnal și zgomot pentru profilurile generate de cele două tehnologii. Am utilizat asimetria (trunchiată) a profilului de densitate a semnalului după eliminarea semnalelor din cele mai înalte și cele mai joase 5 % ale distribuției ca măsură a raportului semnal-zgomot al unui profil. Skewness este o măsură a asimetriei unei distribuții, iar un skewness pozitiv indică faptul că coada din partea dreaptă este mai lungă, ceea ce implică un raport semnal/zgomot bun. În aproape toate cazurile, un profil ChIP-seq are o skewness mai mare decât profilul ChIP-chip corespunzător din aceeași condiție biologică (figura 2 și fișierul suplimentar 1: tabelul S4). Observăm că diferența de skewness depinde de factorul IP, care s-ar putea datora calității diferite a anticorpilor și prevalenței modificărilor histonice sau a evenimentelor de legare. Aceeași concluzie poate fi trasă chiar dacă a fost utilizată o dimensiune diferită a binului (Fișier suplimentar 2: Figura S4). Rezultatele noastre au confirmat observația generală conform căreia ChIP-seq produce de obicei un profil de semnal mai distinctiv decât ChIP-chip.

În continuare, am caracterizat regiunile de îmbogățire din cadrul fiecărui profil ChIP. Pentru a efectua o comparație corectă, am dori să folosim un algoritm care să efectueze apelarea vârfurilor pe datele ChIP-seq și ChIP-chip folosind aceleași criterii. În prezent, mulți algoritmi de apelare a vârfurilor utilizați în mod obișnuit sunt concepuți special pentru analiza datelor ChIP-chip sau ChIP-seq, dar nu pentru ambele. Pentru a depăși această limitare, am identificat vârfurile atât din profilurile ChIP-chip, cât și din cele ChIP-seq, utilizând aceeași euristică de scanare a genomului (a se vedea secțiunea Metode). Rezultatele noastre indică faptul că putem descoperi aproape întotdeauna un număr mai mare de vârfuri și vârfuri mai înguste utilizând date generate de ChIP-seq în comparație cu ChIP-chip atunci când analizăm aceeași probă biologică, iar această concluzie este consecventă indiferent de rigurozitatea criteriilor de identificare utilizate (Figura 2 și Fișierul suplimentar 2: Figura S5). În practică, putem probabil să identificăm un număr și mai mare de vârfuri înguste în datele ChIP-seq dacă utilizăm în mod explicit informațiile specifice șirului în cadrul procedurii de apelare a vârfurilor (pe lângă faptul că deplasăm doar fiecare citire spre capătul său 5′ cu un număr constant de perechi de baze), astfel încât analiza actuală oferă o limită inferioară a eficienței ChIP-seq în comparație cu ChIP-chip. Luate împreună, rezultatele noastre demonstrează că ChIP-seq oferă o rezoluție spațială și un raport semnal-zgomot mai mare.

Reproductibilitatea semnalului la nivel de genom în cadrul și între tehnologii

În continuare, am estimat reproductibilitatea între profilurile ChIP-chip și/sau ChIP-seq la nivel de genom (bini de 1 kb). Pentru a evita prejudecățile datorate diferențelor în ceea ce privește acoperirea genomică și cartografierea secvenței (figura 1e), am exclus regiunile genomice care nu conțin nicio sondă de microarray și regiunile cu o variabilitate neobișnuit de mare între mai multe profiluri INPUT-seq. Coeficientul de corelație Pearson, r, a fost utilizat ca măsură a corelației, deoarece este mai sensibil decât coeficientul de corelație Spearman pentru compararea cozii a două distribuții de semnal, ceea ce este deosebit de important în analiza profilurilor de semnal de îmbogățire ChIP. Corelația între perechile de replici ChIP-chip și între perechile de replici ChIP-seq este, în general, ridicată (mediana r = 0,85 și, respectiv, 0,82), ceea ce indică faptul că ambele tehnologii pot produce rezultate reproductibile. Așa cum era de așteptat, corelația între platforme între perechile de replici ale profilurilor ChIP-chip și ChIP-seq sunt mai modeste (r median = 0,41; Fișier suplimentar 1: Tabelul S5). Se poate ajunge la concluzii similare chiar dacă folosim dimensiuni diferite ale binelor pentru calcularea corelației între profiluri (Fișier suplimentar 2: Figura S6). Un grafic de dispersie reprezentativ care compară fiecare pereche de tehnologii este prezentat în figura 3b-d. Observăm, de asemenea, o corelație pozitivă între skewness și reproductibilitatea interprofil (Fișier suplimentar 2: Figura S7), sugerând că anticorpii mai sensibili pot produce profiluri mai consistente între cele două tehnologii.

Construcția profilului de semnal mediu la TSS și TES

Construcția profilurilor medii de semnal ChIP în jurul unor caracteristici genomice importante, cum ar fi TSS și TES, este o modalitate obișnuită de a vizualiza îmbogățirea semnalului în jurul acestor caracteristici. Prin urmare, am investigat reproductibilitatea profilelor medii TSS și TES (2 kb în sus și 2 kb în jos) pentru fiecare pereche de profile ChIP replicate (Fișier suplimentar 2: Figura S8). Profilurile medii ale celor mai multe perechi de replici sunt foarte consistente. Cu toate acestea, există câteva perechi care sunt semnificativ diferite, în special profilurile de H3K27Me3 și H3K9Me3 atât în stadiul E-16-20 h, cât și în stadiul E-20-24 h (Fișier suplimentar 2: Figurile S8c și S8g). Fără o validare externă, este imposibil de determinat dacă profilurile medii de semnal generate de ChIP-chip sau ChIP-seq sunt mai precise. Cu toate acestea, două linii de dovezi ne-au determinat să credem că este mai probabil ca profilurile medii de semnal generate de ChIP-chip să fie mai precise. În primul rând, toate cele trei replici ChIP-chip la aceste momente de timp au avut profile medii foarte consistente. În al doilea rând, profilurile medii ale semnalului ChIP-seq în aceste condiții biologice au semănat cu tendința de variație a conținutului de GC la TSS și TES (Figura 1c). Corelațiile neobișnuit de scăzute dintre conținutul de GC și profilurile INPUT-seq de la E-16-20 h și E-20-24 h (Figura 1b și Fișierul suplimentar 2: Figura S2b) ne-au determinat să emitem ipoteza că discrepanța observată se datorează unei reprezentări eronate a variației conținutului de GC de către profilurile INPUT-seq respective. Atât H3K27Me3, cât și H3K9Me3 sunt mărci represive care sunt, de obicei, epuizate la TSS-uri și TES-uri, astfel încât orice variație în sustragerea fondului este probabil mult mai pronunțată decât alte mărci histonice care au o îmbogățire puternică a semnalului la aceste caracteristici genomice. Pentru a ne testa ipoteza, am înlocuit fundalul INPUT-seq corespunzător cu INPUT-seq din eșantionul AdultFemale, deoarece acesta are cea mai mare corelație cu variația conținutului de GC. După înlocuire, profilurile medii de semnal generate de ChIP-seq și ChIP-chip în aceste două stadii de dezvoltare sunt în concordanță (Figura 4 și Fișierul suplimentar 2: Figura S9). Acest rezultat este surprinzător, deoarece arată că utilizarea unor INPUT-seq diferite ca și control negativ al aceluiași profil ChIP-seq poate duce la o interpretare substanțial diferită a datelor.

Efectul utilizării unor profiluri de intrare diferite în normalizarea datelor ChIP-seq

După ce am observat impactul INPUT-seq în construirea profilurilor medii TSS și TES, ne-am întrebat dacă utilizarea unor profiluri INPUT-seq diferite pentru normalizarea fondului afectează în mod semnificativ rezultatele de apelare a vârfurilor ChIP-seq. Am folosit SPP pentru a apela vârfurile pentru 10 dintre probele noastre ChIP-seq (CBP, H3K9Ac, H3K9Me3, H3K27Ac, H3K27Me3 la E16-20 h și E20-24 h), unde fiecare profil ChIP a fost normalizat față de patru INPUT-seq diferite ca fundal (intrarea din punctul de timp corespunzător, AdultFemale, AdultMale și E-4-8 h). Aceste profiluri INPUT-seq au fost alese deoarece au o adâncime de secvențiere diferită și o corelație cu conținutul de GC (Figura 1b). O comparație a numărului de vârfuri și a lățimii mediane a vârfurilor este prezentată în figura 5. Am observat o diferență mare în ceea ce privește numărul de vârfuri apelate pentru orice eșantion ChIP-seq atunci când s-a utilizat INPUT-seq diferit ca fundal. În cazul extrem (E-16-24 h, H3K9Me3 ChIP), numărul de vârfuri poate varia de la zero la aproape 40 000 la un FDR de 5 % (Figura 5a). În general, au fost detectate mai multe vârfuri semnificative din punct de vedere statistic (FDR < 0,05) atunci când s-au normalizat față de o probă de ADN de intrare secvențiată în profunzime (AdultMale și E-4-8 h în acest experiment), deși magnitudinea absolută a diferenței variază în funcție de seturile de date ChIP. Diferența în ceea ce privește numărul de vârfuri indică probabil o diferență în ceea ce privește puterea de detecție. Pentru fiecare eșantion ChIP, am calculat proporția de suprapunere între fiecare pereche de seturi de vârfuri generate de patru medii diferite de ADN de intrare (adică șase comparații pe eșantion ChIP). Am constatat că proporția medie de suprapunere în ceea ce privește setul de vârfuri mai mic este de aproximativ 95 %, ceea ce indică faptul că diferențele în ceea ce privește numărul de vârfuri detectate se datorează probabil puterii diferite de a apela vârfurile mai slabe. Am observat că vârfurile puternice (și anume, cele cu FDR de detecție scăzută) au fost mai probabil detectate în diferite seturi de vârfuri (a se vedea fișierul suplimentar 2: Figura S10 pentru un exemplu). Lățimea mediană a vârfurilor detectate este, de asemenea, afectată de utilizarea diferitelor INPUT-seq ca fundal (Figura 5b). Această analiză a arătat că normalizarea folosind INPUT-seq diferite poate avea un impact semnificativ și subapreciat asupra apelării vârfurilor.

Evaluarea variației datorate utilizării unor apelatori de vârfuri diferiți

O altă sursă importantă de variație în analiza profilurilor ChIP-chip și ChIP-seq provine din utilizarea unor algoritmi de analiză diferiți. Un număr mare de instrumente de analiză ChIP-chip și ChIP-seq disponibile publicului au fost dezvoltate până în prezent , și toate acestea utilizează metode diferite pentru deplasarea etichetelor, normalizarea profilului, netezirea, identificarea vârfurilor și calcularea ratei de descoperire falsă. Prin urmare, nu este prea surprinzător să constatăm că diferite dispozitive de apelare a vârfurilor pot genera rezultate destul de diferite în ceea ce privește identificarea situsului de legare, în special atunci când este vorba de vârfuri cu semnale slabe . Utilizând compendiul nostru de seturi de date ChIP-chip și ChIP-seq, am putut evalua în ce măsură variația în identificarea vârfurilor poate fi atribuită utilizării unei tehnologii de profilare diferite și utilizării unor apelatori de vârfuri diferiți. În acest studiu, am analizat profilurile noastre ChIP-chip utilizând doi apelatori de vârfuri: MA2C și Splitter și am analizat profilurile noastre ChIP-seq utilizând alți doi apelatori de vârfuri: MACS și SPP (a se vedea fișierul suplimentar 1: Tabelul S8). Aceste peak callers au fost alese deoarece sunt utilizate pe scară largă, sunt disponibile în mod public și, în general, prezintă performanțe bune în studiile comparative anterioare . Am calculat suprapunerea primelor 1.000 de vârfuri din patru dintre factori (CBP, H3K4Me1, H3K4Me3 și H3K27Me3) în mai multe etape de dezvoltare. Cei patru factori IP au fost aleși deoarece erau profiluri reprezentative care conțineau vârfuri largi (CBP și H3K27Me3) și vârfuri înguste (H3K4Me1 și H3K4Me3). Aici, prezentăm doar rezultatele comparării primelor 1.000 de vârfuri, deoarece acesta este un număr rezonabil din punct de vedere biologic de site-uri de îmbogățire cu grad ridicat de încredere în aceste profiluri. Concluzia generală a acestei analize este robustă în raport cu o varietate de praguri de apelare a vârfurilor (Fișier suplimentar 2: Figura S11). Concordanța între două seturi de vârfuri a fost măsurată prin proporția medie de vârfuri care se suprapun. După cum se arată în figura 6, comparațiile bazate pe profilurile de H3K4Me1 și H3K4Me3 au dat rezultatele așteptate, în care concordanța intraplatformă este mai mare decât concordanța între platforme (adică, seturile de vârfuri generate de doi apelanți de vârfuri pe același profil sunt mai concordante decât seturile de vârfuri generate de doi apelanți de vârfuri pe două profiluri). Cu toate acestea, concordanța intraplatformă poate fi la fel de scăzută ca și concordanța interplatformă atunci când se analizează profilurile H3K27Me3 și CBP, ceea ce implică faptul că variația algoritmilor de apelare a vârfurilor poate fi la fel de mare ca și utilizarea unor tehnologii de profilare diferite pentru anumiți factori IP. Observația conform căreia algoritmii actuali de apelare a vârfurilor produc rezultate mai puțin concordante pentru profilurile ChIP cu domenii largi (CBP și H3K27Me3) decât cele cu vârfuri ascuțite (H3K4Me1 și H3K4Me3) poate sugera că aceștia sunt mai puțin consecvenți în identificarea regiunilor de îmbogățire largă, ceea ce poate fi un subiect interesant pentru investigații suplimentare.