Un pacchetto integrato per l’analisi dei dati di metilazione del DNA bisolfita con mappatura sensibile all’indelebile

Un pacchetto facile da usare, autorun per le analisi di metilazione del DNA

Per completare più convenientemente l’analisi dei dati di metilazione del DNA, abbiamo confezionato tutte le funzioni in un pacchetto facile da usare, autorun per l’analisi di metilazione del DNA. La figura 1 mostra le caratteristiche principali di BatMeth2: 1) BatMeth2 ha prestazioni di allineamento efficienti e accurate. 2) BatMeth2 può calcolare il livello di metilazione del DNA (ML) di singoli siti di citosina o di qualsiasi regione funzionale, come interi cromosomi, regioni geniche, elementi trasponibili (TE), ecc. 3) Dopo l’integrazione di diversi algoritmi statistici, BatMeth2 può eseguire analisi differenziali di metilazione del DNA per qualsiasi regione, qualsiasi numero di campioni di input e requisiti dell’utente. 4) Integrando la visualizzazione dei dati BS-Seq (distribuzione della metilazione del DNA su cromosomi e geni) e l’annotazione della metilazione differenziale, BatMeth2 può visualizzare i dati di metilazione del DNA più chiaramente. Durante l’esecuzione dello strumento BatMeth2, viene generato un rapporto html per le statistiche del campione. I dettagli del report html del campione sono mostrati in http://htmlpreview.github.io/?https://github.com/GuoliangLi-HZAU/BatMeth2/blob/master/BatMeth2-Report/batmeth2.html.

BatMeth2 ha migliori prestazioni di mappatura su dati BS-Seq simulati

Abbiamo prima valutato tutti gli allineatori utilizzando set di dati simulati (senza indel) costituiti da letture con 75 coppie di basi (bp), 100 bp e 150 bp e con diversi tassi di conversione del bisolfito (che vanno da 0 a 100% con passo 10%). Questi set di dati sono stati simulati dal genoma umano (UCSC hg19) utilizzando FASTX-mutate-tools, wgsim (v0.3.0) e il simulatore in SAMtools (v1.1), che consente lo 0,03% di indel, un tasso di errore di base dell’1% nell’intero genoma e un massimo di due mismatch per lettura. Abbiamo mappato le letture simulate sul genoma di riferimento, consentendo al massimo due mismatch. Poiché le posizioni originali delle letture simulate erano note, abbiamo potuto valutare la precisione di tutti i programmi confrontando i loro risultati di mappatura con le posizioni originali.

Per confrontare le prestazioni dei diversi software, una lettura di sequenziamento con indelebili è stata considerata correttamente mappata se le seguenti condizioni erano vere: 1) la lettura è stata univocamente mappata allo stesso filamento da cui è stata simulata e la qualità della mappatura era maggiore di 0; 2) la posizione iniziale riportata della lettura allineata era entro dieci coppie di basi della posizione iniziale originale della lettura simulata; 3) i risultati della mappatura avevano indelebili o mismatch simili alla lettura simulata. Se una qualsiasi di queste condizioni è stata violata, la lettura è stata considerata erroneamente mappata. Poiché BatMeth2 permette un gap nella regione del seme, può trovare posizioni del seme che incorporano indelebili con alta precisione e può evitare posizioni non corrispondenti, che causerebbero il disallineamento delle letture che incorporano indelebili. I risultati in Fig. 2 mostrano che BatMeth2 ha raggiunto il maggior numero di letture correttamente allineate e il minor numero di letture erroneamente allineate in tutti i set di dati di prova a diversi tassi di conversione del bisolfito.

In breve, i risultati dei set di dati wgsim-simulati indel-aberranti mostrano che BatMeth2 ha prestazioni migliori (1~2% meglio del secondo miglior allineatore) rispetto agli altri metodi quando si allineano letture BS generali simulate contenenti una miscela di mismatches e indels. Possiamo vedere che con l’aumento del tasso di conversione BS, la precisione di allineamento di tutti i software si riduce. In queste diverse condizioni, BatMeth2 si comporta meglio.

BatMeth2 ha prestazioni di mappatura migliori su dati BS-Seq reali

