PAMLNeutrality tests

Dr. Robert Kofler

Acknowledgments:

Many thanks go to Dr. Carolin Kosiol who generously shared her PAML course slides, which very hugely inspiring for thislecture.

Introduction

PAML is a widely used and cited tool written by Ziheng Yang

Book: Molecular Evolution - A Statistical Approach; Ziheng Yang, 2014, Oxford University Press

What can you do with PAML?

• Identify positive selectionion• Compare and test phylogenetic trees• Estimate species divergence• Reconstruct ancestral sequences

Identification of positive selection

This is of high interest for many biologists as fixations of advantagous mutations are responsible for phenotypic evolution,evolutionary innovations and species divergence.

Therefore it is necessary to distinguish positive selection from neutral evolution (mutation and drift)

Required input?

• Aligned sequences - mostly multiple sequence alignment• Phylogenetic trees for the aligned sequences

Equally important. What does PAML not do?

• align sequences• gene prediction; for example codeml assumes that sequences are pre-aligned exons, with first nucleotide being first

position of codon; introns are assumed to be already removed• generate phylogenetic trees

Obtaining PAML

http://abacus.gene.ucl.ac.uk/software/paml.html

1

Sub-programs

• baseml: maximum likelyhood analysis of nucleotide sequences• basemlg: implements continous gamma model of Yang (1993)• codeml: maximum likelyhood analysis of codons or amino acids• evolver: simulations of sequences• pamp: parsimony based analysis for estimating substitution rates of Yang and Kumar (1996)• yn00: estimates synonymous and non-synonymous substitution rates between two sequences according to Yang and

Nielsen (2000) → yn00• mcmctree:a bayesian approach (Markov chain monte carlo) for estimating phylogenetic trees• chi2: computes critical chi2 and p values for likelihood ratio tests.

Let’s go to bin directory and have a look# use the path that fits your system# I used paml 4.8cd /Users/robertkofler/programs/paml4.8/binls

## baseml## basemlg## chi2## codeml## evolver## infinitesites## mcmctree## pamp## yn00

Example data sets

PAML comes with a nice set of examples, that allow to explore features and input files of the software but also to showcommon sources of misinterpretation of the data. All examples are from scientific publications# use the path that fits your system# I used paml 4.8cd /Users/robertkofler/programs/paml4.8/examplesls

## 9s.trees## CladeModelCD## DatingSoftBound## HIVNSsites## MHC.Swanson2002MBE## MouseLemurs## README.txt## TipDate.FluH1## TipDate.HIV2## YN00abglobin.result## abglobin.aa## abglobin.nuc## abglobin.trees## dNdSGene1.nuc## horai.nuc## horai.trees## lysin## lysozyme## mtCDNA## mtCDNAape## mtprim9.nuc

I picked some:

2

• HIVNSites: data from HIV-1 to explore which sites are highly conserved and which under positive selection• lysin: sperm lysin from 25 abalone species to test hypothesis that amino acids at the surface are positively selected• MouseLemur: mtDNA of mouse lemurs to estimate divergence dates• TipDate: alignment of 33SIV/HIV-2 sequences

Basic concepts

Universal genetic code

3

Comparing two protein coding DNA sequences

• ds (KS): number of synonymous substitutions per synonymous site• dN (KA): number of non-synonyomous substitutions per non-synonymous site• ω = dN/ds• ω ≈ 1 neutral evolution• ω < 1 purifying selection• ω > 1 positive selection

ω based tests are quite powerful in population genetics for identifying positive selection. They are actually more powerfulthan most other tests since they rely on fewer assumptions.

Yang, Z.; Bielawski, J. P. (2000). “Statistical methods for detecting molecular adaptation”. Trends in ecology & evolution(Personal edition). 15 (12): 496–503

How to count syn. vs non-syn. sites?

4

Counting of synonymous and non-synonymous sites is done for each codon of sequence. For example given the codon TTTthan 1 substitution is synonymous and 8 are non-synonymous (ie. aa-changing)

• S = 1/9 ∗ 3 = 0.333• N = 8/9 ∗ 3 = 2.666

