En integreret pakke til bisulfit DNA-methyleringsdataanalyse med Indel-følsom kortlægning

En brugervenlig, autorun-pakke til DNA-methyleringsanalyser

For at gøre DNA-methyleringsdataanalyser mere bekvemme har vi pakket alle funktionerne i en brugervenlig, autorun-pakke til DNA-methyleringsanalyse. Figur 1 viser de vigtigste funktioner i BatMeth2: 1) BatMeth2 har en effektiv og præcis tilpasningsevne. 2) BatMeth2 kan beregne DNA-methyleringsniveauet (ML) for individuelle cytosin-steder eller enhver funktionel region, såsom hele kromosomer, genregioner, transposable elementer (TE’er) osv. 3) Efter integration af forskellige statistiske algoritmer kan BatMeth2 udføre differentiel DNA-methyleringsanalyse for enhver region, ethvert antal inputprøver og brugerkrav. 4) Ved at integrere BS-Seq datavisualisering (DNA-methyleringsfordeling på kromosomer og gener) og differentiel methyleringsannotation kan BatMeth2 visualisere DNA-methyleringsdataene mere klart. Under udførelsen af BatMeth2-værktøjet genereres en html-rapport for prøvens statistik. Eksempel html-rapportdetaljer vises i http://htmlpreview.github.io/?https://github.com/GuoliangLi-HZAU/BatMeth2/blob/master/BatMeth2-Report/batmeth2.html.

BatMeth2 har bedre kortlægningsydelse på simulerede BS-Seq-data

Vi evaluerede først alle alignere ved hjælp af simulerede datasæt (uden indels) bestående af reads med 75 basepar (bp), 100 bp og 150 bp og med forskellige bisulfitkonverteringsrater (fra 0 til 100 % med trin 10 %). Disse datasæt blev simuleret ud fra det menneskelige genom (UCSC hg19) ved hjælp af FASTX-mutate-tools , wgsim (v0.3.0) og simulatoren i SAMtools (v1.1) , som tillader 0,03 % indels, en basefejlprocent på 1 % i hele genomet og højst to mismatches pr. læsning. Vi kortlagde de simulerede læsninger til referencegenomet, idet vi tillod højst to mismatches. Da de originale positioner for de simulerede læsninger var kendt, kunne vi evaluere nøjagtigheden af alle programmerne ved at sammenligne deres kortlægningsresultater med de originale positioner.

For at sammenligne de forskellige programmers præstationer blev en sekventeringslæsning med indels betragtet som korrekt kortlagt, hvis følgende betingelser var opfyldt: 1) læsningen blev entydigt kortlagt til den samme streng, som den blev simuleret fra, og kortlægningskvaliteten var større end 0; 2) den rapporterede startposition for den justerede læsning var inden for ti basepar af den oprindelige startposition for den simulerede læsning; 3) kortlægningsresultaterne havde lignende indels eller mismatches som den simulerede læsning. Hvis en af disse betingelser ikke blev overholdt, blev læsningen betragtet som forkert kortlagt. Da BatMeth2 tillader ét hul i seed-regionen, kan den finde seed-placeringer med indels med stor nøjagtighed og kan undgå mismatchede placeringer, hvilket ville medføre, at læsninger med indels bliver fejljusteret. Resultaterne i fig. 2 viser, at BatMeth2 opnåede det største antal korrekt justerede læsninger og det laveste antal forkert justerede læsninger i alle testdatasæt ved forskellige bisulfitkonverteringshastigheder.

Kort sagt viser resultaterne fra wgsim-simulerede indel-aberrant-datasæt, at BatMeth2 har bedre præstationer (1~2 % bedre end den næstbedste aligner) end de andre metoder, når der alignes generelle simulerede BS-reads indeholdende en blanding af mismatches og indels. Vi kan se, at med den øgede BS-konverteringshastighed reduceres tilpasningsnøjagtigheden for alle programmerne. Under disse forskellige betingelser klarer BatMeth2 sig bedre.

BatMeth2 har bedre kortlægningsydelse på rigtige BS-Seq-data

