Drzewa filogenetyczne i zegar molekularny

Drzewa filogenetyczne i zegar molekularny

Drzewa filogenetyczne, które tworzyliśmy dotychczas pokazywały relacje filogenetyczne pomiędzy badanymi organizmami i ich stopień podobieństwa, jednak nie pozwalały na odniesienie poszczególnych wydarzeń, takich jak rozdzielanie się kladów, do czasu w którym mogły się wydarzyć.

Drzewa posiadające informacje na temat długości gałęzi (filogramy) pozwalają szacować względne odległości między taksonami ale aby dopasować je do ram czasowych potrzebujemy sposobu ich ,,wykalibrowania'' tak aby można było znaleźć relacje pomiędzy liczbą mutacji a czasem.

Kalibracja drzewa może przebiegać na przynajmniej dwa sposoby. Jeśli znamy tempo ewolucji badanych sekwencji to można dopasować skalę czasową na podstawie liczby mutacji na poszczególnych gałęziach. Można też podejść do problemu z drugiej strony. Znając przybliżony czas w którym doszło do rozdzielenia się kladów, na przykład opierając się na danych opartych na skamieniałościach, można dopasować go do odpowiedniego węzła na drzewie i na tej podstawie wyliczyć związek między czasem a liczbą mutacji. Tutaj pokażę prosty przykład tego drugiego podejścia.

W praktyce trzeba pamiętać, że nie tylko różne rodzaje sekwencji ewoluują z różną prędkością, ale także ten sam typ sekwencji może mieć różne tempo ewolucji u różnych grup organizmów. Oznacza to, że różne części naszego drzewa mogą mieć różne tempo ewolucji. Zatem lepiej jest wyznaczyć kilka ,,markerów'' - w przypadku pierwszego podejścia wyznaczyć tempo ewolucji dla różnych kladów, dla drugiego oznaczyć czas dla kilku węzłów.

Do programów pozwalających na użycie zegara molekularnego do drzew filogenetycznych należą obsługiwany z linii komend r8s, nawiązujący do niego, napisany w Pythonie pyr8s, z którego można korzystać używając linii komend lub interfejsu graficznego, oraz ,,okienkowy'' i dużo bardziej złożony BEAST. Na naszych zajęciach będziemy używali pyr8s.

Pyr8s

Strona domowa programu pyr8s znajduje się pod adresem https://github.com/iTaxoTools/pyr8s.

Program pobieramy używając narzędzia git:

git clone git@github.com:iTaxoTools/pyr8s.git

Można też go pobrać wybierając na stronie zielony przycisk <> Code a następnie Download ZIP.

Ewentualnie możemy wykorzystać wget:

wget https://github.com/iTaxoTools/pyr8s/archive/refs/heads/main.zip
unzip main.zip
mv pyr8s-main pyr8s

Następnie instalujemy program:

cd pyr8s
pip install .

Uruchamianie programu bez instalacji jest opisane na stronie projektu.

W pobranym katalogu pyr8s znajdziemy kilka plików i katalogów:

$ls

build          MANIFEST.in     README.md                              tests
launcher.py    pyproject.toml  requirements.txt
launcher.spec  pyr8s           Sanderson_r8s_1_7_original_manual.pdf
LICENSE.txt    pyr8s.egg-info  setup.py

Plik Sanderson_r8s_1_7_original_manual.pdf zawiera manual programu r8s, na którym został oparty pyr8s, można zatem znaleźć wiele informacji na temat użytkowania programu, my skupimy się na podstawach.

Katalog tests zawiera przykładowe pliki wsadowe programu. Nasz przykładowy plik będzie oparty na jednym z nich (legacy_1), ale nic nie stoi na przeszkodzie, żeby wypróbować też inne.

LICENCE i README to pliki, których nazwy jasno wskazują na ich znaczenie. Pozostałe raczej nie będą nas bezpośrednio interesować.

Utwórz katalog roboczy, np. zegar a w nim plik rosliny.set (przedłużenie pliku nie ma specjalnego znaczenia).