Given a nucleotide sequence of a protein (eg “TTT TTT TTT TTT”) the number of syn. or non-syn. sites is just the sum ofall codons (e.g. S = 4 ∗ 0.333 and N = 4 ∗ 2.66 =)

Calculating ω

Now, we assume an alignment of two protein sequences with a length of 100bp. Using the method shown before we countN = 75 and S = 25. Based on the alignment we observe 5 non-synonymous substitutions and 5 synonymous ones.

ds = 5/25 = 0.2 dn = 5/75 = 0.066 ω = 0.066/0.2 = 0.333

What does it mean when ω < 1?

5

Sources of bias and corrections

Transition, transversion rate difference

Transitions (66.6%) are much more frequent than transversions (33.3%). At the third position of the codon transitions aremore likely to be synonymous than transversions, which can lead to a bias of the synonymous length estimates.

6

7

Codon usage bias

Many codons actually code for the same amino acid, which is also called degeneracy of genetic code (or redundancy). Forexample Alanin (A), here all four codons GCU, GCC, GCA, and GCG would be equally suitable, and in an unbiased codeall should occur equally often in genes (25%). In practice this is not the case for some organism. For example E. coli, S.cerevisae, D. melanogaster, A. thaliana have a codon usage bias, whereas humans have not.

Investigation of genes showed that codon usage bias leads to a reduced number of synonymous sites (thus the opposite effectof the transition/transversion bias).

Multiple mutations in one codon

Frequently two or more substitutions are found in a codon. The problem is that the number of synonymous and non-synonymoussubstitutions depends on the order of the substitions, for example CCT and CAG

• Route 1: CCT(Pro) ↔ CAT (His) ↔ CAG (Gln): S = 0 and N = 2• Route 2: CCT(Pro) ↔ CCG (Pro) ↔ CAG (Gln): S = 1 and N = 1

So what to do? Usually the average over all possible routes is used: S = 0.5 and N = 1.5

Factors to consider for building a model

• transition/transversion rate ratio κ• Biased codo usage for codon j: πj

• nonsynonymous to synonymous rate ratio: ω

For example, the following factors need to be considered to model rates to CTG:

• CTC(Leu) ⇒ CTG (Leu): π• TTG(Leu) ⇒ CTG (Leu): κπ• GTG(Val) ⇒ CTG (Leu): ωπ• CCG(Pro) ⇒ CTG (Leu): κωπ

Based on these factors and given the data PAML computes the most likely values for these factors, ie. the parameters thatbest explain the data.

Branch models vs sites models

The problem is that for an entire protein ω is rarely larger than 1 for all sites along a branch. So frequently site modelsare used, where ω is allowed to vary among codons. Thus some sites could be positively selected while others are not.Alternatively ω may vary among branches, when for example positive selection is strong for a species exposed to a rapidlychanging environment. This are called branch models.

8

Figure 1:

9

10

Literature:

• branch model: Yang, Z. 1998. Likelihood ratio tests for detecting positive selection and application to primate lysozymeevolution. Molecular Biology and Evolution 15:568-573

• site model in HIV-1 Nielsen, R., and Z. Yang. 1998. Likelihood models for detecting positively selected amino acid sitesand applications to the HIV-1 envelope gene. Genetics 148:929-936.

Using PAML

HIVNSites; Sites-model with codeml

We use the example HIVNSites; Investigates HIV-1 (env region) during the infection (from year 3-7) of a single patient. Ittests the hypothesis that some regions could be positively selected during the infection.

Literature: Nielsen, R., and Z. Yang. 1998. Likelihood models for detecting positively selected amino acid sites and applicationsto the HIV-1 envelope gene. Genetics 148:929-936.cd /Users/robertkofler/programs/paml4.8/examples/HIVNSsitesls

## HIVenvSweden.trees## HIVenvSweden.txt## HIVenvSweden4s.trees## HIVenvSweden4s.txt## README.txt## codeml.ctl

Input of sequence data: the PHYLIP format

head /Users/robertkofler/programs/paml4.8/examples/HIVNSsites/HIVenvSweden.txt