Per testare le prestazioni di BatMeth2 su set di dati BS-Seq reali, abbiamo scaricato set di dati BS-Seq paired-end ed estratto a caso 1 milione di letture paired-end 2 × 90 bp da SRA SRR847318, 1 milione di 2 × 101 bp paired-end legge da SRA SRR1035722 e 1 milione di 2 × 125 bp paired-end legge da SRA SRR3503136 per scopi di valutazione. Poiché questi set di dati provengono da linee cellulari o tessuti sani, ci si aspetta che contengano un basso numero di variazioni strutturali. Quindi, abbiamo allineato questi dati reali utilizzando le letture single-end dai set di dati paired-end e valutato i tassi di mappatura concordanti e discordanti dagli allineamenti accoppiati per stimare i tassi di allineamento corretto ed errato. Poiché la dimensione dell’inserto delle letture paired-end era di circa 500 bp, una coppia di letture partner potrebbe essere considerata concordante se fosse mappata entro una distanza nominale di 500 bp; altrimenti, una coppia di letture partner potrebbe essere considerata discordante. Simile ai nostri risultati con i dati simulati, BatMeth2 ha riportato più allineamenti concordanti e meno discordanti sui set di dati reali su un’ampia gamma di punteggi di qualità della mappa, come mostrato in Fig. 3.

Inoltre, la tabella 1 mostra i tempi di esecuzione relativi dei programmi. BatMeth2 con le impostazioni predefinite è stato più veloce della maggior parte degli allineatori pubblicati ed era paragonabile a BWA-meth e BatMeth. Bismark2 (con Bowtie2 come metodo di mappatura fondamentale), BS Seeker2 e BSmooth richiedono tempi di esecuzione più lunghi.

Chiamata di metilazione del DNA

Per valutare l’accuratezza della chiamata di metilazione del DNA tra i diversi software, abbiamo scaricato 450 K di dati di chip di perline dalla linea cellulare IMR90 da ENCODE (Encyclopedia of DNA Elements). Abbiamo anche scaricato i dati di sequenziamento del bisolfito dell’intero genoma (WGBS-Seq) della linea cellulare IMR90 da ENCODE (42.6 Gbases). Per ogni software, abbiamo allineato le letture WGBS-Seq e calcolato il livello di metilazione del DNA. Poi, abbiamo confrontato i risultati con le ML degli stessi siti nei dati 450 K Bead Chip. Quando la differenza tra la ML del DNA dai dati WGBS-Seq dal software e quella dal 450 K Bead Chip era inferiore a 0,2, il risultato della chiamata è stato definito come corretto; altrimenti, è stato considerato errato.

I risultati sono mostrati nella tabella 2. La sovrapposizione tra i risultati corretti di tutti i software è mostrata nel file aggiuntivo 1: Figura S2. Possiamo vedere che BatMeth2 e Biscuit hanno prestazioni simili, che sono migliori di quelle degli altri software. In conclusione, BatMeth2 migliora l’accuratezza sia dell’allineamento BS-read che della chiamata ML del DNA.

BatMeth2 allinea le letture BS consentendo al contempo gli indel di lunghezza variabile

Il cancro contiene una percentuale notevolmente più alta di indel rispetto alle cellule sane. Pertanto, per verificare se BatMeth2 può allineare le letture BS con indelebili di diversa lunghezza, abbiamo scaricato i dati WGBS (75 Gbases) e 450 K Bead Chip da HepG2 (carcinoma epatocellulare, una linea cellulare di cancro) da ENCODE. Abbiamo controllato la distribuzione della lunghezza indel nelle letture dopo l’allineamento dei dati HepG2 WGBS-Seq. File aggiuntivo 1: Figura S3A mostra che le lunghezze degli indel rilevati sono stati principalmente distribuiti nella gamma 1 bp ~ 5 bp, e l’indel più lungo era di 40 bp di lunghezza. Secondo le nostre statistiche, il 2,3% delle letture di allineamento conteneva indelebili. Da questi risultati, sappiamo che BatMeth2 può allineare letture con indel di diverse lunghezze.

