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

Ogólnym celem następnego eksperymentu jest wykorzystanie gęstości pozycji zmapowanych odczytów z danych sekwencjonowania immunoprecypitacji chromatyny do oszacowania a posteriori średniej gęstości odczytów w obrębie genomu. Osiąga się to poprzez wstępne przetwarzanie zmapowanych odczytów ChIP-seq na blokowe profile gęstości, w których ta sama liczba odczytów przypada na niezakładające się koszyki (bins) o szerokości 200 par zasad.

W drugim kroku sąsiednie koszyki o tej samej gęstości są łączone w większy blok, a następnie rekurencyjnie oblicza się posteriorowe średnie gęstości każdego bloku w kontekście wszystkich otaczających go bloków, wykorzystując model bayesowski z filtrami przednimi i wstecznymi. Liczba odczytów dla bloku jest modelowana za pomocą rozkładu Poissona z parametrem theta, który przyjmuje rozkład gamma a priori z parametrami alpha i beta. Następnie ocenia się istotność posteriorowych estymacji średniej gęstości każdego bloku na podstawie tego, czy przekraczają one 90. kwantyl w odniesieniu do gęstości tła kontrolnego wejściowego, aby wygenerować końcowe wyniki wzbogaconych segmentów genomu; otrzymane wyniki ilustrują proces przejścia od surowych odczytów sekwencjonowania, przez posteriorowe estymacje średniej gęstości odczytów, aż po wzbogacone wyspy w danych ChIP-seq podczas analizy BCP.

Ponadto wyniki pokazują, że BCP przewyższa konkurencyjne narzędzie cer. Główną zaletą tej techniki nad istniejącymi metodami, takimi jak CER, jest to, że BCP wykorzystuje najnowsze postępy w zakresie ukrytych modeli markerowych, dzięki czemu lepiej charakteryzuje niuanse analizy danych chipsy niż poprzednie metody heurystyczne. Metoda ta może pomóc w rozwiązaniu kluczowych pytań w dziedzinie epigenomiki, takich jak rola modyfikacji histo, poprzez charakteryzację ich wzorców wzbogacenia w całym genomie.

Chociaż opisana metoda analizy pacjenta może dostarczyć wglądu w analizę danych ChIP-seq, podstawowy schemat może zostać zastosowany również do analizy innych danych sekwencjonowania następnej generacji, takich jak identyfikacja regionów o zróżnicowanym stopniu metylacji w danych sekwencjonowania bisulfitowego, nowych loci transkrypcyjnych w RNA-Seq, wariacji liczby kopii lub dowolnej liczby danych z matryc mikromacierzy (microarray tiling data). Wizualna demonstracja tej metody jest kluczowa dla jasnego zrozumienia metodologii oraz jej zalet. Zalety teoretyczne są zawarte w oprogramowaniu.

Wszystkie zaprezentowane tutaj kroki proceduralne zostały zawarte w jednym pliku wykonywalnym w pakiecie oprogramowania BCP, który jest dostępny do pobrania w tym filmie. Kroki wykonywane przez program są opisane w celu uruchomienia oprogramowania. Wymagane są trzy parametry.

Plik zawierający unikalnie zmapowane odczyty z próbki ChIP oraz analogiczny plik dla odczytów kontrolnych (input), a także nazwa pliku wyjściowego w celu przygotowania plików wejściowych do analizy BCP. Najpierw należy zmapować krótkie odczyty uzyskane z sekwencjonowania do odpowiedniego genomu referencyjnego, korzystając z preferowanego oprogramowania do mapowania krótkich odczytów. Zmapowane lokalizacje powinny zostać przekonwertowane na sześciokolumnowy format BED (browser extensible data) – plik tekstowy z tabulatorami, gdzie każda linia odpowiada jednemu zmapowanemu odczytowi i zawiera: zmapowany chromosom, pozycję startową, pozycję końcową, nazwę odczytu, wynik (score) oraz nić.

Rozszerz lokalizacje mapy chipa i wejść do określonej długości fragmentu. Na przykład rozmiar fragmentu docelowego podczas trawienia enzymatycznego lub sonikacji DNA, zazwyczaj około 200 par zasad. Następnie liczba fragmentów jest agregowana w sąsiednich przedziałach.

Domyślnie rozmiar koszyka (bin size) jest ustawiony na szacowaną długość fragmentu wynoszącą 200 par zasad. Wszelkie możliwe punkty zmiany w zestawie koszyków o identycznej liczbie odczytów najprawdopodobniej wystąpią na zewnętrznych granicach. W związku z tym jest mało prawdopodobne, aby punkt zmiany pojawił się na wewnętrznej granicy między dwoma koszykami o tej samej liczbie odczytów.

W związku z tym należy grupować sąsiednie koszyki o identycznej liczbie odczytów na koszyk w jeden blok. Po przygotowaniu plików wejściowych wywołaj szacowanie BCP, wpisując polecenie widoczne w dolnej części ekranu. Gęstość odczytów każdego bloku jest modelowana jako rozkład Poissona z parametrem średniej theta, który odpowiada mieszaninie rozkładów gamma z parametrami alfa i beta oraz prawdopodobieństwem a priori wystąpienia punktu zmiany w dowolnym bloku.

Warunki brzegowe P dla każdego bloku w ten sposób skutecznie tworzą nieskończony z ukrytym stanem model Markowa lub HMM. Hiperparametry alfa, beta i P są szacowane przy użyciu maksymalnego prawdopodobieństwa a posteriori. Estymaty bayesowskie są obliczane wprost dla każdego bloku theta sub T jako wartość oczekiwana theta sub T przy zadanym why sub T; bardziej tradycyjne, lecz czasochłonne filtry forward i backward, często stosowane w HMM, zostają zastąpione bardziej wydobną obliczeniowo przybliżoną mieszanką o ograniczonej złożoności w celu oszacowania średnich a posteriori theta hat sub T. Wynikowe średnie a posteriori zostaną wygładzone do przybliżonego profilu stałego odcinkowo, zatem bloki z identycznymi theta hat sub T powinny zostać dalej połączone w grupy z zaktualizowanymi współrzędnymi granic.

