A Novel bayésienne changement de point Algorithme pour l'échelle du génome Analyse des données ChIPseq divers types de

Summary

Notre point de changement bayésien (BCP) algorithme s'appuie sur l'état de l'art des progrès de la modélisation ruptures par modèles de Markov cachés et les applique à immunoprécipitation de la chromatine séquençage (ChIPseq) l'analyse des données. BCP se comporte bien dans les types de données à la fois larges et ponctuées, mais il excelle dans l'identification précise robustes, reproductibles îles de l'enrichissement d'histone diffuse.

Abstract

ChIPseq est une technique largement utilisée pour étudier les interactions protéine-ADN. Des profils de densité de lecture sont générés à l'aide de la prochaine séquence d'ADN lié à une protéine et l'alignement du lit court à un génome de référence. Régions enrichies sont révélés comme des pics, qui diffèrent souvent de façon spectaculaire dans la forme, en fonction de la protéine cible ^1. Par exemple, des facteurs de transcription se lient souvent dans un site et la séquence de manière spécifique et ont tendance à produire des pics ponctuées, alors que les modifications des histones sont plus omniprésents et sont caractérisés par de larges îles diffuses d'enrichissement ^2. De manière fiable l'identification de ces régions a été au centre de notre travail.

Algorithmes pour l'analyse des données ChIPseq ont eu recours à diverses méthodes heuristiques, à partir de ^3-5 pour plus rigoureuses modèles statistiques, p.ex. modèles de Markov cachés (HMM) ^6-8. Nous avons cherché une solution qui réduit au minimum la nécessité de difficile à définir, ad hoc paramètres qui sont souventcompromettre la résolution et de réduire l'utilisation intuitive de l'outil. En ce qui concerne les méthodes basées sur HMM, nous avons cherché à limiter les procédures d'estimation de paramètres et simples, les classifications à états finis qui sont souvent utilisés.

En outre, l'analyse des données ChIPseq conventionnelle implique catégorisation de la lecture devrait profils de densité que soit ponctuée diffuse ou suivie par l'application ultérieure de l'outil approprié. Nous avons également pour but de remplacer la nécessité pour ces deux modèles distincts avec un seul modèle plus polyvalent, qui peut répondre avec compétence l'ensemble du spectre des types de données.

Pour atteindre ces objectifs, nous avons d'abord construit un cadre statistique qui, naturellement, modelés ChIPseq structures de données en utilisant une avance fine pointe de la HMM ^9, qui utilise uniquement des formules explicites, une innovation cruciale pour ses avantages de performance. Des modèles plus sophistiqués, puis heuristiques, notre HMM accueille infinies états cachés par unModèle bayésien. Nous l'avons appliquée à l'identification des points de changement de la densité raisonnables lire, ce qui a encore définir des segments de l'enrichissement. Notre analyse a révélé la façon dont notre Changement de point de Bayes (BCP) algorithme a une complexité de calcul réduite, témoigne d'un moment de l'exécution abrégée et empreinte mémoire. L'algorithme PCA a été appliquée avec succès à la fois de pointe et ponctuée d'identification île diffuse avec une précision robuste et limités paramètres définis par l'utilisateur. Cette illustré à la fois sa polyvalence et sa facilité d'utilisation. Par conséquent, nous croyons qu'il peut être mis en œuvre facilement à travers de larges gammes de types de données et les utilisateurs finaux d'une façon qui est facilement comparés et contrastés, ce qui en fait un excellent outil pour ChIPseq analyse de données qui peuvent aider à la collaboration et à la corroboration entre les groupes de recherche. Ici, nous démontrons l'application de la BCP au facteur de transcription existante et les données épigénétiques ^{10,11 12} à illustrer son utilité.

Protocol

1. Préparation des fichiers d'entrée pour l'analyse PCA