In seguito, abbiamo testato l’effetto del rilevamento degli indel sulla chiamata della metilazione del DNA. Per BatMeth2, abbiamo eseguito due opzioni sui dati HepG2: con e senza rilevamento degli indel (cioè, impostare il parametro -I in BatMeth2). Abbiamo anche eseguito Bismark sui dati WGBS-Seq da HepG2 come riferimento per la chiamata di metilazione del DNA con rilevazione degli indelebili, perché Bismark non ha una funzione di chiamata degli indelebili. Abbiamo confrontato la chiamata di metilazione del DNA in BatMeth2 e Bismark con la chiamata dai dati 450 K Bead Chip. I risultati sono mostrati nel file aggiuntivo 1: Figura S3B, dove “BatMeth2-noIndel” corrisponde a BatMeth2 senza rilevazione di indel. Possiamo vedere che, in assenza di rilevamento indel, il risultato di BatMeth2 era solo leggermente migliore di quello di Bismark (con Bowtie1 come metodo di mappatura fondamentale). Il risultato di BatMeth2 con il rilevamento degli indel era significativamente migliore. Inoltre, possiamo vedere che BatMeth2 può rilevare più siti di metilazione del DNA rispetto a BatMeth2-noIndel e Bismark (Bowtie 1). Per capire perché la performance di BatMeth2 con il rilevamento degli indel è migliore, abbiamo definito i siti di metilazione chiamati da BatMeth2 come Risultato A, mentre i siti di metilazione chiamati da BatMeth2-noIndel e Bismark sono stati definiti come Risultato B. Poi, abbiamo lasciato che mclA fosse i siti di metilazione che appaiono nel risultato A ma non nel risultato B. Abbiamo osservato che mclA includeva 23.853 siti di metilazione del DNA e 15.048 (63%) dei 23.853 siti coperti dagli allineamenti delle letture indel chiamate da BatMeth2 con rilevamento indel (vedi file aggiuntivo 1: Figura S3C). Inoltre, abbiamo trovato che i tassi di indel nel risultato A e nel risultato B erano solo 5 e 0%, rispettivamente. Quindi, abbiamo concluso che il rilevamento accurato degli indel può migliorare la metilazione del DNA.

Visualizzazione dei dati di metilazione del DNA

BatMeth2 fornisce strumenti per visualizzare i dati di metilazione. Per illustrare le caratteristiche di visualizzazione di BatMeth2, abbiamo scaricato (1) 117 Gbases di letture single-end dalla linea cellulare umana H9, (2) 105.2 Gbases di letture single-end dalla linea cellulare umana IMR90 e (3) 12.6 Gbases di letture paired-end dal riso wild-type. In primo luogo, BatMeth2 può visualizzare la densità di metilazione della citosina a livello cromosomico. I punti in Fig. 4a rappresentano una finestra scorrevole di 100 kb con un passo di 50 kb. Per consentire la visualizzazione della ML in singoli siti CpG o non-CpG in un browser del genoma, forniamo anche file in formato bed e bigWig (Fig. 4b). Confrontando con la densità di geni e TE, abbiamo osservato che la ML era correlata con la densità TE ed era anticorrelata con la densità del gene (Fig. 4c). Questa tendenza è stata precedentemente osservata nel riso.

In secondo luogo, BatMeth2 può visualizzare le ML dei geni. Più precisamente, BatMeth2 può visualizzare le MLs 2 kb a monte del gene, al sito di inizio della trascrizione (TSS), nel corpo del gene, al sito finale della trascrizione (TES) e 2 kb a valle del corpo del gene. Confrontando le regioni a monte, il corpo e a valle, la Fig. 5a mostra che il DNA ML del corpo del gene è superiore a quello della regione del promotore. Confrontando tutte e cinque le regioni, c’è ovviamente una valle nella regione TSS (Fig. 5b). BatMeth2 può anche calcolare i profili ML intorno agli introni, esoni, regioni intergeniche e TEs (File aggiuntivo 1: Figura S4). Inoltre, BatMeth2 può fornire una mappa di calore di più geni per regione genica per un comodo confronto delle MLs geniche complessive di diversi campioni (Fig. 5c).

In terzo luogo, BatMeth2 può visualizzare la distribuzione della metilazione del DNA. Additional file 1: Figura S5A mostra le distribuzioni di metilazione del DNA nelle linee cellulari H9 e IMR90. Nella figura, il DNA ML è suddiviso in cinque categorie: metilato (M: > 80%), intermedio tra parzialmente metilato e metilato (Mh: 60-80%), parzialmente metilato (H: 40-60%), intermedio tra non metilato e parzialmente metilato (hU: 20-40%), e non metilato (U: < 20%). Come mostrato nel file aggiuntivo 1: Figura S5A, il ML era più alto nella linea cellulare H9 nella categoria M che nella linea cellulare IMR90, soprattutto nel contesto CpG. Nel contesto della sequenza CH, la metilazione CpG è la forma predominante, ma una frazione significativa di citosine metilate si trova nei siti CpA, mentre la ML è inferiore al 40%, in particolare nella linea cellulare H9 (Additional file 1: Figura S5B).

