Integroitu paketti bisulfiitti-DNA-metylaatiodatan analyysiin indel-sensitiivisellä kartoituksella

Helppokäyttöinen autorun-paketti DNA-metylaatioanalyysejä varten

Toteuttaaksemme DNA-metylaatiodatan analyysin helpommin, olemme koonneet kaikki toiminnot helppokäyttöiseen autorun-pakettiin DNA-metylaatiomääritysten analysointia varten. Kuvassa 1 esitetään BatMeth2:n tärkeimmät ominaisuudet: 1) BatMeth2:lla on tehokas ja tarkka kohdistussuoritus. 2) BatMeth2 voi laskea yksittäisten sytosiinipaikkojen tai minkä tahansa toiminnallisen alueen, kuten kokonaisten kromosomien, geenialueiden, transponoituvien elementtien (TE) jne. DNA-metylaatiotason (ML). 3) Erilaisten tilastollisten algoritmien integroinnin jälkeen BatMeth2 voi suorittaa differentiaalisen DNA-metylaatioanalyysin mille tahansa alueelle, mille tahansa määrälle syöttönäytteitä ja käyttäjän vaatimuksille. 4) Integroimalla BS-Seq-datan visualisoinnin (DNA-metylaatiojakauma kromosomeissa ja geeneissä) ja differentiaalisen metylaation annotaation BatMeth2 voi visualisoida DNA-metylaatiotiedot selkeämmin. BatMeth2-työkalun suorittamisen aikana luodaan html-raportti näytteen tilastoista. Näytteen html-raportin yksityiskohdat on esitetty kohdassa http://htmlpreview.github.io/?https://github.com/GuoliangLi-HZAU/BatMeth2/blob/master/BatMeth2-Report/batmeth2.html.

BatMeth2:lla on parempi kartoitussuorituskyky simuloidussa BS-Seq-datassa

Arvioimme ensin kaikkia kohdistimia käyttäen simuloituja datasettejä (ilman indelejä), jotka koostuivat 75 emäsparin (bp), 100 bp:n ja 150 bp:n mittaisista lukemista erilaisilla bisulfiittikonversiolukemisprosentteihin suhteutetuilla lukumäärillä (vaihteluväli oli välillä 0-100 % askeleella 10 %). Nämä tietokokonaisuudet simuloitiin ihmisen genomista (UCSC hg19) käyttäen FASTX-mutate-tools -työkaluja, wgsimiä (v0.3.0) ja SAMtoolsissa (v1.1) olevaa simulaattoria, joka sallii 0,03 % indelejä, 1 %:n emäsvirheprosentin koko genomissa ja korkeintaan kaksi mismatchia lukua kohti. Kartoitimme simuloidut lukemat referenssigenomiin sallimalla enintään kaksi epäsuhtaa. Koska simuloitujen lukujen alkuperäiset sijainnit olivat tiedossa, pystyimme arvioimaan kaikkien ohjelmien tarkkuutta vertaamalla niiden kartoitustuloksia alkuperäisiin sijainteihin.

Vertaillaksemme eri ohjelmistojen suorituskykyä, indeleitä sisältävä sekvensointiluku katsottiin oikein kartoitetuksi, jos seuraavat ehdot täyttyivät: 1) lukema kartoitettiin yksikäsitteisesti samalle säikeelle, josta se simuloitiin, ja kartoituksen laatu oli suurempi kuin 0; 2) kohdistetun lukeman raportoitu alkuasento oli enintään kymmenen emäsparin sisällä simuloidun lukeman alkuperäisestä alkuasennosta; 3) kartoitustulokset sisälsivät samankaltaisia indeleitä tai epäsuhtaisuuksia kuin simuloidussa lukemassa. Jos jokin näistä ehdoista ei täyttynyt, lukua pidettiin väärin kartoitettuna. Koska BatMeth2 sallii yhden aukon siemenalueella, se pystyy löytämään indelejä sisältävät siemenpaikat suurella tarkkuudella ja välttämään epäsopivat paikat, jotka aiheuttaisivat indelejä sisältävien lukujen väärän kohdistuksen. Kuvassa 2 esitetyt tulokset osoittavat, että BatMeth2 saavutti suurimman määrän oikein kohdistettuja lukuja ja pienimmän määrän virheellisesti kohdistettuja lukuja kaikissa testitietoaineistoissa eri bisulfiittimuunnosnopeuksilla.