1. Alignez le lit court produite à partir du séquençage pistes (ChIP et les bibliothèques d'entrée) au génome de référence appropriée en utilisant le logiciel préféré courte alignement lecture. Les emplacements mappées devraient être convertis en des 6 données de la colonne navigateur extensibles (LIT) format ¹³ (UCSC Genome Browser, http://genome.ucsc.edu/ ), une ligne délimité par des tabulations par mappé lecture indiquant le chromosome mappé, la position de départ (0-based), position finale (semi-ouvert), lire le nom, le score (en option), et le volet.

2a. Lire diffuser profils: ChIP Prétraitement Lire densités pour la détection des îles enrichi en données diffuses

1. Étendre la puce et emplacements d'entrée mappées à une longueur prédéterminée fragment, c'est à dire. la taille du fragment cible lors de la digestion enzymatique ou sonication de l'ADN, généralement autour de 200 pb. Compte fragments sont ensuite agrégationted dans les bacs adjacents. Par défaut, la taille bin est réglé sur la longueur du fragment estimé à 200 pb.
2. Tout changement de points possible dans un ensemble de bacs identiques avec chefs lecture sera très probablement à l'automne des limites les plus externes. Par conséquent, il est improbable qu'un point de changement se produira à une limite interne entre deux bacs avec les comtes de lecture mêmes. Alors, les poubelles groupe adjacents identiques, avec des lectures par bac, en un seul bloc, c'est à dire. Format bedGraph ^13.

2b. Ponctuées profils lus: ChIP Prétraitement et fichiers LIT entrée pour la détection de pics dans les données ponctuées

1. Agrégat lectures chevauchantes pour ChIP brin positif et négatif lit séparément. Les brins des densités lues doivent former un profil bimodal de pics positifs et négatifs. Choisissez des paires plus / moins des sommets les plus enrichis et utiliser la distance entre leurs sommets comme une estimation de la longueur des fragments de la bibliothèque.
2. Déplacer la puce et l'entrée se lit la moitié du fragment length du centre et de recalculer la densité de lecture de l', plus décalée et fusionné et moins brin lit. Cette méthode d'estimation de la longueur du fragment a été adoptée à partir de Zhang, et al. ^3. Positions dont la numération de fusion identiques doivent être regroupés en blocs, semblables à l'étape 2a.2.

3. Estimer la moyenne a posteriori Lire la densité de chaque bloc à l'aide de notre rapprochement BCMIX

1. Lire la densité de chaque bloc est modélisée comme une distribution de Poisson, Pois (θ _t), avec une moyenne de paramètre suivant un mélange de distributions Gamma, Γ (α, β), et une probabilité a priori d'un point de changement se produisant à une limite de bloc de p. climatisé Pois (θ _t) sur G (α, β) rend effectivement le modèle HMM un état ​​infini. Estimer l'hyper-paramètres, α, β et p, en utilisant au maximum la probabilité a posteriori.
2. Explicitement calculer les estimations de Bayes pourchaque bloc, θ _t, E (θ _t | γ _Z). Remplacez la consommation plus traditionnelle, mais le temps vers l'avant et vers l'arrière filtres souvent utilisés dans les HMMs, avec le rapprochement Mélange complexité de calcul plus efficacement borné à estimer les moyennes a posteriori, θ _c. Les moyens résultant postérieures sera "lissé" dans un profil approximatif constante par morceaux afin blocs identiques, θ _c, devrait être bloqué avec nouveau contour coordonnées.

4a. Lire diffuser profils: Post-process moyennes a posteriori dans les segments de l'enrichissement diffus

1. Utiliser le numéro d'entrée de lit par chaque nouveau bloc θ _c en tant que le taux d'arrière-plan, Pois (λ _a) et déterminer l'enrichissement à l'aide d'un test d'hypothèse simple basé sur la moyenne si ChIP postérieure, θ _c, δ dépasse un certain seuil. Le 90 ^{e <}/ Sup> quantile est le défaut d est approprié dans la plupart des cas.
2. Fusionner adjacente θ _c blocs qui dépassent l'enrichissement en une seule région et le rapport de fusionner les coordonnées au format LIT simple. Alternativement, on peut signaler le _c θ pour chaque bloc au format bedGraph de préserver les détails en haute résolution des estimations de la densité de lecture.

4b. Ponctuée Lire profils: Post-process moyennes a posteriori en pointe candidats

1. Définir la fréquence de fond, Pois (λ _a), comme étant la moyenne de tous les chiffres de lecture (γ ₂₎ et identifier tous les blocs qui dépassent le seuil, d. Depuis pics ponctuées devraient être plus largement enrichie, le δ par défaut est réglée sur le 99 ^e quantile de Pois (λ _a).
2. Régler le bloc avec le _c maximal θ comme le sommet pic candidat et jouxtent blocs adjacents qui partagent un même den lecturesité (± 1 lire comptent pour permettre une légère variation). Cette région attenante est défini comme un site de liaison candidat.
3. Calculer λ _2, comme les comptes moyens de lecture dans le site candidat ChIP contraignant et test d'hypothèse dans ce contexte d'entrée en fonction ont fait l'hypothèse nulle, H _0, c'est que λ ₁ ≥ λ ₂ et rejeter H ₀ basé sur un seuil de p-value. Pics candidats de sortie au format LIT.

Representative Results

BCP excelle à identifier les régions d'enrichissement large dans les données de modification des histones. Comme point de référence, nous avons déjà comparé nos résultats à ceux de SICER ^3, un outil existant qui a réalisé une excellente performance. Pour mieux illustrer les avantages de BCP, nous avons examiné une modification des histones qui avait été bien étudié pour établir une base pour évaluer les taux de réussite. Dans cet esprit, nous avons ensuite analysé H3K36me3, car il a été démontré à associer fortement avec les organismes de gènes activement transcrits (Figure 1). En revanche, H3K36me3 a également été montré pour être exclusive mutuelle H3K27me3 marques répressives. Nous avons en outre à effet de levier de ces relations connues pour illustrer les avantages de performance de PCA dans l'exactitude des appels île en déterminant la fraction de chevauchement avec des associations connues et dissociations, en corrélation et effet anti-corrélation. Ici, nous avons étayer davantage les avantages de la PCA à l'aide des exemples supplémentairesde haute performance.

Notre travail précédent a démontré une tendance à la taille de l'île beaucoup plus grande dans le BCP, de 23,9 à 25,8 kb, que SICER, de 2,7 à 10,7 kb; plus grandes îles sont plus en ligne avec les attentes classiques des grandes îles diffuses de H3K36me3 enrichissement (PLoS Comp Bio, soumis). Bien sûr, les plus grandes îles ne suffisent pas à indiquer la précision. Donc, nous avons déterminé combien se chevauchent ces régions avaient des gènes connus et une comparaison avec le degré de chevauchement avec l'espace intergénique, une indication du taux de faux positifs (FPR). Couverture gène dans le BCP varient de 0,492 à 0,497 par rapport à 0,276 à 0,437 dans SICER sans impact sévèrement le FPR; intergénique gamme de chevauchement de 0,89 à 0,90 et de 0,85 à 0,98 dans le PCA et SICER, respectivement. Ici, nous présentons une région représentant supplémentaire affichant la relation étroite entre les limites de l'enrichissement et le gène corps-en distinguant clairement actif et réprimered transcription (Figure 1). Cela confirme notre affirmation selon laquelle BCP maintient le chevauchement élevé de gènes actifs par H3K36me3 îles dont les limites sont étroitement liés aux gènes des organismes sans augmenter le degré de chevauchement de faux positifs avec l'espace intergénique, gènes avec transcription refoulé, ou la marque H3K27me3 répressif.

Alors que l'évaluation de la reproductibilité de la BCP-île appelle dans deux ensembles de données répétées, nous avons remarqué BCP ne souffre pas d'une forte dépendance de la profondeur de la couverture lecture dans l'algorithme de concurrence, SICER. Nous fournissons une preuve supplémentaire de la robustesse et la reproductibilité de la BCP en examinant d'autres régions distinctes montrant les limites île cohérentes en dépit de la profondeur de la couverture réduite (simulée par échantillonnage lit à partir de l'ensemble complet de données) (Figure 2).

Pour bien démontrer la polyvalence de la BCP, nous avons obtenu un large éventail de données de modification des histones, y compris la marque ponctuées H3K27ac, H3K9ac et H3K4me3, et la marque diffuse, H3K9me3, en plus de H3K27me3 et H3K36me3. Nous avons analysé ces ensembles de données à l'aide des paramètres par défaut pour les PCA et SICER (figure 3). Ces marques représentent un large éventail de lire les profils de densité et de nous permettre de se concentrer sur une région qui illustre un grand nombre des fonctions les plus couramment associés avec eux. Au centre se trouve H3K36me3 enrichissement au niveau du gène PXDN marquage transcription active. Tomber on pouvait s'y attendre sur le site de début de transcription sont les autres ponctuations, les marques actives, H3K27ac, H3K9ac et H3K4me3. Juste en aval de PXDN est réprimée espace intergénique marquée par H3K27me3 enrichissement. Sur le flanc opposé se trouve un gène H3K27me3 refoulé. Déplacer un pas de plus sur la chromatine sont réduits au silence, comme indiqué par la présence de H3K9me3 enrichissement qui semble indiquer taire de SNTG2 et MYT1L, peut-être dans un sens moins transitoire, puis H3K27me3 répression. Cette région englobe la majorité des phénomènes encontrés dans ChIPseq des modifications des histones et illustre comment la nature dynamique de la BCP peut identifier à la fois l'acétylation ponctuées et H3K4me3 marques tandis que dans le même temps la distinction grandes îles contiguës de H3K27me3 et H3K9me3 répression et H3K36me3 transcription active. Pour rappel, BCP peut faire une telle toutes ces analyses simplement les réglages par défaut et, comme l'a démontré, encore produire des résultats de qualité, indépendamment du type de données. L'algorithme est également rapide et efficace de la mémoire et, par conséquent, offre une utilité pratique convaincantes.

Figure 1. Diffuser lire les profils de densité des modifications des histones. H3K27me3 (en haut) et H3K36me3 (en bas) illustrent les grandes îles, d'enrichissement diffus fortement associés avec des organismes de gènes (cases vertes). H3K27me3 en corrélation avec les gènes réprimés et l'espace intergénique et anticorrelates avec tr activementles organismes de gènes anscribed. L'inverse est vrai pour H3K36me3. Les données sont visualisés dans le navigateur du génome UCSC ( http://genome.ucsc.edu ).

Figure 2. BCP est robuste et reproductible. Île appelle H3K36me3 en deux répétitions et à l'échantillonnage des profondeurs de 30%, 50 et 70 de la réplique 1 jeu de données complet ont été analysés avec BCP. La réplique seconde, avec une couverture en lecture nettement plus faible, produite appels île similaires et le degré de chevauchement a été fortement conservée quelle que soit l'échantillonnage pourcentage. En outre, la précision îles sont restées comme on le voit dans l'alignement près des frontières avec des annotations RefSeq corps de gènes.

Figure 3. BCP est une versaalgorithme de tuile qui peut être appliquée à tous les types de données histone modifications. BCP et SICER ont été utilisés pour analyser la gamme des types de données, de marques de ponctuations comme H3K27ac, H3K9ac et H3K4me3, de diffuser des marques comme H3K36me3, H3K27me3 et H3K9me3. En utilisant les paramètres par défaut pour les deux algorithmes, les îles BCP capturer la densité enrichi quelle que soit leur ampleur tout en SICER fragmente souvent des régions en plusieurs sous-îles. Même dans le cas très large et diffuse des H3K9me3, BCP a des performances raisonnables.

Discussion

Nous avons cherché à développer un modèle d'analyse de données permettant d'identifier ChIPseq deux structures ponctuées et diffuse des données aussi bien. Jusqu'à présent, les régions de l'enrichissement, notamment des régions diffuses, qui reflètent les attentes présupposé de taille grande île, ont été difficiles à identifier. Pour résoudre ces problèmes, nous avons utilisé les plus récentes avancées en matière de technologie HMM, qui possèdent de nombreux avantages par rapport aux modèles heuristiques existantes et HMMs moins innovants.

Notre modèle permet l'utilisation d'un cadre bayésien avec des formules explicites. Il s'agit d'une distinction cruciale de HMM autres, en ce qu'elle nous permet de calculer les moyennes a posteriori, l'. Devrait lire la densité de chaque segment, avec des calculs simples, plutôt que de s'appuyer sur des simulations longues et coûteuses telles que calcul de Markov méthodes Monte Carlo par chaînes Par conséquent, nos temps de calcul et les besoins en mémoire sont considérablement réduits. Utilisation des clusters de calcul haute performance wie cœur double, 2,0 nœuds Ghz avec 2 Go de mémoire 64 bits pour analyser ~ 23 millions H3K27me3 lit ou ~ 21 millions H3K36me3 lit, BCP a pris moins d'une heure pour l'analyse du génome entier par rapport à plusieurs heures à quelques jours nécessaires à d'autres méthodes. Ces gains de temps peuvent être atteints avec seulement la modeste de 2 Go de mémoire.

En outre, nos conditions de modéliser les différents moyens de chaque segment, soit. Pois (θ), lors d'une distribution Gamma continu. Essentiellement, cela permet d'infinies états possibles pour chaque segment. BCP peut fournir plus de classifications binaires simples de enrichie par rapport au fond et préserve les grandeurs de densité de lecture pour chaque segment par le biais des moyens de sortie postérieures.

Nous utilisons aussi de l'algorithme BCMIX de l'efficacité informatique. Cela permet une recherche quasi exhaustive pour changer des points entre l'enrichissement et le fond de toutes les positions génomiques possibles. Ceci fournit une résolution accrue ne confined par des définitions arbitraires fenêtre, avec peu d'impact sur l'exécution ou demandes de mémoire.

Tout cela est réalisé sans perturber la précision, à la fois en théorie, puisque le modèle est statistiquement rigoureuse et ses résultats convergent pour l'estimateur bayésien, aussi bien dans la pratique, comme nous l'avons démontré ici. La couverture gène de nos résultats suggèrent H3K36me3 les appels île sont très précis, sans empiéter sur l'espace intergénique connu mutuellement exclu ou H3K27me3 enrichissement. Les résultats sont remarquablement reproductible et robuste et a montré peu dépendante de la profondeur de couverture, appelant îles similaires avec une couverture élevée du gène et le FPR faible en dépit des profondeurs d'échantillonnage aussi bas que 30%. BCP a été largement utilisé, sans aucun réglage des paramètres par défaut, pour analyser un large éventail de la modification des histones et de données de facteurs de transcription ChIPseq et obtenu de bons résultats dans tous les cas. Nous espérons que grâce à sa grande précision, la robustesse et la reproductibilité, BCP servira comme un moyen efficaceoutil d'analyse de données, la collaboration et la corroboration dans l'avenir.

Materials

Name	Company	Catalog Number	Comments
Linux-based workstation