For at teste BatMeth2’s ydeevne på rigtige BS-Seq-datasæt downloadede vi parvis afsluttede BS-Seq-datasæt og udtog tilfældigt 1 million 2 × 90 bp parvis afsluttede reads fra SRA SRR847318, 1 million 2 × 101 bp paired-end-reads fra SRA SRR1035722 og 1 million 2 × 125 bp paired-end-reads fra SRA SRR3503136 til evalueringsformål. Da disse datasæt er fra sunde cellelinjer eller væv, forventes de at indeholde et lavt antal strukturelle variationer. Derfor tilpassede vi disse reelle data ved hjælp af single-end-reads fra de parrede datasæt og evaluerede de overenstemmende og uoverenstemmende kortlægningsrater fra de parrede tilpasninger for at estimere de korrekte og ukorrekte tilpasningsrater. Da indsætningsstørrelsen af de parvise endeaflæsninger var ca. 500 bp, kunne et par partneraflæsninger betragtes som konkordante, hvis de blev kortlagt inden for en nominel afstand på 500 bp; ellers kunne et par partneraflæsninger betragtes som diskordante. I lighed med vores resultater med de simulerede data rapporterede BatMeth2 flere samstemmende og færre diskordante alignments på de virkelige datasæt over et stort område af kortkvalitetsscorer, som vist i fig. 3.

Dertil kommer, at tabel 1 viser de relative køretider for programmerne. BatMeth2 med standardindstillingerne kørte hurtigere end de fleste af de offentliggjorte alignere og var sammenlignelig med BWA-meth og BatMeth. Bismark2 (med Bowtie2 som den grundlæggende kortlægningsmetode), BS Seeker2 og BSmooth kræver længere køretider.

DNA-methyleringsopkald

For at evaluere nøjagtigheden af DNA-methyleringsopkald blandt forskellige software hentede vi 450 K perlechipdata fra IMR90-cellinjen fra ENCODE (Encyclopedia of DNA Elements). Vi downloadede også helgenombisulfit-sekventeringsdata (WGBS-Seq) fra IMR90-cellinjen fra ENCODE (42,6 Gbaser). For hver software tilpassede vi WGBS-Seq-læsningerne og beregnede DNA-methyleringsniveauet. Derefter sammenlignede vi resultaterne med ML’erne på de samme steder i 450 K Bead Chip-dataene. Når forskellen mellem DNA ML fra WGBS-Seq-dataene fra softwaren og den fra 450 K Bead Chip var mindre end 0,2, blev kaldningsresultatet defineret som korrekt; ellers blev det betragtet som ukorrekt.

Resultaterne er vist i tabel 2. Overlapningen mellem de korrekte resultater fra alle programmerne er vist i Yderligere fil 1: Figur S2. Vi kan se, at BatMeth2 og Biscuit har lignende præstationer, som er bedre end den øvrige softwares præstationer. Konklusionen er, at BatMeth2 forbedrer nøjagtigheden af både BS-read alignment og DNA ML calling.

BatMeth2 justerer BS-reads og tillader samtidig indels af variabel længde

Kræft indeholder en markant højere andel af indels end sunde celler gør. For at verificere, om BatMeth2 kan tilpasse BS-reads med indels af forskellig længde, downloadede vi derfor WGBS-data (75 Gbaser) og 450 K Bead Chip-data fra HepG2 (lever hepatocellulært carcinom, en kræftcellelinje) fra ENCODE. Vi kontrollerede indel-længdefordelingen i læsningerne efter alignering af HepG2 WGBS-Seq-data fra HepG2. Yderligere fil 1: Figur S3A viser, at længderne af de detekterede indels hovedsageligt var fordelt i intervallet 1 bp~ 5 bp, og den længste indel var 40 bp lang. Ifølge vores statistikker indeholdt 2,3 % af alignment-readsene indels. Ud fra disse resultater ved vi, at BatMeth2 kan tilpasse læsninger med indels af forskellig længde.

