Integrovaný balíček pro bisulfitovou analýzu metylačních dat DNA s mapováním citlivým na indel

Snadno použitelný, autorunový balíček pro metylační analýzy DNA

Pro pohodlnější dokončení metylační analýzy DNA jsme všechny funkce zařadili do snadno použitelného, autorunového balíčku pro metylační analýzy DNA. Obrázek 1 ukazuje hlavní funkce balíčku BatMeth2: 1) BatMeth2 má efektivní a přesný výkon při zarovnávání. 2) BatMeth2 dokáže vypočítat úroveň metylace DNA (ML) jednotlivých cytosinových míst nebo libovolných funkčních oblastí, jako jsou celé chromozomy, genové oblasti, transponovatelné elementy (TE) atd. 3) Po integraci různých statistických algoritmů může BatMeth2 provádět diferenciální analýzu metylace DNA pro libovolnou oblast, libovolný počet vstupních vzorků a požadavky uživatele. 4) Díky integraci vizualizace dat BS-Seq (rozložení metylace DNA na chromozomech a genech) a anotace diferenciální metylace může BatMeth2 přehledněji vizualizovat data o metylaci DNA. Během provádění nástroje BatMeth2 je generován html report pro statistiku vzorku. Podrobnosti o ukázkovém html reportu jsou uvedeny v http://htmlpreview.github.io/?https://github.com/GuoliangLi-HZAU/BatMeth2/blob/master/BatMeth2-Report/batmeth2.html.

BatMeth2 má lepší mapovací výkon na simulovaných datech BS-Seq

Nejprve jsme všechny alignnery vyhodnotili pomocí simulovaných souborů dat (bez indelů) sestávajících ze čtení o 75 párech bází (bp), 100 bp a 150 bp a s různou mírou bisulfitové konverze (od 0 do 100 % s krokem 10 %). Tyto soubory dat byly simulovány z lidského genomu (UCSC hg19) pomocí FASTX-mutate-tools , wgsim (v0.3.0) a simulátoru v SAMtools (v1.1) , který umožňuje 0,03 % indelů, 1 % chybovost bází v celém genomu a maximálně dvě neshody na čtení. Simulovaná čtení jsme mapovali na referenční genom, přičemž jsme připustili maximálně dvě neshody. Protože původní pozice simulovaných čtení byly známy, mohli jsme vyhodnotit přesnost všech programů porovnáním jejich mapovacích výstupů s původními pozicemi.

Pro porovnání výkonnosti jednotlivých programů bylo sekvenační čtení s indely považováno za správně namapované, pokud platily následující podmínky: 1) čtení bylo jednoznačně namapováno na stejné vlákno, ze kterého bylo simulováno, a kvalita mapování byla větší než 0; 2) uváděná počáteční pozice zarovnaného čtení byla do deseti párů bází původní počáteční pozice simulovaného čtení; 3) výsledky mapování měly podobné indely nebo neshody jako simulované čtení. Pokud byla některá z těchto podmínek porušena, bylo čtení považováno za chybně namapované. Vzhledem k tomu, že BatMeth2 umožňuje jednu mezeru ve výsevní oblasti, dokáže najít výsevní místa obsahující indely s vysokou přesností a může se vyhnout neshodným místům, která by způsobila chybné zarovnání čtení obsahujících indely. Výsledky na obr. 2 ukazují, že BatMeth2 dosáhl největšího počtu správně zarovnaných čtení a nejnižšího počtu nesprávně zarovnaných čtení ve všech testovacích sadách dat při různých rychlostech bisulfitové konverze.

Ve zkratce, výsledky z datových sad wgsim simulovaných indely a aberacemi ukazují, že BatMeth2 má lepší výkon (o 1~2 % lepší než druhý nejlepší aligner) než ostatní metody při zarovnávání obecných simulovaných BS čtení obsahujících směs neshod a indelů. Vidíme, že se zvýšenou mírou konverze BS se snižuje přesnost zarovnání všech softwarů. Za těchto rozdílných podmínek dosahuje BatMeth2 lepších výsledků.

BatMeth2 má lepší výkonnost při mapování na skutečných datech BS-Seq

