Krzysztof Spalik   
1.II.2004    
  

Wybór właściwego modelu, czyli diabeł tkwi w szczegółach

Pominęliśmy jak dotąd bardzo ważny etap - w jaki sposób oblicza się cząst-kowe prawdopodobieństwa składające się na wartość funkcji wiarygodności dla określonego drzewa. Aby to obliczyć, musimy przyjąć pewne założenia o substytucji nukleotydów, czyli model ewolucji (neutralnej) na poziomie molekularnym. Szczegółowa prezentacja poszczególnych modeli wykracza poza ramy tego omówienia. Warto jednak wiedzieć, że podstawą każdego modelu jest macierz tempa substytucji.
Najbardziej ogólną postać tej macierzy przedstawia ryc. 4. Tempo substytucji jednego nukleotydu przez drugi zależy od średniego tempa substytucji ľ, stałej dla każdego typu podstawienia oraz częstości podstawianego nukleotydu. Ponieważ nukleotydów jest cztery, w najbardziej ogólnym modelu substytucji mamy 12 typów podstawień (nie licząc podstawień synonimicznych). Zwykle jednak zakłada się pełną odwracalność ewolucji (model GTR - general time reversible). Macierz jest zatem symetryczna wzdłuż przekątnej, a tym samym otrzymujemy sześć typów substytucji. Niemalże wszystkie modele substytucji DNA są specjalnymi przypadkami modelu GTR (ryc. 5).
Na przykład jeśli w modelu GTR ograniczymy liczbę typów przekształceń z sześciu do trzech: transwersji (podstawienia puryny przez pirymidynę i odwrotnie) dwóch typów tranzycji (podstawienia jednej puryny przez drugą i jednej pirymidyny przez drugą), to otrzymamy model Tamury i Nei. Jeśli będziemy rozpatrywać tylko tranzycje i transwersje, to otrzymamy model Hasegawa-Kishino-Yano z 1985 r. albo model Felsensteina z 1984 r. Jeśli te modele uprościmy, zakładając tylko jeden rodzaj podstawienia, to otrzymamy model Felsensteina z 1981 r. A jeśli w tym modelu założymy dodatkowo, że częstości występowania wszystkich nukleotydów są równe, to dojdziemy do najprostszego modelu Jukesa-Cantora. Modele substytucji używane są nie tylko w metodzie największej wiarygodności, ale i w metodach odległościowych.

Ryc. 4. Macierz tempa substytucji nukleotydów (Swofford i in. 1996).
Rzędy i kolumny od lewego górnego rogu dotyczą kolejno adeniny, cytozyny. guaniny i tyminy. Parametry: m - bezpośrednie tempo podstawień; a-l - stałe dla każdego typu podstawienia (razem 12); p - frekwencja danego nukleotydu.

Ryc. 5. Zależności między najpowszechniej używanymi modelami substytucji nukleotydów (Swofford i in. 1996).
Modele: GTR - general time reversible; TrN - Tamura i Nei; HKY85 - Hasegawa-Kishino-Yano 1985; F84 - Felsenstein 1984; SYM - model Zharkikha 1994; K3ST - trójparametryczny model Kimury; K2P - dwuparametryczny model Kimury; JC - Jukes-Cantor.

Macierz substytucji, którą przedstawiliśmy, służy do obliczenia macierzy prawdopodobieństwa zmian jednego nukleotydu w drugi. Ta właśnie macierz jest podstawą wyliczania funkcji wiarygodności.

Ryc. 6. Zmiany kształtu rozkładu ? w zależności od wartości parametru ? (Swofford i in. 1996).

W rozważanym dotychczas modelu założyliśmy, że ewoluują wszystkie miejsca w sekwencji oraz że zmieniają się w takim samym tempie. Założenie takie jest oczywiście błędne. Jeśli to założenie nie jest spełnione, metoda największej wiarygodności jest niespójna, czyli ma taką samą wadę, jak metoda parsymonii. Kiedy bowiem część miejsc pozostaje niezmiennych, to funkcja wiarygodności niedoszacowuje liczbę wielokrotnych podstawień. Tym samym źle wyliczona jest długość gałęzi i błędna filogeneza zostaje oceniona najwyżej. Aby temu zapobiec, można oszacować, jaka część miejsc w sekwencji jest silnie konserwowana i nie przyjmuje żadnych zmian. Można także określić rozkład tempa ewolucji w sekwencji.
Zwykle przyjmuje się, że rozkład ten przybiera postać tzw. rozkładu ?. Rozkład ? jest charakteryzowany przez współczynnik kształtu, określany zwykle literą ?. Na ryc. 6 przedstawiono rozkład częstości tempa podstawiania dla różnych wartości parametru ?. Kiedy ? jest niskie, np. 0,5, zauważymy, że najwięcej jest miejsc, które ewoluują wolno, tzn. tempo ich podstawień (na osi x) jest bliskie zeru. Są jednak także nieliczne miejsca, które ewoluują szybko. Im wyższa wartość ?, tym bardziej ujednolica się tempo ewolucji.
Na przykład dla wartości tego współczynnika równej 200, wszystkie pozycje są podstawiane w tempie zbliżonym do 1.

Wybór możliwości jest zatem bardzo duży - kilka podstawowych typów modeli, każdy z możliwością oszacowania (lub nie) miejsc niezmiennych oraz zróżnicowaniem rozkładu tempa podstawień (z różnym parametrem kształtu rozkładu). Na szacowanie filogenezy wpływają również częstości nukleotydów i tempa poszczególnych typów podstawień. Który model wybrać? Odpowiedź nie jest prosta.
Nie należy się kierować samą wartością funkcji wiarygodności. Drzewa uzyskane za pomocą bardziej złożonych modeli, czyli z większą liczbą parametrów, zawsze mają wyższą wartość tej funkcji niż drzewa bazujące na prostszych modelach. Bardziej złożone modele (z większą liczbą stopni swobody) są jednak wrażliwsze na błąd próby. Innym ograniczeniem jest czas obliczeń. Metoda największej wiarygodności jest najbardziej złożoną obliczeniowo metodą szacowania filogenezy. A im więcej parametrów, tym więcej obliczeń do wykonania.

Podobnie jak w wypadku innych modeli, można stosować do ich porównania test wiarygodności albo np. kryterium informacyjne Akaike'go. Trzeba zaznaczyć, że są to środki pomocnicze - badanie ewolucji nie poddaje się w pełni statystyce, dotyczy bowiem odtwarzania przeszłości.

Literatura

Felsenstein, J. 1978. Cases in which parsimony and compatibility methods will be positively misleading. Systematic Zoology 27, 401-410.
Hendy, M.D. & Penny, D. 1989. A framework for the quantitative study of evolutionary trees. Systematic Zoology 38, 297-309.
Swofford, D.L., Olsen, G.J., Waddell, P.J., & Hillis, D.M. 1996. Phylogenetic inference. In D. M. Hillis, C. Moritz, B.K. Mable (ed.), Molecular systematics. 2nd ed. 407-514. Sinauer Associates, Sunderland.