An easy-to-use, autorun package for DNA methylation analyses
An easy-to-use, autorun package for DNA methylation analyses
Aby wygodniej przeprowadzić analizę danych metylacji DNA, spakowaliśmy wszystkie funkcje w łatwy w użyciu, autorun package for DNA methylation analysis. Rysunek 1 przedstawia główne cechy BatMeth2: 1) BatMeth2 ma wydajną i dokładną wydajność wyrównywania. 2) BatMeth2 może obliczać poziom metylacji DNA (ML) poszczególnych miejsc cytozyny lub dowolnych regionów funkcjonalnych, takich jak całe chromosomy, regiony genów, elementy transpozycyjne (TE) itp. 3) Po zintegrowaniu różnych algorytmów statystycznych BatMeth2 może wykonywać analizę różnicową metylacji DNA dla dowolnego regionu, dowolnej liczby próbek wejściowych i wymagań użytkownika. 4) Poprzez integrację wizualizacji danych BS-Seq (dystrybucja metylacji DNA na chromosomach i genach) oraz anotacji metylacji różnicowej, BatMeth2 może wizualizować dane metylacji DNA w sposób bardziej przejrzysty. Podczas wykonywania narzędzia BatMeth2 generowany jest raport html dotyczący statystyk próbki. Szczegóły przykładowego raportu html są pokazane w http://htmlpreview.github.io/?https://github.com/GuoliangLi-HZAU/BatMeth2/blob/master/BatMeth2-Report/batmeth2.html.
Przebieg pracy narzędzia BatMeth2. Dwie duże strzałki oznaczają pliki wejściowe lub wyjściowe
BatMeth2 ma lepszą wydajność mapowania na symulowanych danych BS-Seq
W pierwszej kolejności oceniliśmy wszystkie alignery przy użyciu symulowanych zestawów danych (bez indeli) składających się z odczytów o długości 75 par zasad (bp), 100 bp i 150 bp oraz z różnymi współczynnikami konwersji bisulfitowej (od 0 do 100% z krokiem 10%). Te zestawy danych były symulowane z genomu ludzkiego (UCSC hg19) przy użyciu FASTX-mutate-tools , wgsim (v0.3.0) i symulatora w SAMtools (v1.1) , który dopuszcza 0,03% indeli, 1% błąd bazowy w całym genomie i maksymalnie dwa niedopasowania na odczyt. Symulowane odczyty mapowaliśmy do genomu referencyjnego, dopuszczając co najwyżej dwa niedopasowania. Ponieważ oryginalne pozycje symulowanych odczytów były znane, mogliśmy ocenić dokładność wszystkich programów poprzez porównanie ich wyników mapowania z oryginalnymi pozycjami.
Ocena wszystkich alignerów BS-Seq przy użyciu symulowanych zestawów danych o różnych długościach odczytów z FASTX i wgsim. Symulowane dane z różnymi współczynnikami konwersji bisulfitu są pokazane w różnych kształtach. Wyniki z różnych alignerów pokazane są za pomocą różnych kolorów symboli. Wyniki w pobliżu lewego górnego rogu w każdym panelu pokazują, że oprogramowanie uzyskało większą liczbę poprawnie zmapowanych odczytów i mniejszą liczbę błędnie zmapowanych odczytów. Wyniki z naszego alignera BatMeth2 są najlepsze w różnych symulowanych zestawach danych bisulfitowych
W skrócie, wyniki z wgsim-simulated indel-aberrant datasets pokazują, że BatMeth2 ma lepszą wydajność (1~2% lepszą niż drugi najlepszy aligner) niż inne metody podczas wyrównywania ogólnych symulowanych odczytów BS zawierających mieszaninę niedopasowań i indeli. Można zauważyć, że wraz ze wzrostem współczynnika konwersji BS, dokładność wyrównania wszystkich programów ulega zmniejszeniu. W tych różnych warunkach BatMeth2 radzi sobie lepiej.
BatMeth2 ma lepszą wydajność mapowania na rzeczywistych danych BS-Seq
Aby przetestować wydajność BatMeth2 na rzeczywistych zbiorach danych BS-Seq, pobraliśmy zbiory danych BS-Seq o sparowanych końcach i losowo wyodrębniliśmy 1 milion odczytów o długości 2 × 90 bp o sparowanych końcach z SRA SRR847318, 1 milion 2 × 101 bp sparowanych odczytów z SRA SRR1035722 i 1 milion 2 × 125 bp sparowanych odczytów z SRA SRR3503136 do celów oceny. Ponieważ te zestawy danych pochodzą ze zdrowych linii komórkowych lub tkanek, oczekuje się, że będą zawierać niewielką liczbę zmian strukturalnych. W związku z tym, wyrównaliśmy te rzeczywiste dane przy użyciu odczytów single-end z zestawów danych parowanych i oceniliśmy wskaźniki zgodnego i niezgodnego mapowania z parowanych wyrównań, aby oszacować wskaźniki poprawnego i niepoprawnego wyrównania. Ponieważ wielkość insertu sparowanych odczytów wynosiła około 500 bp, para partnerskich odczytów mogła być uznana za zgodną, jeśli została zmapowana w nominalnej odległości 500 bp; w przeciwnym razie para partnerskich odczytów mogła być uznana za niezgodną. Podobnie jak w przypadku danych symulowanych, BatMeth2 odnotował więcej zgodnych i mniej niezgodnych dopasowań w rzeczywistych zbiorach danych w szerokim zakresie wyników jakości map, jak pokazano na Rys. 3.
Concordance and discordance rates of alignments on real paired-end reads from different aligners. Skumulowana liczba zgodnych i niezgodnych dopasowań od wysokiej do niskiej jakości mapowania dla rzeczywistych odczytów sekwencjonowania bisulfitowego. Jest tylko jeden punkt dla BSmap i alignerów opartych na bowtie oddzielnie, ponieważ te alignery nie mają wyniku jakości map. Bismark-bowtie2L15 oznacza wyrównanie bowtie2 o długości ziarna 15
W tabeli 1 pokazano ponadto względne czasy działania programów. BatMeth2 z domyślnymi ustawieniami działał szybciej niż większość opublikowanych alignerów i był porównywalny z BWA-meth i BatMeth. Bismark2 (z Bowtie2 jako podstawową metodą mapowania), BS Seeker2 i BSmooth wymagają dłuższych czasów działania.
Wywoływanie metylacji DNA
Aby ocenić dokładność wywoływania metylacji DNA wśród różnych programów, pobraliśmy z ENCODE (Encyclopedia of DNA Elements) dane 450 K bead chip z linii komórkowej IMR90. Pobraliśmy również dane z sekwencjonowania bisulfitowego całego genomu (WGBS-Seq) linii komórkowej IMR90 z ENCODE (42.6 Gbases). Dla każdego programu wyrównaliśmy odczyty WGBS-Seq i obliczyliśmy poziom metylacji DNA. Następnie porównaliśmy wyniki z ML w tych samych miejscach w danych 450 K Bead Chip. Gdy różnica między ML DNA z danych WGBS-Seq przez oprogramowanie a ML z danych 450 K Bead Chip była mniejsza niż 0,2, wynik wywołania określano jako prawidłowy; w przeciwnym razie uznawano go za nieprawidłowy.
Wyniki przedstawiono w Tabeli 2. Nakładanie się prawidłowych wyników wszystkich programów przedstawiono w pliku dodatkowym 1: Figure S2. Widzimy, że BatMeth2 i Biscuit mają podobne wyniki, które są lepsze niż w przypadku pozostałych programów. Podsumowując, BatMeth2 poprawia dokładność zarówno wyrównywania odczytów BS, jak i wywoływania ML DNA.
BatMeth2 wyrównuje odczyty BS, pozwalając jednocześnie na występowanie indeli o zmiennej długości
Choroby nowotworowe zawierają znacznie więcej indeli niż zdrowe komórki. Dlatego, aby sprawdzić, czy BatMeth2 może wyrównywać odczyty BS z wcięciami o różnej długości, pobraliśmy z ENCODE dane WGBS (75 Gbaz) i 450 K Bead Chip z HepG2 (rak wątrobowokomórkowy, nowotworowa linia komórkowa). Sprawdziliśmy rozkład długości indeli w odczytach po wyrównaniu danych HepG2 WGBS-Seq. Plik dodatkowy 1: Figura S3A pokazuje, że długości wykrytych indeli były głównie dystrybuowane w zakresie 1 bp~ 5 bp, a najdłuższy indel miał długość 40 bp. Według naszych statystyk, 2,3% odczytów alignmentu zawierało indele. Na podstawie tych wyników wiemy, że BatMeth2 może wyrównywać odczyty z indelami o różnej długości.
Następnie sprawdziliśmy wpływ wykrywania indeli na wywoływanie metylacji DNA. Dla BatMeth2 uruchomiliśmy dwie opcje na danych HepG2: z i bez detekcji indeli (tj. ustawiliśmy parametr -I w BatMeth2). Uruchomiliśmy również Bismark na danych WGBS-Seq z HepG2 jako referencję dla wywoływania metylacji DNA z detekcją indeli, ponieważ Bismark nie posiada funkcji wywoływania indeli. Porównaliśmy wywołanie metylacji DNA w BatMeth2 i Bismark z wywołaniem z danych 450 K Bead Chip. Wyniki są pokazane w pliku dodatkowym 1: Figure S3B, gdzie „BatMeth2-noIndel” odpowiada BatMeth2 bez wykrywania indeli. Widzimy, że przy braku detekcji indeli wynik BatMeth2 był tylko nieznacznie lepszy od wyniku Bismarka (z Bowtie1 jako podstawową metodą mapowania). Wynik BatMeth2 z detekcją indeli był znacznie lepszy. Ponadto można zauważyć, że BatMeth2 może wykryć więcej miejsc metylacji DNA niż BatMeth2-noIndel i Bismark (Bowtie 1). Aby zrozumieć, dlaczego wydajność BatMeth2 w wykrywaniu indeli jest lepsza, zdefiniowaliśmy miejsca metylacji wywołane przez BatMeth2 jako Wynik A, podczas gdy miejsca metylacji wywołane przez BatMeth2-noIndel i Bismark zostały zdefiniowane jako Wynik B. Zaobserwowaliśmy, że mclA obejmuje 23 853 miejsca metylacji DNA i 15 048 (63%) z 23 853 miejsc objętych wyrównaniem odczytów indelowych wywołanych przez BatMeth2 z detekcją indeli (patrz plik dodatkowy 1: Rysunek S3C). Ponadto stwierdziliśmy, że wskaźniki indel w Wyniku A i Wyniku B wynosiły odpowiednio tylko 5 i 0%. Stąd wniosek, że dokładne wykrywanie indeli może poprawić wywoływanie metylacji DNA.
Wizualizacja danych metylacji DNA
BatMeth2 dostarcza narzędzi do wizualizacji danych metylacji. Aby zilustrować funkcje wizualizacyjne BatMeth2, pobraliśmy (1) 117 Gbaz odczytów single-end z ludzkiej linii komórkowej H9, (2) 105,2 Gbaz odczytów single-end z ludzkiej linii komórkowej IMR90 i (3) 12,6 Gbaz sparowanych odczytów z dzikiego ryżu. Po pierwsze, BatMeth2 może wizualizować gęstość metylacji cytozyny na poziomie chromosomów. Kropki na Rys. 4a reprezentują przesuwane okno o długości 100 kb z krokiem 50 kb. Aby umożliwić przeglądanie ML w poszczególnych miejscach CpG lub nie-CpG w przeglądarce genomu, udostępniamy również pliki w formatach bed i bigWig (Rys. 4b). Porównując z gęstością genów i TE, zaobserwowaliśmy, że ML była skorelowana z gęstością TE i antykorelowana z gęstością genów (Rys. 4c). Tendencja ta została wcześniej zaobserwowana u ryżu .
Wizualizacja poziomów metylacji w skali chromosomowej. a Gęstość metylo-cytozyny w ludzkim chromosomie 10. Kropki reprezentują poziomy metylacji w przesuwanych oknach 100Kb z krokiem 50Kb. Czerwone kropki odnoszą się do poziomów metylacji w nici plusowej, a niebieskie do poziomów metylacji w nici minusowej. b Przykład rozkładu poziomów metylacji DNA i różnie zmetylowanych regionów (DMRs) pomiędzy liniami komórkowymi H9 i IMR90 w ludzkim chromosomie 10. c Gęstość genów, elementów transpozonowych (TEs) i poziom metylacji DNA w całym genomie ryżu. Panel A to wyniki wygenerowane z Batmeth2. Panel B to wyniki wizualizacji z przeglądarki UCSC, z plikami BED z Batmeth2
Po drugie, BatMeth2 może wizualizować ML genów. Dokładniej, BatMeth2 może wizualizować MLs 2 kb przed genem, w miejscu startu transkrypcji (TSS), w ciele genu, w miejscu końca transkrypcji (TES) i 2 kb za ciałem genu. Porównując regiony upstream, body i downstream, na Rys. 5a widać, że ML DNA w korpusie genu jest wyższe niż w regionie promotorowym. Porównując wszystkie pięć regionów, widać, że w regionie TSS znajduje się dolina (Rys. 5b). BatMeth2 może również obliczać profile ML wokół intronów, eksonów, regionów intergenicznych i TE (plik dodatkowy 1: Rysunek S4). Dodatkowo, BatMeth2 może dostarczyć mapę cieplną wielu genów według regionów genowych dla wygodnego porównania ogólnych ML genów w różnych próbkach (Rys. 5c).
Wizualizacja metylacji DNA w różnych kontekstach. a Poziomy metylacji DNA w regionach 2Kb powyżej genów, ciał genów, 2Kb poniżej ciał genów. b Profil agregacji poziomów metylacji DNA w obrębie genów. c Mapa cieplna wszystkich genów w regionach 2Kb powyżej genów, ciał genów, 2Kb poniżej ciał genów
Po trzecie, BatMeth2 może wizualizować dystrybucję metylacji DNA. Plik dodatkowy 1: Figura S5A pokazuje rozkłady metylacji DNA w liniach komórkowych H9 i IMR90. Na rysunku, DNA ML jest podzielony na pięć kategorii: metylowany (M: > 80%), pośredni między częściowo metylowanym i metylowanym (Mh: 60-80%), częściowo metylowany (H: 40-60%), pośredni między niemetylowanym i częściowo metylowanym (hU: 20-40%) i niemetylowany (U: < 20%). Jak pokazano w pliku dodatkowym 1: Figura S5A, ML był wyższy w linii komórkowej H9 w kategorii M niż w linii komórkowej IMR90, szczególnie w kontekście CpG. W kontekście sekwencji CH, metylacja CpG jest formą dominującą, ale znaczna część metylowanych cytozyn znajduje się w miejscach CpA, podczas gdy ML wynosi mniej niż 40%, szczególnie w linii komórkowej H9 (plik dodatkowy 1: Rysunek S5B).
Po czwarte, BatMeth2 może analizować korelację pomiędzy poziomem ekspresji genu a ML DNA promotora genu. Zilustrowaliśmy tę cechę wykorzystując linie komórkowe H9 i IMR90. Poziomy ekspresji genów w H9 lub IMR90 zostały podzielone na różne kategorie. Jak pokazano w pliku dodatkowym 1: Figura S5C, wysoko wyrażone geny wykazywały niższe ML w ich regionach promotorowych. Ponadto, podzieliliśmy MLs promotorów genów na pięć kategorii. Wynik w pliku dodatkowym 1: Figura S5D pokazuje, że geny, których promotory miały wyższe wartości ML, wykazywały niższe poziomy ekspresji. Znana jest ujemna korelacja pomiędzy ekspresją genów ssaków a metylacją DNA promotora. Analiza ta dodatkowo wskazuje na dokładność BatMeth2.
Znajdowanie różnie metylowanych cytozyn i regionów (DMCs/DMRs)
Identyfikacja różnie metylowanych cytozyn (DMCs) i różnie metylowanych regionów (DMRs) jest jednym z głównych celów w analizie danych metylacyjnych. Chociaż badacze są czasami zainteresowani korelacją pojedynczych miejsc cytozyny z fenotypem, DMRs są bardzo ważnymi cechami .
Wcześniejsze badania BS-Seq profilowały komórki bez zbierania replik. Dla takich zestawów danych użyliśmy dokładnego testu Fishera, aby wyróżnić różnie metylowane cytozyny (DMCs). Dla zbiorów danych BS-Seq z replikacjami, najbardziej naturalnym modelem statystycznym do określania DMCs jest rozkład beta-binomialny. Wiemy, że istnieje wiele programów umożliwiających analizę danych dotyczących różnicowej metylacji DNA, takich jak methylKit (program do analizy różnicowej wymagający replikacji biologicznych) i Methy-Pipe (program do analizy różnicowej bez duplikacji biologicznych). Nie jest jednak dostępny żaden kompleksowy pakiet obejmujący zarówno mapowanie, jak i analizę różnicową metylacji. Stworzyliśmy więc pakiet, który integruje mapowanie z analizą różnicową. Aby ułatwić identyfikację DMR z danych bisulfitowych bez replik, zintegrowaliśmy dokładny test Fishera w celu wykonania testu hipotezy. Gdy próbka ma dwie lub więcej replik, używamy rozkładu beta-binomialnego do przeprowadzenia analizy metylacji różnicowej. Dostarczamy również pliki bed lub bigWig dla listy DMR. DMR mogą być wizualizowane w przeglądarce genomu (Rys. 4b) z wygenerowanymi plikami bed lub bigWig.
Jako ilustrację, Rys. 6a pokazuje liczbę DMCs i regionów w linii komórkowej IMR90 i w linii komórkowej H9, wykrytych przez BatMeth2 (wartość p< 0.05, meth.diff > = 0.6). BatMeth2 może wizualizować, czy CpGs i DMCs są wzbogacone w niektórych regionach, takich jak regiony genów, CDS, intronów, intergenów, UTR, TE, LTR, LINE i SINE. Rysunek 6b wizualizuje proporcje DMC w różnych regionach genomu. Poza regionami intergenicznymi, nie zaobserwowaliśmy wzbogacenia DMC w żadnych regionach.
Analiza metylacji różnicowej. a Wyniki analizy różnie zmetylowanych regionów (DMRs), różnie zmetylowanych genów (DMGs) i różnie zmetylowanych promotorów (DMPs) pomiędzy liniami komórkowymi H9 i IMR90. b Adnotacja różnie metylowanych cytozyn (DMC) względem różnych właściwości genomowych i elementów powtarzających się. c DMPs zawierają specyficzne dla H9 lub IMR90 indele (pomarańczowe) zajmują znaczną część we wszystkich DMPs (DNA Methylation differential Promoters)
Znaczna część różnie metylowanych promotorów (DMPs) zawiera indele
Wiemy, że indele i metylacja DNA odgrywają ważną rolę w rozwoju tkanek i chorobach. Tutaj badamy związek między różnie metylowanymi promotorami (DMPs) i indelami. Badanie to przeprowadziliśmy wykorzystując odczyty BS-Seq w liniach komórkowych IMR90 i H9. Najpierw wyrównaliśmy odczyty BS-Seq przy użyciu BatMeth2, a następnie wywołaliśmy indeksy przy użyciu narzędzi BisSNP i GATK. Następnie zdefiniowaliśmy indele występujące tylko w H9 lub IMR90 jako indele specyficzne dla linii komórkowej.
Potem wykryliśmy 1384 DMP pomiędzy H9 i IMR90 za pomocą BatMeth2 (wartość p< 0,05, meth.diff > = 0,6). W sumie 236 (17%) spośród wszystkich powyższych DMP zawiera indele, jak pokazano na Rys. 6c. Krótko mówiąc, znaczna część DMP zawiera indele. Dlatego dokładne wyrównanie odczytów BS-Seq w pobliżu tych indeli jest bardzo ważne dla badań i eksploracji metylacji DNA.