Umieść w nim poniższą zawartość:

[ * Blok trees, obowiązkowy ]

begin trees;

tree PLANTS = (Marchantia:0.033817,(Lycopodium:0.040281,((Equisetum:0.048533,(Osmunda:0.033640,Asplenium:0.036526):0.000425):0.011806,((((Cycas:0.009460,Zamia:0.018847):0.005021,Ginkgo:0.014702):1.687e-86,((Pinus:0.021500,(Podocarpac:0.015649,Taxus:0.021081):0.006473):0.002448,(Ephedra:0.029965,(Welwitsch:0.011298,Gnetum:0.014165):0.006883):0.016663):0.006309):0.010855,((Nymphaea:0.016835,(((((Saururus:0.019902,Chloranth:0.020151):1.687e-86,((Araceae:0.020003,(Palmae:0.006005,Oryza:0.031555):0.002933):0.007654,Acorus:0.038488):0.007844):1.777e-83,(Calycanth:0.013524,Lauraceae:0.035902):0.004656):1.687e-86,((Magnolia:0.015119,Drimys:0.010172):0.005117,(Ranunculus:0.029027,((Nelumbo:0.006180,Platanus:0.002347):0.003958,(Buxaceae:0.013294,((Pisum:0.035675,(Fagus:0.009848,Carya:0.008236):0.001459):0.001994,(Ericaceae:0.019136,Solanaceae:0.041396):0.002619):1.687e-86):0.004803):1.687e-86):0.006457):0.002918):0.007348,Austrobail:0.019265):1.687e-86):1.687e-86,Amborella:0.019263):0.003527):0.021625):0.012469):0.019372);
end;

[* Blok z komendami ]

begin rates;

[ * Informacje na temat drzewa]

blformat nsites=952 lengths=persite;

[ * Komenda, która pozwala programowi radzić sobie z gałęziami o zerowej długości ]

collapse;

[* Definicje wybranych elementów drzewa ]

mrca LAND_PLANTS marchantia pisum;
mrca GYMNOSPERMS Ginkgo Gnetum;
mrca ANGIOSPERMS amborella pisum;

[ * Ustawianie sztywnych wartości czasowych dla węzłów ]

fixage taxon=LAND_PLANTS age=450;
constrain taxon=GYMNOSPERMS min_age=310;
constrain taxon=ANGIOSPERMS min_age=200 max_age=245;

[ * Wyliczanie czasów przy użyciu zadeklarowanej metody i algorytmu ]

divtime method=np algorithm=pl;

[ * Czytelna prezentacja wyników ]

showage;

[ * Generowanie chronogramu (tryb tekstowy)- gałęzie proporcjonalne do czasu ]

describe plot=chronogram;

[ * Generowanie filogramu (tryb tekstowy)- gałęzie proporcjonalne do liczby zmian ]

describe plot=phylogram;

[ * Generowanie drzewa w formacie nexus ]

describe plot=tree_description;

[ * Sprawdzamy, czy nie ma innych rozwiązań ]

set num_time_guesses=3;


end;

Na początku pliku znajduje się deklaracja, że jest to plik w formacie NEXUS. Dalej znajduje się komentarz - jak widać ograniczony jest parą nawiasów klamrowych.

Polecenia begin, po którym następuje nazwa oraz end ograniczają blok kodu. W tym przypadku jest to blok trees, niezbędny do działania programu (co jest zrozumiałe). Zauważ, ze komendy kończą się średnikiem.

W bloku znajduje się drzewo (komenda tree), zapisane w formacie newick, które będzie miało nazwę rosliny. Zawartość drzewa została skopiowana z pliku legacy_2. Trzeba tu zaznaczyć, że pyr8s nie wylicza drzew, ale korzysta z już gotowych.

Dopisz do pliku (poniżej komendy end;) blok z poleceniami dla programu r8s:

[* Blok z komendami ]

begin rates;

[ * Informacje na temat drzewa]

