Subskrypcja JoVE jest wymagana do oglądania tego materiału. Zaloguj się lub rozpocznij bezpłatny okres próbny.

Artykuł metodologiczny

Nowy bajesowski algorytm detekcji punktów zmiany w analizie genomowej zróżnicowanych typów danych ChIP-seq

10.7K wyświetleń

DOI:

10.3791/4273

10 grudnia 2012

W tym artykule

Podsumowanie

Nasz algorytm Bayesian Change Point (BCP) opiera się na najnowocześniejszych osiągnięciach w modelowaniu punktów zmiany przy użyciu ukrytych modeli Markowa i stosuje je w analizie danych z sekwencjonowania immunoprecypitacji chromatyny (ChIPseq). BCP sprawdza się dobrze zarówno w przypadku danych o charakterze szerokim, jak i punktowym, lecz szczególnie wyróżnia się w dokładnej identyfikacji trwałych i powtarzalnych wysp rozproszonego wzbogacenia histonów.

Streszczenie

ChIPseq jest powszechnie stosowaną techniką badania oddziaływań białko-DNA. Profile gęstości odczytów generuje się poprzez sekwencjonowanie nowej generacji DNA związanego z białkami i dopasowanie krótkich odczytów do genomu referencyjnego. Regiony wzbogacone objawiają się jako piki, których kształt często drastycznie różni się w zależności od białka docelowego1. Na przykład czynniki transkrypcyjne często wiążą się w sposób specyficzny dla miejsca i sekwencji i mają tendencję do tworzenia pików punktowych, podczas gdy modyfikacje histonów są bardziej rozproszone i charakteryzują się szerokimi, rozmytymi wyspami wzbogacenia2. Niezawodna identyfikacja tych regionów była głównym celem naszej pracy.

Algorytmy analizujące dane ChIP-seq wykorzystują różnorodne metodologie, od heurystyki3-5 po bardziej rygorystyczne modele statystyczne, np. ukryte modele Markowa (HMMs)6-8. Poszukiwaliśmy rozwiązania, które zminimalizowałoby konieczność stosowania trudnych do zdefiniowania parametrów ad hoc, które często obniżają rozdzielczość i ograniczają intuicyjność narzędzia. W odniesieniu do metod opartych na HMM dążyliśmy do ograniczenia procedur szacowania parametrów oraz prostych, skończonych klasyfikacji stanów, które są często wykorzystywane.

Dodatkowo konwencjonalna analiza danych ChIPseq obejmuje kategoryzację oczekiwanych profili gęstości odczytów jako punktowe lub rozproszone, a następnie zastosowanie odpowiedniego narzędzia. Naszym dalszym celem było zastąpienie potrzeby stosowania tych dwóch odrębnych modeli jednym, bardziej wszechstronnym modelem, który potrafi obsłużyć cały zakres typów danych.

Aby zrealizować te cele, skonstruowaliśmy najpierw ramy statystyczne, które naturalnie modelowały struktury danych ChIPseq, wykorzystując nowoczesne osiągnięcia w zakresie HMM9, opierające się wyłącznie na formułach jawnych – innowację tę uznano za kluczową dla osiągnięcia przewagi wydajnościowej. Nasz model HMM, bardziej zaawansowany niż modele heurystyczne, dopuszcza nieskończoną liczbę stanów ukrytych dzięki zastosowaniu modelu bayesowskiego. Wykorzystaliśmy go do identyfikacji prawdopodobnych punktów zmiany gęstości odczytów, które dodatkowo definiują segmenty wzbogacenia. Nasza analiza wykazała, że algorytm bayesowskich punktów zmiany (BCP) charakteryzuje się zredukowaną złożonością obliczeniową, co potwierdzono skróconym czasem uruchomienia oraz mniejszym zapotrzebowaniem na pamięć. Algorytm BCP został z powodzeniem zastosowany zarówno do identyfikacji punktowych szczytów, jak i rozproszonych wysp, wykazując wysoką dokładność przy ograniczonej liczbie parametrów definiowanych przez użytkownika. Zademonstrowało to zarówno jego wszechstronność, jak i łatwość obsługi. W związku z tym uważamy, że może on zostać łatwo wdrożony dla szerokiego zakresu typów danych i użytkowników końcowych w sposób umożliwiający łatwe porównywanie i zestawianie wyników, co czyni go doskonałym narzędziem do analizy danych ChIPseq, wspierającym współpracę i weryfikację wyników między grupami badawczymi. W niniejszej pracy prezentujemy zastosowanie BCP do istniejących danych dotyczących czynników transkrypcyjnych10,11 oraz danych epigenetycznych12, aby zilustrować jego użyteczność.