Dernæst testede vi effekten af indel-detektion på DNA-methyleringsopkald. For BatMeth2 kørte vi to muligheder på HepG2-dataene: med og uden indel-detektion (dvs. indstillet -I-parameter i BatMeth2). Vi kørte også Bismark på WGBS-Seq-dataene fra HepG2 som en reference for DNA-methyleringsopkald med indel-detektion, fordi Bismark ikke har en funktion til opkald af indels. Vi sammenlignede kaldelsen af DNA-methylering i BatMeth2 og Bismark med kaldelsen fra 450 K Bead Chip-dataene. Resultaterne er vist i Additional file 1: Figur S3B, hvor “BatMeth2-noIndel” svarer til BatMeth2 uden indel-detektion. Vi kan se, at i fravær af indel-detektion var resultatet af BatMeth2 kun en smule bedre end resultatet af Bismark (med Bowtie1 som den grundlæggende kortlægningsmetode). Resultatet af BatMeth2 med indel-detektion var betydeligt bedre. Desuden kan vi se, at BatMeth2 kan påvise flere DNA-methyleringssteder end BatMeth2-noIndel og Bismark (Bowtie 1). For at forstå, hvorfor BatMeth2’s præstation med indel-detektion er bedre, definerede vi de methyleringssteder, der blev fundet af BatMeth2, som resultat A, mens de methyleringssteder, der blev fundet af BatMeth2-noIndel og Bismark, blev defineret som resultat B. Derefter lod vi mclA være de methyleringssteder, der optræder i Resultat A, men ikke i Resultat B. Vi observerede, at mclA omfattede 23.853 DNA-methyleringssteder og 15.048 (63 %) af de 23.853 steder, der var dækket af alignementerne af indel-reads, der blev kaldt af BatMeth2 med indel-detektion (se Yderligere fil 1: Figur S3C). Desuden fandt vi, at indel-raterne i Resultat A og Resultat B kun var henholdsvis 5 og 0 %. Derfor konkluderede vi, at nøjagtig indel-detektion kan forbedre DNA-methyleringskaldet.

Visualisering af DNA-methyleringsdata

BatMeth2 indeholder værktøjer til visualisering af methyleringsdata. For at illustrere visualiseringsfunktionerne i BatMeth2 downloadede vi (1) 117 Gbase af single-end reads fra den menneskelige H9-cellinje, (2) 105.2 Gbase af single-end reads fra den menneskelige IMR90-cellinje og (3) 12.6 Gbase af paired-end reads fra wild-type ris. For det første kan BatMeth2 visualisere cytosinmethyleringstætheden på kromosomniveau. Prikkerne i fig. 4a repræsenterer et glidende vindue på 100 kb med et trin på 50 kb. For at muliggøre visning af ML på individuelle CpG- eller ikke-CpG-steder i en genombrowser, leverer vi også filer i bed- og bigWig-formater (Fig. 4b). Ved at sammenligne med tætheden af gener og TE’er observerede vi, at ML var korreleret med TE-tætheden og var antikorreleret med gentætheden (Fig. 4c). Denne tendens er tidligere blevet observeret i ris .

For det andet kan BatMeth2 visualisere genernes ML’er. Mere præcist kan BatMeth2 visualisere ML’erne 2 kb opstrøms for genet, ved transkriptionens startsted (TSS), i genlegemet, ved transkriptionens slutsted (TES) og 2 kb nedstrøms for genlegemet. Ved sammenligning af regionerne opstrøms, kroppen og nedstrøms viser fig. 5a, at DNA-ML’en i genkroppen er højere end i promotorregionen. Sammenligner man alle fem regioner, er der tydeligvis en dal i TSS-regionen (fig. 5b). BatMeth2 kan også beregne ML-profiler omkring introner, exoner, intergeniske regioner og TE’er (Additional file 1: Figur S4). Derudover kan BatMeth2 give et varmekort over flere gener efter genregion til bekvem sammenligning af de samlede gen-ML’er for forskellige prøver (Fig. 5c).

For det tredje kan BatMeth2 visualisere fordelingen af DNA-methylering. Additional file 1: Figur S5A viser DNA-methyleringsfordelingerne i H9- og IMR90-cellelinjerne. I figuren er DNA ML opdelt i fem kategorier: methyleret (M: > 80 %), mellem delvist methyleret og methyleret (Mh: 60-80 %), delvist methyleret (H: 40-60 %), mellem ikke-methyleret og delvist methyleret (hU: 20-40 %) og ikke-methyleret (U: < 20 %). Som vist i Additional file 1: Figur S5A, var ML højere i H9-cellinjen i M-kategorien end i IMR90-cellinjen, især i CpG-konteksten. I CH-sekvenskonteksten er CpG-methylering den fremherskende form, men en betydelig andel af methylerede cytosiner findes på CpA-steder, mens ML er mindre end 40 %, især i H9-cellinjen (Yderligere fil 1: Figur S5B).

