# 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](https://github.com/iTaxoTools/pyr8s).

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ść:

```
#NEXUS 

[ * 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;

```

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=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;
```

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''.

![Opis węzłów](https://1220929090-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LWbPsK9zO6CYBs0XQug%2F-LWbSZEjjHNxnoMXUUv7%2F-LWbSbcJAxElUVM8xSgT%2Ffigtree-01.png?generation=1547924056964802\&alt=media)

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:

![Wiek węzłów](https://1220929090-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LWbPsK9zO6CYBs0XQug%2F-LWbSZEjjHNxnoMXUUv7%2F-LWbSbcOTU51ndiSQzsc%2Ffigtree-02.png?generation=1547924056565978\&alt=media)

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.

![Drzewo](https://1220929090-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2F-LWbPsK9zO6CYBs0XQug%2Fuploads%2Fgit-blob-d4784457c32123e21b0b427198c46c1ef737da36%2Ffigtree-03.png?alt=media)

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.