blformat nsites=952 lengths=persite;

[ * Komenda, która pozwala programowi radzić sobie z gałęziami o zerowej długości ]

collapse;

[* Definicje wybranych elementów drzewa ]

mrca LAND_PLANTS marchantia pisum;
mrca GYMNOSPERMS Ginkgo Gnetum;
mrca ANGIOSPERMS amborella pisum;

[ * Ustawianie sztywnych wartości czasowych dla węzłów ]

fixage taxon=LAND_PLANTS age=450;
constrain taxon=GYMNOSPERMS min_age=310;
constrain taxon=ANGIOSPERMS min_age=200 max_age=245;

[ * Wyliczanie czasów przy użyciu zadeklarowanej metody i algorytmu ]

divtime method=lf algorithm=tn;

[ * Czytelna prezentacja wyników ]

showage;

[ * Generowanie chronogramu (tryb tekstowy)- gałęzie proporcjonalne do czasu ]

describe plot=chronogram;

[ * Generowanie filogramu (tryb tekstowy)- gałęzie proporcjonalne do liczby zmian ]

describe plot=phylogram;

[ * Generowanie drzewa w formacie nexus ]

describe plot=tree_description;

[ * Sprawdzamy, czy nie ma innych rozwiązań ]

set num_time_guesses=3;
divtime method=lf algorithm=tn;

end;

Komenda blformat opisuje drzewo, które badamy. nsites to liczba miejsc w sekwencjach, które zostały użyte przy tworzeniu drzewa. lengths informuje o znaczeniu liczb określających długość gałęzi: opcja persite oznacza liczbę mutacji na miejsce w sekwencji, jeśli oznaczałyby całkowitą liczbę mutacji to użylibyśmy opcji total.

Komenda collapse, nie wchodząc w szczegóły pozwala radzić sobie programowi w sytuacjach, w których gałęzie mają zerową długość.

Następnie definiujemy/nazywamy, używając komendy mrca dwa wewnętrzne węzły. Nazwa komendy to skrót od recent common ancestor, czyli ostatni wspólny przodek. Podając nazwy dwóch taksonów, określamy, że podana nazwa odnosi się do węzła w którym rozdzieliły się ich drogi ewolucyjne.

W naszym przykładzie oznaczamy węzły ostatnich wspólnych przodków (badanych) roślin lądowych (LAND_PLANTS), nagonasiennych (GYMNOSPERMS) oraz okrytonasiennych (ANGIOSPERMS).

Dalej, wykorzystujemy utworzone wcześniej nazwy do umiejscowienia ich w czasie, co pozwala na ,,wyskalowanie'' drzewa. Komenda fixage ustawia ,,na sztywno'' czas, natomiast constrain pozwala określić zakres czasowy w którym doszło do rozdzielenia się linii ewolucyjnych, poprzez wyznaczenie jednej (maksymalną lub minimalną) lub obu granic. W naszym pliku czas jest określony w milionach lat.

Komenda divtime wylicza czasy dla pozostałych węzłów przy użyciu wybranej metody i algorytmu. Nie będziemy ich analizować. Wyniki zostają wypisane w sposób mało czytelny, dlatego dalej używamy komendy skowage, która je prezentuje w sposób bardziej przejrzysty.

Następnie trzykrotnie zostaje wywołana komenda descriptions z różnymi wartościami parametru plot, dzięki której otrzymujemy wyniki w formie drzew w różnych formatach. Dwa pierwsze drzewa prezentowane są w trybie wizualno-tekstowym: chronogram, w którym długości gałęzi odpowiadają czasowi (dlatego są wyrównane do prawej krawędzi) oraz filogram, posiadający gałęzie o długościach proporcjonalnych do liczby mutacji. Trzecie drzewo jest zapisane w formacie nexus.