For det fjerde kan BatMeth2 analysere sammenhængen mellem genekspressionsniveauet og genpromotor-DNA ML. Vi illustrerede denne funktion ved hjælp af H9- og IMR90-cellelinjerne. Ekspressionsniveauerne for generne i H9 eller IMR90 blev inddelt i forskellige kategorier. Som vist i Yderligere fil 1: Figur S5C, udviste de højt udtrykte gener lavere ML’er i deres promotorregioner. Desuden inddelte vi ML’erne for genpromotorer i fem kategorier. Resultatet i Additional file 1: Figur S5D viser, at gener med promotorer med højere ML-værdier udviste lavere ekspressionsniveauer. Den negative korrelation mellem ekspressionen af pattedyrsgener og DNA-methylering af promotorer er kendt . Denne analyse viser yderligere nøjagtigheden af BatMeth2.

Finding af differentielt methylerede cytosiner og regioner (DMCs/DMRs)

Identifikation af differentielt methylerede cytosiner (DMCs) og differentielt methylerede regioner (DMRs) er et af de vigtigste mål i methyleringsdataanalyse. Selv om forskere lejlighedsvis er interesseret i at korrelere enkelte cytosinsteder til en fænotype , er DMR’er meget vigtige funktioner .

Første BS-Seq-undersøgelser profilerede celler uden indsamling af replikater. For sådanne datasæt brugte vi Fisher’s eksakte test til at skelne differentielt methylerede cytosiner (DMC’er). For BS-Seq-datasæt med replikater er den mest naturlige statistiske model til at kalde DMC’er beta-binomialfordeling . Vi ved, at en række softwareprogrammer kan udføre differentiel DNA-methyleringsdataanalyse, såsom methylKit (et differentielt analyseprogram, der kræver biologiske replikater) og Methy-Pipe (et differentielt analyseprogram uden biologisk replikation). Der findes imidlertid ingen omfattende pakke, der omfatter både kortlægning og differentiel methyleringsanalyse. Vi har derfor udviklet en pakke, der integrerer kortlægning med differentiel analyse. For at lette identifikationen af DMR’er fra bisulfitdata uden replikater integrerede vi Fisher’s exact-test for at udføre en hypotesetest. Når en prøve har to eller flere replikater, bruger vi beta-binomialfordelingen til at udføre differentiel methyleringsanalyse. Vi leverer også bed- eller bigWig-filer til listen over DMR’er. DMR’erne kan visualiseres i en genombrowser (Fig. 4b) med de genererede bed- eller bigWig-filer.

Som illustration viser Fig. 6a antallet af DMC’er og regioner i IMR90-cellinjen og i H9-cellinjen, som detekteret af BatMeth2 (p-værdi< 0,05, meth.diff > = 0,6). BatMeth2 kan visualisere, om CpG’er og DMC’er er beriget i visse regioner, såsom gen-, CDS-, intron-, intergeniske, UTR-, TE-, LTR-, LINE- og SINE-regioner. Figur 6b visualiserer andelene af DMC’er i forskellige genomiske regioner. Bortset fra de intergeniske regioner observerede vi ikke DMC-berigelse i nogen regioner.

En væsentlig andel af differentielt methylerede promotorer (DMP’er) indeholder indels

Vi ved, at indels og DNA-methylering spiller en vigtig rolle i vævsudvikling og sygdomme . Her undersøger vi forholdet mellem differentielt methylerede promotorer (DMP’er) og indels. Vi udførte denne undersøgelse ved hjælp af BS-Seq-læsninger i IMR90- og H9-cellelinjer. Vi justerede først BS-Seq-læsningerne ved hjælp af BatMeth2; derefter blev indels kaldt ved hjælp af BisSNP- og GATK-værktøjer. Efterfølgende definerede vi de indels, der kun forekommer i H9 eller IMR90, som cellelinjespecifikke indels.

Dernæst påviste vi 1384 DMP’er mellem H9 og IMR90 ved hjælp af BatMeth2 (p-værdi< 0,05, meth.diff > = 0,6). I alt 236 (17 %) blandt alle ovenstående DMP’er indeholder indels, som vist i fig. 6c. Kort sagt indeholder en betydelig del af DMP’erne indels. Derfor er en nøjagtig tilpasning af BS-Seq-reads i nærheden af disse indels meget vigtig for forskning og udforskning af DNA-methylering.