Lyhyesti sanottuna wgsim-simuloitujen indel-aberrantti-datasettien tulokset osoittavat, että BatMeth2:lla on parempi suorituskyky (1~2 % parempi kuin toiseksi parhaimmalla alignerilla) kuin muilla menetelmillä, kun se linjaa yleisiä simuloituja biosulfiittilukulukulukulauseita (BS reads), jotka sisältävät sekoituksia epäsuhtaisuudesta (mismatches) ja indeleistä. Näemme, että kun BS-muunnosnopeus kasvaa, kaikkien ohjelmistojen kohdistustarkkuus heikkenee. Näissä erilaisissa olosuhteissa BatMeth2 toimii paremmin.

BatMeth2:lla on parempi kartoitussuorituskyky todellisessa BS-Seq-datassa

Testataksemme BatMeth2:n suorituskykyä todellisissa BS-Seq-tietoaineistoissa, latasimme parittaisia BS-Seq-tietoaineistoja ja poimimme satunnaisotannalla miljoona 2 × 90 bp:n parittaista lukulukua SRA SRR847318:sta, 1 miljoona 2 × 101 bp pareittain päättyvää lukua SRA SRR1035722:sta ja 1 miljoona 2 × 125 bp pareittain päättyvää lukua SRA SRR3503136:sta arviointia varten. Koska nämä tietokokonaisuudet ovat peräisin terveistä solulinjoista tai kudoksista, niiden odotetaan sisältävän vain vähän rakenteellisia variaatioita. Näin ollen kohdistimme nämä todelliset tiedot käyttämällä pareittain kerättyjen tietokokonaisuuksien yksipäätteisiä lukuja ja arvioimme pareittaisista kohdistuksista saadut yhteneväiset ja epäyhteneväiset kartoitusmäärät arvioidaksemme oikeiden ja virheellisten kohdistusten osuudet. Koska parilukujen inserttikoko oli noin 500 bp, parilukuparia voitiin pitää yhtenevänä, jos ne karttuivat 500 bp:n nimellisellä etäisyydellä; muussa tapauksessa parilukuparia voitiin pitää epäyhtenevänä. Samoin kuin simuloidulla datalla saamamme tulokset, BatMeth2 raportoi enemmän yhteneviä ja vähemmän epäyhteneviä kohdistuksia todellisilla dataseteillä laajalla kartan laatupisteiden vaihteluvälillä, kuten kuvassa 3 on esitetty.

Taulukossa 1 on lisäksi esitetty ohjelmien suhteelliset suoritusajat. BatMeth2 toimi oletusasetuksilla nopeammin kuin useimmat julkaistut kohdistimet ja oli verrattavissa BWA-methiin ja BatMethiin. Bismark2 (jossa Bowtie2 on perustavanlaatuinen kartoitusmenetelmä), BS Seeker2 ja BSmooth vaativat pidemmät suoritusajat.

DNA-metylaation kutsuminen

Arvioidaksemme DNA-metylaation kutsumisen tarkkuutta eri ohjelmistojen välillä latasimme 450 K:n bead chip -datan IMR90-solulinjalta ENCODE:stä (ENCODE = Encyclopedia of DNA Elements). Latasimme myös IMR90-solulinjan koko genomin bisulfiittisekvensointitiedot (WGBS-Seq) ENCODE:sta (42,6 Gbases). Kunkin ohjelmiston osalta kohdistimme WGBS-Seq-lukemat ja laskimme DNA-metylaatiotason. Sitten verrattiin tuloksia 450 K Bead Chip -datan samojen kohtien ML-arvoihin. Kun ohjelmiston WGBS-Seq-datasta saadun DNA ML:n ja 450 K Bead Chip -aineistosta saadun ML:n välinen ero oli alle 0,2, kutsutulos määriteltiin oikeaksi; muussa tapauksessa sitä pidettiin virheellisenä.

