Genaue Vorhersage der zelltypspezifischen Bindung von Transkriptionsfaktoren

Im Rahmen des ENCODE-DREAM-Wettbewerbs wurde eine große Anzahl von Ansätzen, die von 40 internationalen Teams entwickelt wurden, anhand von 13 zelltypspezifischen ChIP-seq-Assays für 12 verschiedene TFs beim Menschen verglichen (Additional file 1: Abbildung S1). Ein Satz von 109 Datensätzen für die gleichen (und zusätzliche) TFs in anderen Zelltypen wurde für das Training bereitgestellt. Die Trainingsdaten umfassten zelltypspezifische DNase-seq-Daten, zelltypspezifische RNA-seq-Daten, genomische Sequenzen und Annotationen sowie In-silico-DNA-Formvorhersagen. Darüber hinaus wurden zelltypspezifische und TF-spezifische ChIP-seq-Daten und abgeleitete Beschriftungen für Trainingschromosomen bereitgestellt, während die Vorhersagen nur für die verbleibenden Chromosomen chr1, chr8 und chr21 ausgewertet wurden, die keine ChIP-seq-Trainingsdaten enthielten. Für 200-bp-Regionen, die um 50 bp verschoben sind, wurden genomweite Vorhersagen über die Wahrscheinlichkeit, dass eine bestimmte Region einen ChIP-seq-Peak überlappt, von den beteiligten Teams angefordert. Die Vorhersagen wurden anhand (i) der Fläche unter der ROC-Kurve (AUC-ROC), (ii) der Fläche unter der Präzisions-Rückruf-Kurve (AUC-PR), (iii) des Rückrufs bei 10 % FDR und (iv) des Rückrufs bei 50 % FDR für jeden der 13 Testdatensätze bewertet. Diese wurden für jeden Datensatz auf der Grundlage des durchschnittlichen, normalisierten Rangs, der für jedes dieser Maße in 10 Bootstrap-Stichproben der herausgehaltenen Chromosomen erzielt wurde, aggregiert, und ein endgültiges Ranking wurde als Durchschnitt dieser Rangstatistiken erhalten (vgl. https://www.synapse.org/#!Synapse:syn6131484/wiki/405275).

Als Ergebnis dieses Rankings belegte der in dieser Arbeit vorgestellte Ansatz (Team „J-Team“) gemeinsam mit dem Ansatz des Teams „Yuanfang Guan“ den ersten Platz.

Im Folgenden untersuchen wir den Einfluss verschiedener Aspekte des vorgeschlagenen Ansatzes auf die endgültige Vorhersageleistung. Erstens untersuchen wir die Auswirkungen verschiedener Gruppen von verwandten Merkmalen (DNase-seq-Daten, Motiv-Scores, RNA-seq-Daten, sequenzbasierte und annotationsbasierte Merkmale) auf die Vorhersageleistung. Zweitens untersuchen wir die Bedeutung des iterativen Trainingsansatzes im Gegensatz zu einem Training auf Basis der ursprünglichen Trainingsdaten. Drittens vergleichen wir die Leistung der Vorhersagen von Klassifikatoren, die auf Trainingsdaten für einzelne Zelltypen trainiert wurden, mit der Leistung der aggregierten Vorhersage, die durch Mittelwertbildung über diese Zelltypen erzielt wurde. Schließlich wenden wir die vorgeschlagene Methode zur Vorhersage der zelltypspezifischen TF-Bindung für 31 TFs in 22 weiteren primären Zelltypen an, was insgesamt 682 Vorhersagespuren ergibt.

Auswirkung der Merkmalssätze auf die Vorhersageleistung

Wir verwenden die Vorhersageleistung, die durch den vorgeschlagenen Ansatz unter Verwendung aller Merkmalssätze (Abschnitt „Merkmale“), des iterativen Trainingsverfahrens (Abschnitt „Iteratives Training“) und der Aggregation über alle Trainingszelltypen (Abschnitt „Vorhersageschema“) erzielt wurde, als Basis für alle weiteren Vergleiche (Abb. 1; „alle Merkmale“). In diesem Manuskript betrachten wir AUC-PR als primäres Leistungsmaß, da AUC-PR aussagekräftiger über die Klassifizierungsleistung bei stark unausgewogenen Klassifizierungsproblemen ist und Recall auf den verschiedenen FDR-Stufen eher instabil ist, da es einzelnen Punkten auf der Präzisions-Recall-Kurve entspricht. Die AUC-PR-Werte werden mit dem R-Paket PRROC berechnet, das auch im ENCODE-DREAM-Wettbewerb verwendet wurde.

Wir stellen fest, dass die durch AUC-PR gemessene Vorhersageleistung zwischen den verschiedenen Transkriptionsfaktoren stark variiert (Abb. 1), mit einem mittleren AUC-PR-Wert von 0,4098. Die beste Vorhersageleistung wird für CTCF, der ein langes und informationsreiches Bindungsmotiv hat, in zwei verschiedenen Zelltypen (IPSC und PC-3) erreicht. Eine überdurchschnittliche Leistung wird auch für FOXA1 und HNF4A in Leberzellen erzielt. Für die meisten anderen TFs finden wir AUC-PR-Werte um 0,4, während wir für NANOG und REST eine eher niedrige Vorhersagegenauigkeit beobachten.

Um den Beitrag ausgewählter Merkmale zur endgültigen Vorhersageleistung zu analysieren, schließen wir systematisch Gruppen verwandter Merkmale aus den Eingabedaten für Training und Vorhersage aus. Als Basiswert messen wir den AUC-PR für den Klassifikator unter Verwendung aller Merkmalsätze. Zusätzlich messen wir den AUC-PR, wenn wir jeden einzelnen Merkmalsatz ausschließen, wobei die Differenz dieser beiden AUC-PR-Werte die Verbesserung quantifiziert, die durch die Einbeziehung des Merkmalsatzes erzielt wird (Abb. 2a).

Wir beobachten die größte Auswirkung für den Satz von Merkmalen, die aus DNase-seq-Daten abgeleitet wurden. Die Verbesserung der AUC-PR durch die Einbeziehung von DNase-seq-Daten schwankt zwischen 0,087 für E2F1 und 0,440 für HNF4A mit einem Median von 0,252.

Merkmale, die auf Motiv-Scores basieren (einschließlich de novo entdeckter Motive und solcher aus Datenbanken), tragen ebenfalls erheblich zur endgültigen Vorhersageleistung bei. Hier beobachten wir große Verbesserungen für einige TFs, nämlich 0,231 für CTCF in IPSC-Zellen, 0,175 für CTCF in PC-3-Zellen und 0,167 für FOXA1. Im Gegensatz dazu beobachten wir einen Rückgang der Vorhersageleistung im Fall von JUND (- 0,080), wenn motivbasierte Merkmale einbezogen werden. Für die übrigen TFs finden wir Verbesserungen der AUC-PR zwischen 0,008 und 0,079. Des Weiteren betrachten wir zwei Untergruppen von Motiven, nämlich alle Motive, die durch De-novo-Motiventdeckung aus den Challenge-Daten gewonnen wurden, und alle Slim/LSlim-Modelle, die die Abhängigkeiten zwischen den Motiven erfassen. Für Motive aus der de-novo-Motiventdeckung finden wir eine Verbesserung für 9 der 13 Datensätze, und für das Slim/LSlim-Modell finden wir eine Verbesserung für 10 der 13 Datensätze. Die absoluten Verbesserungen (Median von 0,011 bzw. 0,006) sind jedoch eher gering, möglicherweise weil (i) Motive, die durch De-novo-Motiventdeckung erhalten wurden, redundant zu denen in Datenbanken sein könnten und (ii) Intra-Motiv-Abhängigkeiten und Heterogenitäten, die durch Slim/LSlim-Modelle erfasst werden, teilweise durch Variationen in den Motiven aus verschiedenen Quellen abgedeckt werden könnten.

Besonders hervorzuheben sind RNA-seq-basierte Merkmale (Median 0.001), annotationsbasierte Merkmale (0,000) und sequenzbasierte Merkmale (0,001) haben fast keinen Einfluss auf die Vorhersageleistung.

Da der Satz von DNase-seq-basierten Merkmalen ziemlich vielfältig ist, einschließlich Merkmalen, die von Fold-Enrichment-Spuren, Peak-Listen oder Variationen zwischen Zelltypen abgeleitet sind, zielen wir darauf ab, den Einfluss verwandter Gruppen dieser Merkmale weiter zu untersuchen. Zu diesem Zweck testen wir weiter, wie die Vorhersageleistung durch das Entfernen spezifischer Gruppen von DNase-seq-Merkmalen (vgl. Zusatzdatei 1: Text S2) aus dem vollständigen Merkmalssatz (Zusatzdatei 1: Abbildung S2) beeinflusst wird. Wir stellen fest, dass keine dieser Merkmalsgruppen allein einen großen Einfluss auf die Vorhersageleistung hat, obwohl allmähliche Unterschiede zu beobachten sind, da die Einbeziehung von Fold-Enrichment-basierten und Peak-basierten Merkmalen einen weitgehend positiven Beitrag leistet, während der Einfluss der anderen Merkmalsgruppen eher unklar ist. Dies könnte auf die umfangreichen Redundanzen und Korrelationen zurückzuführen sein, die zwischen diesen verschiedenen Gruppen bestehen, wodurch der Verlust einer einzelnen Merkmalsgruppe weitgehend kompensiert werden kann.

Daher testen wir zusätzlich ein Szenario, bei dem das Weglassen aller DNase-seq-basierten Merkmale (d. h. die Daten hinter dem Geigenplot „DNase-seq“ in Abb. 2a) als Basisfall betrachtet wird und nur eine der spezifischen Gruppen zu diesem reduzierten Merkmalssatz hinzugefügt wird (Abb. 2b). Zunächst einmal ist festzustellen, dass alle Merkmalsgruppen einen positiven Beitrag zur Gesamtvorhersageleistung leisten. Den größten Beitrag leistet die Gruppe „Faltenanreicherung“, aber auch verwandte Gruppen wie „Long Range“, die im Wesentlichen einen Mittelwert über breitere Fenster der Faltenanreicherungsspur bildet, und „Peak-based“, die Peaks verwendet, die ursprünglich auf der Grundlage der DNase-Seq-Abdeckung aufgerufen wurden. Den geringsten Beitrag finden wir in der Gruppe „Variation“, die die Variation bzw. Erhaltung des DNase-seq-Signals zwischen den Zelltypen misst. Da der Beitrag jeder einzelnen Gruppe von Merkmalen positiv ist, betrachten wir im Folgenden den kompletten Satz DNase-seq-basierter Merkmale.

Nachdem wir festgestellt haben, dass DNase-seq-basierte und motivbasierte Merkmale einen großen Einfluss auf die Vorhersageleistung haben, haben wir auch die Vorhersageleistung des vorgeschlagenen Ansatzes getestet, indem wir nur Merkmale verwendet haben, die auf DNase-seq-Daten bzw. TF-Motiven basieren. Alle anderen Merkmale, d. h. RNA-seq-basierte Merkmale, annotationsbasierte Merkmale und auf Rohsequenz basierende Merkmale, wurden ausgeschlossen. Wir stellen fest (Abb. 1), dass Klassifikatoren, die ausschließlich motivbasierte Merkmale verwenden, für einige TFs (CTCF und in gewissem Maße E2F1 und GABPA) bereits eine angemessene Vorhersageleistung erbringen, während wir für die übrigen TFs AUC-PR-Werte unter 0,12 beobachten. Dies kann durch die große Anzahl von falsch-positiven Vorhersagen erklärt werden, die typischerweise von Ansätzen erzeugt werden, die ausschließlich Motivinformationen verwenden, was nur im Fall von langen, spezifischen Motiven vermieden werden kann, wie es bei CTCF der Fall ist.

Klassifikatoren, die nur DNase-seq-basierte Merkmale verwenden, liefern für viele der untersuchten TFs eine bemerkenswerte Leistung (Abb. 1), die nur für die beiden CTCF-Datensätze niedriger ist als für den motivbasierten Klassifikator. Für einige Datensätze (insbesondere JUND, aber auch EGR1, MAX) stellen wir sogar fest, dass ein Klassifikator, der nur auf DNase-seq-Daten basiert, den Klassifikator, der alle Merkmale verwendet, übertrifft.

Im Fall von JUND kann die Leistungssteigerung bei Vernachlässigung aller Nicht-DNase-Merkmale wahrscheinlich auf eine starke Anpassung der Klassifikatorparameter an entweder zelltypspezifische Bindungsmotive oder zelltypspezifische Co-Bindung mit anderen TFs zurückgeführt werden, da JUND der einzige Datensatz mit einer verbesserten Leistung ist, wenn man, wie oben beschrieben, motivbasierte Merkmale ausschließt. Für alle drei TFs finden wir eine Verbesserung der Vorhersageleistung, wenn die Klassifikatorparameter auf den Trainingschromosomen des Testzelltyps trainiert werden (Fall „innerhalb des Zelltyps“; Zusatzdatei 1: Abbildung S3).

Da DNase-seq-basierte und motivbasierte Merkmale die primären Merkmalsätze zu sein scheinen, die die Vorhersageleistung beeinflussen, untersuchen wir schließlich die Vorhersageleistung eines Klassifikators, der nur diese beiden Merkmalsätze verwendet. Wir stellen fest, dass die Vorhersageleistung, die nur DNase-seq-basierte und motivbasierte Merkmale verwendet, weitgehend identisch mit der des Klassifikators ist, der alle Merkmale verwendet (Abb. 1), wobei wir den größten Verlust an AUC-PR für TAF1 (0,017) und den größten Gewinn an AUC-PR für NANOG (0,007) beobachten. Wir stellen ein ähnliches Verhalten für den Fall innerhalb des Zelltyps fest (Zusatzdatei 1: Abbildung S3). Da die ausgelassenen Merkmalsätze alle RNA-seq-basierten Merkmale enthalten, hat dies auch zur Folge, dass ein zelltypspezifischer Assay (nämlich DNase-seq) für die Vorhersage der TF-Bindung ausreicht, was den Bereich der Zelltypen mit leicht verfügbaren experimentellen Daten, auf die der vorgeschlagene Ansatz angewendet werden kann, erweitert.

Iteratives Training verbessert die Vorhersageleistung

Als zweiten Schlüsselaspekt des vorgeschlagenen Ansatzes untersuchen wir die Auswirkungen des iterativen Trainingsverfahrens auf die endgültige Vorhersageleistung. Zu diesem Zweck vergleichen wir für jeden TF die AUC-PR-Werte, die sich aus der Mittelung der Vorhersagen aller fünf Klassifikatoren ergeben, die sich aus dem iterativen Trainingsverfahren für alle Trainingszelltypen ergeben, mit den AUC-PR-Werten, die sich aus der Mittelung der anfänglich trainierten Klassifikatoren für alle Trainingszelltypen ergeben, d. h. Klassifikatoren, die nur auf den anfänglichen Trainingsdaten trainiert wurden (Abschnitt „Anfängliche Trainingsdaten“).

Für 11 der 13 Testdatensätze beobachten wir eine Verbesserung der Vorhersageleistung durch das iterative Trainingsverfahren (Abb. 3). Die größten Verbesserungen werden für E2F1 (0,114), FOXA2 (0,085), NANOG (0,08), FOXA1 (0,063) und MAX (0,061) erzielt. Darunter befinden sich TFs, für die wir eine gute Leistung bei der ausschließlichen Verwendung von DNase-seq-basierten Merkmalen (E2F1, MAX) beobachtet haben, und TFs, für die die Kombination mit motivbasierten Merkmalen vorteilhaft war (FOXA1, FOXA2, NANOG), was darauf hindeutet, dass die zusätzlichen negativen Regionen, die in den Iterationen 2 bis 5 hinzugefügt wurden, keine Verzerrung in Richtung eines dieser beiden Merkmalstypen bewirken. Für vier dieser fünf TFs wurden nur ein (FOXA2, NANOG, FOXA1) oder zwei (E2F1) Trainingszelltypen bereitgestellt, und die Variation zwischen den verschiedenen Klassifikatoren aus dem iterativen Training kann dazu beitragen, eine Überanpassung zu vermeiden. Im Gegensatz dazu stellen wir einen Leistungsabfall für JUND (0,041) und auch TAF1 (0,01) fest, der durch eine stärkere Betonung der zelltypspezifischen Bindungsregionen in den nachfolgenden Iterationen des iterativen Trainingsverfahrens verursacht werden könnte. Diese Hypothese wird auch durch die Beobachtung gestützt, dass das iterative Trainingsverfahren immer dann zu einer Steigerung der Vorhersageleistung führt, wenn Klassifikatorparameter auf den Trainingschromosomen des Testzelltyps trainiert werden (Additional file 1: Abbildung S4).

Die Mittelung der Vorhersagen verbessert die zufällige Auswahl der Zelltypen

Für 9 der 12 betrachteten TFs, werden Daten für mehr als einen Trainingszelltyp mit den Challenge-Daten bereitgestellt. Eine zentrale Frage könnte daher die Wahl des Zelltyps sein, der für das Training und anschließend für die Vorhersagen für den Testzelltyp verwendet wird. Die einzigen zelltypspezifischen experimentellen Daten, die für diese Wahl zur Verfügung stehen, sind jedoch DNase-seq- und RNA-seq-Daten, während die Ähnlichkeit der Zelltypen von der betrachteten TF abhängen könnte. Tatsächlich haben Ähnlichkeitsmaße, die aus DNase-seq-Daten (z. B. Jaccard-Koeffizienten überlappender DNase-seq-Peaks, Korrelation von Profilen) oder aus RNA-seq-Daten (z. B., Korrelation von TPM-Werten) haben sich in Vorstudien zu den Trainingszelltypen als nicht aussagekräftig in Bezug auf die Ähnlichkeit der TF-Bindungsregionen erwiesen.

Daher betrachten wir die Wahl des Trainingszelltyps als latente Variable und mitteln über die von den jeweiligen Klassifikatoren generierten Vorhersagen (siehe Abschnitt „Vorhersageschema“). Da die Etiketten der Testzelltypen nach der Challenge zur Verfügung gestellt wurden, können wir nun die Auswirkungen dieser Wahl auf die Vorhersageleistung bewerten und auch die Vorhersageleistung von Klassifikatoren testen, die auf einzelnen Zelltypen trainiert wurden (Abb. 4). 4).