Pro otestování výkonnosti BatMeth2 na skutečných sadách dat BS-Seq jsme stáhli sady dat BS-Seq s párovým koncem a náhodně extrahovali 1 milion 2 × 90 bp čtení s párovým koncem ze SRA SRR847318, 1 milion 2 × 101 bp párových čtení ze SRA SRR1035722 a 1 milion 2 × 125 bp párových čtení ze SRA SRR3503136 pro účely hodnocení. Protože tyto soubory dat pocházejí ze zdravých buněčných linií nebo tkání, očekává se, že budou obsahovat nízký počet strukturních variací. Proto jsme tato skutečná data zarovnali pomocí jednokoncových čtení ze sad dat s párovým koncem a vyhodnotili jsme míru souhlasného a nesouhlasného mapování z párových zarovnání, abychom odhadli míru správného a nesprávného zarovnání. Protože velikost inzertu čtení s párovým koncem byla přibližně 500 bp, dvojice partnerských čtení mohla být považována za konkordantní, pokud byla mapována v nominální vzdálenosti 500 bp; v opačném případě mohla být dvojice partnerských čtení považována za diskordantní. Podobně jako naše výsledky se simulovanými daty zaznamenal BatMeth2 více konkordantních a méně diskordantních zarovnání na reálných souborech dat ve velkém rozsahu skóre kvality mapování, jak ukazuje obr. 3.

V tabulce 1 jsou navíc uvedeny relativní doby běhu programů. BatMeth2 s výchozím nastavením běžel rychleji než většina publikovaných zarovnávačů a byl srovnatelný s BWA-meth a BatMeth. Bismark2 (s Bowtie2 jako základní mapovací metodou), BS Seeker2 a BSmooth vyžadují delší dobu běhu.

Volání metylace DNA

Pro vyhodnocení přesnosti volání metylace DNA mezi různými softwary jsme z ENCODE (Encyclopedia of DNA Elements) stáhli data 450 K kuličkových čipů z buněčné linie IMR90. Stáhli jsme také data celogenomového bisulfitového sekvenování (WGBS-Seq) buněčné linie IMR90 z ENCODE (42,6 Gbází). Pro každý software jsme zarovnali čtení WGBS-Seq a vypočítali úroveň metylace DNA. Poté jsme výsledky porovnali s ML na stejných místech v datech 450 K Bead Chip. Pokud byl rozdíl mezi ML DNA z dat WGBS-Seq podle softwaru a ML z 450 K Bead Chip menší než 0,2, byl výsledek volání definován jako správný; v opačném případě byl považován za nesprávný.

Výsledky jsou uvedeny v tabulce 2. Překrývání správných výsledků všech softwarů je uvedeno v doplňkovém souboru 1: Obrázek S2. Vidíme, že BatMeth2 a Biscuit mají podobné výsledky, které jsou lepší než výsledky ostatních softwarů. Závěrem lze říci, že BatMeth2 zlepšuje přesnost zarovnání čtení BS i volání ML DNA.

BatMeth2 zarovnává čtení BS a zároveň umožňuje indely s proměnlivou délkou

Rakovina obsahuje výrazně vyšší podíl indelů než zdravé buňky. Proto jsme pro ověření, zda BatMeth2 dokáže zarovnat BS čtení s indely různé délky, stáhli z ENCODE data WGBS (75 Gbases) a 450 K Bead Chip data z HepG2 (jaterní hepatocelulární karcinom, nádorová buněčná linie). Po zarovnání dat HepG2 WGBS-Seq jsme zkontrolovali rozložení délky indelů ve čteních. Doplňkový soubor 1: Obrázek S3A ukazuje, že délky zjištěných indelů byly rozloženy převážně v rozmezí 1 bp ~ 5 bp a nejdelší indel měl délku 40 bp. Podle našich statistik obsahovalo indely 2,3 % zarovnávacích čtení. Z těchto výsledků vyplývá, že BatMeth2 dokáže zarovnat čtení s indely různých délek.