Tulokset on esitetty taulukossa 2. Kaikkien ohjelmistojen oikeiden tulosten päällekkäisyys on esitetty lisätiedostossa 1: Kuva S2. Näemme, että BatMeth2:lla ja Biscuitilla on samanlaiset suoritukset, jotka ovat parempia kuin muiden ohjelmistojen. Yhteenvetona voidaan todeta, että BatMeth2 parantaa sekä BS-lukujen kohdistamisen että DNA ML -kutsun tarkkuutta.

BatMeth2 kohdistaa BS-lukuja sallien samalla vaihtelevan pituiset indelit

Syöpä sisältää huomattavan paljon suurempaa indelien osuutta kuin terveet solut. Sen vuoksi tarkistaaksemme, voiko BatMeth2 kohdistaa BS-lukemia, joissa on eripituisia indellejä, latasimme ENCODE:stä WGBS-datan (75 Gbases) ja 450 K Bead Chip -datan HepG2:sta (maksan hepatosellulaarinen karsinooma, syöpäsolulinja). Tarkistimme indelien pituusjakauman lukemissa HepG2:n WGBS-Seq-datan kohdistamisen jälkeen. Lisätiedosto 1: Kuva S3A osoittaa, että havaittujen indelien pituudet jakautuivat pääasiassa 1 bp~ 5 bp:n välille, ja pisin indel oli 40 bp pitkä. Tilastojemme mukaan 2,3 prosenttia alignointilukemista sisälsi indellejä. Näiden tulosten perusteella tiedämme, että BatMeth2 pystyy kohdistamaan lukuja, joissa on eripituisia indellejä.

Seuraavaksi testasimme indelien havaitsemisen vaikutusta DNA-metylaation kutsumiseen. BatMeth2:n osalta ajoimme HepG2-datalla kaksi vaihtoehtoa: indeleiden havaitsemisen kanssa ja ilman sitä (eli asetimme -I-parametrin BatMeth2:ssa). Ajoimme myös Bismarkia HepG2:n WGBS-Seq-datalla vertailuarvona DNA-metylaation kutsumiselle indelien havaitsemisen kanssa, koska Bismarkissa ei ole indelien kutsumistoimintoa. Vertasimme DNA-metylaation kutsumista BatMeth2:ssa ja Bismarkissa 450 K Bead Chip -datasta saatuun kutsumiseen. Tulokset on esitetty lisätiedostossa 1: Kuva S3B, jossa ”BatMeth2-noIndel” vastaa BatMeth2:ta ilman indel-tunnistusta. Näemme, että ilman indel-tunnistusta BatMeth2:n tulos oli vain hieman parempi kuin Bismarkin tulos (Bowtie1:n ollessa peruskartoitusmenetelmä). BatMeth2:n tulos indelien havaitsemisen kanssa oli huomattavasti parempi. Lisäksi voidaan havaita, että BatMeth2 pystyy havaitsemaan enemmän DNA-metylaatiokohtia kuin BatMeth2-noIndel ja Bismark (Bowtie 1). Ymmärtääksemme, miksi BatMeth2:n suorituskyky indelien havaitsemisen kanssa on parempi, määrittelimme BatMeth2:n havaitsemat metylaatiokohdat tulokseksi A, kun taas BatMeth2-noIndelin ja Bismarkin havaitsemat metylaatiokohdat määriteltiin tulokseksi B. Havaitsimme, että mclA sisälsi 23 853 DNA:n metylaatiokohdetta ja 15 048 (63 %) niistä 23 853 kohdasta, jotka BatMeth2:n indel-tunnistuksella kutsumien indel-lukujen kohdistukset kattoivat (ks. lisätiedosto 1: kuva S3C). Lisäksi havaitsimme, että indel-asteet tuloksessa A olivat vain 5 % ja tuloksessa B 0 %. Näin ollen päädyimme siihen, että tarkka indelien havaitseminen voi parantaa DNA-metylaation kutsumista.

DNA-metylaatiotietojen visualisointi