Für alle Testdatensätze mit mehreren verfügbaren Trainingszelltypen stellen wir fest, dass die gemittelte Vorhersage AUC-PR-Werte ergibt, die über dem Median der AUC-PR-Werte liegen, die für einzelne Trainingszelltypen erzielt wurden. Diese Verbesserung ist besonders ausgeprägt für REST, GABPA und MAX.

Um weiter zu untersuchen, ob die Mittelwertbildung über Klassifikatoren für einzelne Zelltypen konservierte Bindungsregionen (d.h. Regionen, die in der Mehrzahl der Zelltypen als „gebunden“ gekennzeichnet sind) gegenüber zelltypspezifischen Bindungsregionen begünstigt, bewerten wir die Vorhersageleistung für solche Regionen auch separat (Zusatzdatei 1: Abbildung S5). Insbesondere betrachten wir eine gebundene Region als konserviert, wenn sie auch in mindestens drei von vier Trainingszelltypen als „gebunden“ gekennzeichnet ist, und wir betrachten eine gebundene Region als zelltypspezifisch, wenn diese Region in höchstens einem von vier Trainingszelltypen als „gebunden“ gekennzeichnet ist. Das erste, was uns in Zusatzdatei 1: Abbildung S5 ist, dass die absoluten AUC-PR-Werte für zelltypspezifische Regionen wesentlich niedriger sind als für konservierte Regionen. Eine Erklärung könnte ein Unterschied in der Klassen(un)gleichheit aufgrund der ausgewählten Untergruppen von Regionen sein. Dieser allgemeine Trend bleibt jedoch bestehen, wenn man die AUC-ROC betrachtet (Zusätzliche Datei 1: Abbildung S6). Zweitens stellen wir fest, dass die Abweichung zwischen Klassifikatoren, die aus verschiedenen Trainingszelltypen gelernt wurden, in den meisten Fällen für die zelltypspezifischen Regionen größer ist als für die konservierten Regionen. Das Verhalten in Bezug auf die absolute Leistung ist ähnlich für die einzelnen Klassifikatoren, ihre mittlere Leistung und die Leistung der Mittelung über Klassifikatoren für einzelne Zelltypen. Wir stellen fest, dass der AUC-PR, der durch die Mittelwertbildung erzielt wird, immer besser ist als die mediane Leistung für einzelne Zelltypen für konservierte Regionen, aber dasselbe gilt auch, wenn zelltypspezifische Regionen für sieben der neun Datensätze mit mehr als einem Trainingszelltyp berücksichtigt werden.