Dále jsme testovali vliv detekce indelů na volání metylace DNA. Pro BatMeth2 jsme na datech HepG2 spustili dvě možnosti: s detekcí indelů a bez ní (tj. nastavení parametru -I v BatMeth2). Na datech WGBS-Seq z HepG2 jsme také spustili Bismark jako referenci pro volání metylace DNA s detekcí indelů, protože Bismark nemá funkci pro volání indelů. Porovnali jsme volání metylace DNA v BatMeth2 a Bismark s voláním z dat 450 K Bead Chip. Výsledky jsou uvedeny v doplňkovém souboru 1: S3B, kde „BatMeth2-noIndel“ odpovídá BatMeth2 bez detekce indelů. Vidíme, že při absenci detekce indelů byl výsledek metody BatMeth2 jen o málo lepší než výsledek metody Bismark (se základní mapovací metodou Bowtie1). Výsledek metody BatMeth2 s detekcí indelů byl výrazně lepší. Dále můžeme vidět, že BatMeth2 dokáže detekovat více míst metylace DNA než BatMeth2-noIndel a Bismark (Bowtie 1). Abychom pochopili, proč je výkon BatMeth2 s detekcí indelů lepší, definovali jsme metylační místa vyvolaná BatMeth2 jako výsledek A, zatímco metylační místa vyvolaná BatMeth2-noIndel a Bismarkem byla definována jako výsledek B. Pak jsme nechali mclA jako metylační místa, která se objevují ve výsledku A, ale ne ve výsledku B. Zjistili jsme, že mclA zahrnuje 23 853 metylačních míst DNA a 15 048 (63 %) z 23 853 míst pokrývají zarovnání čtení indelů vyvolaných BatMeth2 s detekcí indelů (viz doplňkový soubor 1: obrázek S3C). Kromě toho jsme zjistili, že podíl indelů ve výsledku A byl pouze 5 % a ve výsledku B 0 %. Proto jsme dospěli k závěru, že přesná detekce indelů může zlepšit volání metylace DNA.

Vizualizace metylačních dat DNA

BatMeth2 poskytuje nástroje pro vizualizaci metylačních dat. Pro ilustraci vizualizačních funkcí BatMeth2 jsme stáhli (1) 117 Gbase jednokoncových čtení z lidské buněčné linie H9, (2) 105,2 Gbases jednokoncových čtení z lidské buněčné linie IMR90 a (3) 12,6 Gbases párových čtení z divokého typu rýže. Za prvé, BatMeth2 dokáže vizualizovat hustotu metylace cytosinu na úrovni chromozomů. Tečky na obr. 4a představují posuvné okno o velikosti 100 kb s krokem 50 kb. Aby bylo možné prohlížet ML na jednotlivých CpG nebo ne-CpG místech v prohlížeči genomu, poskytujeme také soubory ve formátech bed a bigWig (obr. 4b). Porovnáním s hustotou genů a TE jsme zjistili, že ML koreluje s hustotou TE a je antikorelován s hustotou genů (obr. 4c). Tato tendence byla již dříve pozorována u rýže .

Druhé, BatMeth2 dokáže vizualizovat ML genů. Přesněji řečeno, BatMeth2 dokáže vizualizovat ML 2 kb před genem, v místě začátku transkripce (TSS), v těle genu, v místě konce transkripce (TES) a 2 kb za tělem genu. Při porovnání oblastí před, v těle a za genem je na obr. 5a vidět, že ML DNA v těle genu je vyšší než ML v oblasti promotoru. Při porovnání všech pěti oblastí je zřejmé, že v oblasti TSS je údolí (obr. 5b). Program BatMeth2 dokáže také vypočítat profily ML v okolí intronů, exonů, intergenních oblastí a TE (Doplňkový soubor 1: Figure S4). Kromě toho může BatMeth2 poskytnout tepelnou mapu více genů podle genových oblastí pro pohodlné porovnání celkových ML genů různých vzorků (obr. 5c).

Zatřetí, BatMeth2 může vizualizovat distribuci metylace DNA. Additional file 1: Obrázek S5A ukazuje rozložení metylace DNA v buněčných liniích H9 a IMR90. Na obrázku je DNA ML rozdělena do pěti kategorií: metylovaná (M: > 80 %), intermediární mezi částečně metylovanou a metylovanou (Mh: 60-80 %), částečně metylovaná (H: 40-60 %), intermediární mezi nemetylovanou a částečně metylovanou (hU: 20-40 %) a nemetylovaná (U: < 20 %). Jak je uvedeno v doplňkovém souboru 1: S5A, ML byl vyšší u buněčné linie H9 v kategorii M než u buněčné linie IMR90, zejména v kontextu CpG. V kontextu sekvence CH je převažující formou metylace CpG, ale významná část metylovaných cytosinů se nachází v místech CpA, zatímco ML je méně než 40 %, zejména u buněčné linie H9 (Additional file 1: Figure S5B).