BatMeth2 tarjoaa työkaluja metylaatiotietojen visualisointiin. BatMeth2:n visualisointiominaisuuksien havainnollistamiseksi latasimme (1) 117 Gt:n edestä yksittäislukemia ihmisen H9-solulinjasta, (2) 105,2 Gt:n edestä yksittäislukemia ihmisen IMR90-solulinjasta ja (3) 12,6 Gt:n edestä pareittain luettuja lukemia villityypin riisistä. BatMeth2:lla voidaan ensinnäkin visualisoida sytosiinimetylaatiotiheys kromosomitasolla. Pisteet kuvassa 4a edustavat 100 kb:n liukuikkunaa, jonka askel on 50 kb. Jotta ML:ää voidaan tarkastella yksittäisissä CpG- tai muissa kuin CpG-kohdissa genomiselaimessa, tarjoamme tiedostot myös bed- ja bigWig-muodossa (kuva 4b). Vertailemalla geenien ja TE:iden tiheyteen havaitsimme, että ML korreloi TE-tiheyden kanssa ja oli antikorreloitunut geenitiheyden kanssa (kuva 4c). Tämä suuntaus on aiemmin havaittu riisissä .

Toisekseen BatMeth2 voi visualisoida geenien ML:t. Tarkemmin sanottuna BatMeth2 voi visualisoida ML:t 2 kb geenistä ylävirtaan, transkription aloituskohdassa (TSS), geenirungossa, transkription lopetuskohdassa (TES) ja 2 kb geenirungosta alavirtaan. Kun verrataan ylävirran, geenirungon ja alavirran alueita, kuvasta 5a nähdään, että geenirungon DNA:n ML on korkeampi kuin promoottorin alueella. Kun verrataan kaikkia viittä aluetta, TSS-alueella on selvästi laakso (kuva 5b). BatMeth2 voi myös laskea ML-profiilit intronien, eksonien, intergeenisten alueiden ja TE:iden ympärillä (lisätiedosto 1: kuva S4). Lisäksi BatMeth2 voi tarjota lämpökartan useista geeneistä geenialueittain eri näytteiden geenien kokonais-ML:ien kätevää vertailua varten (kuva 5c) (kuva 5c).

Kolmanneksi BatMeth2:lla pystytään havainnollistamaan DNA-metylaatioiden jakauma. Lisätiedosto 1: Kuvassa S5A esitetään DNA-metylaatiojakaumat H9- ja IMR90-solulinjoissa. Kuvassa DNA ML on jaettu viiteen luokkaan: metyloitunut (M: > 80 %), osittain metyloituneen ja metyloituneen välimuoto (Mh: 60-80 %), osittain metyloitunut (H: 40-60 %), metyloitumattoman ja osittain metyloituneen välimuoto (hU: 20-40 %) ja metyloitumaton (U: < 20 %). Kuten lisätiedostosta 1 käy ilmi: Kuva S5A, ML oli korkeampi H9-solulinjassa M-luokassa kuin IMR90-solulinjassa, erityisesti CpG-kontekstissa. CH-sekvenssikontekstissa CpG-metylaatio on vallitseva muoto, mutta merkittävä osa metyloituneista sytosiinista löytyy CpA-kohdista, kun taas ML on alle 40 %, erityisesti H9-solulinjassa (Additional file 1: Figure S5B).

Neljänneksi BatMeth2:lla voidaan analysoida geenin ilmentymistason ja geenin promoottori-DNA:n ML:n välistä korrelaatiota. Havainnollistimme tätä ominaisuutta käyttämällä H9- ja IMR90-solulinjoja. H9- tai IMR90-geenien ilmentymistasot jaettiin eri luokkiin. Kuten lisätiedostossa 1 on esitetty: Kuva S5C, voimakkaasti ilmentyneillä geeneillä oli alhaisemmat ML:t promoottorialueillaan. Lisäksi jaoimme geenien promoottorien ML:t viiteen luokkaan. Tulos lisätiedostossa 1: Kuva S5D osoittaa, että geeneillä, joiden promoottoreilla oli korkeammat ML-arvot, oli alhaisemmat ekspressiotasot. Nisäkkäiden geenien ilmentymisen ja promoottorin DNA-metylaation välinen negatiivinen korrelaatio tunnetaan . Tämä analyysi osoittaa edelleen BatMeth2:n tarkkuutta.

Differentiaalisesti metyloitujen sytosiinien ja alueiden (DMC:t/DMR:t)