## 13 273#### U68496 GTAGTAATTAGATCTGAAAACTTCTCGAACAATGCTAAAACCATAATAGTACAGCTAAATAAATCTGTAGAAATTAATTGTACAAGACCCAACAACAATACAAGAAGAAGTATACATTTTGGACCAGGGAAAGCATTTTATGCAGGAGAAATAATAGGAGATATAAGACAAGCATATTGTACCCTTAATGGAACAGAATGGAATAACACTTTAAAACAGGTAGCTGAAAAATTAAGAGAACAATTTATTAAAACAATAGTTTTTAATCAATCC## U68497 ATAGTAATTAGATCTGAAAACTTCTCGAACAATGCTAAAACCATAATAGTACAGCTAAATAAATCTGTAGAAATTAATTGTACAAGACCCAACAACAATACAAGAAGAAGTATACATTTCGGACCAGGGAAAGCATTTTATGCAGGAGAAATAATAGGAGATATAAGACAAGCATATTGTACTCTTAATGGAGCAGAATGGAATAACACTGTAAAACAGGTAGCTGCAAAATTAAGAGAACAATTTAATAAAACAATAATCTTTAATCAATCC## U68498 GTAGTAATTAGATCTGAAAACTTCACGAACAATGCTAAAACCATAATAGTACAGCTAAATAAATCTGTAGAAATTAATTGTGTAAGACCCGGCAACAATACAAGAAGAAGTATACATATAGGACCAGGGAGAGCATATTATACAGGAGAAGTAATAGGAGATATAAGACAAGCACATTGTAACCTTAGTAGAACAGACTGGAATAAAACTTTAAAACAGGTAGCTGAAAAATTAAGAGAACAATTTAATACAACAATAGTCTTTAATCAATCC## U68499 ATAGTAATTAGATCTGAAAACTTCACGAACAATGCTAAAACCATAATAGTACAGCTAAATAAATCTGTAGAAATTAAGTGTGAAAGACCCAACAACAATACAAGAAAAAGTGTACATATAGGACCAGGGAGAGCATATTATACAGGAGAAATAATAGGAGATATAAGACAAGCACATTGTAACCTTAGTGGAACAGAATGGAGGGAAACTTTAAAACAGGTAGCTGAAAAATTAAGAGAACAATTTAATAAAACAATAGTCTTTAATCAATCC## U68500 ATAGTAATTAGATCTGAAAACTTCACGAACAATGCTAAAACCATAATAGTACATCTAAATGAATCTGTAGAAATTATTTGTGAAAGACCCAACAACAATACAAGAAAAAGTGTACATATGGGACCAGGGAGAGCATATTACACAGGAGAAATAATAGGAGATATAAGACAAGCACATTGTAACATTAGTAGAACAAATTGGACGGAAACTTTAAAACAGGTAGCTGAAAAATTAAGAGAACAATTTAATAAAACAATAGTCTTTAATCAATCC## U68501 GTAGTAATTAGATCTGAAAACTTCACGAACAATGCTAAGACCATAATAGTACAGCTAAATAAACCTGTAAAAATTAATTGTACAAGACCCAACAACAATGCAAAAATAAGAATACATATAGGACCAGGGAGACCATTTTATACAGCAGGAGAAATAGGAAATATAAGACAAGCACATTGTAACCTTAGTAGAACAGACTGGAATAACACTTTAAAACTGGTAGCTGAAAAATTAAGAGAACAATTTAATAAAACAATAGTCTTTAATCAATCC## U68502 GTAGTAATTAGATCTGAAAACTTCACGAACAATGCTAAGACCATAATAGTACAGCTAAATAACTCTGTAGCAATTAAGTGTGAAAGACCCAACAACAATACAAGAAAAAGTATACCTATAGGACCAGGGAGAGCCTTTTATACAACAGGAGACATAGGAGATATAAGACAAGCACATTGTAACCTTAGTAGAAAAGACTGGAATGACACTTTAAGACAGGTAGTTGGAAAGTTAAGAGAACAATTTGGAAGAACAATAATCTTTAATCAATCC## U68503 ATAGTAATTAGATCTGAAAACTTCACGAACAATGCTAAAACCATAATAGTACAGCTAAAGGAACCTGTAGACATTACTTGTGAAAGACCCAGCAACAATACAAGAAAAAGTATACATATAGGACCAGGAAAAGCATTTTATGCAACAGGAGAAATAGGAGATATAAGACGAGCACATTGTAACCTTAATAGAACAGCATGGAATAAAACTTTAAAACAGGTAGTTGAAAAATTAAGAGAACAATTTAAGAAAACAATAACCTTTAACCAATCC