Za čtvrté, BatMeth2 může analyzovat korelaci mezi úrovní genové exprese a ML DNA promotoru genu. Tuto funkci jsme ilustrovali pomocí buněčných linií H9 a IMR90. Úrovně exprese genů v H9 nebo IMR90 byly rozděleny do různých kategorií. Jak je uvedeno v doplňkovém souboru 1: S5C, vysoce exprimované geny vykazovaly nižší ML ve svých promotorových oblastech. Dále jsme MLs promotorů genů rozdělili do pěti kategorií. Výsledek v Doplňkovém souboru 1: S5D ukazuje, že geny s promotory s vyššími hodnotami ML vykazovaly nižší úrovně exprese. Negativní korelace mezi expresí savčích genů a promotorovou metylací DNA je známá . Tato analýza dále naznačuje přesnost metody BatMeth2.

Nalezení diferenciálně metylovaných cytosinů a oblastí (DMCs/DMRs)

Identifikace diferenciálně metylovaných cytosinů (DMCs) a diferenciálně metylovaných oblastí (DMRs) je jedním z hlavních cílů analýzy metylačních dat. Ačkoli se výzkumníci občas zajímají o korelaci jednotlivých cytosinových míst s fenotypem , DMR jsou velmi důležitými rysy .

Počáteční studie BS-Seq profilovaly buňky bez sběru replikátů. U takových souborů dat jsme použili Fisherův přesný test k rozlišení diferenciálně metylovaných cytosinů (DMC). Pro soubory dat BS-Seq s replikáty je nejpřirozenějším statistickým modelem pro označení DMCs beta-binomické rozdělení . Víme, že analýzu dat diferenciální metylace DNA může provádět řada softwarových programů, například methylKit (program pro diferenciální analýzu, který vyžaduje biologické replikáty) a Methy-Pipe (program pro diferenciální analýzu bez biologických replikátů). Není však k dispozici žádný komplexní balík zahrnující jak mapování, tak diferenciální metylační analýzu. Proto jsme vyvinuli balíček, který integruje mapování s diferenční analýzou. Abychom usnadnili identifikaci DMR z bisulfitových dat bez replikátů, integrovali jsme Fisherův přesný test k provedení testu hypotéz. Pokud má vzorek dvě nebo více replikátů, používáme k provedení diferenciální metylační analýzy beta-binomické rozdělení. Poskytujeme také soubory bed nebo bigWig pro seznam DMR. DMR lze vizualizovat v prohlížeči genomu (obr. 4b) pomocí vygenerovaných souborů bed nebo bigWig.

Pro ilustraci jsou na obr. 6a uvedeny počty DMC a oblastí v buněčné linii IMR90 a v buněčné linii H9, jak byly zjištěny pomocí BatMeth2 (hodnota p< 0,05, meth.diff > = 0,6). BatMeth2 dokáže zobrazit, zda jsou CpGs a DMCs obohaceny v některých oblastech, jako jsou genové, CDS, intronové, intergenové, UTR, TE, LTR, LINE a SINE oblasti. Obrázek 6b vizualizuje podíly DMCs v různých genomových oblastech. Kromě intergenových oblastí jsme obohacení DMC nepozorovali v žádné oblasti.

Podstatný podíl diferenciálně metylovaných promotorů (DMP) obsahuje indely

Víme, že indely a metylace DNA hrají důležitou roli ve vývoji tkání a onemocnění . Zde zkoumáme vztah mezi diferenciálně metylovanými promotory (DMP) a indely. Tuto studii jsme provedli pomocí čtení BS-Seq v buněčných liniích IMR90 a H9. Čtení BS-Seq jsme nejprve zarovnali pomocí nástroje BatMeth2; poté jsme indely vyvolali pomocí nástrojů BisSNP a GATK. Následně jsme indely, které se vyskytují pouze v H9 nebo IMR90, definovali jako indely specifické pro buněčné linie.

Poté jsme pomocí BatMeth2 detekovali 1384 DMP mezi H9 a IMR90 (hodnota p< 0,05, meth.diff > = 0,6). Celkem 236 (17 %) ze všech výše uvedených DMP obsahuje indely, jak ukazuje obr. 6c. Stručně řečeno, značná část DMP obsahuje indely. Proto je přesné zarovnání čtení BS-Seq v blízkosti těchto indelů velmi důležité pro výzkum a zkoumání metylace DNA.

.