Protokół

1. Przygotowanie plików wejściowych do analizy BCP

  1. Krótkie odczyty uzyskane z sekwencjonowania (biblioteki ChIP oraz biblioteki kontrolne/input) należy dopasować do odpowiedniego genomu referencyjnego przy użyciu preferowanego oprogramowania do mapowania krótkich odczytów. Zmapowane lokalizacje powinny zostać przekonwertowane do 6-kolumnowego formatu browser extensible data (BED)13 (przeglądarka genomu UCSC, http://genome.ucsc.edu/), gdzie każdy wiersz oddzielony tabulatorem odpowiada jednemu zmapowanemu odczytowi i wskazuje zmapowany chromosom, pozycję startową (liczoną od 0), pozycję końcową (półotwartą), nazwę odczytu, wynik (opcjonalnie) oraz nić.

2a. Rozproszone profile odczytów: wstępne przetwarzanie gęstości odczytów ChIP w celu wykrycia wzbogaconych wysp w danych rozproszonych

  1. Rozszerz zmapowane lokalizacje ChIP oraz input do określonej długości fragmentu, t. j.wielkość fragmentów docelowych podczas trawienia enzymatycznego lub sonikacji DNA, zazwyczaj około 200 bp. Następnie liczba fragmentów jest agregowana w sąsiednich koszykach (bins). Domyślnie rozmiar koszyka jest ustawiony na szacowaną długość fragmentu wynoszącą 200 bp.
  2. Wszelkie możliwe punkty zmiany w zbiorze koszyków o identycznej liczbie odczytów najprawdopodobniej znajdą się na ich zewnętrznych granicach. W związku z tym mało prawdopodobne jest, aby punkt zmiany wystąpił na wewnętrznej granicy między dwoma koszykami o tej samej liczbie odczytów. Należy zatem pogrupować sąsiednie koszyki o identycznej liczbie odczytów na koszyk w jeden blok. t. j.format bedGraph13.

2b. Profile odczytów punktowych: Wstępne przetwarzanie plików BED dla ChIP i Input w celu wykrywania pików w danych punktowych

  1. Zagreguj nakładające się odczyty dla nici plus i minus z ChIP osobno. Gęstości odczytów specyficznych dla nici powinny tworzyć bimodalny profil pików plus i minus. Wybierz pary pików plus/minus o największym wzbogaceniu i wykorzystaj odległość między ich wierzchołkami jako szacunek długości fragmentu biblioteki.
  2. Przesuń odczyty ChIP i input o połowę długości fragmentu w stronę centrum i ponownie oblicz gęstość odczytów dla przesuniętych i scalonych odczytów nici plus i minus. Metodologia szacowania długości fragmentu została zaadaptowana z Zhang, et al.3. Pozycje o identycznej liczbie scalonych odczytów powinny zostać pogrupowane w bloki, analogicznie do kroku 2a.2.

3. Oszacowanie a posteriori średniej gęstości odczytów każdego bloku przy użyciu naszej aproksymacji BCMIX

  1. Gęstość odczytów każdego bloku jest modelowana jako rozkład Poissona, Pois(θt), z parametrem średniej opartym na mieszaninie rozkładów Gamma, Γ(α,β), oraz prawdopodobieństwem a priori wystąpienia punktu zmiany na dowolnej granicy bloków równym p. Warunkowanie Pois(θt) względem G(α,β) w praktyce czyni ten model ukrytym modelem Markowa (HMM) o nieskończonej liczbie stanów. Hiperparametry α, β oraz p należy oszacować za pomocą maksymalnej wiarygodności a posteriori.
  2. Oblicz jawnie estymatory Bayesa dla każdego bloku, θt, jako E(θtZ). Zastąp tradycyjne, lecz czasochłonne filtry forward i backward, często stosowane w HMM, bardziej wydobyczną obliczeniowo aproksymacją Bounded Complexity Mixture w celu oszacowania średnich a posteriori, θc. Wynikowe średnie a posteriori zostaną „wygładzone” do przybliżonego profilu stałego odcinkowo, zatem bloki o identycznych wartościach θc powinny zostać ponownie pogrupowane z aktualizacją współrzędnych granic.

4a. Rozproszone profile odczytów: przetwarzanie końcowe średnich a posteriori na segmenty rozproszonego wzbogacenia

  1. Przyjmij liczbę odczytów wejściowych dla każdego nowego bloku θc jako tempo tła, Pois(λa), i określ stopień wzbogacenia za pomocą prostego testu hipotez, sprawdzając, czy a posteriori średnia wartość ChIP, θc, przekracza określony próg δ. Domyślną wartością δ jest kwantyl 90ty, który jest odpowiedni w większości przypadków.
  2. Połącz sąsiednie bloki θc, które wykazują wzbogacenie, w jeden region i podaj współrzędne złączeń w prostym formacie BED. Alternatywnie można podać wartość θc dla każdego bloku w formacie bedGraph, aby zachować szczegóły wysokiej rozdzielczości szacunków gęstości odczytów.

4b. Punktowe profile odczytów: postprocessing średnich a posteriori w celu wyłonienia kandydatów na piki

  1. Zdefiniuj stopę tła, Pois(λa), jako średnią ze wszystkich liczb odczytów (γ2) i zidentyfikuj wszystkie bloki, które przekraczają próg d. Ponieważ oczekuje się, że piki punktowe będą znacznie bardziej wzbogacone, domyślna wartość δ została ustawiona na 99.th kwantyl Pois(λa).
  2. Ustaw blok o maksymalnej wartości θc jako wierzchołek kandydującego piku i dołącz sąsiednie bloki o podobnej gęstości odczytów (±1 odczyt, aby uwzględnić niewielkie wariacje). Ten połączony obszar definiuje się jako kandydujące miejsce wiązania.
  3. Oblicz λ2 jako średnią liczbę odczytów w kandydującym miejscu wiązania ChIP i przeprowadź test hipotezy względem tła wejściowego, gdzie hipotezą zerową, H0, jest λ1 λ2, a następnie odrzuć H0 na podstawie progu wartości p. Wyeksportuj kandydujące piki w formacie BED.

Dostęp ograniczony. Zaloguj się lub rozpocznij wersję próbną, aby wyświetlić tę treść.

Wyniki

BCP doskonale identyfikuje regiony szerokiego wzbogacenia w danych dotyczących modyfikacji histonów. Jako punkt odniesienia porównaliśmy wcześniej nasze wyniki z wynikami uzyskanymi za pomocą SICER3, istniejącego narzędzia, które wykazało wysoką skuteczność. Aby jak najlepiej zilustrować zalety BCP, zbadaliśmy dobrze poznaną modyfikację histonową, co pozwoliło stworzyć podstawę do oceny wskaźników sukcesu. Mając to na uwadze, przeanalizowaliśmy H3K36me3, ponieważ wykazano, że silnie wiąże się ona z aktywnie tr...

Dostęp ograniczony. Zaloguj się lub rozpocznij wersję próbną, aby wyświetlić tę treść.

Dyskusja

Postawiliśmy sobie za cel opracowanie modelu analizy danych ChIP-seq, który byłby w stanie równie skutecznie identyfikować struktury danych punktowe i rozproszone. Do tej pory trudne było wykrywanie regionów wzbogacenia, w szczególności regionów rozproszonych, które odzwierciedlają założenie o dużej wielkości wysp. Aby rozwiązać te problemy, wykorzystaliśmy najnowsze osiągnięcia w technologii HMM, które oferują wiele zalet w porównaniu z istniejącymi modelami heurystycznymi oraz mniej innowacyjnymi modelami HMM.

Dostęp ograniczony. Zaloguj się lub rozpocznij wersję próbną, aby wyświetlić tę treść.

Oświadczenia

Nie zadeklarowano żadnych konfliktów interesów.

Podziękowania

Nagroda fundacji STARR (MQZ), grant NIH ES017166 (MQZ), grant NSF DMS0906593 (HX).

Dostęp ograniczony. Zaloguj się lub rozpocznij wersję próbną, aby wyświetlić tę treść.

Materiały

Lista materiałów użytych w tym artykule
NazwaFirmaNumer katalogowyKomentarze
Stacja robocza z systemem Linux

Bibliografia

  1. Park, P. J. ChIP-seq: advantages and challenges of a maturing technology. Nat. Rev. Genet. 10, 669-680 (2009).
  2. Barski, A., et al. High-resolution profiling of histone methylations in the human genome. Cell. 129, 823-837 (2007).
  3. Zhang, Y., et al. Model-based Analysis of ChIP-Seq (MACS). Genome Biol. 9, R137(2008).
  4. Zang, C., et al. A clustering approach for identification of enriched domains from histone modification ChIP-Seq data. Bioinformatics. 25, 1952-1958 (2009).
  5. Jothi, R., Cuddapah, S., Barski, A., Cui, K., Zhao, K. Genome-wide identification of in vivo protein-DNA binding sites from ChIP-Seq data. Nucleic Acids Res. 36, 5221-5231 (2008).
  6. Qin, Z. S., et al. HPeak: an HMM-based algorithm for defining read-enriched regions in ChIP-Seq data. BMC Bioinformatics. 11, 369(2010).
  7. Song, Q., Smith, A. D. Identifying dispersed epigenomic domains from ChIP-Seq data. Bioinformatics. 27, 870-871 (2011).
  8. Spyrou, C., Stark, R., Lynch, A. G., Tavaré, S. BayesPeak: Bayesian analysis of ChIP-seq data. BMC Bioinformatics. 10, 299(2009).
  9. Lai, T., Xing, H. A simple Bayesian approach to multiple change-points. Statistica Sinica. , (2011).
  10. Robertson, G., et al. Genome-wide profiles of STAT1 DNA association using chromatin immunoprecipitation and massively parallel sequencing. Nat. Methods. 4, 651-657 (2007).
  11. Stitzel, M. L., et al. Global epigenomic analysis of primary human pancreatic islets provides insights into type 2 diabetes susceptibility loci. Cell Metab. 12, 443-455 (2010).
  12. Bernstein, B. E., et al. The NIH Roadmap Epigenomics Mapping Consortium. Nat. Biotechnol. 28, 1045-1048 (2010).
  13. Karolchik, D., et al. The UCSC Table Browser data retrieval tool. Nucleic Acids Res. 32, 493-496 (2004).
  14. Matys, V., et al. TRANSFAC: transcriptional regulation, from patterns to profiles. Nucleic Acids Res. 31, 374-378 (2003).
  15. Portales-Casamar, E., et al. JASPAR 2010: the greatly expanded open-access database of transcription factor binding profiles. Nucleic Acids Res. 38, D105-D110 (2010).

Dostęp ograniczony. Zaloguj się lub rozpocznij wersję próbną, aby wyświetlić tę treść.

Przedruki i uprawnienia

Tagi

Analiza danych ChIP seqUkryte modele MarkowaWzbogacenie w skali ca ego genomuIdentyfikacja modyfikacji histonowychWi zanie czynnik w transkrypcyjnychModelowanie rozk adu Poissonarednia g sto a posterioriWykrywanie region w wzbogaconychSekwencjonowanie nowej generacji