Daher können wir argumentieren, dass die Mittelwertbildung über die zelltypspezifischen Klassifikatoren im Allgemeinen genauere Vorhersagen liefert, als durch eine uninformierte Wahl eines spezifischen Trainingszelltyps erreicht werden würde.

Wir stellen jedoch auch fest, dass für fast alle Testdatensätze mit mehreren Trainingszelltypen (die einzige Ausnahme ist CTCF für den Zelltyp PC-3) die beste Vorhersageleistung, die für einen der einzelnen Trainingszelltypen erzielt wurde, in einigen Fällen erhebliche Verbesserungen gegenüber dem vorgeschlagenen Mittelwertverfahren gebracht hätte. Bemerkenswert ist, dass die Varianz der AUC-PR zwischen den verschiedenen Trainingszelltypen für JUND besonders ausgeprägt ist, was die frühere Hypothese stützt, dass einige Merkmale, z. B. Bindungsmotive oder Co-Bindung von TFs, für JUND sehr zelltypspezifisch sind. Im Allgemeinen würde die Ableitung informativer Maße für die TF-spezifische Zelltyp-Ähnlichkeit auf der Grundlage von zelltypspezifischen Assays und vorläufigen Vorhersagen von Bindungsstellen wahrscheinlich zu einer weiteren Steigerung der Leistung von computergestützten Ansätzen für die Vorhersage zelltypspezifischer TF-Bindung führen.

