1 Biomatematyka Dr Wioleta Drobik-CzwarnoŁańcuchy Markowa Biomatematyka Dr Wioleta Drobik-Czwarno
2 Polecane Łańcuchy Markowa wizualnie:
3 Procesy Markowa Ciąg doświadczeń w których wynik kolejnego eksperymentu zależy tylko od wyniku poprzedniego Oznacza to, że dochodząc do każdego stanu, łańcuch „zapomina”, skąd przyszedł, a prawdopodobieństwa przejścia w następnym ruchu zależą tylko od położenia bieżącego Wyniki poszczególnych prób są stanami łańcucha Ciąg zmiennych losowych o wartościach całkowitych nazwiemy łańcuchem Markowa jeśli: jeśli nie zależy od numeru próby (t) to taki łańcuch nazwiemy jednorodnym
4 Modele z czasem dyskretnymCzym różnią się łańcuchy Markowa od równań rekurencyjnych i modeli z czasem dyskretnym? Łańcuchy Markowa: można przewidzieć kolejny stan tylko z pewnym prawdopodobieństwem to w jakim stanie aktualnie znajduje się model też wyznaczamy z pewnym prawdopodobieństwem w dalszych etapach łańcucha Łańcuch Markowa jest to procesem stochastycznym gdyż zmiana stanów ma charakter losowy
5 Łańcuchy Markowa Jednorodny łańcuch Markowa - przykłady:rzut monetą – dwa możliwe wyniki, w każdej próbie, orzeł lub reszka, a prawdopodobieństwo otrzymania wynosi zawsze 0,5 niezależnie od tego jakie wyniki otrzymamy w przeszłości – jednorodny łańcuch Markowa rosyjska ruletka, w której używamy sześciostrzałowego rewolweru z jedną kulą. Za każdym razem kręcimy bębenkiem, więc prawdopodobieństwo, że pistolet wystrzeli jest zawsze jednakowe. Niejednorodny łańcuch Markowa – przykład: Rosyjska ruletka, w której używamy sześciostrzałowego rewolweru z jedną kulą, ale bębenkiem kręcimy tylko raz
6 Rzut monetą Jest prostym procesem stochastycznym, dwa możliwe wyniki O lub R Graf Macierz przejścia – prawdopodobieństwo przejścia z jednego stanu do drugiego Przyszły stan orzeł reszka S0 = stan początkowy orzeł orzeł reszka reszka Suma wierszy jest zawsze równa 1 Obecny stan
7 Jak znaleźć prawdopodobieństwa przejścia?Prawdopodobieństwa przejścia łańcucha Markowa można znaleźć doświadczalnie, zliczając przypadki przejść między stanami Znaleziony łańcuch może być użyty do symulacji badanego zjawiska
8 Co się stanie w przyszłości?Rozkład stacjonarny Rozkład prawdopodobieństw do którego dąży łańcuch Markowa Odpowiednio długo symulowany łańcuch Markowa zbiega do swojego stanu stacjonarnego (jeżeli taki istnieje) To, czy taki rozkład istnieje i jest jednoznacznie wyznaczony, zależy m.in. od rodzaju stanów występujących w danym łańcuchu
9 Jak znaleźć rozkład stacjonarny?Sposób 1: Oblicz Pt, gdzie t jest wysoką potęgą i zwróć jeden z wierszy otrzymanej macierzy Sposób 2: MCMC (Marcov Chain Monte Carlo) Symulujemy wstępnie dużą liczbę kroków łańcucha tak by zbiegł do rozkładu stacjonarnego Dla określonego punktu przez N iteracji zliczamy ilości stanów jakie przyjął łańcuch Za prawdopodobieństwo przyjęcia stanu i przyjmujemy:
10 Pchła Pchła skacze między podłogą, kotem, psem a człowiekiem Graf:Za każdym razem wybiera miejsce docelowe z takim samym prawdopodobieństwem (1/3) Pchła ginie po wskoczeniu na człowieka Graf: Źródło:
11 Macierz przejścia Z definicji łańcucha Markowa wynika, że prawdopodobieństwo przejścia ze stanu i do j jest zawsze jednakowe. Oznaczmy je pij Reprezentacja macierzowa łańcucha Markowa 1 – podłoga 2 – pies 3 – kot – człowiek
12 Zadania Wyznacz macierz przejścia i graf dla sytuacji, w której, pchła skacze po kocie, podłodze i człowieku. Skacząc na człowieka ginie.
13 Modele Markowa w biologiiWykorzystanie w biologii molekularnej: Ewolucja molekularna i konstrukcja drzew filogenetycznych Poszukiwanie charakterystycznych wzorów sekwencji np. sekwencje regulatorowe Analiza struktury białek I wiele, wiele innych
14 Modele ewolucji sekwencjiMakak CACATG Szympans CAAATG Człowiek CACATG Po co tworzymy modele ewolucji sekwencji? Jak bardzo analizowane sekwencje różnią się od siebie? Często brak jest sekwencji przodka Oszacowanie liczby mutacji musi opierać się na różnicach pomiędzy współczesnymi sekwencjami Konieczne jest uwzględnienie wielokrotnych mutacji tej samej pozycji Najprostszy sposób: dopasowanie sekwencji oraz wyznaczenie odsetka pozycji na których obserwujemy różnice
15 Modelowanie sekwencji DNAPrzykład: Jakie będzie prawdopodobieństwo, że kolejnym nukleotydem (xi) będzie A,C,G lub T jeżeli poprzednim (xi-1) było G? Przykładowe prawdopodobieństwa przejścia: P(xi = A | xi-1 = G) = 0,16 P(xi = C | xi-1 = G) = 0,34 P(xi = G | xi-1 = G) = 0,38 P(xi = T | xi-1 = G) = 0,12 Zakończenie sekwencji (E) pozwala na uwzględnienie długości sekwencji w modelu
16 Model Jukesa-Cantora (JC)Model ewolucji sekwencji opracowany przez Jukesa i Cantora w 1969. Opisuje pojedynczą pozycje w dopasowaniu pary sekwencji DNA, na której może znajdować się jeden z nukelotydów A,C,G lub T Stosowany dla niekodujących fragmentów DNA Założenia: Zajście mutacji w jednym miejscu nie wpływa na pozostałe miejsca Jednakowe częstości nukleotydów (f=0,25) Jednakowe prawdopodobieństwo każdej substytucji
17 Model Jukesa-Cantora (JC)Graf: Prawdopodobieństwo, że na pojedynczej pozycji obserwujemy nukleotyd A, C, G lub T w czasie t (t)=(PA, PC,PG,PT) Macierz przejścia: r – tempo mutacji (liczba mutacji po jednym kroku podzielona przez długość fragmentu DNA)
18 Ukryte modele Markowa (HMM)Skrót HMM pochodzi od angielskiego Hidden Markov Models Są one rozszerzeniem definicji łańcucha Markowa. Modelują proces stochastyczny, którego pewne właściwości nie są znane Przyjmujemy, że udaje się obserwować symbole (liczby, znaki) emitowane przez układ, a nie jego stany wewnętrzne Definiuje się stany ukryte oraz stany początkowe i końcowe Prawdopodobieństwo wystąpienia znaku będzie zależeć od: Znaków występujących na poprzedzających do pozycjach Ukrytego stanu w jakim znajduje się model
19 Dekodowanie Chcąc poznać najbardziej prawdopodobne rozmieszczenie poszczególnych stanów w sekwencji musimy przeprowadzić dekodowanie modelu Algorytm Viterbiego Identyfikacja najbardziej prawdopodobnej ścieżki przejścia przez model, czyli sekwencji przejścia przez model stanów ukrytych, której prawdopodobieństwo jest największe Dla każdej pozycji wyznaczamy prawdopodobieństwo jej przynależności do jednego ze stanów
20 Dekodowanie Algorytmy sufiksowy i prefiksowy (ang. forward and backward algorithms) Rozważamy wszystkie możliwe ścieżki w modelu w celu wyznaczenia dla każdej pozycji przynależności do określonego stanu Dla pewnych pozycji w sekwencji prawdopodobieństwo jest bliskie 1, natomiast dla pozostałych zbliża się do 0 Dla pozostałych pozycji prawdopodobieństwo będzie mniej jednoznaczne – mogą one znaleźć się na krańcach jednoznacznie przypisanych fragmentów lub w ich środku
21 Nieuczciwe kasyno Możliwe są dwa stany: kostka uczciwa oraz nieuczciwaKasyno ma dwa rodzaje kości: Uczciwa – symetryczna, prawdopodobieństwo wyrzucenia określonej liczby oczek jest równe zawsze 1/6 Nieuczciwa – Prawdopodobieństwo wyrzucenia 6 jest równe 1/2, natomiast pozostałych oczek 1/10 Możliwe są dwa stany: kostka uczciwa oraz nieuczciwa Układ może zmieniać swój stan z pewnym prawdopodobieństwem, ale samego stanu nie jesteśmy w stanie zaobserwować. Obserwujemy jedynie wyniki rzutu kostką.
22 Wyspy CpG Wyspa CpG to odcinek DNA o długości co najmniej 200 pz, charakteryzujący się zawartością par CG powyżej 50% W dinukleotydach CpG, poza wyspami CpG, C często ulega metylacji co z dużym prawdopodobieństwem prowadzi do deaminacji i mutacji w T
23 Wyspy CpG Występują w 60% sekwencji promotorowych, gdzie metylacja cytozyny mogłaby wyciszać ekspresję genów na poziomie transkrypcji Ponieważ wyspy CpG występują w pobliżu genów często są chronione przed metylacją dzięki specyficznym białkom wiążącym się z DNA lub pozycjonowaniu nukleosomów, co blokuje dostęp metylotransferaz Rozpoznanie wysp CpG może być wykorzystane jako wskazówka dla znalezienia ciekawych biologicznie fragmentów genomu
24 Wyspy CpG a HMM Czy krótka sekwencja DNA pochodzi z wyspy CpG?Zaproponowany model oparty na HMM: Model ma osiem stanów: A+, C+, G+, T+, A−, C−, G−, T−. Stan z plusem np. A+ oznacza , że znajdujemy się w rejonie wyspy CpG i czytamy A Stan z minusem, np. A− ,oznacza, że znajdujemy się poza rejonem wyspy CpG i czytamy A Emitowane symbole: będąc w stanie emitujemy A ξ (gdzie ξ ∈ {−, +} z prawdopodobieństwem 1, każdy inny symbol jest emitowany w z prawdopodobieństwem 0
25 Wyspy CpG a HMM Wymagane jest posiadane danych w których stany są znane, na ich podstawie wyznaczamy prawdopodobieństwa przejścia oddzielnie dla wysp CpG oraz pozostałych fragmentów sekwencji Ile razy nukleotyd i pojawia się za nukleotydem j w sekwencji oznaczonej + (wyspa CpG) Analogiczne obliczenia prowadzimy dla sekwencji oznaczonej - Wzory:
26 Wyspy CpG a HMM W obszarze wysp CpG zaobserwujemy większe prawdopodobieństwo przejścia w C i G niż poza nią Przykład: Posiadając krótką sekwencje x obliczymy p(x) dla każdego łańcucha Markowa szansę p(x|+) i p(x|-)
27 Wyspy CpG Iloraz szans:Jeżeli suma ilorazów szans będzie większa od 0 wnioskujemy, że x pochodzi z wyspy CpG Suma Ilorazów szans dla sekwencji CGCG Wzory:
28 Zadanie: Czy sekwencja CCACG pochodzi z wyspy CpG?
29 Literatura Foryś U Modelowanie matematyczne w biologii i medycynie. Higgs P., Atwood T Bioinformatyka i ewolucja molekularna. PWN. a/2013/08/29/O_polowaniu_na_pchle_i_czekaniu/ Tiuryn J Wstęp do obliczeniowej biologii molekularnej. Wykłady