BLOK 3: TVORBA  STROMŮ ZE SEKVENCÍ DNA

Programů na konstrukci stromů ze sekvencí DNA je celá řada. Zřejmě nejvšestrannějším programem pro konstrukci stromů z DNA sekvencí je PAUP.

PAUP ( $85  - http://paup.csit.fsu.edu/)

Program pro práci především se sekvencemi DNA. Umožňuje analýzu distančními metodami (neighbor-joining, minimum evolution, metoda nejmenších čtverců), metodou maximální parsimonie a maximální věrohodnosti (heuristic search, star-decomposition a quartet-puzzling) s mnoha obměnami a kombinacemi. Provádí testy (test vhodného nukleotidového složení sekvencí, Kishino-Hasegava )a výpočty různých koeficientů. Vstupem je alignment ve formátu Nexus.

Jak získat alignment ve formátu nexus

Editoru BioEdit zadejte File → export → sequence alignment. Zadejte název souboru a typ souboru Paup/nexus. Soubor s alignmentem obsahuje tzv. DATA block, který začíná begin data; , končí end; a obsahuje alignment s potřebnými informacemi.

begin data;

 dimensions ntax=2 nchar=40;

 format datatype=dna interleave gap=-;

  matrix

     R183  AAGGAAGCACACTTCGGTCA TAGATTAAGCCATGCAAGTG

     R186  AAGGAAGCACACTTCGGTGA TTGATTAAGCCATGCAAGTG ;

  end;

Práce v PAUPu (platí pro Windows verzi)

•          Alignment otevřeme v PAUPu volbou File → Open a ovládáme vše z příkazové řádky. Pracuje se ve třech módech (distanční, parsimonní, věrohodnostní). Mód se změní příkazem set criterion=distance/parsimony/likelihood. Je třeba dávat pozor, ve kterém módu se nacházíte. Seznam použitelných příkazů získáte, zadáte-li otazník. Nápovědu pro použití příkazu získáte, zadáte-li jméno příkazu mezeru a otazník (např. set ?). PAUP má podrobný manuál.

•          Další možností je vepsat přímo do souboru sled příkazů, které se hned po zpuštění automaticky provedou. Je to možné udělat buď v PAUPu pokud soubor zpustíme v režimu edit (dialogové okno Open, vlevo dole) nebo v textovém editoru jako je Notepad. Příkazy pro konstrukci stromu se uvádějí v rámci PAUP bloku, který začíná begin Paup;  a končí end;. Příklady takových PAUP bloků jsou uvedeny níže.

•          Kromě DATA a PAUP existují další bloky (ASSUMPTION, CODON, TREE...) více v manuálu.

PAUP block na konstrukci a bootstrap stromu metodou  Neighbor-joining

begin Paup;                                                                                                                 z příkazového řádku netřeba zadávat

set criterion=distance;                                                                                                přepnutí na distanční mód

dset distance=k2p;                                                                                                     nastavení distančního koeficientu na k2p, možno měnit i další parametry modelu

nj;                                                                                                                                konstrukce stromu metodou neighbor-joining

savetrees file=stromnj1.tre brlens=yes;                                                                     uloží strom do souboru strom1 s délkami větví

bootstrap nrep=1000 search=nj keepall=yes treefile=stromnj2-bstrees.tre;           bootstrap stromu, při kterém budou do stromu zařazeny i nody s nejvyšší podporou, i když je méně než 50%. Stromy z jednotlivých replikátů bootstrapingu ukládá do souboru stromnj2-bstrees.tre                                             

savetrees file=strom2 from=1 to=1000 savebootp=nodelabels;                              uloží strom do souboru strom2 s hodnotami bootstrapu uvedených na nodech.

 end;

PAUP block na konstrukci a bootstrap stromu metodou  maximální parsimonie

begin Paup;

set criterion=parsimony;

hsearch addseq=random nrep=10 swap=TBR;                                               hledání počátečního stromu metodou heuristic search s náhodným přidáváním taxonů a provedení swapu (hledání lepšího stromu                                                                                                                                                                   v okolí přehazováním větví) to vše s deseti opakováními a výběrem nejlepšího stromu ze všech

showtrees;                                                                                                        ukáže strom.

savetrees file=strom1mp.tre brlens=yes;

bootstrap search=heuristic nrep=1000 keepall=yes strommp2-bstrees.tre;  bootstrap, parametry hledání zůstávají od příkazu hsearch. Parametry hsearch možno měnit za lomítkem – např.  addseq=asis swap=TBR

savetrees file=strom2 from=1 to=10000 savebootp=nodelabels;                   hodnotu za to= je u parsimonie vhodné nastavit vyšší než počet replikátů bootstrapu, neboť hrozí, že parsimonie nalezne více stejně parsimoních stromů než je počet bootstrapových replikátů.

end;

Zobrazení, editování a tisk stromů

Program Treeview (freeware – http://taxonomy.zoology.gla.ac.uk/rod/rod.html)

Velmi jednoduchá obsluha. Umožňuje různá znázornění stromu (ikonky pod hlavním menu), zakořenění outgroupy (tree -> define otgroup a tree -> root with outgroup ) a export do grafických programů ve formátu wmf (tree -> print trees -> picture).

Práce s gapy

Pokud je nevyhodíte rovnou při úpravě alignmentu je možné s nimi různě zacházet při analýze.

1.         V programu PAUP je u maximální parsimonie je možno zvolit:

1.         Považování gapů za chybějící informace příkaz:  pset  gapmode=missing

2.         Považování gapů za pátý stav sekvence příkaz: pset  gapmode=newstate

Nejlepší je, jako ostatně vždy, vyzkoušet všechny možnosti a přesvědčit se jestli to mění výsledek nebo ne.

Zahrnutí rozdílné substituční rychlosti v různých pozicích do analýzy metodou maximální věrohodnosti, dvě možnosti:

1. Počítat s gama (Γ) rozložením této rychlosti a integrace věrohodnosti každé pozice podle tohoto rozložení. V praxi se rozdělí toto rozložení do několika (stačí i čtyři) diskrétních kategorií a věrohodnost pro každou pozici se získá součtem věrohodností pro každou kategorii vážených pravděpodobností, že do této kategorie pozice patří. Je potřeba znát koeficient α funkce Γ. Některé programy jej umí určit (Treepuzzle, PAUP).

2. Rozdělit přímo konkrétní pozice alignmentu do kategorií a určit u nich relativní substituční rychlost. Treepuzzle to umí. Je možné zadat rozdělení alignmentu do kategorií a jejich rychlosti i v Phylipu (tato metoda „site specific rates“ se moc neosvědčuje, do analýzy vstupuje velké množství parametrů a tím i větší chyba).

Určování hodnot parametrů modelů, dvě možnosti:

1.         Pomocí maximální věrohodnosti. Vyberou se takové hodnoty parametrů, které maximalizují věrohodnost dané topologie. Jako počáteční topologie se používá strom vytvořený metodou neighbor-joining nebo maximální parsimonie bez jakýchkoliv parametrů. Většinou to stačí. Přesnější je provést iteraci, tj.  vytvořit nový strom na základě určených hodnot parametrů a použít jej pro určení nových, přesnějších, hodnot parametrů. Opakovat to dokud se budou hodnoty parametrů měnit. Tento postup se používá častěji (program Modeltest, PAUP – příkaz Lscores, Treepuzzle).

2.         Určit je odjinud. Např.  parametr α funkce Γ lze odhadnout proložením negativně binomického rozdělení do rozdělení počtu změn pro jednotlivé pozice určených maximální parsimonií.

Výběr vhodného modelu (distanční metody a maximální věrohodnost)

S přibývajícími parametry model lépe vystihuje substituční procesy, ale roste rozptyl výběru výsledných hodnot. To vnáší do výsledku chybu náhody, která je větší při větším počtu OTU a kratších sekvencích. Tato chyba hraje větší roli u distančních metod než u maximální věrohodnosti. Je potřeba vybrat optimální model, ne ten nejsložitější.

Distanční metody 

Používat co nejjednodušší modely, raději bez gama korekce na rozdílné substituční rychlosti v různých pozicích (i když je to teoreticky možné). Dvě doporučení lidí co tomu rozumí.

Nei a Kumar:

1.         Jestliže je (d) vypočítané Jukes-Cantorovým modelem asi 0,05 nebo méně, používat Jukes-Cantorův model nebo p distanci.

2.         Jestliže je (d) větší,  je velký rozdíl v počtu tranzicí a tranzverzí (R>5) a délka sekvencí je veliká  je možno použít Kimura 2 parametrovou distanci nebo gama korigovanou distanci. Pokud jsou sekvence krátké a je mnoho OTU je lepší opět použít p distanci.

Swoford a kol. dodávají:

1.         Složitější modely je možné používat, jestliže je počet pozic větší než 2000.

2.         Výhodná je LogDet distance, protože má malý rozptyl, je hodně obecná a odolná vůči vlivu různého nukleotidového složení.

LogDet distance je citlivá na rozdíly v substituční rychlosti různých pozic a nedá se na tento vliv korigovat. Při  jejím použití je potřeba vyhodit z alignmentu alespoň všechny invariabilní pozice (PAUP- příkaz: Exclude constant). Jako vždy je při výběru vhodného modelu nejlepší vyzkoušet více možností a dívat se co se děje. Pokud dávají dva modely stejné výsledky je lepší použít ten jednodušší.

Pozor: Vše, co bylo řečeno, platí pro konstrukce stromů. Pokud chcete vypočítat délky větví u určitého stromu, použijte nejsložitější model, ten je nejpřesnější.

Mnoho dobrých fylogenetiků (Peter Foster) používá i u distančních metod složité modely. Velmi oblíbená je ML(maximum likelihood) distance. 

Maximální věrohodnost

V rámci maximální věrohodnosti je možné rozhodovat, jestli složitější model dává signifikantně lepší výsledek pomocí likelihood ratio testu (LRT).

δ=2(ln L1-lnL0 )

lnL1….věrohodnost stromu podle složitějšího modelu

lnL0….věrohodnost stromu podle jednoduššího modelu (nulová hypotéza)

Hodnota statistiky δ je vždy větší než 0. Její rozložení lze obecně určit simulací. Ve speciálním případě (pokud je jednodušší model obsažen ve složitějším modelu) má tato statistika zhruba rozložení χ2 se stupni volnosti odpovídajícími rozdílu v počtu volných parametrů mezi modely.

Modeltest 3.7 (freeware - http://darwin.uvigo.es/software/modeltest.html)

Tento test je inkorporován v programu Modeltest, kde se provádí hierarchicky od nejjednodušších modelů ke složitějším. Postup výběru modelu je následující:

•          Určit parametry pro všechny evoluční modely, které podporuje PAUP. Nejprve spusťte v PAUPu soubor s vašim alignmentem. Poté spusťte v tomtéž PAUPu soubor modelblockPAUPb10, který je uložen pod Modeltest3.7 folder/paupblock. Program určí pro každý model hodnoty parametrů, které maximalizují likelihood NJ stromu. Výstupem bude soubor .scores, ve kterém jsou tyto hodnoty uloženy.

•          Soubor .scores spusťte v programu Modeltest3.7.win.exe uloženém v adresáři Modeltest3.7 folder/bin. Můžete to provést následovně. Umístěte soubor .scores do stejného adresáře jako Modeltest3.7.win.exe. V Total Commanderu přetáhněte myší soubor Modeltest3.7.win.exe do příkazové řádky, udělejte mezeru za ni napište znaménko < a přetáhněte tam soubor .scrores a uzavřete znaménkem >. Udělejte mezeru a za ni napište vámi zvolený název výstupního souboru (např. nejmodel.txt). Příkazová řádka vypadá tedy asi takto:  

Modeltest3.7.win.exe <alignment.scores> nejmodel.txt

Po zpuštění programu klávesou ENTER program vytvoří soubor s vámi zvoleným jménem.

Výstupní soubor obsahuje hodnoty likelihoodu stromu pro všechny modely a hodnoty signifikance rozdílů mezi likelihoody modelů tak, jak byly hierarchicky srovnávány. Na konec program vybere nejvhodnější modely podle dvou kriterií: výsledku LRT (likelihood ratio testu) a AIC (Akaike information criterion). Pro vybraný model napíše příkazový řádek pro PAUP, ve kterém se nadefinuje PAUPu tento model. Příklad takového řádku:

Lset  Base=(0.4009 0.2060 0.1059)  Nst=2  TRatio=3.3424  Rates=gamma  Shape=0.4337  Pinvar=0;