Erstellung einer Sammlung von zelltypspezifischen TF-Bindungsspuren

Nachdem wir festgestellt haben, dass ein einziger experimenteller Assay, nämlich DNase-seq, für die Vorhersage zelltypspezifischer TF-Bindung mit modernster Genauigkeit ausreicht, können wir nun die Klassifikatoren, die wir für die Trainingszelltypen und TFs erhalten haben, für Vorhersagen für weitere Zelltypen verwenden. Zu diesem Zweck verwenden wir die Klassifikatoren, die nur DNase-seq-basierte und motivbasierte Merkmale berücksichtigen, aber weder RNA-seq-basierte Merkmale, annotationsbasierte Merkmale noch Merkmale, die auf der Rohsequenz basieren, was sich als vergleichbar mit der Vorhersageleistung des vollständigen Modells erwiesen hat (vgl. Abb. 1, Abschnitt „Auswirkung von Merkmalssätzen auf die Vorhersageleistung“). Zu diesem Zweck luden wir DNase-seq-Daten für eine Sammlung von primären Zelltypen und Geweben herunter (siehe Abschnitt „Daten“), verarbeiteten diese auf die gleiche Weise wie die ursprünglichen Challenge-Daten und extrahierten anschließend DNase-seq-abhängige Merkmale (Abschnitt „Merkmale“). Wir haben dann die trainierten Klassifikatoren für alle 31 TFs, die in der Challenge berücksichtigt wurden, auf diese 22 DNase-seq-Merkmalssätze angewandt, um insgesamt 682 Vorhersagespuren zu erhalten.

