3 MODELE STANU
Pomysł połączenia predykcji w kompletną strukturę genu polega na tym, że „gramatyczne” ograniczenia mogą wykluczyć niektóre błędne złożenia eksonów. Gramatyczna struktura problemu została podkreślona przez Davida Searlsa (Searls, 1992; Dong i Searls, 1994), który zaproponował również wykorzystanie metod gramatyk formalnych z informatyki i lingwistyki. Programowanie dynamiczne można często wygodnie opisać za pomocą pewnego rodzaju automatu skończonego stanu (Searls i Murphy, 1995; Durbin et al., 1998). Model może mieć stan dla początku translacji (S), jeden dla miejsc donorowych (D), jeden dla miejsc akceptorowych (A) i jeden dla zakończenia translacji (T). Za każdym razem, gdy następuje przejście z jednego stanu do drugiego, dodawany jest wynik (lub kara). W przypadku przejścia ze stanu donora do stanu akceptora, wynik intronu jest dodawany do całkowitego wyniku, i tak dalej. Na rysunku 11.2 pokazany jest diagram stanów dla powyższego prostego algorytmu programowania dynamicznego. Dla każdej zmiennej w algorytmie istnieje odpowiadający jej stan o tej samej nazwie, a ponadto potrzebny jest stan początkowy i końcowy.
Rysunek 11.2. Automat stanów skończonych odpowiadający prostemu algorytmowi DP.
Zaletą takiego sformułowania jest to, że programowanie dynamiczne dla znalezienia maksymalnego wyniku (lub minimalnej kary) jest bardziej ogólnego typu, a zatem dodawanie nowych stanów lub nowych przejść jest łatwe. Na przykład narysowanie diagramu stanów dla bardziej ogólnego algorytmu programowania dynamicznego, dopuszczającego dowolną liczbę genów, a także genów częściowych, jest proste (rysunek 11.3), natomiast zapisanie go jest kłopotliwe. Podobnie diagram stanów dla naszkicowanego powyżej algorytmu frame-aware przedstawiony jest na rysunku 11.4.
Rysunek 11.3. Model z rysunku 11.2 z dodanymi przejściami, które umożliwiają przewidywanie dowolnej liczby genów i genów częściowych, w których sekwencja zaczyna się lub kończy w środku eksonu lub intronu.
Rysunek 11.4. Model zapewniający spójność ramki w całym genie. Podobnie jak na dwóch poprzednich rysunkach, linie kropkowane odpowiadają regionom intergenicznym, przerywane – intronom, a pełne – regionom kodującym (eksonom).
Jeżeli stosowaną punktacją są logiczne prawdopodobieństwa lub logiczne szanse, to automat skończonego stanu jest w istocie ukrytym modelem Markowa (HMM), a te zostały ostatnio wprowadzone do poszukiwania genów przez kilka grup. Jedyną zasadniczą różnicą w stosunku do schematów programowania dynamicznego omówionych powyżej jest to, że modele te są w pełni probabilistyczne, co ma pewne zalety. Jedną z tych zalet jest to, że problem ważenia jest łatwiejszy.
VEIL (Henderson et al., 1997) jest zastosowaniem HMM do problemu znajdowania genów. W tym modelu wszystkie sensory są HMM-ami. Moduł eksonów jest w zasadzie niejednorodnym łańcuchem Markowa pierwszego rzędu, co zostało opisane powyżej. Jest to naturalny porządek dla implementacji w HMM, ponieważ wtedy każde z warunkowych prawdopodobieństw niejednorodnego łańcucha Markowa odpowiada prawdopodobieństwu przejścia z jednego stanu do drugiego w HMM. W modelu pierwszego rzędu nie da się uniknąć kodonów stop w ramce odczytu, ale w VEIL w sprytny sposób dodaje się kilka dodatkowych stanów, dzięki czemu prawdopodobieństwo wystąpienia kodonu stop wynosi zero. W podobny sposób tworzone są czujniki miejsc splice’owania. Poszczególne moduły są następnie łączone zasadniczo tak, jak na rysunku 11.2 (tzn. nie jest wymuszana spójność ramek). Połączony model jest jednym wielkim HMM, a wszystkie przejścia mają przypisane prawdopodobieństwa. Prawdopodobieństwa te mogą być oszacowane na podstawie zestawu danych treningowych metodą największej wiarygodności. W przypadku łączenia modeli sprowadza się to zasadniczo do zliczania wystąpień różnych typów przejść w zbiorze danych. W związku z tym, niejawne ważenie poszczególnych sensorów nie jest tak naprawdę problemem.
Mimo, że sposób znajdowania optymalnej struktury genu jest podobny w duchu do powyższego programowania dynamicznego, w praktyce wygląda zupełnie inaczej. Dzieje się tak dlatego, że programowanie dynamiczne odbywa się na poziomie poszczególnych stanów we wszystkich podmodelach, a takich stanów w VEIL jest ponad 200. Ponieważ model jest w pełni probabilistyczny, można obliczyć prawdopodobieństwo dowolnej sekwencji stanów dla danej sekwencji DNA. Taka sekwencja stanów (zwana ścieżką) decyduje o przyporządkowaniu eksonów i intronów. Jeśli ścieżka przechodzi przez model egzonu, to ta część sekwencji jest oznaczana jako egzon; jeśli przechodzi przez model intronu, to jest oznaczana jako intron, itd. Algorytm programowania dynamicznego, który nazywany jest algorytmem Viterbiego, znajduje najbardziej prawdopodobną ścieżkę przez model dla danej sekwencji i na tej podstawie wyprowadzana jest przewidywana struktura genu (patrz Rabiner (1989) dla ogólnego wprowadzenia do HMM).
Ten model probabilistyczny ma tę zaletę, że rozwiązuje problem ważenia poszczególnych sensorów. Można wykazać, że estymacja parametrów metodą największego prawdopodobieństwa jest optymalna, jeśli istnieje wystarczająca ilość danych treningowych i jeśli statystyczny charakter genów może być opisany przez taki model. Słabym elementem VEIL jest model eksonów pierwszego rzędu, który prawdopodobnie nie jest w stanie uchwycić statystyki regionów kodujących, a większość innych metod wykorzystuje modele czwartego lub piątego rzędu.
Obecnie rozwijana jest wyszukiwarka genów oparta na HMM o nazwie HMMgene. Podstawowa metoda jest taka sama jak VEIL, ale zawiera kilka rozszerzeń standardowej metodologii HMM, które zostały opisane przez Krogha (1997). Jednym z najważniejszych jest modelowanie regionów kodujących za pomocą niejednorodnego łańcucha Markowa czwartego rzędu zamiast łańcucha pierwszego rzędu. Odbywa się to poprzez niemal trywialne rozszerzenie standardowego formalizmu HMM, które pozwala na dowolnego rzędu łańcuch Markowa w stanie modelu, podczas gdy standardowy HMM ma prosty bezwarunkowy rozkład prawdopodobieństwa na czterech bazach (odpowiadający 0. rzędowi). Model jest świadomy ramek i może przewidywać dowolną liczbę genów i genów częściowych, więc ogólna struktura modelu jest taka jak na rysunku 11.4 z dodanymi przejściami, aby umożliwić rozpoczęcie i zakończenie w intronach, jak na rysunku 11.3.
Jak już wspomniano, metoda estymacji z maksymalnym prawdopodobieństwem działa dobrze, jeśli struktura modelu może opisać prawdziwą statystykę genów. Jest to założenie bardzo wyidealizowane, dlatego HMMgene stosuje inną metodę estymacji parametrów zwaną warunkową maksymalną wiarygodnością (Juang i Rabiner, 1991; Krogh, 1994). Luźno mówiąc, maksymalne prawdopodobieństwo maksymalizuje prawdopodobieństwo sekwencji DNA w zbiorze treningowym, podczas gdy warunkowe maksymalne prawdopodobieństwo maksymalizuje prawdopodobieństwo struktur genowych tych sekwencji, czyli to, co nas interesuje. Tego typu optymalizacja jest koncepcyjnie podobna do tej stosowanej w GeneParser, gdzie również optymalizowana jest dokładność predykcji. HMMgene wykorzystuje również do predykcji struktury genu algorytm programowania dynamicznego, inny niż algorytm Viterbiego. Wszystkie te metody przyczyniły się do wysokiej wydajności HMMgene.
Genie jest kolejnym przykładem probabilistycznego modelu stanów, który nazywany jest uogólnionym HMM (Kulp i in., 1996; Reese i in., 1997). Rysunek 11.4 jest w rzeczywistości strukturą stanu Genie, a zarówno ten rysunek jak i rysunek 11.2 są w zasadzie skopiowane z Kulp et al. (1996). W Genie, czujniki sygnału (miejsca splotu) i treści (potencjał kodowania) są sieciami neuronowymi, a wyjście tych sieci jest interpretowane jako prawdopodobieństwo. Interpretacja ta wymaga estymacji dodatkowych parametrów prawdopodobieństwa, które działają jak wagi na sensorach. Tak więc, mimo że model jest sformułowany jako model probabilistyczny, problem wag nadal pojawia się w ukryciu. Algorytm predykcji jest niemal identyczny z algorytmem programowania dynamicznego z ostatniego rozdziału. Wersja Genie uwzględnia również podobieństwa baz danych jako część sensora eksonów (Kulp i in., 1997).
Istnieją dwie główne zalety uogólnionych HMM w porównaniu ze standardowymi HMM. Po pierwsze, poszczególne sensory mogą być dowolnego typu, np. sieci neuronowe, podczas gdy w standardowym HMM są one ograniczone przez ramy HMM. Po drugie, rozkład długości (np. regionów kodowania) może być brany pod uwagę w sposób jawny, podczas gdy naturalnym rozkładem długości dla HMM jest rozkład geometryczny, który maleje wykładniczo wraz z długością. Możliwe jest jednak dość zaawansowane modelowanie długości w HMM, jeśli używa się kilku stanów (Durbin et al., 1998). Z drugiej strony, zaletą systemu takiego jak HMMgene jest to, że jest to jeden zintegrowany model, który może być optymalizowany jednocześnie dla maksymalnej dokładności predykcji.
Innym programem do wyszukiwania genów opartym na uogólnionym HMM jest GENSCAN (Burge i Karlin, 1997). Główne różnice pomiędzy strukturą stanu GENSCAN a strukturą Genie lub HMMgene polegają na tym, że GENSCAN modeluje sekwencję w obu kierunkach jednocześnie. W wielu wyszukiwarkach genów, takich jak te opisane powyżej, geny są najpierw przewidywane na jednej nici, a następnie na drugiej. Jednoczesne modelowanie obu nici zostało przeprowadzone z dużym sukcesem w GeneMark, a podobna metoda jest zaimplementowana w GENSCAN. Jedną z zalet (i być może główną) jest to, że taka konstrukcja pozwala uniknąć przewidywania nakładających się genów na dwie nici, które przypuszczalnie są bardzo rzadkie w ludzkim genomie. GENSCAN modeluje dowolną liczbę genów i genów częściowych, podobnie jak HMMgene. Sensory w GENSCAN są podobne do tych używanych w HMMgene. Na przykład, czujnik kodujący jest niejednorodnym łańcuchem Markowa piątego rzędu. Sensory sygnałowe są w zasadzie zależnymi od pozycji macierzami wagowymi, a więc są również bardzo podobne do tych z HMMgene, ale w modelach miejsc splice’owania występują bardziej zaawansowane funkcje. GENSCAN modeluje również promotory oraz 5′ i 3′ UTR.
.