In quarto luogo, BatMeth2 può analizzare la correlazione tra il livello di espressione genica e la ML del DNA promotore del gene. Abbiamo illustrato questa caratteristica usando le linee cellulari H9 e IMR90. I livelli di espressione dei geni in H9 o IMR90 sono stati divisi in diverse categorie. Come mostrato nel file aggiuntivo 1: Figura S5C, i geni altamente espressi hanno esibito MLs inferiori nelle loro regioni promotrici. Inoltre, abbiamo diviso le MLs dei promotori dei geni in cinque categorie. Il risultato nel file aggiuntivo 1: Figura S5D mostra che i geni con promotori che hanno valori ML più alti hanno esibito livelli di espressione più bassi. La correlazione negativa tra l’espressione dei geni dei mammiferi e la metilazione del DNA del promotore è nota. Questa analisi indica ulteriormente la precisione di BatMeth2.

Ricercare le citochine e le regioni metilate differentemente (DMCs/DMRs)

L’identificazione delle citochine metilate differentemente (DMCs) e delle regioni metilate differentemente (DMRs) è uno degli obiettivi principali nell’analisi dei dati di metilazione. Anche se i ricercatori sono occasionalmente interessati a correlare singoli siti di citosina a un fenotipo, le DMR sono caratteristiche molto importanti.

I primi studi BS-Seq hanno profilato le cellule senza raccogliere repliche. Per tali set di dati, abbiamo usato il test esatto di Fisher per discernere le citochine metilate in modo diverso (DMC). Per BS-Seq set di dati con repliche, il modello statistico più naturale per chiamare DMCs è beta-binomiale distribuzione. Sappiamo che un certo numero di programmi software in grado di eseguire l’analisi differenziale dei dati di metilazione del DNA, come methylKit (un programma di analisi differenziale che richiede repliche biologiche) e Methy-Pipe (un programma di analisi differenziale senza duplicazione biologica). Tuttavia, non è disponibile alcun pacchetto completo che includa sia la mappatura che l’analisi differenziale della metilazione. Così, abbiamo sviluppato un pacchetto che integra la mappatura con l’analisi differenziale. Per facilitare l’identificazione di DMR da dati bisolfiti senza repliche, abbiamo integrato il test esatto di Fisher per eseguire un test di ipotesi. Quando un campione ha due o più repliche, usiamo la distribuzione beta-binomiale per eseguire l’analisi di metilazione differenziale. Forniamo anche file bed o bigWig per l’elenco dei DMR. I DMR possono essere visualizzati in un browser del genoma (Fig. 4b) con i file bed o bigWig generati.

Come illustrazione, la Fig. 6a mostra il numero di DMC e le regioni nella linea cellulare IMR90 e nella linea cellulare H9, come rilevato da BatMeth2 (valore p< 0.05, meth.diff > = 0.6). BatMeth2 può visualizzare se CpGs e DMCs sono arricchiti in alcune regioni, come il gene, CDS, introne, intergenico, UTR, TE, LTR, LINE e SINE regioni. Figura 6b visualizza le proporzioni di DMCs in diverse regioni genomiche. A parte le regioni intergeniche, non abbiamo osservato arricchimento DMC in nessuna regione.

Una proporzione sostanziale di promotori differenzialmente metilati (DMPs) contengono indel

Sappiamo che indel e metilazione del DNA svolgono un ruolo importante nello sviluppo dei tessuti e nelle malattie. Qui, esaminiamo la relazione tra i promotori differenzialmente metilati (DMPs) e gli indel. Abbiamo eseguito questo studio utilizzando le letture BS-Seq nelle linee cellulari IMR90 e H9. Abbiamo prima allineato le letture BS-Seq usando BatMeth2; poi, gli indel sono stati chiamati usando gli strumenti BisSNP e GATK. Successivamente, abbiamo definito gli indel che si verificano solo in H9 o IMR90 come indel specifici della linea cellulare.

Poi, abbiamo rilevato 1384 DMPs tra H9 e IMR90 da BatMeth2 (p value< 0.05, meth.diff > = 0.6). Un totale di 236 (17%) tra tutti i DMPs di cui sopra contengono indel, come mostrato in Fig. 6c. In breve, una parte sostanziale dei DMPs contiene indel. Pertanto, l’allineamento accurato di BS-Seq legge vicino a questi indel è molto importante per la ricerca e l’esplorazione della metilazione del DNA.