Für die ausgewählten Zelltypen (Zusatzdatei 1: Tabelle S5) sind nur wenige zelltyp- und TF-spezifische ChIP-seq-Daten verfügbar (Zusatzdatei 1: Tabelle S6). Einerseits bedeutet dies, dass die vorhergesagten TF-Bindungsspuren wertvolle, neue Informationen für die untersuchte Sammlung von 31 TFs liefern. Andererseits bietet dies die Möglichkeit, ein Benchmarking und eine Überprüfung der Richtigkeit der Vorhersagen für die Untergruppe dieser TFs und Zelltypen durchzuführen, für die entsprechende ChIP-seq-Daten vorliegen. Für das Benchmarking erhalten wir zusätzlich die „entspannten“ und (wo verfügbar) „konservativen“ Peak-Dateien von ENCODE und leiten die zugehörigen Labels („gebunden“, „ungebunden“, „mehrdeutig“) gemäß dem für die ENCODE-DREAM-Herausforderung vorgeschlagenen Verfahren ab.

Für CTCF mit verfügbaren ChIP-seq-Peaks für mehrere Zelltypen finden wir im Allgemeinen eine Vorhersageleistung, die mit der Leistung vergleichbar ist, die für die Herausforderungsdaten beobachtet wurde (vgl. Additional file 1: Tabelle S4). Für diese Zelltypen liegen die AUC-PR-Werte (Zusätzliche Datei 1: Tabelle S7) zwischen 0,7720 und 0,8197, wenn konservative und entspannte Peaks verfügbar sind und wenn die Spender zwischen den DNase-seq- und ChIP-seq-Experimenten übereinstimmen, während die Leistung bei nicht übereinstimmenden Spendern (0,7322) und im Falle fehlender konservativer Peaks (0,7270) etwas geringer ist. Für JUN, MAX und MYC sind aufgrund fehlender Replikate nur entspannte Peaks aus ENCODE verfügbar. Hier finden wir AUC-PR-Werte von 0,6310 für JUN, was wesentlich höher ist als bei den Challenge-Daten; 0,4004 für MAX, was etwas niedriger ist als bei den Challenge-Daten; und 0,1989 für MYC, das nicht zu den Test-TFs in der Challenge gehörte, aber in der Leaderboard-Runde eine wesentlich bessere Leistung erzielte.