• line 1: number of species followed by the sequence lenght• next lines: species-name (sequence name) and the sequence (separator either newline or two space)

Input of the tree:

head /Users/robertkofler/programs/paml4.8/examples/HIVNSsites/HIVenvSweden.trees

## 13 1#### (U68496: 0.023746, U68497: 0.079178, ((U68498: 0.050407, (U68499: 0.028429, U68500: 0.099890): 0.083460): 0.023647, (U68501: 0.131681, (U68502: 0.191150, (U68503: 0.193336, ((((U68504: 0.018984, U68506: 0.070370): 0.055603, U68505: 0.030301): 0.036688, U68507: 0.061120): 0.034773, U68508: 0.198346): 0.050591): 0.064013): 0.033695): 0.048101): 0.152813);###### // end of file. The rest are notes.########

11

• line 1: species and trees• line 2: the tree similar to the Newick format

the control file:

All parameters are set with the control file (*.ctl)cat /Users/robertkofler/programs/paml4.8/examples/HIVNSsites/codeml.ctl

## seqfile = HIVenvSweden.txt * sequence data file name## treefile = HIVenvSweden.trees * tree structure file name#### outfile = mlc * main result file name## noisy = 3 * 0,1,2,3,9: how much rubbish on the screen## verbose = 0 * 1: detailed output, 0: concise output## runmode = 0 * 0: user tree; 1: semi-automatic; 2: automatic## * 3: StepwiseAddition; (4,5):PerturbationNNI; -2: pairwise#### seqtype = 1 * 1:codons; 2:AAs; 3:codons-->AAs## CodonFreq = 2 * 0:1/61 each, 1:F1X4, 2:F3X4, 3:codon table## clock = 0 * 0: no clock, unrooted tree, 1: clock, rooted tree## aaDist = 0 * 0:equal, +:geometric; -:linear, {1-5:G1974,Miyata,c,p,v}## model = 0#### NSsites = 0 1 2## * 0:one w; 1:NearlyNeutral; 2:PositiveSelection; 3:discrete;## * 4:freqs; 5:gamma;6:2gamma;7:beta;8:beta&w;9:beta&gamma;10:3normal## icode = 0 * 0:standard genetic code; 1:mammalian mt; 2-10:see below## Mgene = 0 * 0:rates, 1:separate; 2:pi, 3:kappa, 4:all#### fix_kappa = 0 * 1: kappa fixed, 0: kappa to be estimated## kappa = .3 * initial or fixed kappa## fix_omega = 0 * 1: omega or omega_1 fixed, 0: estimate## omega = 1.3 * initial or fixed omega, for codons or codon-based AAs## ncatG = 10 * # of categories in the dG or AdG models of rates#### getSE = 0 * 0: don't want them, 1: want S.E.s of estimates## RateAncestor = 0 * (0,1,2): rates (alpha>0) or ancestral states (1 or 2)#### Small_Diff = .45e-6## cleandata = 1 * remove sites with ambiguity data (1:yes, 0:no)?## fix_blength = 0 * 0: ignore, -1: random, 1: initial, 2: fixed

Most important: * model - set the branch models * Nssites - sites models (branch models must be set to model=0); severalsites models can be specified at once, separated by space

12

Figure 2:

13

Most widely approaches for detecting positive selection are to compare models M1a with M2a or M7 with M8.

Run codeml

