Een eenvoudig te gebruiken, autorun pakket voor DNA methylatie analyses
Om de DNA methylatie data analyse handiger te maken, hebben we alle functies verpakt in een eenvoudig te gebruiken, autorun pakket voor DNA methylatie analyse. Figuur 1 toont de belangrijkste kenmerken van BatMeth2: 1) BatMeth2 heeft efficiënte en nauwkeurige uitlijningsprestaties. 2) BatMeth2 kan het DNA-methyleringsniveau (ML) berekenen van individuele cytosineplaatsen of functionele regio’s, zoals hele chromosomen, genregio’s, transposable elementen (TE’s), enz. 3) Na de integratie van verschillende statistische algoritmen kan BatMeth2 differentiële DNA-methyleringsanalyses uitvoeren voor elke regio, elk aantal inputmonsters en gebruikersvereisten. 4) Door het integreren van BS-Seq data visualisatie (DNA methylatie distributie op chromosomen en genen) en differentiële methylatie annotatie, kan BatMeth2 de DNA methylatie data duidelijker visualiseren. Tijdens de uitvoering van de BatMeth2 tool, wordt een html rapport gegenereerd voor de statistieken van het monster. Details van het html-verslag worden getoond in http://htmlpreview.github.io/?https://github.com/GuoliangLi-HZAU/BatMeth2/blob/master/BatMeth2-Report/batmeth2.html.
De workflow van BatMeth2. De twee grote pijlen betekenen invoer- of uitvoerbestanden
BatMeth2 heeft betere mapping-prestaties op gesimuleerde BS-Seq-gegevens
We hebben eerst alle aligners geëvalueerd met gesimuleerde datasets (zonder indels) bestaande uit leest met 75 basenparen (bp), 100 bp en 150 bp en met verschillende bisulfietconversiepercentages (variërend van 0 tot 100% met stap 10%). Deze datasets werden gesimuleerd van het menselijk genoom (UCSC hg19) met behulp van FASTX-mutate-tools , wgsim (v0.3.0) en de simulator in SAMtools (v1.1) , die 0,03% indels, een 1% base error rate in het hele genoom en een maximum van twee mismatches per lees mogelijk maakt. We brachten de gesimuleerde data in kaart met het referentie genoom, waarbij we maximaal twee mismatches toestonden. Omdat de originele posities van de gesimuleerde lezingen bekend waren, konden we de nauwkeurigheid van alle programma’s evalueren door hun mapping outputs te vergelijken met de originele posities.
Om de prestaties van de verschillende software te vergelijken, werd een sequencing lees met indels beschouwd als correct in kaart gebracht als de volgende voorwaarden waar waren: 1) de lees was uniek in kaart gebracht op dezelfde streng als waarvandaan was gesimuleerd en de mapping kwaliteit was groter dan 0; 2) de gerapporteerde startpositie van de uitgelijnde lees was binnen tien basenparen van de oorspronkelijke startpositie van de gesimuleerde lees; 3) de mapping resultaten hadden soortgelijke indels of mismatches als de gesimuleerde lees. Als aan een van deze voorwaarden niet werd voldaan, werd de aflezing als foutief in kaart gebracht beschouwd. Omdat BatMeth2 één hiaat in de zaadregio toestaat, kan het zaadlocaties met indels met hoge nauwkeurigheid vinden en kan het mismatches vermijden, die zouden leiden tot foutieve uitlijning van lezingen met indels. De resultaten in Fig. 2 laten zien dat BatMeth2 het grootste aantal correct uitgelijnde reads en het laagste aantal incorrect uitgelijnde reads heeft bereikt in alle testdatasets bij verschillende bisulfietconversiesnelheden.
Evaluatie van alle BS-Seq aligners met behulp van gesimuleerde datasets met verschillende leeslengtes van FASTX en wgsim. Gesimuleerde gegevens met verschillende bisulfietconversies worden in verschillende vormen weergegeven. Resultaten van verschillende aligners worden getoond met verschillende kleuren van de symbolen. De resultaten bij de linkerbovenhoek in elk paneel laten zien dat de software meer correct gemapte leest en het lagere aantal incorrect gemapte leest bereikte. De resultaten van onze aligner BatMeth2 zijn de beste in de verschillende gesimuleerde bisulfiet datasets
In het kort, de resultaten van wgsim-gesimuleerde indel-aberrant datasets laten zien dat BatMeth2 betere prestaties (1 ~ 2% beter dan de tweede top aligner) dan de andere methoden heeft bij het afstemmen van algemene gesimuleerde BS leest die een mengsel van mismatches en indels. We kunnen zien dat met de toegenomen BS conversie, de uitlijn nauwkeurigheid van alle software wordt verminderd. Onder deze verschillende omstandigheden presteert BatMeth2 beter.
BatMeth2 heeft betere mapping-prestaties op echte BS-Seq-gegevens
Om de prestaties van BatMeth2 op echte BS-Seq-datasets te testen, hebben we gepaarde-end BS-Seq-datasets gedownload en willekeurig 1 miljoen 2 × 90 bp gepaarde-end reads van SRA SRR847318 geëxtraheerd, 1 miljoen 2 × 101 bp gepaarde-end leest van SRA SRR1035722 en 1 miljoen 2 × 125 bp gepaarde-end leest van SRA SRR3503136 voor evaluatiedoeleinden. Omdat deze datasets afkomstig zijn van gezonde cellijnen of weefsels, wordt verwacht dat ze een laag aantal structurele variaties bevatten. Daarom hebben we deze echte data uitgelijnd met single-end reads van de gepaarde datasets en hebben we de concordante en discordante mapping rates van de gepaarde uitlijningen geëvalueerd om de correcte en incorrecte uitlijningspercentages te schatten. Omdat de insert grootte van de gepaarde-end leest ongeveer 500 bp was, kon een paar van partner leest worden beschouwd als concordant als ze werden gemapped binnen een nominale afstand van 500 bp; anders kon een paar van partner leest worden beschouwd als discordant. Vergelijkbaar met onze resultaten met de gesimuleerde data, rapporteerde BatMeth2 meer concordante en minder discordante alignments op de echte datasets over een groot bereik van kaart kwaliteit scores, zoals getoond in Fig. 3.
Concordance en discordance rates van alignments op echte paired-end reads van verschillende aligners. Cumulatieve telling van concordante en discordante uitlijningen van hoge naar lage mapping kwaliteit voor echte bisulfiet sequencing leest. Er is slechts één punt voor BSmap en de aligners op basis van bowtie apart, omdat deze aligners geen kaartkwaliteitsscore hebben. Bismark-bowtie2L15 betekent bowtie2 alignment met zaadlengte 15
Bovendien toont tabel 1 de relatieve runtimes van de programma’s. BatMeth2 met de standaard instellingen liep sneller dan de meeste van de gepubliceerde aligners en was vergelijkbaar met BWA-meth en BatMeth. Bismark2 (met Bowtie2 als de fundamentele mapping-methode), BS Seeker2 en BSmooth hebben langere runtijden nodig.
DNA methylation calling
Om de nauwkeurigheid van DNA-methylation calling tussen verschillende software te evalueren, hebben we 450 K bead chip gegevens van de IMR90 cellijn gedownload van ENCODE (Encyclopedia of DNA Elements). We hebben ook gedownload hele-genoom bisulfiet sequencing (WGBS-Seq) gegevens van de IMR90 cellijn van ENCODE (42,6 Gbases). Voor elke software hebben we de WGBS-Seq reads uitgelijnd en het niveau van DNA-methylering berekend. Vervolgens hebben we de resultaten vergeleken met de MLs op dezelfde plaatsen in de 450 K Bead Chip gegevens. Wanneer het verschil tussen de DNA ML van de WGBS-Seq-gegevens door de software en die van de 450 K Bead Chip minder dan 0,2 was, werd het aanroepende resultaat als correct gedefinieerd; anders werd het als incorrect beschouwd.
De resultaten zijn weergegeven in tabel 2. De overlap tussen de correcte resultaten van alle software wordt getoond in Extra bestand 1: Figuur S2. We kunnen zien dat BatMeth2 en Biscuit vergelijkbare prestaties hebben, die beter zijn dan die van de andere software. BatMeth2 verbetert dus de nauwkeurigheid van zowel BS-lees alignment als DNA ML calling.
BatMeth2 alignments BS-lees while allowing for variable-length indels
Cancer bevat een aanzienlijk hoger aandeel indels dan gezonde cellen doen. Daarom, om na te gaan of BatMeth2 BS leest kan uitlijnen met indels van verschillende lengtes, hebben we WGBS data (75 Gbases) en 450 K Bead Chip data van HepG2 (lever hepatocellulair carcinoom, een kanker cellijn) gedownload van ENCODE. We controleerden de indel lengteverdeling in de leest na de uitlijning van HepG2 WGBS-Seq gegevens. Extra bestand 1: Figuur S3A toont aan dat de lengtes van de gedetecteerde indels werden voornamelijk verdeeld in de 1 bp ~ 5 bp bereik, en de langste indel was 40 bp in lengte. Volgens onze statistieken, 2,3% van de alignment leest bevatte indels. Uit deze resultaten weten we dat BatMeth2 kan uitlijnen leest met indels van verschillende lengtes.
Volgende, testten we het effect van indel detectie op DNA methylering roepen. Voor BatMeth2, draaiden we twee opties op de HepG2 gegevens: met en zonder indel detectie (dat wil zeggen, stel -I parameter in BatMeth2). We hebben ook Bismark uitgevoerd op de WGBS-Seq gegevens van HepG2 als referentie voor DNA-methylering met indel-detectie, omdat Bismark geen indel-oproepfunctie heeft. We vergeleken het oproepen van DNA-methylering in BatMeth2 en Bismark met het oproepen uit de 450 K Bead Chip-gegevens. De resultaten zijn te zien in Additional file 1: Figuur S3B, waarbij “BatMeth2-noIndel” overeenkomt met BatMeth2 zonder indel-detectie. We kunnen zien dat, bij afwezigheid van indel-detectie, het resultaat van BatMeth2 slechts iets beter was dan dat van Bismark (met Bowtie1 als de fundamentele mapping-methode). Het resultaat van BatMeth2 met indeldetectie was aanzienlijk beter. Bovendien kunnen we zien dat BatMeth2 meer DNA-methyleringsplaatsen kan detecteren dan BatMeth2-noIndel en Bismark (Bowtie 1). Om te begrijpen waarom de prestatie van BatMeth2 met indel-detectie beter is, hebben we de methyleringsplaatsen die door BatMeth2 worden opgeroepen gedefinieerd als Resultaat A, terwijl de methyleringsplaatsen die door BatMeth2-noIndel en Bismark worden opgeroepen gedefinieerd zijn als Resultaat B. Vervolgens laten we mclA de methylatie sites zijn die voorkomen in Resultaat A maar niet in Resultaat B. We zagen dat mclA 23.853 DNA methylatie sites omvatte en 15.048 (63%) van de 23.853 sites die gedekt werden door de uitlijningen van indel reads opgeroepen door BatMeth2 met indel detectie (zie Additional file 1: Figuur S3C). Bovendien vonden we dat de indel percentages in Resultaat A en Resultaat B slechts 5 en 0% waren, respectievelijk. Vandaar dat we concludeerden dat nauwkeurige indel detectie kan verbeteren DNA methylatie calling.
Visualisatie van DNA methylatie data
BatMeth2 biedt tools om de methylatie data te visualiseren. Ter illustratie van de visualisatie functies van BatMeth2, we gedownload (1) 117 Gbases van single-end leest van de menselijke H9 cellijn, (2) 105,2 Gbases van single-end leest van de menselijke IMR90 cellijn en (3) 12,6 Gbases van gepaarde-end leest van wild-type rijst. Ten eerste, BatMeth2 kan visualiseren cytosine methylatie dichtheid op het chromosoom niveau. De stippen in Fig. 4a vertegenwoordigen een glijdend venster van 100 kb met een stap van 50 kb. Om het bekijken van de ML op individuele CpG of niet-CpG sites in een genoom browser mogelijk te maken, bieden wij ook bestanden in bed en bigWig formaten (Fig. 4b). Door vergelijking met de dichtheid van genen en TEs, zagen we dat het ML gecorreleerd was met de TE dichtheid en anticorreleerde met de gen dichtheid (Fig. 4c). Deze tendens is eerder waargenomen in rijst .
Visualisatie van de methyleringsniveaus in chromosoomschaal. a De methyl-cytosine dichtheid in humaan chromosoom 10. De stippen vertegenwoordigen de methylering niveaus in glijdende vensters van 100Kb met een stap van 50Kb. De rode stippen verwijzen naar de methylering niveaus in de plus streng, en de blauwe stippen verwijzen naar de methylering niveaus in de min-streng. b Een voorbeeld over de distributies van de DNA-methylering niveaus en differentieel gemethyleerd regio’s (DMRs) tussen H9 en IMR90 cellijnen in het menselijk chromosoom 10. c De dichtheid van genen, transposon elementen (TEs) en het niveau van DNA-methylering in het hele rijst-genoom. Paneel A is de resultaten gegenereerd uit Batmeth2. Paneel B is de visualisatie resultaten van UCSC browser, met de BED-bestanden van Batmeth2
Ten tweede, BatMeth2 kan visualiseren de MLs van genen. Meer precies kan BatMeth2 de ML’s 2 kb stroomopwaarts van het gen, op de transcriptiestartplaats (TSS), in het genlichaam, op de transcriptie-eindplaats (TES) en 2 kb stroomafwaarts van het genlichaam visualiseren. Vergelijking van de upstream-, body- en downstream-regio’s toont in fig. 5a dat het DNA ML van het genlichaam hoger is dan dat in de promotorregio. Bij vergelijking van alle vijf regio’s is er duidelijk een dal in de TSS-regio (Fig. 5b). BatMeth2 kan ook de ML profielen rond introns, exons, intergene regio’s en TE’s berekenen (Additional file 1: Figuur S4). Bovendien kan BatMeth2 een heat map van meerdere genen per gen regio voor een handige vergelijking van de totale gen MLs van verschillende monsters (Fig. 5c).
Visualisatie van DNA-methylatie onder verschillende contexten. a De DNA-methyleringsniveaus in 2Kb regio’s stroomopwaarts van genen, genlichamen, 2Kb stroomafwaarts van genlichamen. b Het aggregatieprofiel van DNA-methyleringsniveaus over genen. c De hittekaart van alle genen in 2Kb regio’s stroomopwaarts van genen, genlichamen, 2Kb stroomafwaarts van genlichamen
Derde, BatMeth2 kan de verdeling van DNA-methylering visualiseren. Additional file 1: Figuur S5A toont de DNA methylering distributies in de H9 en IMR90 cellijnen. In de figuur is het DNA ML verdeeld in vijf categorieën: gemethyleerd (M: > 80%), intermediair tussen gedeeltelijk gemethyleerd en gemethyleerd (Mh: 60-80%), gedeeltelijk gemethyleerd (H: 40-60%), intermediair tussen niet-gemethyleerd en gedeeltelijk gemethyleerd (hU: 20-40%), en niet-gemethyleerd (U: < 20%). Zoals te zien is in Additional file 1: Figuur S5A, was de ML hoger in de H9 cellijn in de M categorie dan in de IMR90 cellijn, vooral in de CpG context. In de CH-sequentiecontext is CpG-methylering de overheersende vorm, maar een significante fractie van gemethyleerde cytosines wordt gevonden op CpA-locaties, terwijl de ML minder dan 40% bedraagt, met name in de H9-cellijn (Additional file 1: figuur S5B).
Vierde, BatMeth2 kan de correlatie analyseren tussen genexpressieniveau en genpromoter-DNA ML. We illustreerden deze functie met behulp van de H9 en IMR90 cellijnen. De expressieniveaus van de genen in H9 of IMR90 werden onderverdeeld in verschillende categorieën. Zoals getoond in Extra bestand 1: Figuur S5C, vertoonden de hoog tot expressie komende genen lagere MLs in hun promotorgebieden. Verder verdeelden we de MLs van de genpromotors in vijf categorieën. Het resultaat in Additional file 1: Figuur S5D laat zien dat genen met promotors met hogere ML-waarden lagere expressieniveaus vertoonden. De negatieve correlatie tussen de expressie van zoogdiergenen en promoter DNA-methylering is bekend. Deze analyse geeft verder de nauwkeurigheid van BatMeth2 aan.
Het vinden van differentieel gemethyleerde cytosines en regio’s (DMCs/DMRs)
Het identificeren van differentieel gemethyleerde cytosines (DMCs) en differentieel gemethyleerde regio’s (DMRs) is een van de belangrijkste doelen bij de analyse van methyleringsgegevens. Hoewel onderzoekers af en toe geïnteresseerd zijn in het correleren van enkele cytosine sites aan een fenotype, zijn DMR’s zeer belangrijke kenmerken.
Early BS-Seq studies profileerden cellen zonder replicaten te verzamelen. Voor dergelijke datasets hebben we Fisher’s exact test gebruikt om differentieel gemethyleerde cytosines (DMC’s) te onderscheiden. Voor BS-Seq datasets met herhalingen, de meest natuurlijke statistische model om DMC’s te bellen is bèta-binomiale verdeling. We weten dat een aantal software programma’s differentiële DNA methylatie data analyse kunnen uitvoeren, zoals methylKit (een differentieel analyse programma dat biologische herhalingen vereist) en Methy-Pipe (een differentieel analyse programma zonder biologische herhalingen). Er is echter geen uitgebreid pakket beschikbaar dat zowel mapping als differentiële methylatie analyse omvat. Daarom hebben wij een pakket ontwikkeld dat mapping integreert met differentiële analyse. Om de identificatie van DMR’s uit bisulfiet data zonder replicaten te vergemakkelijken, hebben we Fisher’s exact test geïntegreerd om een hypothesetoets uit te voeren. Wanneer een monster twee of meer replicaten heeft, gebruiken we de bèta-bbinomiale verdeling om differentiële methyleringsanalyse uit te voeren. We leveren ook bed of bigWig bestanden voor de lijst van DMRs. De DMR’s kunnen worden gevisualiseerd in een genoom browser (Fig. 4b) met de gegenereerde bed of bigWig files.
Ter illustratie, Fig. 6a toont de aantallen DMC’s en regio’s in de IMR90 cellijn en in de H9 cellijn, zoals gedetecteerd door BatMeth2 (p waarde< 0.05, meth.diff > = 0.6). BatMeth2 kan visualiseren of CpG’s en DMC’s verrijkt zijn in bepaalde regio’s, zoals gen-, CDS-, intron-, intergenic-, UTR-, TE-, LTR-, LINE- en SINE-regio’s. Figuur 6b visualiseert de verhoudingen van DMC’s in verschillende genomische regio’s. Afgezien van de intergene regio’s hebben we in geen enkele regio een verrijking met DMC’s waargenomen.
Differentiële methyleringsanalyse. a Analyseresultaten van differentieel gemethyleerde regio’s (DMR’s), differentieel gemethyleerde genen (DMG’s) en differentieel gemethyleerde promotors (DMP’s) tussen H9- en IMR90-cellijnen. b Annotatie van differentieel gemethyleerde cytosines (DMC) tegen verschillende genomische eigenschappen en herhalingselementen. c DMPs bevatten H9 of IMR90 specifieke-indels (oranje) bezetten een substantieel deel in de alle DMPs (DNA Methylation differential Promoters)
Een substantieel deel van differentieel gemethyleerde promoters (DMPs) bevatten indels
We weten dat indels en DNA methylering een belangrijke rol spelen in weefselontwikkeling en ziekten . Hier onderzoeken we de relatie tussen differentieel gemethyleerde promotors (DMPs) en indels. We hebben deze studie uitgevoerd met de BS-Seq gegevens in IMR90 en H9 cellijnen. We hebben eerst de BS-Seq gegevens uitgelijnd met BatMeth2; daarna werden de indels opgeroepen met BisSNP en GATK tools. Vervolgens definieerden we de indels die alleen in H9 of IMR90 voorkomen als cellijn-specifieke indels.
Daarna ontdekten we 1384 DMPs tussen H9 en IMR90 door BatMeth2 (p waarde< 0.05, meth.diff > = 0.6). In totaal 236 (17%) van alle bovengenoemde DMP’s bevatten indels, zoals te zien is in fig. 6c. Kortom, een aanzienlijk deel van de DMP’s bevat indels. Daarom is nauwkeurige uitlijning van BS-Seq leest in de buurt van deze indels zeer belangrijk voor onderzoek en exploratie van DNA-methylering.