Die 682 genomweiten Vorhersagespuren sind immer noch recht groß (ca. 880 MB pro Spur) und beanspruchen daher beträchtlichen Speicherplatz, der dem typischen Nutzer möglicherweise nicht zur Verfügung steht, während die Mehrzahl der Regionen wahrscheinlich nicht durch die interessierende TF gebunden ist. Daher verdichten wir diese Vorhersagen zu Listen mit vorhergesagten Peaks im narrowPeak-Format, indem wir zusammenhängende Abschnitte mit hoher Bindungswahrscheinlichkeit verbinden und einen Schwellenwert von 0,6 (entspannt) bzw. 0,8 (konservativ) auf die maximale Wahrscheinlichkeit anwenden, die in einem vorhergesagten „Peak“ beobachtet wird. Wir stellen diese Peak-Dateien unter https://www.synapse.org/#!Synapse:syn11526239(doi:10.7303/syn11526239) zum Download bereit.

Um einen Eindruck von der Qualität der vorhergesagten Peaks zu erhalten, berechnen wir außerdem Jaccard-Koeffizienten auf der Grundlage von Peak-Überlappungen (berechnet mit dem R-Paket GenomicRanges) zwischen den vorhergesagten Peak-Dateien und denen der entsprechenden, verfügbaren ChIP-seq-Peaks (Additional file 1: Table S9, S11) und stellen fest, dass diese weitgehend mit der vorherigen Bewertung auf der Grundlage der abgeleiteten Bezeichnungen übereinstimmen.

