Ein einfach zu bedienendes Autorun-Paket für DNA-Methylierungsanalysen
Um die Analyse von DNA-Methylierungsdaten komfortabler zu gestalten, haben wir alle Funktionen in ein einfach zu bedienendes Autorun-Paket für DNA-Methylierungsanalysen gepackt. Abbildung 1 zeigt die Hauptmerkmale von BatMeth2: 1) BatMeth2 hat eine effiziente und genaue Alignment-Leistung. 2) BatMeth2 kann den DNA-Methylierungsgrad (ML) einzelner Cytosinstellen oder beliebiger funktioneller Regionen, wie ganzer Chromosomen, Genregionen, transponierbarer Elemente (TEs) usw., berechnen. 3) Nach der Integration verschiedener statistischer Algorithmen kann BatMeth2 differenzielle DNA-Methylierungsanalysen für jede Region, jede Anzahl von Eingangsproben und Benutzeranforderungen durchführen. 4) Durch die Integration der BS-Seq-Datenvisualisierung (DNA-Methylierungsverteilung auf Chromosomen und Genen) und der differenziellen Methylierungsannotation kann BatMeth2 die DNA-Methylierungsdaten deutlicher visualisieren. Während der Ausführung des BatMeth2-Tools wird ein HTML-Bericht für die Statistik der Probe erstellt. Die Details des html-Berichts sind in http://htmlpreview.github.io/?https://github.com/GuoliangLi-HZAU/BatMeth2/blob/master/BatMeth2-Report/batmeth2.html dargestellt.
Der Arbeitsablauf von BatMeth2. Die beiden großen Pfeile bedeuten Eingabe- bzw. Ausgabedateien
BatMeth2 hat eine bessere Mapping-Leistung bei simulierten BS-Seq-Daten
Wir haben zunächst alle Aligner anhand simulierter Datensätze (ohne Indels) bewertet, die aus Reads mit 75 Basenpaaren (bp), 100 bp und 150 bp und mit unterschiedlichen Bisulfit-Umwandlungsraten (von 0 bis 100 % mit einem Schritt von 10 %) bestanden. Diese Datensätze wurden anhand des menschlichen Genoms (UCSC hg19) mit FASTX-mutate-tools , wgsim (v0.3.0) und dem Simulator in SAMtools (v1.1) simuliert, der 0,03 % Indels, eine Basenfehlerrate von 1 % im gesamten Genom und maximal zwei Mismatches pro Read zulässt. Wir haben die simulierten Reads auf das Referenzgenom gemappt und dabei maximal zwei Mismatches zugelassen. Da die ursprünglichen Positionen der simulierten Reads bekannt waren, konnten wir die Genauigkeit aller Programme bewerten, indem wir ihre Mapping-Ergebnisse mit den ursprünglichen Positionen verglichen.
Um die Leistungen der verschiedenen Programme zu vergleichen, wurde ein Sequenzierungs-Read mit Indels als korrekt gemappt betrachtet, wenn die folgenden Bedingungen erfüllt waren: 1) der Read wurde eindeutig auf denselben Strang gemappt, von dem aus er simuliert wurde, und die Mapping-Qualität war größer als 0; 2) die gemeldete Startposition des alignierten Reads lag innerhalb von zehn Basenpaaren der ursprünglichen Startposition des simulierten Reads; 3) die Mapping-Ergebnisse wiesen ähnliche Indels oder Mismatches auf wie der simulierte Read. Wenn eine dieser Bedingungen nicht erfüllt war, wurde der Read als falsch gemappt betrachtet. Da BatMeth2 eine Lücke in der Seed-Region zulässt, kann es Seed-Positionen, die Indels enthalten, mit hoher Genauigkeit finden und nicht übereinstimmende Positionen vermeiden, die dazu führen würden, dass Reads, die Indels enthalten, falsch ausgerichtet werden. Die Ergebnisse in Abb. 2 zeigen, dass BatMeth2 in allen Testdatensätzen bei unterschiedlichen Bisulfit-Konversionsraten die größte Anzahl korrekt ausgerichteter Reads und die geringste Anzahl falsch ausgerichteter Reads erzielte.
Evaluierung aller BS-Seq-Aligner anhand simulierter Datensätze mit unterschiedlichen Leselängen aus FASTX und wgsim. Simulierte Daten mit unterschiedlichen Bisulfit-Konversionsraten sind in verschiedenen Formen dargestellt. Die Ergebnisse der verschiedenen Aligner sind durch unterschiedliche Farben der Symbole dargestellt. Die Ergebnisse in der linken oberen Ecke jedes Feldes zeigen, dass die Software eine größere Anzahl korrekt zugeordneter Reads und eine geringere Anzahl falsch zugeordneter Reads erzielt hat. Die Ergebnisse unseres Aligners BatMeth2 sind die besten in den verschiedenen simulierten Bisulfit-Datensätzen
Kurz gesagt zeigen die Ergebnisse von wgsim-simulierten Indel-Aberrant-Datensätzen, dass BatMeth2 eine bessere Leistung (1~2% besser als der zweitbeste Aligner) als die anderen Methoden hat, wenn es um das Alignment von allgemeinen simulierten BS-Reads geht, die eine Mischung aus Mismatches und Indels enthalten. Wir können sehen, dass mit der erhöhten BS-Konvertierungsrate die Alignment-Genauigkeit aller Software reduziert wird. Unter diesen unterschiedlichen Bedingungen schneidet BatMeth2 besser ab.
BatMeth2 hat eine bessere Mapping-Leistung auf realen BS-Seq-Daten
Um die Leistung von BatMeth2 auf realen BS-Seq-Datensätzen zu testen, luden wir Paired-End-BS-Seq-Datensätze herunter und extrahierten zufällig 1 Million 2 × 90 bp Paired-End-Reads aus SRA SRR847318, 1 Million 2 × 101 bp Paired-End-Reads aus SRA SRR1035722 und 1 Million 2 × 125 bp Paired-End-Reads aus SRA SRR3503136 zu Evaluierungszwecken. Da diese Datensätze von gesunden Zelllinien oder Geweben stammen, ist davon auszugehen, dass sie nur eine geringe Anzahl von Strukturvariationen enthalten. Daher haben wir diese realen Daten mit Single-End-Reads aus den Paired-End-Datensätzen ausgerichtet und die übereinstimmenden und nicht übereinstimmenden Zuordnungsraten aus den gepaarten Ausrichtungen ausgewertet, um die korrekten und falschen Ausrichtungsraten zu schätzen. Da die Insertgröße der Paired-End-Reads etwa 500 bp betrug, konnte ein Paar von Partner-Reads als übereinstimmend angesehen werden, wenn sie innerhalb eines nominalen Abstands von 500 bp gemappt wurden; andernfalls konnte ein Paar von Partner-Reads als nicht übereinstimmend angesehen werden. Ähnlich wie bei unseren Ergebnissen mit den simulierten Daten meldete BatMeth2 bei den realen Datensätzen über einen großen Bereich von Kartenqualitätswerten mehr konkordante und weniger diskordante Ausrichtungen, wie in Abb. 3 dargestellt.
Konkordanz- und Diskordanzraten von Ausrichtungen bei realen Paired-End-Reads von verschiedenen Alignern. Kumulative Anzahl der übereinstimmenden und nicht übereinstimmenden Alignments von hoher bis niedriger Mapping-Qualität für echte Bisulfit-Sequenzierungs-Reads. Für BSmap und die auf Bowtie basierenden Aligner gibt es jeweils nur einen Punkt, da diese Aligner keine Bewertung der Abbildungsqualität haben. Bismark-Bowtie2L15 bedeutet Bowtie2-Alignment mit Seed-Länge 15
In Tabelle 1 sind außerdem die relativen Laufzeiten der Programme angegeben. BatMeth2 lief mit den Standardeinstellungen schneller als die meisten der veröffentlichten Aligner und war mit BWA-meth und BatMeth vergleichbar. Bismark2 (mit Bowtie2 als grundlegende Mapping-Methode), BS Seeker2 und BSmooth benötigen längere Laufzeiten.
DNA-Methylierungs-Calling
Um die Genauigkeit des DNA-Methylierungs-Calling verschiedener Software zu bewerten, luden wir 450 K Bead-Chip-Daten der IMR90-Zelllinie von ENCODE (Encyclopedia of DNA Elements) herunter. Außerdem luden wir Ganzgenom-Bisulfit-Sequenzierungsdaten (WGBS-Seq) der IMR90-Zelllinie von ENCODE herunter (42,6 Gbases). Für jede Software richteten wir die WGBS-Seq-Daten aus und berechneten den Grad der DNA-Methylierung. Dann verglichen wir die Ergebnisse mit den MLs an denselben Stellen in den 450 K Bead Chip-Daten. Wenn die Differenz zwischen der DNA-ML aus den WGBS-Seq-Daten durch die Software und der aus den 450 K Bead Chip-Daten weniger als 0,2 betrug, wurde das Calling-Ergebnis als korrekt definiert; andernfalls wurde es als falsch angesehen.
Die Ergebnisse sind in Tabelle 2 dargestellt. Die Überlappung der korrekten Ergebnisse aller Softwareprogramme ist in Zusatzdatei 1 dargestellt: Abbildung S2. Es ist zu erkennen, dass BatMeth2 und Biscuit ähnliche Leistungen aufweisen, die besser sind als die der anderen Software. Zusammenfassend lässt sich sagen, dass BatMeth2 die Genauigkeit sowohl des BS-Read-Alignments als auch des DNA ML-Calling verbessert.
BatMeth2 aligniert BS-Reads unter Berücksichtigung von Indels variabler Länge
Krebs enthält einen deutlich höheren Anteil an Indels als gesunde Zellen. Um zu überprüfen, ob BatMeth2 BS-Reads mit Indels unterschiedlicher Länge alignieren kann, haben wir WGBS-Daten (75 Gbases) und 450 K Bead Chip-Daten von HepG2 (Leber hepatozelluläres Karzinom, eine Krebszelllinie) von ENCODE heruntergeladen. Nach dem Alignment der HepG2 WGBS-Seq-Daten überprüften wir die Verteilung der Indel-Längen in den Reads. Zusätzliche Datei 1: Abbildung S3A zeigt, dass die Längen der entdeckten Indels hauptsächlich im Bereich von 1 bp bis 5 bp verteilt waren, und das längste Indel war 40 bp lang. Nach unserer Statistik enthielten 2,3 % der Alignment-Reads Indels. Anhand dieser Ergebnisse wissen wir, dass BatMeth2 Reads mit Indels unterschiedlicher Länge alignieren kann.
Als Nächstes testeten wir die Auswirkungen der Indel-Erkennung auf das DNA-Methylierungs-Calling. Für BatMeth2 ließen wir zwei Optionen auf den HepG2-Daten laufen: mit und ohne Indel-Erkennung (d. h., wir setzten den Parameter -I in BatMeth2). Wir ließen auch Bismark auf den WGBS-Seq-Daten von HepG2 als Referenz für das DNA-Methylierungs-Calling mit Indel-Erkennung laufen, da Bismark keine Indel-Calling-Funktion hat. Wir verglichen das Callen von DNA-Methylierung in BatMeth2 und Bismark mit dem Callen aus den 450 K Bead Chip-Daten. Die Ergebnisse sind in Zusatzdatei 1 dargestellt: Abbildung S3B, wobei „BatMeth2-noIndel“ BatMeth2 ohne Indel-Erkennung entspricht. Es ist zu erkennen, dass das Ergebnis von BatMeth2 ohne Indel-Erkennung nur geringfügig besser ist als das von Bismark (mit Bowtie1 als grundlegende Kartierungsmethode). Das Ergebnis von BatMeth2 mit Indel-Erkennung war deutlich besser. Außerdem können wir sehen, dass BatMeth2 mehr DNA-Methylierungsstellen erkennen kann als BatMeth2-noIndel und Bismark (Bowtie 1). Um zu verstehen, warum die Leistung von BatMeth2 bei der Indel-Erkennung besser ist, haben wir die von BatMeth2 gefundenen Methylierungsstellen als Ergebnis A definiert, während die von BatMeth2-noIndel und Bismark gefundenen Methylierungsstellen als Ergebnis B definiert wurden. Dann ließen wir mclA die Methylierungsstellen sein, die in Ergebnis A, aber nicht in Ergebnis B auftauchen. Wir beobachteten, dass mclA 23.853 DNA-Methylierungsstellen und 15.048 (63 %) der 23.853 Stellen umfasste, die von den Alignments der Indel-Reads abgedeckt wurden, die von BatMeth2 mit Indel-Erkennung aufgerufen wurden (siehe Zusatzdatei 1: Abbildung S3C). Darüber hinaus stellten wir fest, dass die Indel-Raten in Ergebnis A und Ergebnis B nur 5 bzw. 0 % betrugen. Daraus schlussfolgerten wir, dass eine genaue Indel-Erkennung die DNA-Methylierungsbestimmung verbessern kann.
Visualisierung von DNA-Methylierungsdaten
BatMeth2 bietet Werkzeuge zur Visualisierung der Methylierungsdaten. Um die Visualisierungsfunktionen von BatMeth2 zu veranschaulichen, haben wir (1) 117 Gbases von Single-End-Reads aus der menschlichen H9-Zelllinie, (2) 105,2 Gbases von Single-End-Reads aus der menschlichen IMR90-Zelllinie und (3) 12,6 Gbases von Paired-End-Reads aus Wildtyp-Reis heruntergeladen. Zunächst kann BatMeth2 die Cytosin-Methylierungsdichte auf Chromosomenebene sichtbar machen. Die Punkte in Abb. 4a stellen ein gleitendes Fenster von 100 kb mit einem Schritt von 50 kb dar. Um die ML an einzelnen CpG- oder Nicht-CpG-Stellen in einem Genombrowser betrachten zu können, stellen wir auch Dateien im bed- und bigWig-Format zur Verfügung (Abb. 4b). Beim Vergleich mit der Dichte von Genen und TEs konnten wir feststellen, dass die ML mit der TE-Dichte korreliert und mit der Gendichte antikorreliert ist (Abb. 4c). Diese Tendenz wurde zuvor bei Reis beobachtet.
Visualisierung der Methylierungsniveaus in der Chromosomenskala. a Die Methyl-Cytosin-Dichte im menschlichen Chromosom 10. Die Punkte stellen die Methylierungsgrade in gleitenden Fenstern von 100Kb mit einem Schritt von 50Kb dar. Die roten Punkte beziehen sich auf die Methylierungsniveaus im Plus-Strang und die blauen Punkte auf die Methylierungsniveaus im Minus-Strang. b Ein Beispiel für die Verteilung der DNA-Methylierungsniveaus und der differentiell methylierten Regionen (DMRs) zwischen H9- und IMR90-Zelllinien im menschlichen Chromosom 10. c Die Dichte von Genen, Transposon-Elementen (TEs) und das Niveau der DNA-Methylierung im gesamten Reisgenom. Tafel A zeigt die mit Batmeth2 erzeugten Ergebnisse. Tafel B zeigt die Visualisierungsergebnisse des UCSC-Browsers mit den BED-Dateien von Batmeth2
Zweitens kann BatMeth2 die MLs von Genen visualisieren. Genauer gesagt kann BatMeth2 die MLs 2 kb stromaufwärts des Gens, an der Transkriptionsstartstelle (TSS), im Genkörper, an der Transkriptionsendstelle (TES) und 2 kb stromabwärts des Genkörpers sichtbar machen. Vergleicht man die stromaufwärts, im Genkörper und stromabwärts gelegenen Regionen, so zeigt Abb. 5a, dass die ML der DNA im Genkörper höher ist als in der Promotorregion. Vergleicht man alle fünf Regionen, so zeigt sich ein Tal in der TSS-Region (Abb. 5b). BatMeth2 kann auch die ML-Profile um Introns, Exons, intergene Regionen und TEs berechnen (Zusatzdatei 1: Abbildung S4). Darüber hinaus kann BatMeth2 eine Heatmap mehrerer Gene nach Genregionen erstellen, um die Gesamt-MLs verschiedener Proben bequem vergleichen zu können (Abb. 5c).
Visualisierung der DNA-Methylierung in verschiedenen Kontexten. a Die DNA-Methylierungsniveaus in 2Kb-Regionen stromaufwärts von Genen, Genkörpern, 2Kb stromabwärts von Genkörpern. b Das Aggregationsprofil der DNA-Methylierungsniveaus über Gene hinweg. c Die Wärmekarte aller Gene in 2Kb-Regionen stromaufwärts von Genen, Genkörpern, 2Kb stromabwärts von Genkörpern
Drittens kann BatMeth2 die Verteilung der DNA-Methylierung visualisieren. Additional file 1: Abbildung S5A zeigt die DNA-Methylierungsverteilungen in den Zelllinien H9 und IMR90. In der Abbildung ist die DNA-ML in fünf Kategorien unterteilt: methyliert (M: > 80%), zwischen teilweise methyliert und methyliert (Mh: 60-80%), teilweise methyliert (H: 40-60%), zwischen nicht methyliert und teilweise methyliert (hU: 20-40%) und nicht methyliert (U: < 20%). Wie in Zusatzdatei 1: Abbildung S5A, war die ML in der H9-Zelllinie in der M-Kategorie höher als in der IMR90-Zelllinie, insbesondere im CpG-Kontext. Im CH-Sequenzkontext ist die CpG-Methylierung die vorherrschende Form, aber ein signifikanter Anteil der methylierten Cytosine findet sich an CpA-Stellen, während die ML weniger als 40 % beträgt, insbesondere in der H9-Zelllinie (Zusätzliche Datei 1: Abbildung S5B).
Viertens kann BatMeth2 die Korrelation zwischen Genexpressionsniveau und Genpromoter-DNA-ML analysieren. Wir haben diese Funktion anhand der Zelllinien H9 und IMR90 veranschaulicht. Die Expressionsniveaus der Gene in H9 oder IMR90 wurden in verschiedene Kategorien unterteilt. Wie in Zusatzdatei 1: Abbildung S5C zu sehen ist, wiesen die hoch exprimierten Gene niedrigere MLs in ihren Promotorregionen auf. Außerdem unterteilten wir die MLs der Genpromotoren in fünf Kategorien. Das Ergebnis in Additional file 1: Abbildung S5D zeigt, dass Gene mit Promotoren mit höheren ML-Werten geringere Expressionswerte aufwiesen. Die negative Korrelation zwischen der Expression von Säugetiergenen und der DNA-Methylierung des Promotors ist bekannt. Diese Analyse ist ein weiterer Hinweis auf die Genauigkeit von BatMeth2.
Finden von differenziell methylierten Cytosinen und Regionen (DMCs/DMRs)
Die Identifizierung von differenziell methylierten Cytosinen (DMCs) und differenziell methylierten Regionen (DMRs) ist eines der wichtigsten Ziele bei der Analyse von Methylierungsdaten. Obwohl Forscher gelegentlich daran interessiert sind, einzelne Cytosinstellen mit einem Phänotyp zu korrelieren, sind DMRs sehr wichtige Merkmale.
In frühen BS-Seq-Studien wurden Zellen profiliert, ohne Wiederholungen zu sammeln. Für solche Datensätze haben wir den exakten Test von Fisher verwendet, um differenziell methylierte Cytosine (DMCs) zu erkennen. Für BS-Seq-Datensätze mit Wiederholungen ist das natürlichste statistische Modell zur Bestimmung von DMCs die Beta-Binomialverteilung. Es ist bekannt, dass eine Reihe von Softwareprogrammen eine differenzielle DNA-Methylierungsdatenanalyse durchführen kann, wie z. B. methylKit (ein differenzielles Analyseprogramm, das biologische Replikate erfordert) und Methy-Pipe (ein differenzielles Analyseprogramm ohne biologische Duplikate). Es gibt jedoch kein umfassendes Paket, das sowohl die Kartierung als auch die differentielle Methylierungsanalyse beinhaltet. Daher haben wir ein Paket entwickelt, das die Kartierung mit der differentiellen Analyse verbindet. Um die Identifizierung von DMRs aus Bisulfit-Daten ohne Replikate zu erleichtern, haben wir den exakten Test von Fisher integriert, um einen Hypothesentest durchzuführen. Wenn eine Probe zwei oder mehr Replikate hat, verwenden wir die Beta-Binomialverteilung, um eine differenzielle Methylierungsanalyse durchzuführen. Wir stellen auch bed- oder bigWig-Dateien für die Liste der DMRs zur Verfügung. Die DMRs können in einem Genombrowser (Abb. 4b) mit den erzeugten bed- oder bigWig-Dateien visualisiert werden.
Zur Veranschaulichung zeigt Abb. 6a die Anzahl der DMCs und Regionen in der IMR90-Zelllinie und in der H9-Zelllinie, wie sie von BatMeth2 erkannt wurden (p-Wert< 0,05, meth.diff > = 0,6). BatMeth2 kann anzeigen, ob CpGs und DMCs in bestimmten Regionen angereichert sind, wie z. B. in Gen-, CDS-, Intron-, intergenen, UTR-, TE-, LTR-, LINE- und SINE-Regionen. Abbildung 6b veranschaulicht die Anteile von DMCs in verschiedenen genomischen Regionen. Abgesehen von den intergenischen Regionen konnten wir in keiner Region eine DMC-Anreicherung beobachten.
Differenzielle Methylierungsanalyse. a Analyseergebnisse der differentiell methylierten Regionen (DMRs), differentiell methylierten Gene (DMGs) und differentiell methylierten Promotoren (DMPs) zwischen den Zelllinien H9 und IMR90. b Annotation von differentiell methylierten Cytosinen (DMC) anhand verschiedener genomischer Eigenschaften und Wiederholungselemente. c DMPs enthalten H9- oder IMR90-spezifische Indels (orange), die einen erheblichen Anteil an allen DMPs ausmachen (DNA Methylation differential Promoters)
Ein erheblicher Anteil der differentiell methylierten Promotoren (DMPs) enthält Indels
Wir wissen, dass Indels und DNA-Methylierung eine wichtige Rolle bei der Gewebeentwicklung und bei Krankheiten spielen. Hier untersuchen wir die Beziehung zwischen differenziell methylierten Promotoren (DMP) und Indels. Für diese Studie verwendeten wir die BS-Seq-Reads der Zelllinien IMR90 und H9. Zunächst richteten wir die BS-Seq-Reads mit BatMeth2 aus; dann wurden die Indels mit BisSNP- und GATK-Tools bestimmt. Anschließend definierten wir die Indels, die nur in H9 oder IMR90 vorkommen, als zelllinienspezifische Indels.
Danach entdeckten wir mit BatMeth2 1384 DMPs zwischen H9 und IMR90 (p-Wert< 0,05, meth.diff > = 0,6). Insgesamt 236 (17 %) aller oben genannten DMPs enthalten Indels, wie in Abb. 6c dargestellt. Kurz gesagt, ein erheblicher Anteil der DMPs enthält Indels. Daher ist ein genaues Alignment von BS-Seq-Reads in der Nähe dieser Indels für die Erforschung der DNA-Methylierung sehr wichtig.