Ett lättanvänt autorun-paket för DNA-metyleringsanalyser
För att göra DNA-metyleringsdataanalysen mer bekväm har vi paketerat alla funktioner i ett lättanvänt autorun-paket för DNA-metyleringsanalys. Figur 1 visar huvudfunktionerna i BatMeth2: 1) BatMeth2 har en effektiv och noggrann anpassningsprestanda. 2) BatMeth2 kan beräkna DNA-metyleringsnivån (ML) för enskilda cytosinplatser eller alla funktionella regioner, t.ex. hela kromosomer, genregioner, transposerbara element (TE) osv. 3) Efter integrering av olika statistiska algoritmer kan BatMeth2 utföra differentiell DNA-metyleringsanalys för vilken region som helst, vilket antal ingångsprover som helst och användarens krav. 4) Genom att integrera BS-Seq-datavisualisering (DNA-metyleringsfördelning på kromosomer och gener) och differentiell metyleringsannotation kan BatMeth2 visualisera DNA-metyleringsdata tydligare. Under utförandet av BatMeth2-verktyget genereras en html-rapport för statistiken för provet. Detaljer om exempel på html-rapport visas i http://htmlpreview.github.io/?https://github.com/GuoliangLi-HZAU/BatMeth2/blob/master/BatMeth2-Report/batmeth2.html.
Arbetsflödet för BatMeth2. De två stora pilarna betyder inmatnings- eller utdatafiler
BatMeth2 har bättre kartläggningsprestanda på simulerade BS-Seq-data
Vi utvärderade först alla aligners med hjälp av simulerade dataset (utan indels) som består av reads med 75 baspar (bp), 100 bp och 150 bp och med olika bisulfitkonverteringshastigheter (från 0 till 100 % med steg 10 %). Dessa dataset simulerades från det mänskliga genomet (UCSC hg19) med hjälp av FASTX-mutate-tools , wgsim (v0.3.0) och simulatorn i SAMtools (v1.1) , som tillåter 0,03 % indels, en felfrekvens på 1 % baser i hela genomet och högst två felmatchningar per läsning. Vi mappade de simulerade läsningarna till referensgenomet och tillät högst två felmatchningar. Eftersom de simulerade läsningarnas ursprungliga positioner var kända kunde vi utvärdera noggrannheten hos alla program genom att jämföra deras mappningsresultat med de ursprungliga positionerna.
För att jämföra de olika programvarornas prestanda ansågs en sekvenseringsavläsning med indels vara korrekt mappad om följande villkor var uppfyllda: 1) Avläsningen mappades unikt till samma sträng som den simulerades från och mappningskvaliteten var större än 0. 2) Den rapporterade startpositionen för den anpassade avläsningen låg inom tio baspar av den simulerade avläsningens ursprungliga startposition. 3) Mappningsresultaten hade liknande indels eller mismatchningar som den simulerade avläsningen. Om något av dessa villkor inte uppfylldes ansågs läsningen vara felaktigt mappad. Eftersom BatMeth2 tillåter en lucka i fröregionen kan den hitta fröplatser som innehåller indels med hög noggrannhet och kan undvika felmatchade platser, vilket skulle leda till att läsningar som innehåller indels blir feljusterade. Resultaten i figur 2 visar att BatMeth2 uppnådde det största antalet korrekt anpassade läsningar och det lägsta antalet felaktigt anpassade läsningar i alla testdataset vid olika bisulfitkonverteringshastigheter.
Utvärdering av alla BS-Seq-justerare med hjälp av simulerade dataset med olika läsningslängder från FASTX och wgsim. Simulerade data med olika bisulfitkonverteringshastigheter visas i olika former. Resultat från olika aligners visas med olika färger på symbolerna. Resultaten nära det övre vänstra hörnet i varje panel visar att programvaran uppnådde ett större antal korrekt mappade läsningar och ett lägre antal felaktigt mappade läsningar. Resultaten från vår alignerare BatMeth2 är bäst i de olika simulerade bisulfit-datasetterna
I korthet visar resultaten från wgsim-simulerade indel-aberrant-datasetterna att BatMeth2 har bättre prestanda (1~2 % bättre än den näst bästa aligneraren) än de andra metoderna när man aligariserar generella simulerade BS-reads som innehåller en blandning av felmatchningar och indels. Vi kan se att med den ökade BS-konverteringsfrekvensen minskar anpassningsnoggrannheten för alla programvarorna. Under dessa olika förhållanden presterar BatMeth2 bättre.
BatMeth2 har bättre mappningsprestanda på riktiga BS-Seq-data
För att testa BatMeth2:s prestanda på riktiga BS-Seq-datasatser laddade vi ner parvisa BS-Seq-datasatser och extraherade slumpmässigt 1 miljon 2 × 90 bp parvisa reads från SRA SRR847318, 1 miljon 2 × 101 bp parade slutresultat från SRA SRR1035722 och 1 miljon 2 × 125 bp parade slutresultat från SRA SRR3503136 för utvärderingsändamål. Eftersom dessa dataset kommer från friska cellinjer eller vävnader förväntas de innehålla ett lågt antal strukturella variationer. Därför anpassade vi dessa riktiga data med hjälp av single-end reads från de parvisa datamängderna och utvärderade de överensstämmande och avvikande kartläggningsfrekvenserna från de parvisa anpassningarna för att uppskatta frekvensen för korrekt och felaktig anpassning. Eftersom insättningsstorleken för de parade ändläsningarna var ungefär 500 bp kunde ett par partnerläsningar anses vara överensstämmande om de mappades inom ett nominellt avstånd på 500 bp, annars kunde ett par partnerläsningar anses vara diskordanta. I likhet med våra resultat med de simulerade uppgifterna rapporterade BatMeth2 fler överensstämmande och färre diskordanta anpassningar på de verkliga datamängderna över ett stort intervall av kvalitetspoäng för kartor, vilket visas i figur 3.
Samstämmighets- och diskordansfrekvenser för anpassningar på verkliga parvisade läsuppgifter från olika anpassningsverktyg. Kumulativt antal concordant och discordant alignments från hög till låg kartläggningskvalitet för riktiga bisulfitsekvenseringsavläsningar. Det finns endast en poäng för BSmap och de aligners som baseras på bowtie separat, eftersom dessa aligners inte har någon poäng för kartkvalitet. Bismark-bowtie2L15 betyder bowtie2-justering med frölängd 15
I tabell 1 visas dessutom de relativa körtiderna för programmen. BatMeth2 med standardinställningarna körde snabbare än de flesta publicerade aligners och var jämförbar med BWA-meth och BatMeth. Bismark2 (med Bowtie2 som grundläggande kartläggningsmetod), BS Seeker2 och BSmooth kräver längre körtider.
DNA-metyleringskallelse
För att utvärdera noggrannheten av DNA-metyleringskallelse bland olika programvaror laddade vi ner 450 K bead-chipdata från cellinjen IMR90 från ENCODE (Encyclopedia of DNA Elements). Vi laddade också ner helgenombisulfitsekvensering (WGBS-Seq) av data från cellinjen IMR90 från ENCODE (42,6 Gbases). För varje programvara anpassade vi WGBS-Seq-avläsningarna och beräknade nivån av DNA-metylering. Därefter jämförde vi resultaten med MLs på samma platser i 450 K Bead Chip-data. När skillnaden mellan DNA ML från WGBS-Seq-data från programvaran och DNA ML från 450 K Bead Chip var mindre än 0,2 definierades callresultatet som korrekt, annars ansågs det vara felaktigt.
Resultaten visas i tabell 2. Överlappningen mellan de korrekta resultaten från alla programvaror visas i tilläggsfil 1: Figur S2. Vi kan se att BatMeth2 och Biscuit har liknande prestationer, som är bättre än de andra programvarornas. Sammanfattningsvis förbättrar BatMeth2 noggrannheten hos både BS-read alignment och DNA ML calling.
BatMeth2 alignar BS-reads samtidigt som det tillåter indels av varierande längd
Cancer innehåller en anmärkningsvärt högre andel indels än vad friska celler gör. För att kontrollera om BatMeth2 kan anpassa BS-reads med indels av olika längd laddade vi därför ner WGBS-data (75 Gbases) och 450 K Bead Chip-data från HepG2 (leverhepatocellulärt karcinom, en cancercellinje) från ENCODE. Vi kontrollerade indel-längdfördelningen i läsningarna efter anpassningen av HepG2 WGBS-Seq-data från HepG2. Ytterligare fil 1: Figur S3A visar att längden på de upptäckta indelerna huvudsakligen fördelades i intervallet 1 bp~ 5 bp, och det längsta indelet var 40 bp långt. Enligt vår statistik innehöll 2,3 % av de avläsningar som gjordes vid anpassningen indels. Från dessa resultat vet vi att BatMeth2 kan anpassa läsningar med indels av olika längd.
Nästan testade vi effekten av indeldetektering på DNA-metylering. För BatMeth2 körde vi två alternativ på HepG2-data: med och utan indeldetektion (dvs. inställd -I-parameter i BatMeth2). Vi körde också Bismark på WGBS-Seq-data från HepG2 som en referens för DNA-metyleringskallning med indeldetektion, eftersom Bismark inte har någon indel-kallningsfunktion. Vi jämförde DNA-metyleringskallelsen i BatMeth2 och Bismark med kallelsen från 450 K Bead Chip-data. Resultaten visas i tilläggsfil 1: Figur S3B, där ”BatMeth2-noIndel” motsvarar BatMeth2 utan indel-detektion. Vi kan se att i avsaknad av indeldetektering var BatMeth2:s resultat endast något bättre än Bismarks (med Bowtie1 som grundläggande kartläggningsmetod). Resultatet av BatMeth2 med indeldetektion var betydligt bättre. Dessutom kan vi se att BatMeth2 kan upptäcka fler DNA-metyleringsställen än BatMeth2-noIndel och Bismark (Bowtie 1). För att förstå varför BatMeth2:s resultat med indeldetektion är bättre definierade vi de metyleringsställen som BatMeth2 identifierade som resultat A, medan de metyleringsställen som BatMeth2-noIndel och Bismark identifierade definierades som resultat B. Därefter lät vi mclA vara de metyleringsställen som förekommer i resultat A men inte i resultat B. Vi observerade att mclA omfattade 23 853 DNA-metyleringsställen och 15 048 (63 %) av de 23 853 ställen som täcktes av anpassningarna av indel-avläsningar som kallades av BatMeth2 med indel-detektion (se tilläggsfil 1: figur S3C). Dessutom fann vi att indelfrekvenserna i resultat A och resultat B endast var 5 respektive 0 %. Därför drog vi slutsatsen att noggrann indeldetektion kan förbättra DNA-metyleringskallandet.
Visualisering av DNA-metyleringsdata
BatMeth2 tillhandahåller verktyg för att visualisera metyleringsdata. För att illustrera visualiseringsfunktionerna i BatMeth2 laddade vi ner (1) 117 Gbase av single-end reads från den mänskliga H9-cellinjen, (2) 105,2 Gbase av single-end reads från den mänskliga IMR90-cellinjen och (3) 12,6 Gbase av paired-end reads från vildtypris. För det första kan BatMeth2 visualisera cytosinmetyleringstätheten på kromosomnivå. Punkterna i figur 4a representerar ett glidfönster på 100 kb med ett steg på 50 kb. För att möjliggöra visning av ML på enskilda CpG- eller icke-CpG-platser i en genomwebbläsare tillhandahåller vi också filer i formaten bed och bigWig (fig. 4b). Genom att jämföra med tätheten av gener och TE:er observerade vi att ML var korrelerad med TE-tätheten och var antikorrelerad med gentätheten (fig. 4c). Denna tendens har tidigare observerats i ris .
Visualisering av metyleringsnivåerna i kromosomskalan. a Metylcytosintätheten i människans kromosom 10. Punkterna representerar metyleringsnivåerna i glidande fönster på 100 Kb med ett steg på 50 Kb. De röda prickarna avser metyleringsnivåerna i plussträngen och de blå prickarna avser metyleringsnivåerna i minussträngen. b Ett exempel på fördelningen av DNA-metyleringsnivåerna och differentiellt metylerade regioner (DMR) mellan cellinjerna H9 och IMR90 i människans kromosom 10. c Tätheten av gener, transposonelement (TEs) och DNA-metyleringsnivån i hela risgenomet. Panel A är resultaten från Batmeth2. Panel B är visualiseringsresultaten från UCSC browser, med BED-filer från Batmeth2
För det andra kan Batmeth2 visualisera genernas MLs. Mer exakt kan BatMeth2 visualisera MLs 2 kb uppströms från genen, vid transkriptionsstartplatsen (TSS), i genkroppen, vid transkriptionsslutplatsen (TES) och 2 kb nedströms från genkroppen. Om man jämför regionerna uppströms, kroppen och nedströms visar figur 5a att DNA ML i genkroppen är högre än i promotorregionen. Om man jämför alla fem regionerna finns det uppenbarligen en dal i TSS-regionen (fig. 5b). BatMeth2 kan också beräkna ML-profiler runt introner, exoner, intergeniska regioner och TEs (Ytterligare fil 1: Figur S4). Dessutom kan BatMeth2 tillhandahålla en värmekarta över flera gener per genregion för en bekväm jämförelse av de övergripande gen-MLs för olika prover (fig. 5c).
Visualisering av DNA-metylering under olika kontexter. a DNA-metyleringsnivåerna i 2Kb-regioner uppströms gener, genkroppar, 2Kb nedströms genkroppar. b Aggregeringsprofilen av DNA-metyleringsnivåer över gener. c Värmekartan över alla gener i 2Kb-regioner uppströms gener, genkroppar, 2Kb nedströms genkroppar
För det tredje kan BatMeth2 visualisera fördelningen av DNA-metylering. Additional file 1: Figur S5A visar fördelningen av DNA-metylering i cellinjerna H9 och IMR90. I figuren är DNA ML uppdelat i fem kategorier: metylerat (M: > 80 %), mellan delvis metylerat och metylerat (Mh: 60-80 %), delvis metylerat (H: 40-60 %), mellan icke metylerat och delvis metylerat (hU: 20-40 %) och icke metylerat (U: < 20 %). Som framgår av tilläggsfil 1: Figur S5A var ML högre i H9-cellinjen i kategorin M än i IMR90-cellinjen, särskilt i CpG-kontexten. I CH-sekvenskontexten är CpG-metylering den dominerande formen, men en betydande andel metylerade cytosiner finns på CpA-platser, medan ML är mindre än 40 %, särskilt i H9-cellinjen (Additional file 1: Figur S5B).
För det fjärde kan BatMeth2 analysera korrelationen mellan genuttrycksnivå och ML i genpromotor-DNA. Vi illustrerade denna funktion med hjälp av cellinjerna H9 och IMR90. Uttrycksnivåerna för generna i H9 eller IMR90 delades in i olika kategorier. Såsom visas i tilläggsfil 1: Figur S5C uppvisade de högt uttryckta generna lägre MLs i sina promotorområden. Dessutom delade vi in MLs för genernas promotorer i fem kategorier. Resultatet i Additional file 1: Figur S5D visar att gener med promotorer med högre ML-värden uppvisade lägre uttrycksnivåer. Den negativa korrelationen mellan uttrycket av däggdjursgener och promotorernas DNA-metylering är känd . Denna analys visar ytterligare på BatMeth2:s noggrannhet.
Finnande av differentiellt metylerade cytosiner och regioner (DMCs/DMRs)
Identifiering av differentiellt metylerade cytosiner (DMCs) och differentiellt metylerade regioner (DMRs) är ett av de viktigaste målen vid analys av metyleringsdata. Även om forskare ibland är intresserade av att korrelera enskilda cytosinplatser med en fenotyp är DMR mycket viktiga egenskaper.
I tidiga BS-Seq-studier profilerade man celler utan att samla in replikat. För sådana dataset använde vi Fishers exakta test för att urskilja differentiellt metylerade cytosiner (DMC). För BS-Seq-dataset med replikat är den naturligaste statistiska modellen för att identifiera DMC:er beta-binomialfördelning . Vi vet att ett antal mjukvaruprogram kan utföra differentiell DNA-metyleringsdataanalys, t.ex. methylKit (ett program för differentiell analys som kräver biologiska replikat) och Methy-Pipe (ett program för differentiell analys utan biologiska replikat). Det finns dock inget heltäckande paket som omfattar både kartläggning och differentiell metyleringsanalys. Därför har vi utvecklat ett paket som integrerar kartläggning med differentiell analys. För att underlätta identifieringen av DMR från bisulfitdata utan replikat integrerade vi Fishers exakta test för att utföra ett hypotestest. När ett prov har två eller flera replikat använder vi beta-binomialfördelningen för att utföra differentiell metyleringsanalys. Vi tillhandahåller också bed- eller bigWig-filer för listan över DMRs. DMR:erna kan visualiseras i en genombrowser (fig. 4b) med de genererade bed- eller bigWig-filerna.
Som illustration visar fig. 6a antalet DMC:er och regioner i IMR90-cellinjen och i H9-cellinjen, som detekterats av BatMeth2 (p-värde< 0,05, meth.diff > = 0,6). BatMeth2 kan visualisera om CpGs och DMCs är berikade i vissa regioner, t.ex. gen-, CDS-, intron-, intergeniska-, UTR-, TE-, LTR-, LINE- och SINE-regioner. Figur 6b visualiserar proportionerna av DMCs i olika genomiska regioner. Bortsett från de intergeniska regionerna observerade vi ingen DMC-anrikning i några regioner.
Differentiell metyleringsanalys. a Analysresultat av differentiellt metylerade regioner (DMR), differentiellt metylerade gener (DMG) och differentiellt metylerade promotorer (DMP) mellan cellinjerna H9 och IMR90. b Annotering av differentiellt metylerade cytosiner (DMC) mot olika genomiska egenskaper och upprepade element. c DMPs innehåller H9- eller IMR90-specifika indels (orange) utgör en betydande andel av alla DMPs (DNA Methylation differential Promoters)
En betydande andel av differentiellt metylerade promotorer (DMPs) innehåller indels
Vi vet att indels och DNA-metylering spelar en viktig roll för vävnadsutveckling och sjukdomar . Här undersöker vi förhållandet mellan differentiellt metylerade promotorer (DMP) och indels. Vi utförde den här studien med hjälp av BS-Seq-läsningar i cellinjerna IMR90 och H9. Vi anpassade först BS-Seq-avläsningarna med hjälp av BatMeth2; därefter identifierades indels med hjälp av BisSNP- och GATK-verktygen. Därefter definierade vi de indels som endast förekommer i H9 eller IMR90 som cellinjespecifika indels.
Därefter upptäckte vi 1384 DMPs mellan H9 och IMR90 med hjälp av BatMeth2 (p-värde< 0,05, meth.diff > = 0,6). Totalt 236 (17 %) av alla DMP ovan innehåller indels, vilket visas i fig. 6c. Kort sagt innehåller en betydande andel av DMP:erna indels. Därför är noggrann anpassning av BS-Seq-avläsningar nära dessa indels mycket viktig för forskning och utforskning av DNA-metylering.