Schließlich erlauben die Daten für CTCF einen Vergleich der Überlappung zwischen vorhergesagten Peak-Listen und experimentell bestimmten Peak-Listen mit den Überlappungen, die für (i) technische Replikate (Zusatzdatei 1: Tabelle S12) und (ii) biologische Replikate (Zusatzdatei 1: Tabelle S10) beobachtet wurden. Wir stellen fest, dass die Überschneidungen zwischen Vorhersagen und IDR-geschwächten Peaks geringer sind als die zwischen IDR-geschwächten Peaks und/oder technischen Replikaten. Für CTCF sind drei unabhängige Experimente für „Vorhautfibroblasten“-Gewebe verfügbar, und wir verwenden zwei unabhängige DNase-seq-Proben für dieses Gewebe für unsere Vorhersage. Beim Vergleich der Jaccard-Koeffizienten in diesen beiden Situationen (vgl. Additional file 1: Tabellen S9, S10) stellen wir fest, dass die Jaccard-Koeffizienten zwischen Vorhersagen und IDR-geschwächten Peaks zwischen 0,568 und 0,693 variieren, während wir Jaccard-Koeffizienten zwischen 0,658 und 0,72 für biologische Replikate beobachten. Auf der Grundlage dieser begrenzten Daten können wir schlussfolgern, dass die rechnerischen Vorhersagen nur geringfügig weniger konsistent sind als die biologischen Replikate, zumindest für CTCF.

