Egy könnyen használható, autorun csomag DNS-metilációs elemzésekhez
A DNS-metilációs adatok elemzésének kényelmesebb elvégzése érdekében az összes funkciót egy könnyen használható, autorun csomagba csomagoltuk a DNS-metilációs elemzéshez. Az 1. ábra mutatja a BatMeth2 fő funkcióit: 1) A BatMeth2 hatékony és pontos igazítási teljesítményt nyújt. 2) A BatMeth2 képes kiszámítani az egyes citozin helyek vagy bármely funkcionális régió, például teljes kromoszómák, génrégiók, transzponálható elemek (TE-k) stb. DNS-metilációs szintjét (ML). 3) A különböző statisztikai algoritmusok integrálása után a BatMeth2 képes differenciális DNS-metilációs elemzést végezni bármilyen régióra, bármilyen számú bemeneti mintára és felhasználói igényre. 4) A BS-Seq adatok vizualizációjának (DNS-metiláció eloszlása a kromoszómákon és a géneken) és a differenciális metilációs annotációnak az integrálásával a BatMeth2 képes a DNS-metilációs adatokat világosabban megjeleníteni. A BatMeth2 eszköz végrehajtása során a minta statisztikáiról html jelentés készül. A html-jelentés mintájának részletei a http://htmlpreview.github.io/?https://github.com/GuoliangLi-HZAU/BatMeth2/blob/master/BatMeth2-Report/batmeth2.html.
A BatMeth2 munkafolyamata. A két nagy nyíl a bemeneti vagy kimeneti fájlokat jelenti
A BatMeth2 jobb leképezési teljesítményt nyújt szimulált BS-Seq adatokon
Először az összes illesztőprogramot szimulált (indelek nélküli) adatkészletekkel értékeltük, amelyek 75 bázispár (bp), 100 bp és 150 bp méretű olvasatokból álltak, különböző biszulfit-konverziós arányokkal (0 és 100% között, 10%-os lépéssel). Ezeket az adatkészleteket a humán genomból (UCSC hg19) szimuláltuk a FASTX-mutate-tools , a wgsim (v0.3.0) és a SAMtools (v1.1) szimulátorának segítségével, amely 0,03%-os indeleket, 1%-os bázishibaarányt engedélyez a teljes genomban és olvasásonként legfeljebb két eltérést. A szimulált olvasatokat leképeztük a referencia genomra, legfeljebb két eltérést megengedve. Mivel a szimulált leolvasások eredeti pozíciói ismertek voltak, az összes program pontosságát értékelni tudtuk azáltal, hogy összehasonlítottuk a leképezési eredményeiket az eredeti pozíciókkal.
A különböző szoftverek teljesítményének összehasonlításához egy indeleket tartalmazó szekvenáló olvasatot akkor tekintettünk helyesen leképezettnek, ha a következő feltételek teljesültek: 1) az olvasatot egyértelműen ugyanarra a szálra képeztük le, amelyről szimuláltuk, és a leképezés minősége nagyobb volt, mint 0; 2) az összehangolt olvasat bejelentett kezdőpozíciója tíz bázispáron belül volt a szimulált olvasat eredeti kezdőpozíciójához képest; 3) a leképezési eredmények hasonló indeleket vagy mismatcheket mutattak, mint a szimulált olvasat. Ha e feltételek bármelyike sérült, a leolvasást hibásan leképezettnek tekintettük. Mivel a BatMeth2 egy hézagot enged meg a magterületben, nagy pontossággal képes megtalálni az indeleket tartalmazó maghelyeket, és el tudja kerülni a hibásan illeszkedő helyeket, amelyek az indeleket tartalmazó olvasatok téves illesztését eredményeznék. A 2. ábrán látható eredmények azt mutatják, hogy a BatMeth2 érte el a legnagyobb számú helyesen igazított olvasatot és a legkisebb számú helytelenül igazított olvasatot az összes tesztadatkészletben a különböző biszulfitkonverziós arányok mellett.
A BS-Seq összes igazítójának értékelése a FASTX és wgsim különböző olvasáshosszúságú szimulált adatkészleteivel. A szimulált adatok különböző biszulfit-konverziós rátákkal különböző alakzatokkal vannak ábrázolva. A különböző igazítók eredményei a szimbólumok különböző színeivel jelennek meg. Az egyes panelek bal felső sarka közelében lévő eredmények azt mutatják, hogy a szoftver több helyesen leképezett olvasatot és a helytelenül leképezett olvasatok alacsonyabb számát érte el. A BatMeth2 alignerünk eredményei a legjobbak a különböző szimulált biszulfitos adathalmazokban
Röviden, a wgsim-szimulált indel-aberrant adathalmazok eredményei azt mutatják, hogy a BatMeth2 jobb teljesítményt nyújt (1~2%-kal jobbat, mint a második legjobb aligner), mint a többi módszer a hibás illesztéseket és indeleket tartalmazó általános szimulált BS olvasatok illesztésekor. Láthatjuk, hogy a BS konverziós arány növekedésével az összes szoftver igazítási pontossága csökken. Ezekben a különböző körülmények között a BatMeth2 jobban teljesít.
A BatMeth2 jobb illesztési teljesítményt nyújt valós BS-Seq-adatokon
A BatMeth2 teljesítményének tesztelésére valós BS-Seq-adatkészleteken párosított végű BS-Seq-adatkészleteket töltöttünk le, és véletlenszerűen 1 millió 2 × 90 bp párosított végű olvasatot vontunk ki az SRA SRR847318-ból, 1 millió 2 × 101 bp páros végű olvasatot az SRA SRR1035722-ből és 1 millió 2 × 125 bp páros végű olvasatot az SRA SRR3503136-ból értékelés céljából. Mivel ezek az adatkészletek egészséges sejtvonalakból vagy szövetekből származnak, várhatóan kevés strukturális variációt tartalmaznak. Ezért ezeket a valós adatokat a páros végű adatkészletekből származó egyvégű leolvasások segítségével igazítottuk, és a helyes és helytelen igazítási arányok becsléséhez kiértékeltük a páros igazításokból származó egyező és nem egyező leképezési arányokat. Mivel a páros végű leolvasások inzertmérete körülbelül 500 bp volt, egy pár partnerolvasás akkor tekinthető konkordánsnak, ha a leképezésük 500 bp névleges távolságon belül történt; ellenkező esetben egy pár partnerolvasás diszkordánsnak tekinthető. A szimulált adatokkal kapott eredményeinkhez hasonlóan a BatMeth2 a valós adatkészleteken a térképminőségi pontszámok széles tartományában több konkordáns és kevesebb diszkordáns igazítást jelentett, amint azt a 3. ábra mutatja.
A különböző igazítókból származó valós páros végű olvasatokon végzett igazítások konkordancia és diszkordancia aránya. Az egyező és nem egyező igazítások kumulatív száma a magas és az alacsony leképezési minőség között valós biszulfit szekvenáló leolvasásokhoz. A BSmap és a bowtie-alapú alignerek esetében külön-külön csak egy pont van, mivel ezek az alignerek nem rendelkeznek térképminőségi pontszámmal. A Bismark-bowtie2L15 a bowtie2 igazítást jelenti 15
maghosszúsággal
Az 1. táblázat mutatja továbbá a programok relatív futási idejét. A BatMeth2 az alapértelmezett beállításokkal gyorsabban futott, mint a legtöbb publikált aligner, és összehasonlítható volt a BWA-meth és a BatMeth programokkal. A Bismark2 (Bowtie2-vel mint alapvető leképezési módszerrel), a BS Seeker2 és a BSmooth hosszabb futási időt igényel.
DNS-metilációs hívás
A DNS-metilációs hívás különböző szoftverek közötti pontosságának értékeléséhez az ENCODE (Encyclopedia of DNA Elements) adatbázisából töltöttük le az IMR90 sejtvonal 450 K bead chip adatait. Az ENCODE-ból letöltöttük az IMR90 sejtvonal teljes génállományú biszulfit-szekvenálási (WGBS-Seq) adatait is (42,6 Gbázist). Minden egyes szoftver esetében összehangoltuk a WGBS-Seq leolvasásokat, és kiszámítottuk a DNS-metiláció szintjét. Ezután összehasonlítottuk az eredményeket a 450 K Bead Chip adatainak azonos helyeken mért ML-értékeivel. Ha a szoftver által a WGBS-Seq adatokból származó DNS ML és a 450 K Bead Chip adataiból származó DNS ML közötti különbség kisebb volt, mint 0,2, a hívási eredményt helyesnek, ellenkező esetben helytelennek tekintettük.
Az eredményeket a 2. táblázat mutatja. Az összes szoftver helyes eredményeinek átfedését az 1. kiegészítő fájl mutatja: S2 ábra. Látható, hogy a BatMeth2 és a Biscuit hasonló teljesítményt nyújt, amely jobb, mint a többi szoftveré. Összefoglalva, a BatMeth2 javítja mind a BS-olvasások összehangolásának, mind a DNS ML megnevezésének pontosságát.
A BatMeth2 összehangolja a BS-olvasásokat, miközben figyelembe veszi a változó hosszúságú indeleket
A rákos sejtek jelentősen nagyobb arányban tartalmaznak indeleket, mint az egészséges sejtek. Ezért annak ellenőrzésére, hogy a BatMeth2 képes-e különböző hosszúságú indeleket tartalmazó BS olvasatok igazítására, letöltöttük az ENCODE-tól a WGBS adatokat (75 Gbázist) és a HepG2 (máj hepatocelluláris karcinóma, egy rákos sejtvonal) 450 K Bead Chip adatait. A HepG2 WGBS-Seq-adatok összehangolása után ellenőriztük az indel-hosszúság eloszlását a leolvasásokban. Kiegészítő fájl 1: S3A ábra azt mutatja, hogy a detektált indelek hossza főként az 1 bp~ 5 bp tartományban oszlott meg, és a leghosszabb indel 40 bp hosszúságú volt. Statisztikáink szerint az illesztési leolvasások 2,3%-a tartalmazott indelt. Ezekből az eredményekből tudjuk, hogy a BatMeth2 képes a különböző hosszúságú indeleket tartalmazó olvasatok igazítására.
A következőkben az indelek detektálásának hatását vizsgáltuk a DNS-metilációs hívásra. A BatMeth2 esetében két opciót futtattunk a HepG2-adatokon: indel-detektálással és anélkül (azaz a BatMeth2 -I paraméterének beállításával). A Bismarkot a HepG2 WGBS-Seq adatain is lefuttattuk referenciaként az indel-detektálással történő DNS-metiláció-meghíváshoz, mivel a Bismark nem rendelkezik indel-meghívó funkcióval. Összehasonlítottuk a BatMeth2 és a Bismark DNS-metilációs hívását a 450 K Bead Chip adatokból származó hívással. Az eredményeket az 1. kiegészítő fájlban mutatjuk be: S3B ábra, ahol a “BatMeth2-noIndel” megfelel a BatMeth2-nek indel-felismerés nélkül. Láthatjuk, hogy indel-detektálás hiányában a BatMeth2 eredménye csak kis mértékben volt jobb, mint a Bismarké (Bowtie1 alapvető térképezési módszerrel). A BatMeth2 eredménye indel detektálással lényegesen jobb volt. Továbbá látható, hogy a BatMeth2 több DNS-metilációs helyet képes detektálni, mint a BatMeth2-noIndel és a Bismark (Bowtie 1). Annak megértéséhez, hogy miért jobb a BatMeth2 teljesítménye az indel detektálással, a BatMeth2 által megnevezett metilációs helyeket A eredményként, míg a BatMeth2-noIndel és Bismark által megnevezett metilációs helyeket B eredményként definiáltuk. Ezután az mclA legyen az A eredményben megjelenő, de a B eredményben nem szereplő metilációs helyek. Megfigyeltük, hogy az mclA 23 853 DNS-metilációs helyet tartalmazott, és a BatMeth2 által indel-detektálással hívott indel olvasatok igazításai által lefedett 23 853 helyből 15 048-at (63%) (lásd Kiegészítő fájl 1: S3C ábra). Ezenkívül azt találtuk, hogy az A és B eredményben az indelek aránya mindössze 5, illetve 0% volt. Ezért arra a következtetésre jutottunk, hogy a pontos indel-detektálás javíthatja a DNS-metilációs hívást.
A DNS-metilációs adatok vizualizálása
A BatMeth2 eszközöket biztosít a metilációs adatok vizualizálásához. A BatMeth2 vizualizációs funkcióinak illusztrálására letöltöttünk (1) 117 Gbázist egyvégű olvasatokból a humán H9 sejtvonalból, (2) 105,2 Gbázist egyvégű olvasatokból a humán IMR90 sejtvonalból és (3) 12,6 Gbázist párosított olvasatokból a vad típusú rizsből. Először is, a BatMeth2 képes a citozin-metilációs sűrűséget kromoszómaszinten megjeleníteni. A 4a. ábrán a pontok egy 100 kb-os csúszóablakot képviselnek 50 kb-os lépéssel. Ahhoz, hogy az ML-t az egyes CpG- vagy nem-CpG-helyeken egy genomböngészőben meg lehessen tekinteni, a fájlokat bed és bigWig formátumban is megadjuk (4b. ábra). A gének és TE-k sűrűségével összehasonlítva megfigyeltük, hogy az ML korrelált a TE-sűrűséggel és antikorrelált a génsűrűséggel (4c. ábra). Ezt a tendenciát korábban rizsben is megfigyeltük .
A metilációs szintek vizualizációja a kromoszómaskálán. a A metil-citozin sűrűsége a humán 10-es kromoszómában. A pontok a metilációs szinteket ábrázolják 100Kb-os csúszóablakokban, 50Kb-os lépéssel. A piros pontok a plusz szál metilációs szintjére, a kék pontok pedig a mínusz szál metilációs szintjére utalnak. b Példa a DNS-metilációs szintek és a differenciálisan metilált régiók (DMR) eloszlására a H9 és IMR90 sejtvonalak között a 10-es humán kromoszómán. c A gének, transzpozon elemek (TE-k) sűrűsége és a DNS-metilációs szint a teljes rizs genomban. Az A panel a Batmeth2 segítségével generált eredményeket mutatja. A B panel az UCSC böngészőből származó vizualizációs eredmények, a Batmeth2
BED fájljaival>A BatMeth2 képes vizualizálni a gének ML-jét. Pontosabban, a BatMeth2 képes vizualizálni az ML-eket a géntől 2 kb-mal feljebb, a transzkripció kezdőhelyénél (TSS), a géntestben, a transzkripció véghelyénél (TES) és a géntesttől 2 kb-mal lejjebb. Az 5a. ábra az upstream, a test és a downstream régiók összehasonlításával azt mutatja, hogy a géntest DNS ML-je magasabb, mint a promóter régióé. Ha mind az öt régiót összehasonlítjuk, a TSS régióban nyilvánvalóan egy völgy található (5b. ábra). A BatMeth2 az intronok, exonok, intergénikus régiók és TE-k körüli ML-profilokat is ki tudja számítani (Additional file 1: S4 ábra). Ezenkívül a BatMeth2 több génről génrégiónként egy hőtérképet is tud nyújtani a különböző minták teljes gén ML-értékeinek kényelmes összehasonlításához (5c. ábra).
A DNS-metiláció vizualizációja különböző kontextusokban. a A DNS-metilációs szintek a gének, géntestek előtti 2Kb-os régiókban, géntestek utáni 2Kb-os régiókban. b A DNS-metilációs szintek aggregációs profilja a gének között. c Az összes gén hőtérképe a gének, géntestek előtti 2Kb-os régiókban, géntestek utáni 2Kb-os régiókban
Harmadszor, a BatMeth2 képes a DNS-metiláció eloszlásának vizualizálására. Additional file 1: S5A ábra a DNS-metiláció eloszlását mutatja a H9 és az IMR90 sejtvonalakban. Az ábrán a DNS ML öt kategóriára van felosztva: metilált (M: > 80%), köztes a részben metilált és a metilált között (Mh: 60-80%), részben metilált (H: 40-60%), köztes a nem metilált és a részben metilált között (hU: 20-40%), és nem metilált (U: < 20%). Amint az 1. kiegészítő fájlban látható: S5A ábra, az ML magasabb volt a H9 sejtvonalban az M kategóriában, mint az IMR90 sejtvonalban, különösen a CpG-kontextusban. A CH szekvencia-kontextusban a CpG-metiláció az uralkodó forma, de a metilált citozinok jelentős része megtalálható a CpA helyeken, míg az ML kevesebb, mint 40%, különösen a H9 sejtvonalban (Additional file 1: S5B ábra).
Negyedszer, a BatMeth2 képes elemezni a génexpressziós szint és a gén promóter DNS ML közötti korrelációt. Ezt a funkciót a H9 és az IMR90 sejtvonalak segítségével szemléltettük. A H9 vagy IMR90 gének expressziós szintjét különböző kategóriákba soroltuk. Amint az az 1. kiegészítő fájlban látható: S5C ábra, a magasan expresszált gének alacsonyabb ML-t mutattak a promóter régióikban. Továbbá a gének promótereinek ML-értékeit öt kategóriába soroltuk. Az eredmény az 1. kiegészítő fájlban látható: S5D ábra azt mutatja, hogy a magasabb ML-értékkel rendelkező promóterekkel rendelkező gének alacsonyabb expressziós szintet mutattak. Az emlős gének expressziója és a promóter DNS-metiláció közötti negatív korreláció ismert . Ez az elemzés tovább jelzi a BatMeth2 pontosságát.
Differenciálisan metilált citozinok és régiók (DMCs/DMRs)
A differenciálisan metilált citozinok (DMCs) és differenciálisan metilált régiók (DMRs) azonosítása a metilációs adatok elemzésének egyik fő célja. Bár a kutatókat esetenként az egyes citozinhelyek fenotípussal való korrelációja érdekli, a DMR-ek nagyon fontos jellemzők.
A korai BS-Seq vizsgálatokban a sejteket ismétlések gyűjtése nélkül profilozták. Az ilyen adathalmazok esetében Fisher egzakt tesztjét használtuk a differenciálisan metilált citozinok (DMC-k) megkülönböztetésére. Ismétlésekkel rendelkező BS-Seq-adatkészletek esetében a DMC-k megnevezésére a legtermészetesebb statisztikai modell a béta-binomiális eloszlás. Tudjuk, hogy számos szoftverprogram képes differenciális DNS-metilációs adatok elemzésére, például a methylKit (biológiai ismétléseket igénylő differenciális elemző program) és a Methy-Pipe (biológiai ismétlések nélküli differenciális elemző program). Nem áll azonban rendelkezésre olyan átfogó csomag, amely a térképezést és a differenciális metilációs elemzést is magában foglalja. Ezért kifejlesztettünk egy olyan csomagot, amely integrálja a térképezést és a differenciális elemzést. A DMR-ek azonosításának megkönnyítése érdekében a biszulfit adatokból, replikátumok nélkül, integráltuk a Fisher-féle egzakt tesztet a hipotézisvizsgálat elvégzéséhez. Ha egy minta két vagy több ismétléssel rendelkezik, a béta-binomiális eloszlást használjuk a differenciális metilációs elemzés elvégzéséhez. A DMR-ek listájához bed vagy bigWig fájlokat is biztosítunk. A DMR-ek a genom-böngészőben (4b. ábra) vizualizálhatók a generált bed- vagy bigWig-fájlokkal.
A 6a. ábra szemléltetésképpen mutatja a DMC-k és régiók számát az IMR90 sejtvonalban és a H9 sejtvonalban, ahogyan azt a BatMeth2 kimutatta (p-érték< 0,05, meth.diff > = 0,6). A BatMeth2 képes láthatóvá tenni, hogy a CpG-k és DMC-k bizonyos régiókban, például gén, CDS, intron, intergenikus, UTR, TE, LTR, LINE és SINE régiókban feldúsulnak-e, vagy sem. A 6b. ábra a DMC-k arányát szemlélteti a különböző genomi régiókban. Az intergenikus régiókon kívül egyetlen régióban sem tapasztaltunk DMC-gazdagodást.
Differenciális metilációs elemzés. a Differenciálisan metilált régiók (DMR), differenciálisan metilált gének (DMG) és differenciálisan metilált promóterek (DMP) elemzési eredményei a H9 és IMR90 sejtvonalak között. b A differenciálisan metilált citozinok (DMC) annotációja a különböző genomiális tulajdonságokkal és ismétlődő elemekkel szemben. c A H9 vagy IMR90 specifikus indeleket tartalmazó DMP-k (narancssárga) jelentős arányt foglalnak el az összes DMP-ben (DNS-metilációs differenciált promóterek)
A differenciálisan metilált promóterek (DMP-k) jelentős hányada tartalmaz indeleket
Tudjuk, hogy az indelek és a DNS-metiláció fontos szerepet játszanak a szövetek fejlődésében és a betegségekben . Itt a differenciálisan metilált promóterek (DMP-k) és az indelek közötti kapcsolatot vizsgáljuk. Ezt a vizsgálatot IMR90 és H9 sejtvonalak BS-Seq leolvasásainak felhasználásával végeztük. Először a BS-Seq leolvasásokat igazítottuk a BatMeth2 segítségével; majd az indeleket BisSNP és GATK eszközökkel hívtuk ki. Ezt követően azokat az indeleket, amelyek csak a H9-ben vagy az IMR90-ben fordulnak elő, sejtvonal-specifikus indeleknek definiáltuk.
Ezután a BatMeth2 segítségével 1384 DMP-t detektáltunk a H9 és az IMR90 között (p érték< 0,05, meth.diff > = 0,6). A fenti DMP-k közül összesen 236 (17%) tartalmaz indelt, amint az a 6c. ábrán látható. Röviden, a DMP-k jelentős része tartalmaz indeleket. Ezért a BS-Seq olvasatok pontos igazítása ezen indelek közelében nagyon fontos a DNS-metiláció kutatása és feltárása szempontjából.