to Introduction Biosystematics - Zool 575
Introduction to Biosystematics Lecture 20 - Model choice & corrected data
Outline
1. Introduction, data correction, simple models 2. Complex models, Among Site Rate Variation 3. Model Choice - Akaike Information Criterion
Derek S. Sikes University of Calgary Zool 575
Four steps - each should be explained in methods
1. Character (data) selection (not too fast, not too slow) Why did you choose these data? 2. Alignment of Data (hypotheses of primary homology) How did you align your data? 3. Analysis selection (choose the best model / method(s)) - data exploration Why did you chose your analysis method? 4. Conduct analysis - strategies etc.
Lack of a reproducible method resulted in three major approaches: (Explosion in 1960s) 1. Phenetics - similarity / distances only, not evolution, not phylogeny 2. Cladistics - phylogeny inferred using characters & parsimony 3. Statistical Phylogenetics - phylogeny inferred using corrected data & best fitting model [both distance & character data]
Lack of a reproducible method resulted in three major approaches: (Explosion in 1960s) 1. Phenetics - similarity / distances only, not evolution, not phylogeny [uncorrected distance data] 2. Cladistics - phylogeny inferred using characters & parsimony [+/- uncorrected character data] 3. Statistical Phylogenetics - phylogeny inferred using corrected data & best fitting model [both distance & character data]
Saturation
Saturation graph as time proceeds DNA distances also increase, to a point of saturation
observed DNA distances p-distance Observable change increases in a linear fashion (x ~ y) for a while Only so much change is observable time Real change continues with time
1
Introduction to Biosystematics - Zool 575
Saturation
Saturation graph as time proceeds DNA distances also increase, to a point of saturation
DNA Distances actual observed Observable change increases in a linear fashion (x ~ y) for a while Only so much change is observable time Real change continues with time
Selection of Molecular characters
Why?
- constraints: for a given gene, some sites essentially do not change (preventing DNA distances from reaching 100%) - even for regions that are variable, typically DNA distances cant go beyond 75% since 1/4 of the changes would be to the same nucleotide - other sites do change: for a given comparison of 2 taxa a variable site might have changed: once: (good) two or more times: (bad) - multiple hits information lost
Selection of Molecular characters
example Species1 Species2 ATGCCTGGACTTATAA ATGCCGGGAGATATAA . ..
Selection of Molecular characters
Example - recent divergence, no saturation Ancestral sequence ATGCCTGGACTTATAA ATGCCTGGACATATAA ATGCCTGGAGATATAA ATGCCGGGAGATATAA . .. ATGCCTGGACTTATAA 1 ATGCCTGGACTTATAA 1 ATGCCTGGACTTATAA 1 . ..
3 changes observed - this is the minimum and is only the ACTUAL number of changes if there have been no multiple hits (or back mutations) i.e. each site changed only once since speciation / divergence
3 changes observed, 3 actual changes
Selection of Molecular characters
Example - ancient divergence, with saturation Ancestral sequence TTGCGTGGACTTATAA TTGCGTGGACATATAT ATGCCTGGAGTTATAA ATGCCGGGAGATATAA . .. ATGCGTGGACTTATTA 4 ATGCCTGGACTTAAAA 7 ATGCCTGGACTTATAA 4 . ..
Saturation
Transitions & Transversions Recall: A, G = purines T, C = pyrimidines
most mutations are transitions [Ts] (replacement of a purine with a purine, or a pyrimidine with a pyrimidine) Transversions [Tv] (purine replaced by a pyrimidine) are less common
3 changes observed, 15 actual changes (2 parallelisms (sites 1 & 5) = homoplasy)
2
Introduction to Biosystematics - Zool 575
Saturation
Another option to deal with saturation Convert DNA data to purines & pyrimidines in matrix use R for purines, Y for pyrimidines Removes all the transitions from the data & leaves only transversions Less extreme than using only amino acids (which may ignore a lot of signal in the DNA data)
Saturation
This is essentially an extreme form of weighting (a method of correcting the data) transitions are weighted = 0 and transversions are weighted = 1 some opt for less extreme weighting: e.g. transitions = 1 transversions = 2 The logic is the less common evolutionary change should be less homoplasious & have more phylogenetic information (signal)
Parsimony
- Variants of parsimony: Fitch (unordered 0->1, 0->2 ) Wagner (ordered states 0->1->2 etc) Dollo (each state originates only once) Transversion Parsimony (recode data into purine & pyrimidines) Generalized parsimony
Weighting
Weighted Parsimony Generalized parsimony - Can be Fitch, unweighted, unordered or.. - Can weight certain characters or weight certain transformations between states eg can weight tranvsersions twice transitions - Computationally costly - No objective method to determine weights
Weighting
Weighted Parsimony - Generalized parsimony - a priori weighting has been criticized as being too subjective - but equal weighting is still weighting eg a dataset with evidence of massive transition bias equal weighting (1:1) is a stronger and less realistic assumption than unequal weighting (2:1)
Weighting
Some (mostly cladists) consider weighting (or correcting) of data in any way to be wrong They argue that the subjectively decided weighting scheme may alter the results In which case the results are not due to the signal in the data but due to the weighting scheme itself The key is whether the weighting scheme is justified
3
Introduction to Biosystematics - Zool 575
Models of Sequence Evolution
- Given that observed distances may underestimate actual distances - Methods were developed to correct the distances / data - to estimate the actual distances - Great variety of models of sequence evolution: - simple (unrealistic, restrictive) - complex (realistic, unrestrictive)
Models of Sequence Evolution
Some features many models deal with: 1. The tendency of one state to change to another (rate) 2. The frequency of the different states (base composition) 3. Variation in rates of change among characters (ASRV) - some sites change a lot, others change little or not at all
Models of Sequence Evolution
The tendency of one state to change to another - Described by a rate matrix - The most complex in PAUP* is the GTR model, (general time reversible) - The simplest is JC69 - Various intermediates: HKY, F81, K2P
General time reversible model (GTR)
Different probabilities for each of 6 substitution types (lset nst = 6) 5 free rate parameters (6th is set to make the total = 1.0)
Jukes and Cantor (1969) Model
Models of Sequence Evolution
Base Composition - can be (in PAUP*): - equal (eg JC69) A=C=T=G - specified (based on what?) - prior analysis - empirical (equal to observed - this is a good estimate of the values estimated by maximum likelihood & much faster) - estimated by ML for DNA there are 4 bases so this involves 3 free parameters
Same probability for each of 6 substitution types (lset nst = 1) no free rate parameters
4
Introduction to Biosystematics - Zool 575
Models of Sequence Evolution
Among Site Rate Variation (ASRV)
- if present in data and not accounted for analysis can be mislead - inconsistent, erroneous - can be accommodated by adding extra parameters to any model - pInvar - proportion of invariable sites [ eg JC+I] - gamma - distribution of ASRV controlled by one parameter, alpha [ eg GTR+I+G]
The shape of the gamma distribution for different values of alpha. If alpha = infinity, rate is equal for all sites (no ASRV)
Low alpha (eg 0.5 and lower) = high ASRV High alpha (eg > 2) = low ASRV
Models of Sequence Evolution
How are these parameters chosen?
- They can be +/- arbitrarily set by the user (not recommended!) - They can be estimated from a near optimal tree and set to help find the optimal tree (searches go faster if these parameters are not being estimated during the search) - They can be estimated from the data during the search (values that maximize the likelihood are chosen) - Estimation from the data is superior to methods for Weighting of Parsimony
Models of Sequence Evolution
You can use models that: Deal with different transition/transversion ratios Deal with unequal base composition Deal with heterogeneity of rates across sites Deal with heterogeneity of the substitution process (different rates across lineages, different rates at different parts of the tree) Increased realism & complexity requires more parameters The more free parameters, the better your model fits your data (good) The more free parameters, the higher the variance of the estimate (bad) Should, therefore, fit a model to your data (more later)
Models of Sequence Evolution
Some commonly used models Jukes and Cantor (JC69): All base compositions equal (0.25 each), rate of change from one base to another is the same Kimura 2-Parameter (K2P): All base compositions equal (0.25 each), different substitution rate for transitions and transversions). Hasegawa-Kishino-Yano (HKY): Like the K2P, but with base composition free to vary. General Time Reversible (GTR): Base composition free to vary, all possible substitutions can differ. All these models can be extended to accommodate invariable sites and site-to-site rate variation.
Models of Sequence Evolution
- All assume: 1) rate of change at a site is independent of prior states (related to # 2) (Markov property) 2) probability of change is constant over time (homogeneity) 3) base composition (frequency of A, C, G, T) is in equilibrium (stationarity), ie it is constant over time 4) rate of change symmetrical is (eg A to T is same as T to A) (time reversible)
5
Introduction to Biosystematics - Zool 575
Models of Sequence Evolution
5) a shared set of branch lengths (a common mechanism) - if branch lengths are going to be estimated with all the data we must assume all the data share (ie evolved at) the same historical rate [coins] - eg if a branch is short ( = low change) for some characters, it is short for all the characters - Parsimony does not make this assumption branch lengths are ignored (no common mechanism) or we might say Parsimony assumes characters do not share a common set of branch lengths
Models of Sequence Evolution
Models can be used to correct both distance data and character data - Distance data can be corrected using models via a method of moments estimator (simple but less precise) or a maximum likelihood estimator (complex but more precise) - Character data are corrected via maximum likelihood estimation (or Bayesian)
Models of Sequence Evolution
For distance data: - distances are the branch lengths - correction to observed (D) to yield estimate of actual (d) - simplest correction: Jukes and Cantor 1969 (JC69) d = -3/4ln (1- 4/3 D) eg D = 0.119, d = 0.130 eg D = 0.00482, d = 0.00483
Models of Sequence Evolution
Saturation graph
DNA Distances actual observed
Small distances need little correction Larger distances need more
time
Models of Sequence Evolution
Models of Sequence Evolution
The Jukes and Cantor 1969 model is the simplest but also the most restrictive and least realistic
All substitutions have an equal probability and base frequencies are equal fewest parameters, most assumptions
6
Introduction to Biosystematics - Zool 575
Models of Sequence Evolution
The assumptions of JC69 are wrong 1. Assumes equal base frequencies - but most taxa have nucleotide biases - eg. Insect mtDNA heavy AT bias (80%+) 2. Assumes rates of change among bases are equal - but we know transitions far exceed transversions due to proofreading enzymes
Models of Sequence Evolution
The assumptions of JC69 are wrong 3. Assumes all sites (characters) change with same probability (no among-site rate variation, ASRV) - structural & functional aspects of molecule constrain change at certain sites eg 3rd codon positions much more variable than 1st and 2nd positions all other models improve on these restrictive & unrealistic assumptions by relaxing them
Models of Sequence Evolution
The assumptions of JC69 are wrong
Models of Sequence Evolution
More complex models usually fit the data better
The expectations of JC69 do not fit the reality
Models of Sequence Evolution
- Models specify many aspects of the evolutionary process e.g. - the probability of change (rates) - the composition of states - Models can be used to generate simulated data (create artificial data) - More typically they are used to estimate the likelihood of observed data - I think (hope) it will be easier to understand models if we start with simulation of data
Simulation of data using a model
- Models can specify anything - from crazy (p(A)=1.0; p(C )=0.0; p(T)=0.0; p(G)=0.0 are A) but no one would do this!
all states
- to more realistic (p(A)=0.25; p(C )=0.25; p(T)=0.25; p(G)=0.25 all states are equally common ) like JC69 - to really realistic (p(A)=0.35; p(C )=0.12; p(T)=0.45; p(G)=0.08 )
7
Introduction to Biosystematics - Zool 575
Simulation of data using a model
1. Start with a tree - you fix the topology on which the program will evolve the data. 2. A starting sequence is generated (or given) for the ancestor AGTCAGTC
Simulation of data using a model
Crazy asymmetrical model: rate parameters specify all mutations change to A, but start with JC69 base frequencies Use program Seq-Gen to simulate data with this model AGTCAGAC AGTCAGTC AGTAAATA AGTCAATA AAAAAATA AATAAGAC AAAAAAAC
AATCAAAA AGTCAAAA AAAAAAAA Eventually the sequences would be all As
3. The program then mutates (evolves) the ancestral sequence up the branches to generate the final observed data
Simulation of data using a model
JC69 model: rate parameters specify all mutations equiprobable; start with equal base frequencies (note: stochasticity will generate small changes from equality - more obvious in smaller sequences) AGTCAGGC AGTCAGTC AGTGACTT AGTCACTT AGTCACGT AGCTACGT Final sequences with similar frequencies of 4 bases ATCGACTA ACTAAGGC ACGAATGC
Choosing a Model
Model choice can be the single most important step in the analysis (after the alignment) Using the wrong model (model misspecification) for the analysis can yield incorrect results = another source of phylogenetic error Perform a model - fitting exercise to find a model that captures the complexity of the data but does not overwhelm the datas ability to estimate its parameters
AATCACGT
Choosing a Model
Increased realism & complexity requires more parameters The more free parameters, the better your model fits your data (good) The more free parameters, the higher the variance of the estimate (bad) However, being too complex is generally preferred over being too simple (esp with Bayesian analyses) also, our most complex models are still often much simpler than the data
NCM - process
NCM - pattern always shortest tree
CM - process & pattern Loose analogy of Over-fitting Not always shortest tree
8
Introduction to Biosystematics - Zool 575
Phylogenetic inference: Model selection
To determine which model fits the data best - to avoid over - and under - fitting: the Akaike Information Criterion Parsimony (MP) vs. Mkv+!
Phylogenetic inference: Model selection
Testing models of evolution - Modeltest Version 3.06 (c) Copyright, 1998-2000 David Posada ([e-mail]) Department of Zoology, Brigham Young University WIDB 574, Provo, UT 84602, USA _______________________________________________________________ Tue Nov 18 11:12:00 2003 Input file: model.scores Input format: Paup matrix file
AIC = -2(-lnL) + 2(k) -lnL for MP = 146.95 -lnL for Mkv!= 322.2 parameters = 1767 parameters = 59 AIC=3828 AIC=762
** Log Likelihood scores ** JC F81 K80 HKY TrNef TrN K81 K81uf TIMef TIM TVMef TVM SYM GTR = = = = = = = = = = = = = = 3260.3015 3141.3474 3173.9417 2995.3235 3153.1194 2986.3464 3137.0684 2987.9604 3116.0508 2977.2561 3105.5444 2978.9819 3084.5220 2966.5007 +I 3260.3015 3141.3474 3173.9417 2995.3235 3153.1194 2986.3464 3137.0684 2987.9604 3116.0508 2977.2561 3105.5444 2978.9819 3084.5220 2966.5007 +G 3268.2341 3147.3276 3179.7322 2995.5459 3157.2998 2987.1013 3143.3892 2991.4077 3121.4861 2979.6755 3108.6196 2984.5586 3087.6189 2969.7056 +I+G 3260.2649 3140.7678 3173.3606 2988.6401 3151.5359 2979.6970 3136.2917 2983.5828 3114.4956 2971.7051 3103.8760 2976.3069 3082.5017 2961.8000
Lower AIC is better - models are rewarded for their fit to the data but penalized for having superfluous parameters
Tuffley & Steel (1997) -lnL of MP = (nsites+nsteps)*ln1/4 MP parameters = nsites x nbranches
Phylogenetic inference: Model selection
** Akaike Information Criterion (AIC) ** Model selected: GTR+I+G -lnL = 2961.8000 AIC = 5943.6001 Base frequencies: freqA = 0.3485 freqC = 0.1311 freqG = 0.1268 freqT = 0.3936 Substitution model: Rate matrix R(a) [A-C] = 16.9162 R(b) [A-G] = 72.5632 R(c) [A-T] = 29.8403 R(d) [C-G] = 0.0000 R(e) [C-T] = 280.2567 R(f) [G-T] = 1.0000 Among-site rate variation Proportion of invariable sites (I) = Variable sites (G) Gamma distribution shape parameter =
Phylogenetic inference: Model selection
Alternative to AIC - Bayes Factors - needed for proper Bayesian model selection (AIC requires maximum likelihoods but Bayesian analyses yield only marginalized likelihoods) - Simple to compute: to compare model 0 (M0 ) to model 1 (M1) the Bayes factor in favor of model 1 over 0 is the ratio of the model likelihoods:
B 10 = harmonic mean of marginalized log-likelihood of model 1 harmonic mean of marginalized log-likelihood of model 0
0.6435 1.9963
Recall that dividing logs is the same as subtracting so
Phylogenetic inference: Model selection
Alternative to AIC - Bayes Factors - model 0 score: -31,634 - model 1 score: -30,962 - B10 = |-31,634| - |-30,962| = 672 - compare to interpretation chart for Bayes factors (from Kass & Raftery,
1995)
Correction doesnt always work
- Random error +/- gone with huge dataset of 124,026 characters - Systematic error evident in ME analysis (tree on right) (even with corrected distance data!) (& loss of data) - 100% branch support values indicate no random error
1-3 3-20 20-150 150+
= not worth more than a bare mention = positive = strong = very strong
Kass, R.E. and Raftery, A.E. (1995). Bayes factors. Journal of the American Statistical Association 90: 773-795.
9
Introduction to Biosystematics - Zool 575
Terms - from lecture & readings
corrected description of statistical phylogenetics Fitch Parsimony Wagner Parsimony Dollo Parsimony Transversion Parsimony Generalized Parsimony Weighting equal weighting rate matrix JC69 K2P HKY GTR base composition Among Site Rate Variation pInvar (I) gamma (G) alpha homogeneous Markov models stationarity time reversible common mechanism no common mechanism method of moments estimator maximum likelihood estimator model misspecification model fitting modeltest
Study questions
What are the three main components to a model of sequence evolution? What does a low alpha value (eg 0.33) of the shape parameter of the gamma distribution mean about your data? What are the three assumptions of the JC69 model that are wrong and why are they wrong? What is good and bad about increasing the number of parameters in a model? Why is equal weighting a form of subjective a priori weighting?
10