First copy the whole folder to a new location than run codeml. It will automatically use the control (ctl) file.cp -rf /Users/robertkofler/programs/paml4.8/examples/HIVNSsites /Users/robertkofler/tmpcd /Users/robertkofler/tmp/HIVNSsites/Users/robertkofler/programs/paml4.8/bin/codeml

Now lets look at the content of the folder. Obviously there are several novel files.cd /Users/robertkofler/tmp/HIVNSsitesls -l

## total 272## -rw-r--r-- 1 robertkofler staff 917 Jan 18 10:43 2NG.dN## -rw-r--r-- 1 robertkofler staff 917 Jan 18 10:43 2NG.dS## -rw-r--r-- 1 robertkofler staff 917 Jan 18 10:43 2NG.t## -rw-r--r--@ 1 robertkofler staff 546 Jan 18 10:40 HIVenvSweden.trees## -rw-r--r--@ 1 robertkofler staff 3881 Jan 18 10:40 HIVenvSweden.txt## -rw-r--r--@ 1 robertkofler staff 46 Jan 18 10:40 HIVenvSweden4s.trees## -rw-r--r--@ 1 robertkofler staff 1149 Jan 18 10:40 HIVenvSweden4s.txt## -rw-r--r--@ 1 robertkofler staff 1006 Jan 18 10:40 README.txt## -rw-r--r--@ 1 robertkofler staff 1734 Jan 18 10:40 codeml.ctl## -rw-r--r-- 1 robertkofler staff 13219 Jan 18 10:43 lnf## -rw-r--r-- 1 robertkofler staff 27840 Jan 18 10:43 mlc## -rw-r--r-- 1 robertkofler staff 15103 Jan 18 10:43 rst## -rw-r--r-- 1 robertkofler staff 535 Jan 18 10:43 rst1## -rw-r--r-- 1 robertkofler staff 36370 Jan 18 10:43 rub

The output is the mlc fileless /Users/robertkofler/tmp/HIVNSsites/mlc

Excercise: interpretation of output

We have the following resources available * Manual: http://abacus.gene.ucl.ac.uk/software/pamlDOC.pdf * FAQ: http://abacus.gene.ucl.ac.uk/software/pamlFAQs.pdf * First publication: http://www.genetics.org/content/148/3/929.long *Second publication: http://www.genetics.org/content/155/1/431.long

Answer the following questions

• How many sequences have been used and what is the number of codons per sequence?• Which codons have been positively selected?• Whats the ω of the positively selected codons?• How would you interpret these data, thinking about the tree and not just two sequences?• Where do we see that we did not accidentially use a branch model?• Do we have a codon usage bias?• What is κ in the three used models?• What is ω in the three used models?• What is the number of non-synonymous sites? What is the number of synonymous sites?• Why is N and S different between the models. Which factors could explain that?• Where to find the log likelihood of a model (lnL)? And where the number of parameters (np)?• We have three models which one is now the best? (a low log likelihood or a high?)

Some additional explanations

NEB vs BEB

• Naive Empirical Bayes (NEB) ignores sampling errors in parameter estimates.

14

• Bayes Empirical Bayes (BEB) accounts for sampling errors by integrating over a prior.

Model selection; test if one model is better than another one

Needs the degrees of freedom (to take the number of estimated parameters into account; in principle always a model withfewer estimated parameters is preferable) and the log likelihood (a higher log likelihood is better),

Model 1 (no selection): np = 26 and lnL = −1114.64 Model 2 (yes selection): np = 28 and lnL = −1106.44

Model 2 estimates more parameters but also has a higher log likelihood. So is the difference significant?

Now we can compute df = 28− 26 = 2 2δL = 2(−1106.44−−1114.64) = 16.4

Comparing this to a chi-squared distribution/Users/robertkofler/programs/paml4.8/bin/chi2

We can can conclude that model 2 (selection) is significnatly better suited to explain the data than model 1 (no selection)

Reminder

Remember, in practice mostly M1a (1) is compared to M2a (2) and M7 to M8.

• M7-M8: Yang, Nielsen, Goldman, Pedersen (2000 Genetics 155:431-449• M1a-M2a: Yang, Wong, Nielsen (2005. Mol Biol Evol).

15