BCP wykorzystuje liczbę odczytów wejściowych na blok jako poziom tła i określa stopień wzbogacenia. Wykorzystuje do tego prosty test hipotez, oparty na sprawdzeniu, czy średnia gęstość w pozycji chipa dla danego bloku przekracza określony próg istotności. Domyślnym progiem jest 90. kwantyl, który jest odpowiedni w większości przypadków.

Następnie BCP łączy sąsiadujące bloki z najwyższą gęstością średnią posterior, które przekraczają poziom wzbogacenia, w jeden region i raportuje złączone współrzędne w przeglądarce. Format danych BCP o rozszerzalnej strukturze doskonale radzi sobie z identyfikacją obszarów szerokiego wzbogacenia w danych dotyczących modyfikacji histonów. Tutaj wyniki BCP są porównywane z wynikami cser, istniejącego narzędzia, które wykazało wysoką wydajność; wcześniejsze prace tego laboratorium nad tym samym tematem, badające trymetylację H3K36, wykazały tendencję do znacznie większego rozmiaru wysp w BCP niż w cser.

Większe wyspy są bardziej zgodne z konwencjonalnym oczekiwaniem dotyczącym szerokich, rozproszonych wysp wzbogacenia o trimetylację H3K36. Większe wyspy nie są jednak jedynym wskaźnikiem dokładności. Dlatego też, znana asocjacja wysp trimetylacji H3K36 z ciałami aktywnie transkrybowanych genów, a także ich wzajemna wykluczalność z wyspami trimetylacji H3K27, została wykorzystana do oceny wydajności BCP i CER w porównaniu z CER. BCP identyfikował większe, ciągłe wyspy, które lepiej obejmowały ciała genów, nie zwiększając przy tym stopnia nakładania się z wyspami trimetylacji H3K27.

BCP utrzymuje wysoki stopień pokrycia aktywnych genów przez wyspy trimetylacji H3K36, z granicami ściśle dopasowanymi do ciał genów, bez zwiększania stopnia fałszywie dodatniego pokrycia z genami w przestrzeniach międzygenowych o stłumionej transkrypcji lub z represywnym markerem trimetylacji H3K27. Podczas oceny powtarzalności wyznaczania wysp przez BCP w dwóch powtórzonych zestawach danych zaobserwowano, że BCP nie wykazuje silnej zależności od głębokości pokrycia odczytami, co różni go od konkurencyjnego algorytmu cser. Dodatkowe dowody na odporność i powtarzalność BCP dostarcza analiza kolejnych odrębnych regionów, wykazująca spójne granice wysp mimo zmniejszonej głębokości pokrycia. Aby w pełni zademonstrować wszechstronność BCP, pozyskano szerokie spektrum danych o modyfikacjach histonów, obejmujące markery punktowe: acetylację H3K27, acetylację H3K9 oraz trimetylację H3K4, a także marker rozproszony: trimetylację H3K9, oprócz trimetylacji H3K27 i trimetylacji H3K36. Zestawy danych te przeanalizowano przy użyciu domyślnych ustawień parametrów zarówno dla BCP, jak i cser.

W centrum znajduje się wzbogacenie trimetylacji H3K36 w genie PXDN, co oznacza aktywną transkrypcję; zgodnie z oczekiwaniami, w miejscu rozpoczęcia transkrypcji znajdują się dodatkowe punktowe markery aktywności: acetylacja H3K27, acetylacja H3K9 oraz trimetylacja H3K4. Bezpośrednio za genem PXDN znajduje się represyjna przestrzeń międzygenowa, charakteryzująca się wzbogaceniem trimetylacji H3K27; po przeciwnej stronie znajduje się gen represyjny z trimetylacją H3K27. Przechodząc o jeden krok dalej.

Nasza wyciszona chromatyna, wskazana przez wzbogacenie o trimetylację H3K9, wydaje się wskazywać na wyciszenie SN TG2 oraz MYT1L, być może w sposób mniej przejściowy niż represja poprzez trimetylację H3K27. Region ten obejmuje większość zjawisk napotkanych podczas analizy ChIP-seq modyfikacji histonowych. Ilustruje on, w jaki sposób dynamiczny charakter BCP pozwala na identyfikację zarówno punktowej acetylacji, jak i znaczników trimetylacji H3K4, przy jednoczesnym rozróżnianiu dużych, ciągłych wysp represji poprzez trimetylację H3K27 i H3K9, a także aktywnej transkrypcji związanej z trimetylacją H3K36.

Wykonanie tego algorytmu zajmuje około 30 minut, w zależności od liczby odczytów i wyników analizy genomu. Procedura ta nie wymaga znaczącej optymalizacji, która jest często niezbędna w przypadku innych metod. Dzięki zastosowaniu BBCP można badać wiele różnych białek docelowych w immunoprecypitacji chromatyny, w tym różne modyfikacje histonów oraz transkrypcyjne czynniki wiążące DNA, aby odpowiedzieć na dodatkowe pytania dotyczące mechanizmów epigenomicznych i regulacji genów.

Po obejrzeniu tego filmu powinni Państwo dobrze rozumieć, w jaki sposób BCP jest wykorzystywany do identyfikacji obszarów wzbogaconych w rozproszone marki histonowe w analizie danych ChIP-seq.