Differentiaalisesti metyloitujen sytosiinien (DMC:t) ja differentiaalisesti metyloitujen alueiden (DMR:t) tunnistaminen on yksi tärkeimmistä tavoitteista metylaatiodatan analysoinnissa. Vaikka tutkijat ovat toisinaan kiinnostuneita korreloimaan yksittäisiä sytosiinikohtia fenotyyppiin , DMR:t ovat erittäin tärkeitä ominaisuuksia .

Varhaisissa BS-Seq-tutkimuksissa soluja profiloitiin keräämättä toistoja. Tällaisten tietokokonaisuuksien osalta käytimme Fisherin tarkkaa testiä erottaaksemme eri tavoin metyloituneet sytosiinit (DMC). BS-Seq-tietoaineistoissa, joissa on toistoja, luonnollisin tilastollinen malli DMC:iden nimeämiseen on beeta-binomijakauma . Tiedämme, että useat ohjelmistot voivat suorittaa DNA:n metylaatiodatan differentiaalisen analyysin, kuten methylKit (differentiaalinen analyysiohjelma, joka vaatii biologisia replikaatteja) ja Methy-Pipe (differentiaalinen analyysiohjelma ilman biologisia replikaatteja). Saatavilla ei kuitenkaan ole kattavaa pakettia, joka sisältäisi sekä kartoituksen että differentiaalisen metylaatioanalyysin. Niinpä kehitimme paketin, joka yhdistää kartoituksen ja differentiaalianalyysin. Helpottaaksemme DMR:ien tunnistamista bisulfiittidatasta, jossa ei ole toistoja, integroimme Fisherin tarkan testin hypoteesitestin suorittamiseksi. Kun näytteessä on kaksi tai useampia toistoja, käytämme beeta-binomijakaumaa differentiaalisen metylaatioanalyysin suorittamiseen. Tarjoamme myös bed- tai bigWig-tiedostoja DMR-luetteloa varten. DMR:t voidaan visualisoida genomiselaimessa (kuva 4b) luotujen bed- tai bigWig-tiedostojen avulla.

Kuvassa 6a on havainnollistettu DMC:iden ja alueiden lukumäärät IMR90-solulinjassa ja H9-solulinjassa BatMeth2:lla havaittuna (p-arvo< 0,05, meth.diff > = 0,6). BatMeth2 voi visualisoida, ovatko CpG:t ja DMC:t rikastuneet tietyillä alueilla, kuten geeni-, CDS-, introni-, intergeenisillä, UTR-, TE-, LTR-, LINE- ja SINE-alueilla. Kuvassa 6b visualisoidaan DMC:iden osuudet eri genomialueilla. Intergeenisiä alueita lukuun ottamatta emme havainneet DMC-rikastumista millään alueella.

Huomattava osa differentiaalisesti metyloituneista promoottoreista (DMP:istä) sisältää indellejä

Tietääksemme tiedämme, että indeleillä ja DNA:n metylaatiolla on merkittävä rooli kudosten kehityksessä ja sairauksissa . Tässä tutkimme eri tavoin metyloituneiden promoottorien (DMP) ja indeleiden välistä suhdetta. Teimme tämän tutkimuksen käyttämällä BS-Seq-lukemia IMR90- ja H9-solulinjoissa. Kohdistimme BS-Seq-lukemat ensin BatMeth2:n avulla, minkä jälkeen indelit kutsuttiin BisSNP- ja GATK-työkaluilla. Tämän jälkeen määrittelimme indelit, jotka esiintyvät vain H9:ssä tai IMR90:ssä, solulinjakohtaisiksi indeleiksi.

Tällöin havaitsimme BatMeth2:lla 1384 DMP:tä H9:n ja IMR90:n välillä (p-arvo< 0,05, met.diff > = 0,6). Kaikkien edellä mainittujen DMP:iden joukossa on yhteensä 236 (17 %) indeliä, kuten kuvasta 6c käy ilmi. Lyhyesti sanottuna huomattava osa DMP:istä sisältää indeleitä. Siksi BS-Seq-lukujen tarkka kohdistaminen näiden indeleiden läheisyydessä on erittäin tärkeää DNA-metylaation tutkimuksessa ja selvittämisessä.