Un pachet ușor de utilizat, autorun pentru analize de metilare a ADN-ului
Pentru a finaliza mai convenabil analiza datelor de metilare a ADN-ului, am împachetat toate funcțiile într-un pachet ușor de utilizat, autorun pentru analiza de metilare a ADN-ului. Figura 1 prezintă principalele caracteristici ale BatMeth2: 1) BatMeth2 are o performanță de aliniere eficientă și precisă. 2) BatMeth2 poate calcula nivelul de metilare a ADN (ML) a situsurilor de citosină individuale sau a oricăror regiuni funcționale, cum ar fi cromozomi întregi, regiuni de gene, elemente transpozabile (TE) etc. 3) După integrarea diferiților algoritmi statistici, BatMeth2 poate efectua analize diferențiale de metilare a ADN-ului pentru orice regiune, orice număr de probe de intrare și cerințele utilizatorului. 4) Prin integrarea vizualizării datelor BS-Seq (distribuția metilației ADN pe cromozomi și gene) și a adnotării metilației diferențiale, BatMeth2 poate vizualiza mai clar datele de metilare a ADN. În timpul executării instrumentului BatMeth2, se generează un raport html pentru statisticile eșantionului. Detaliile eșantionului de raport html sunt prezentate în http://htmlpreview.github.io/?https://github.com/GuoliangLi-HZAU/BatMeth2/blob/master/BatMeth2-Report/batmeth2.html.
Figura 1
Funcția de lucru a BatMeth2. Cele două săgeți mari semnifică fișierele de intrare sau de ieșire
BatMeth2 are performanțe mai bune de cartografiere pe date BS-Seq simulate
Am evaluat mai întâi toate alinierele folosind seturi de date simulate (fără indels) constând în lecturi cu 75 de perechi de baze (pb), 100 pb și 150 pb și cu diferite rate de conversie a bisulfitului (variind de la 0 la 100% cu pas de 10%). Aceste seturi de date au fost simulate din genomul uman (UCSC hg19) folosind FASTX-mutate-tools , wgsim (v0.3.0) și simulatorul din SAMtools (v1.1) , care permite 0,03% indels, o rată de eroare de bază de 1% în întregul genom și un maxim de două nepotriviri pe citire. Am cartografiat citirile simulate la genomul de referință, permițând cel mult două nepotriviri. Deoarece pozițiile originale ale citirilor simulate erau cunoscute, am putut evalua acuratețea tuturor programelor comparând rezultatele cartografierii lor cu pozițiile originale.
Pentru a compara performanțele diferitelor programe informatice, o citire de secvențiere cu indels a fost considerată corect cartografiată dacă erau îndeplinite următoarele condiții: 1) citirea a fost cartografiată în mod unic pe aceeași filieră de pe care a fost simulată și calitatea cartografierii a fost mai mare decât 0; 2) poziția de pornire raportată a citirii aliniate a fost în termen de zece perechi de baze față de poziția de pornire originală a citirii simulate; 3) rezultatele cartografierii au avut indels sau neconcordanțe similare cu cele ale citirii simulate. În cazul în care oricare dintre aceste condiții a fost încălcată, citirea a fost considerată greșit cartografiată. Deoarece BatMeth2 permite un singur decalaj în regiunea de însămânțare, acesta poate găsi locații de însămânțare care includ indeli cu o precizie ridicată și poate evita locațiile nepotrivite, ceea ce ar face ca lecturile care includ indeli să fie aliniate greșit. Rezultatele din Fig. 2 arată că BatMeth2 a obținut cel mai mare număr de citiri aliniate corect și cel mai mic număr de citiri aliniate incorect în toate seturile de date de testare la diferite rate de conversie a bisulfitului.
Evaluarea tuturor alinierilor BS-Seq utilizând seturi de date simulate cu diferite lungimi de citire din FASTX și wgsim. Datele simulate cu diferite rate de conversie a bisulfitului sunt prezentate în forme diferite. Rezultatele de la diferite aliniere sunt prezentate cu culori diferite ale simbolurilor. Rezultatele din apropierea colțului din stânga sus în fiecare panou arată că software-ul a obținut un număr mai mare de citiri corect cartografiate și un număr mai mic de citiri cartografiate incorect. Rezultatele de la alignerul nostru BatMeth2 sunt cele mai bune în diferitele seturi de date de bisulfite simulate
În rezumat, rezultatele din seturile de date indel-aberrante simulate de wgsim arată că BatMeth2 are performanțe mai bune (cu 1~2% mai bune decât al doilea aligner de top) decât celelalte metode atunci când aliniază lecturi BS simulate generale care conțin un amestec de nepotriviri și indeli. Putem observa că, odată cu creșterea ratei de conversie BS, precizia de aliniere a tuturor programelor software este redusă. În aceste condiții diferite, BatMeth2 se comportă mai bine.
BatMeth2 are o performanță mai bună de cartografiere pe date BS-Seq reale
Pentru a testa performanța BatMeth2 pe seturi de date BS-Seq reale, am descărcat seturi de date BS-Seq împerecheate și am extras la întâmplare 1 milion de citiri împerecheate de 2 × 90 bp din SRA SRR847318, 1 milion de citiri de 2 × 101 bp împerecheate din SRA SRR1035722 și 1 milion de citiri de 2 × 125 bp împerecheate din SRA SRR3503136 în scopul evaluării. Deoarece aceste seturi de date provin din linii celulare sau țesuturi sănătoase, se așteaptă ca acestea să conțină un număr redus de variații structurale. Prin urmare, am aliniat aceste date reale utilizând lecturi cu un singur capăt din seturile de date împerecheate și am evaluat ratele de cartografiere concordante și discordante din alinierile împerecheate pentru a estima ratele de aliniere corectă și incorectă. Având în vedere că dimensiunea de inserție a lecturilor de capăt împerecheate a fost de aproximativ 500 pb, o pereche de lecturi partenere ar putea fi considerată concordantă dacă au fost mapate la o distanță nominală de 500 pb; în caz contrar, o pereche de lecturi partenere ar putea fi considerată discordantă. Similar cu rezultatele noastre cu datele simulate, BatMeth2 a raportat mai multe alinieri concordante și mai puține discordante pe seturile de date reale pe o gamă largă de scoruri de calitate a hărților, după cum se arată în Fig. 3.
Ratele de concordanță și discordanță ale alinierilor pe lecturi reale împerecheate de la diferiți alinieri. Numerele cumulate ale alinierilor concordante și discordante de la o calitate de cartografiere ridicată la una scăzută pentru lecturi reale de secvențiere cu bisulfit. Există un singur punct pentru BSmap și pentru aliniatoarele bazate pe bowtie separat, deoarece aceste aliniatoare nu au un scor de calitate a hărților. Bismark-bowtie2L15 înseamnă alinierea bowtie2 cu lungimea semințelor 15
În plus, tabelul 1 prezintă timpii de execuție relativi ai programelor. BatMeth2 cu setările implicite a rulat mai repede decât majoritatea aliniatoarelor publicate și a fost comparabil cu BWA-meth și BatMeth. Bismark2 (cu Bowtie2 ca metodă fundamentală de cartografiere), BS Seeker2 și BSmooth necesită timpi de execuție mai lungi.
DNA methylation calling
Pentru a evalua acuratețea apelării metilației ADN între diferite programe software, am descărcat 450 K date de cipuri de mărgele din linia celulară IMR90 de la ENCODE (Encyclopedia of DNA Elements). De asemenea, am descărcat datele de secvențiere a bisulfitului întregului genom (WGBS-Seq) din linia celulară IMR90 de la ENCODE (42,6 Gbase). Pentru fiecare software, am aliniat citirile WGBS-Seq și am calculat nivelul de metilare a ADN-ului. Apoi, am comparat rezultatele cu ML-urile la aceleași situri din datele Bead Chip 450 K. Atunci când diferența dintre ML de ADN din datele WGBS-Seq de către software și cea din 450 K Bead Chip a fost mai mică de 0,2, rezultatul apelării a fost definit ca fiind corect; în caz contrar, a fost considerat incorect.
Rezultatele sunt prezentate în tabelul 2. Suprapunerea între rezultatele corecte ale tuturor programelor software este prezentată în fișierul suplimentar 1: Figura S2. Se poate observa că BatMeth2 și Biscuit au performanțe similare, care sunt mai bune decât cele ale celorlalte software-uri. În concluzie, BatMeth2 îmbunătățește acuratețea atât a alinierii citirilor BS, cât și a apelării ADN ML.
BatMeth2 aliniază citirile BS, permițând în același timp indelurile de lungime variabilă
Cancerul conține o proporție considerabil mai mare de indeluri decât celulele sănătoase. Prin urmare, pentru a verifica dacă BatMeth2 poate alinia lecturi BS cu indeli de diferite lungimi, am descărcat de la ENCODE date WGBS (75 Gbase) și date Bead Chip 450 K de la HepG2 (carcinom hepatocelular hepatic, o linie celulară canceroasă). Am verificat distribuția lungimii indelilor în lecturi după alinierea datelor HepG2 WGBS-Seq. Fișier suplimentar 1: Figura S3A arată că lungimile indelurilor detectate au fost distribuite în principal în intervalul 1 bp~ 5 bp, iar cel mai lung indel a avut o lungime de 40 bp. Conform statisticilor noastre, 2,3 % din citirile de aliniere conțineau indeli. Din aceste rezultate, știm că BatMeth2 poate alinia lecturi cu indeluri de diferite lungimi.
În continuare, am testat efectul detectării indelurilor asupra apelării metilării ADN. Pentru BatMeth2, am rulat două opțiuni pe datele HepG2: cu și fără detectarea indelurilor (adică, am setat parametrul -I în BatMeth2). De asemenea, am rulat Bismark pe datele WGBS-Seq din HepG2 ca referință pentru apelarea metilației ADN cu detectarea indelurilor, deoarece Bismark nu are o funcție de apelare a indelurilor. Am comparat apelarea metilației ADN în BatMeth2 și Bismark cu apelarea din datele 450 K Bead Chip. Rezultatele sunt prezentate în fișierul suplimentar 1: Figura S3B, unde „BatMeth2-noIndel” corespunde BatMeth2 fără detectarea indelurilor. Putem observa că, în absența detectării indelurilor, rezultatul BatMeth2 a fost doar puțin mai bun decât cel al Bismark (cu Bowtie1 ca metodă fundamentală de cartografiere). Rezultatul BatMeth2 cu detectarea indelurilor a fost semnificativ mai bun. În plus, putem observa că BatMeth2 poate detecta mai multe situsuri de metilare a ADN-ului decât BatMeth2-noIndel și Bismark (Bowtie 1). Pentru a înțelege de ce performanța BatMeth2 cu detectarea indelurilor este mai bună, am definit situsurile de metilare apelate de BatMeth2 drept Rezultatul A, în timp ce situsurile de metilare apelate de BatMeth2-noIndel și Bismark au fost definite drept Rezultatul B. Apoi, am lăsat mclA să fie situsurile de metilare care apar în Rezultatul A, dar nu și în Rezultatul B. Am observat că mclA a inclus 23 853 de situsuri de metilare a ADN-ului și 15 048 (63 %) din cele 23 853 de situsuri acoperite de alinierile de citiri indel numite de BatMeth2 cu detectarea indelurilor (a se vedea Fișierul suplimentar 1: Figura S3C). În plus, am constatat că ratele de indel în Rezultatul A și Rezultatul B au fost de numai 5 și, respectiv, 0 %. Prin urmare, am concluzionat că detectarea precisă a indelurilor poate îmbunătăți apelarea metilării ADN.
Vizualizarea datelor de metilare a ADN
BatMeth2 oferă instrumente pentru a vizualiza datele de metilare. Pentru a ilustra caracteristicile de vizualizare ale BatMeth2, am descărcat (1) 117 Gbase de citiri single-end de la linia celulară umană H9, (2) 105,2 Gbase de citiri single-end de la linia celulară umană IMR90 și (3) 12,6 Gbase de citiri paired-end de la orez de tip sălbatic. În primul rând, BatMeth2 poate vizualiza densitatea de metilare a citosinei la nivelul cromozomului. Punctele din Fig. 4a reprezintă o fereastră glisantă de 100 kb cu un pas de 50 kb. Pentru a permite vizualizarea ML la situri CpG sau non-CpG individuale într-un browser de genom, furnizăm, de asemenea, fișiere în format bed și bigWig (Fig. 4b). Prin comparație cu densitatea genelor și a TE, am observat că ML a fost corelată cu densitatea TE și a fost anticorrelată cu densitatea genelor (Fig. 4c). Această tendință a fost observată anterior la orez .
Vizualizarea nivelurilor de metilare în scara cromozomului. a Densitatea de metil-citosină în cromozomul uman 10. Punctele reprezintă nivelurile de metilare în ferestre glisante de 100Kb cu un pas de 50Kb. Punctele roșii se referă la nivelurile de metilare în șirul plus, iar punctele albastre se referă la nivelurile de metilare în șirul minus. b Un exemplu privind distribuția nivelurilor de metilare a ADN-ului și a regiunilor diferențiat metilate (DMR) între liniile celulare H9 și IMR90 în cromozomul uman 10. c Densitatea genelor, a elementelor transpozonice (TE) și nivelul de metilare a ADN-ului în întregul genom al orezului. Panoul A reprezintă rezultatele generate de Batmeth2. Panoul B reprezintă rezultatele vizualizării din browserul UCSC, cu fișierele BED de la Batmeth2
În al doilea rând, BatMeth2 poate vizualiza ML-urile genelor. Mai precis, BatMeth2 poate vizualiza ML-urile la 2 kb în amonte de genă, la locul de început al transcripției (TSS), în corpul genei, la locul de sfârșit al transcripției (TES) și la 2 kb în aval de corpul genei. Comparând regiunile din amonte, din corp și din aval, Fig. 5a arată că ML ADN din corpul genei este mai mare decât cel din regiunea promotoare. Comparând toate cele cinci regiuni, există în mod evident o vale în regiunea TSS (Fig. 5b). BatMeth2 poate calcula, de asemenea, profilurile ML în jurul intronilor, exonilor, regiunilor intergenice și TE-urilor (Fișier suplimentar 1: Figura S4). În plus, BatMeth2 poate furniza o hartă termică a mai multor gene în funcție de regiunea genetică pentru o comparație convenabilă a ML globale ale genelor din diferite eșantioane (Fig. 5c).
Vizualizare a metilației ADN în diferite contexte. a Nivelurile de metilare a ADN-ului în regiuni de 2Kb în amonte de gene, corpuri de gene, 2Kb în aval de corpuri de gene. b Profilul de agregare a nivelurilor de metilare a ADN-ului între gene. c Harta termică a tuturor genelor în regiuni de 2Kb în amonte de gene, corpuri de gene, 2Kb în aval de corpuri de gene
În al treilea rând, BatMeth2 poate vizualiza distribuția metilării ADN-ului. Fișierul suplimentar 1: Figura S5A prezintă distribuțiile de metilare a ADN-ului în liniile celulare H9 și IMR90. În figură, ADN ML este împărțit în cinci categorii: metilat (M: > 80 %), intermediar între parțial metilat și metilat (Mh: 60-80 %), parțial metilat (H: 40-60 %), intermediar între nemetilat și parțial metilat (hU: 20-40 %) și nemetilat (U: < 20 %). După cum se arată în fișierul suplimentar 1: Figura S5A, ML a fost mai mare în linia celulară H9 din categoria M decât în linia celulară IMR90, în special în contextul CpG. În contextul secvenței CH, metilarea CpG este forma predominantă, dar o fracțiune semnificativă de citozine metilate se găsește în situsurile CpA, în timp ce ML este mai mică de 40 %, în special în linia celulară H9 (Fișier suplimentar 1: Figura S5B).
În al patrulea rând, BatMeth2 poate analiza corelația dintre nivelul de expresie a genelor și ML a ADN-ului promotor al genelor. Am ilustrat această caracteristică folosind liniile celulare H9 și IMR90. Nivelurile de expresie ale genelor din H9 sau IMR90 au fost împărțite în diferite categorii. După cum se arată în fișierul suplimentar 1: Figura S5C, genele foarte bine exprimate au prezentat ML-uri mai mici în regiunile lor promotoare. În plus, am împărțit ML-urile promotorilor genelor în cinci categorii. Rezultatul din fișierul suplimentar 1: Figura S5D arată că genele cu promotori care au valori ML mai mari au prezentat niveluri de expresie mai scăzute. Corelația negativă dintre expresia genelor de mamifere și metilarea ADN-ului promotor este cunoscută . Această analiză indică în continuare acuratețea BatMeth2.
Descoperirea citosinelor și regiunilor diferențial metilate (DMCs/DMRs)
Identificarea citosinelor diferențial metilate (DMCs) și a regiunilor diferențial metilate (DMRs) este unul dintre obiectivele majore în analiza datelor de metilare. Deși cercetătorii sunt ocazional interesați de corelarea situsurilor de citosine unice cu un fenotip , DMR-urile sunt caracteristici foarte importante .
Primele studii BS-Seq au profilat celulele fără a colecta replici. Pentru astfel de seturi de date, am folosit testul exact al lui Fisher pentru a discerne citosinele metilate diferențiat (DMC). Pentru seturile de date BS-Seq cu replici, cel mai natural model statistic pentru a numi DMC-urile este distribuția beta-binomială . Știm că o serie de programe software pot efectua analiza diferențială a datelor de metilare a ADN-ului, cum ar fi methylKit (un program de analiză diferențială care necesită replici biologice) și Methy-Pipe (un program de analiză diferențială fără replici biologice). Cu toate acestea, nu este disponibil niciun pachet cuprinzător care să includă atât cartografierea, cât și analiza diferențială a metilației. Astfel, am dezvoltat un pachet care integrează cartografierea cu analiza diferențială. Pentru a facilita identificarea DMR-urilor din datele de bisulfit fără duplicate, am integrat testul exact al lui Fisher pentru a efectua un test de ipoteză. Atunci când un eșantion are două sau mai multe replici, folosim distribuția beta-binomială pentru a efectua analiza de metilare diferențială. De asemenea, furnizăm fișiere bed sau bigWig pentru lista de DMR-uri. DMR-urile pot fi vizualizate într-un browser de genom (Fig. 4b) cu ajutorul fișierelor bed sau bigWig generate.
Ca o ilustrare, Fig. 6a arată numărul de DMR-uri și regiuni în linia celulară IMR90 și în linia celulară H9, așa cum au fost detectate de BatMeth2 (valoare p< 0,05, meth.diff > = 0,6). BatMeth2 poate vizualiza dacă CpG-urile și DMC-urile sunt îmbogățite în anumite regiuni, cum ar fi regiunile genetice, CDS, intron, intergenice, UTR, TE, LTR, LINE și SINE. Figura 6b vizualizează proporțiile de DMC-uri în diferite regiuni genomice. În afară de regiunile intergenice, nu am observat îmbogățirea DMC în nicio regiune.
Analiză de metilare diferențiată. a Rezultatele analizei regiunilor diferențiat metilate (DMR), a genelor diferențiat metilate (DMG) și a promotorilor diferențiat metilați (DMP) între liniile celulare H9 și IMR90. b Adnotarea citosinelor diferențiat-metilate (DMC) în raport cu diferite proprietăți genomice și elemente repetate. c DMP-urile conțin indeli specifici H9 sau IMR90 (portocaliu) ocupă o proporție substanțială în toate DMP-urile (promotori diferențiați de metilare a ADN-ului)
O proporție substanțială de promotori diferențiat metilați (DMP-uri) conțin indeli
Știm că indelii și metilarea ADN-ului joacă un rol important în dezvoltarea țesuturilor și în boli . Aici, examinăm relația dintre promotorii diferențiat metilați (DMPs) și indels. Am realizat acest studiu folosind citirile BS-Seq în liniile celulare IMR90 și H9. Am aliniat mai întâi citirile BS-Seq folosind BatMeth2; apoi, indelurile au fost apelate folosind instrumentele BisSNP și GATK. Ulterior, am definit indelurile care apar numai în H9 sau IMR90 ca indeluri specifice liniei celulare.
Apoi, am detectat 1384 DMP-uri între H9 și IMR90 prin BatMeth2 (valoare p< 0,05, meth.diff > = 0,6). Un total de 236 (17%) dintre toate DMP-urile de mai sus conțin indels, așa cum se arată în Fig. 6c. Pe scurt, o proporție substanțială de DMP-uri conțin indeli. Prin urmare, alinierea exactă a lecturilor BS-Seq în apropierea acestor indels este foarte importantă pentru cercetarea și explorarea metilării ADN-ului.
.