Auf der Grundlage der vorhergesagten Peak-Listen können wir auch die vorhergesagten Bindungseigenschaften der verschiedenen TFs zwischen den Zelltypen vergleichen. Zunächst untersuchen wir die Anzahl der vorhergesagten Peaks pro TF und Zelltyp (Zusatzdatei 1: Abbildung S7). Wir finden eine eindeutige Gruppe von sehr häufig vorkommenden TFs (CTCF, GATA3, SPI1, CEBPB, FOXA1, FOXA2, MAX), die typischerweise auch eine große Anzahl von Peaks in den Trainingsdaten aufweisen. Unter diesen TFs finden wir zelltypspezifische Muster, die vom ubiquitär vorkommenden CTCF bis zu einer stark variierenden Häufigkeit von GATA3 reichen. Die übrigen TFs weisen eine wesentlich geringere Anzahl vorhergesagter Peaks mit ähnlichen Mustern auf, z. B. für ATF7/ARID3A/NANOG oder EP300/TEAD4/JUND, wobei letztere Gruppe Berichten zufolge in distalen Enhancern mitbindet. Als nächstes untersuchen wir die Stabilität der Peak-Vorhersagen, d. h. die Jaccard-Koeffizienten der Peaks, die für jeden der TFs in verschiedenen Zelltypen vorhergesagt wurden (Additional file 1: Abbildung S8). Auch hier finden wir erhebliche Unterschiede zwischen den TFs, wobei GABPA, CTCF und REST im Median Jaccard-Koeffizienten von über 0,7 aufweisen. Bemerkenswert ist, dass CTCF eine der TFs mit der größten Anzahl vorhergesagter Peaks war (Median 37 455), während wir für REST (Median 3 364) und GABPA (Median 5 430) eine Größenordnung weniger vorhergesagte Peaks beobachtet haben. Am anderen Ende der Skala finden wir indirekt bindende TFs wie EP300, oder TFs, die sehr spezifisch für Zelltypen sind, die in unseren Daten unterrepräsentiert sind, wie NANOG (Stammzellen) und HNF4A (Leber, Niere, Darm). Schließlich untersuchen wir die gemeinsame Bindung von TFs, indem wir den durchschnittlichen Jaccard-Koeffizienten über Zelltypen hinweg für jedes TF-Paar berechnen (Additional file 1: Abbildung S9). Hier beobachten wir verschiedene Gruppen von gemeinsam auftretenden TFs wie CTCF/ZNF143 oder FOXA1/FOXA2, von denen bekannt ist, dass sie in vivo interagieren. Darüber hinaus finden wir eine größere Gruppe von TFs mit erheblichen Überschneidungen zwischen ihren vorhergesagten Peaks, darunter YY1, MAX, CREB1, MYC, E2F6, E2F1 und TAF1. Da TAF1 (TATA-box binding protein associated factor 1) mit der Transkriptionsinitiierung an der TATA-Box in Verbindung gebracht wird, könnte eine Erklärung darin liegen, dass die Bindungsstellen dieser TFs an Kernpromotoren angereichert sind. Tatsächlich wurde die Bindung an proximale Promotoren für MYC/MAX, CREB1, YY1 und E2F-Faktoren nachgewiesen.

Vereinfachte Catchitt-Implementierung bringt konkurrenzfähige Leistung

Wir vergleichen schließlich Catchitt, die vereinfachte Implementierung des iterativen Trainingsansatzes, der Chromatin-Zugänglichkeit und Motiv-Scores kombiniert, mit der Challenge-Implementierung, die DNase-seq-basierte und motivbasierte Merkmale für den Fall innerhalb des Zelltyps verwendet. Zu diesem Zweck wählen wir fünf Kombinationen von Zelltyp und Transkriptionsfaktor aus, die den Bereich der in der Herausforderung beobachteten Leistungswerte abdecken. Konkret betrachten wir NANOG und TAF1, die die niedrigsten AUC-PR-Werte (vgl. Additional file 1: Abbildung S3) für die Challenge-Implementierung erzielten; CTCF in IPSC-Zellen, die den größten AUC-PR-Wert erzielten; und FOXA1 und HNF4A, die mittlere AUC-PR-Werte erzielten, aber erheblich vom iterativen Training profitierten (vgl. Additional file 1: Abbildung S4). Wir fassen die Ergebnisse dieses Vergleichs in Additional file 1: Tabelle S13 zusammen. Trotz einer etwa zehnfachen Reduzierung der Anzahl der berücksichtigten Motive und weiterer Vereinfachungen (Abschnitt „Catchitt: eine optimierte Open-Source-Implementierung“) liefert Catchitt immer noch konkurrenzfähige AUC-PR-Werte. Ein Ranking der Catchitt-Ergebnisse innerhalb der ursprünglichen Ergebnisse der Challenge zeigt, dass die mit den Catchitt-Ergebnissen erzielte Leistung nur zwei Ränge unter der Challenge-Implementierung liegt, die DNase-seq-basierte und motivbasierte Merkmale verwendet. Wie zuvor stellen wir eine erhebliche Verbesserung der Vorhersageleistung aufgrund des iterativen Trainingsverfahrens fest.