Istnieje możliwość, że dla danego zestawu danych istnieje więcej niż jedno optymalne rozwiązanie. Dlatego na końcu sprawdzamy trzykrotnie (można ustawić inną wartość) obliczenia, za każdym razem program (nie wchodząc w szczegóły) używając innych wewnętrznych wartości startowych. Jeśli otrzymujemy dla każdej serii wyliczeń otrzymujemy inne wyniki, oznacza to, że istnieją różne optymalne rozwiązania.

Zapisz plik i uruchom, obejrzyj wyniki na ekranie i spróbuj dopasować je do kolejnych poleceń. Wynik warto też zapisać w pliku:

pyr8s rosliny.set > rosliny.out

Jak widać program wczytał drzewo i zakończył działanie. Teraz skopiuj z pliku wynikowego (lub z terminala) drzewo wynikowe w formacie newick i zapisz je w osobnym pliku (np. rosliny-drzewo.newick). Na początku dodaj tree rosliny =.

tree rosliny = (Marchantia:450.0,(Lycopodium:410.78210385570924,((Equisetum:329.47918630858095,Osmunda:329.47918630858095,Asplenium:329.47918630858095):55.7563560412508,(((Cycas:274.97325536020344,Zamia:274.97325536020344):35.02821286834137,Ginkgo:310.0014682285448,((Pinus:264.6701687673707,(Podocarpac:203.0680348291571,Taxus:203.0680348291571):61.60213393821357):16.010086935898187,(Ephedra:178.72746475802242,(Welwitsch:118.30319962228805,Gnetum:118.30319962228805):60.42426513573437):101.95279094524645):29.32121252527594)GYMNOSPERMS:29.60655820058639,(Nymphaea:244.9959805814245,(Saururus:197.7567107520732,Chloranth:197.7567107520732,((Araceae:124.94255524908729,(Palmae:108.1784605504707,Oryza:108.1784605504707):16.76409469861659):37.556244753081074,Acorus:162.49880000216837):35.25791074990482,(Calycanth:169.95575255556145,Lauraceae:169.95575255556145):27.800958196511743,((Magnolia:133.15916544718178,Drimys:133.15916544718178):45.623197584123574,(Ranunculus:144.5527306406023,(Nelumbo:95.33856419094742,Platanus:95.33856419094742):49.21416644965488,(Buxaceae:119.0911854038647,(Pisum:109.7240316974227,(Fagus:101.75051974022425,Carya:101.75051974022425):7.973511957198454):9.367153706441997,(Ericaceae:110.75160222156798,Solanaceae:110.75160222156798):8.339583182296721):25.461545236737607):34.229632390703046):18.974347720767838):47.239269829351315,Austrobail:244.9959805814245,Amborella:244.9959805814245)ANGIOSPERMS:94.6120458477067):45.627515920700546):25.54656150587749):39.21789614429076)LAND_PLANTS;

Zauważ, że zawiera ono opisy węzłów, które im nadaliśmy a także wyliczone czasy. Brak natomiast wartości odpowiadających długości gałęzi w odniesieniu do liczby mutacji.

Otwórz plik w programie FigTree. Przy pytaniu o nazwę etykiet dla węzłów wpisz ,,Opis''.

Jak zwykle dopasuj czcionki itp. aby drzewo było bardziej czytelne. Następnie zaznacz zakładkę ,,Node Labels'' wyświetlaną wartość na ,,Node ages''. Po dopasowaniu parametrów czcionki i liczby miejsc po przecinku powinieneś uzyskać mniej więcej taki obraz:

Jak widać, na drzewie widoczny jest wyliczony (i ustawiony przez nas) wiek dla poszczególnych węzłów, oznaczający miliony lat, które upłynęły od momentu rozdzielnia się poszczególnych linii ewolucyjnych.

Ustawiając ,,opisy'' w ,,Branch Labels'', otrzymasz drzewo dodatkowo z etykietkami, które przypisaliśmy do węzłów.

Powyższy przykład pokazał podstawy używania programu pyr8s. Studiując dołączone do programu pliki z przykładami a zwłaszcza tutorial można dowiedzieć się znacznie więcej na temat jego możliwości, do czego zachęcam.

Last updated