An Improved Codon Modeling Approach for Accurate Estimation of the Mutation Bias

0.85 MB·English

by  Thibault Latrille, Nicolas Lartillot

An Improved Codon Modeling Approach for Accurate
Estimation of the Mutation Bias
Thibault Latrille *,1,2 and Nicolas Lartillot 1
1CNRS,LaboratoiredeBiom(cid:2)etrieetBiologie(cid:2)EvolutiveUMR,Universit(cid:2)e deLyon,Universit(cid:2)e Lyon1,Villeurbanne,France
2(cid:2)EcoleNormaleSup(cid:2)erieuredeLyon,Universit(cid:2)edeLyon,Universit(cid:2)eLyon1,Lyon,France
*Correspondingauthor:E-mail:[email protected].
Associateeditor:TalPupko
Abstract
Phylogenetic codon models are routinely used to characterize selective regimes in coding sequences. Their parametric D
o
design, however, is still a matter of debate, in particular concerning the question of how to account for differing w
n
nucleotide frequencies and substitution rates. This problem relates to the fact that nucleotide composition in lo
a
protein-coding sequences is the result of the interactions between mutation and selection. In particular, because of de
d
thestructureofthegeneticcode,thenucleotidecompositiondiffersbetweenthethreecodingpositions,withthethird fro
positionshowingamoreextremecomposition.Yet,phylogeneticcodonmodelsdonotcorrectlycapturethisphenom- m
h
enonandinsteadpredictthatthenucleotidecompositionshouldbethesameforallthreepositions.Alternatively,some ttp
s
models allow for different nucleotide rates at the three positions, an approach conflating the effects of mutation and ://a
selection on nucleotide composition. In practice, it results in inaccurate estimation of the strength of selection. c
a
d
Conceptually, the problem comes from the fact that phylogenetic codon models do not correctly capture the fixation e
m
biasactingagainstthemutationalpressureatthemutation–selectionequilibrium.Toaddressthisproblemandtomore ic
.o
accuratelyidentifymutationratesandselectionstrength,wepresentanimprovedcodonmodelingapproachwherethe u
p
fixationrateisnotseenasascalar,butasatensor.Thisapproachgivesanaccuraterepresentationofhowmutationand .co
m
selectionopposeeachotheratequilibriumandyieldsareliableestimateofthemutationalprocess,whiledisentangling /m
themeanfixationprobabilitiesprevailingindifferentmutationaldirections. b
e
/a
Keywords:codonmodels,phylogenetics,nucleotidebias,mutation–selectionmodels. rtic
le
/3
9
/2
Introduction by invoking a parameter x, acting multiplicatively on non- /ms
a
synonymous substitutions rates (Goldman and Yang 1994; c
0
Phylogeneticcodonmodelsarenowroutinelyusedinmany 0
MuseandGaut1994). Usingaparametric modelautomati- 5
domainsofbioinformaticsandmolecularevolutionarystud- /6
callycorrectsforthemultiplicityissuescreatedbythecom- A5
0
ies. Oneoftheir main applications has been to characterize 3
plexstructureofthegenetic codeandby unevenmutation r5
the genes, sites (Nielsen and Yang 1998; Yang et al. 2005; ratesbetweennucleotides.Asaresult,xcapturesthenet,or ti05 b
Murrell et al. 2012), or lineages (Zhang et al. 2005; cy
aggregate,effectofselectiononnonsynonymousmutations, l g
Kosakovsky Pond et al. 2011) having experienced positive alsocalledd =d (DosReis2015;SpielmanandWilke2015). eue
N S s
selection(Murrelletal.2015;Enardetal.2016).Moregener- In reality, the selective effects associated with nonsynon- t o
n
ally, these models highlight the respective contributions of ymous mutations depends on the context (site-specificity) 1
5
mutation,selection,geneticdrift(Teufeletal.2018),andbi- andtheaminoacidsinvolvedinthetransition(Kosioletal. F
e
ased gene conversion (Kosiol and Anisimova 2019; Pouyet 2007). Attempts at an explicit modeling of these complex bru
andGilbert2021),andthecausesoftheirvariationbetween selectivelandscapeshavealsobeendone,leadingtomecha- ary
genes (Zhang and Yang 2015) or across species (Seo et al. nistic codon models, based on the mutation–selection for- 20
2
2004;Popadinetal.2007;LartillotandPoujol2011). malism (Halpern and Bruno 1998). These models, further 2
Conceptually, codon models take advantage of the fact developedinmultipleinferenceframeworks(Rodrigueetal.
thatsynonymousandnonsynonymoussubstitutionsaredif- 2010; Tamuri et al. 2012), sometimes using empirically in-
ferentially impacted by selection. Assuming synonymous formed fitness landscapes (Bloom 2014), could have many
mutations are neutral, the synonymous substitution rate is interesting applications, such as inferring the distribution of
equal to the underlying mutation rate (Kimura 1983). fitness effects (Tamuriet al. 2012)or detecting genes under
Nonsynonymous substitutions, on the other hand, reflect adaptation(RodrigueandLartillot2017;Rodrigueetal.2021),
the combined effect of mutation and selection (Ohta or even phylogenetic inference (Ren et al. 2005). However,
1995). Phenomenological codon models formalize this idea they are computationally complex and potentially sensitive
(cid:2)TheAuthor(s)2022.PublishedbyOxfordUniversityPressonbehalfoftheSocietyforMolecularBiologyandEvolution.
ThisisanOpenAccessarticledistributedunderthetermsoftheCreativeCommonsAttributionLicense(https://creativecommons.
org/licenses/by/4.0/),whichpermitsunrestrictedreuse,distribution,andreproductioninanymedium,providedtheoriginalworkis Open Access
properlycited.
Mol.Biol.Evol.39(2):msac005 doi:10.1093/molbev/msac005 AdvanceAccesspublicationJanuary11,2022 1
MBE
LatrilleandLartillot . doi:10.1093/molbev/msac005
totheviolationoftheirassumptionsaboutthefitnessland- confoundthetestsfordetectingpositiveselectionand,more
scape (such as site independence). For these reasons, phe- generally,theestimationofx(Galtieretal.2009;Ratnakumar
nomenological codon models remain an attractive, etal.2010;Lartillotetal.2013;Figuetetal.2014;Bol(cid:2)ıvaretal.
potentiallymorerobust,althoughstillperfectibleapproach. 2019).EvenintheabsenceofgBGC,however,unevenmuta-
Theparametricdesignofphenomenologicalcodonmod- tionrates varyingacross speciescanhaveanimportant im-
els,relyingonasingleaggregateparameterx(orsite-specific pactontheestimationofthestrengthofselection(Gu(cid:2)eguen
x),raisesthequestionwhethertheyaccuratelyestimatethe
andDuret2018).Allthissuggeststhat,evenbeforeintroduc-
underlyingselectiveandmutationalprocess.First,simulations
inggBGCincodonmodels,correctlyformalizingtheinterplay
underamutation–selectionformalismhaveshownthatthe
between mutation and selection in current codon models
strengthofselectionisestimatedreliablybyphenomenolog-
would be an important first step, which is the focus of this
icalcodonmodels(SpielmanandWilke2015).Morespecif-
manuscript.
ically, the model originally proposed by Muse and Gaut
Inthisdirection,thekeypointthatneedstobecorrectly
(1994), hereafter called MG, gives an accurate estimate of D
formalized is the following. If the nucleotide’s realized fre- o
the underlying x. However, several observations suggest w
quenciesaretheresult ofacompromisebetweenmutation n
thatthemutationalprocessisnotaccuratelyestimated.For lo
andselection,thenthisimpliesthatthestrengthofselection ad
instance, in their simplest form (Goldman and Yang 1994; e
MuseandGaut1994),codonmodelspredictthatthenucle- iFsonrointstthaencsea,mifethbeetmweuetantiaolnl npurcolceeostsidiesAorT-abmiaisneod,atchidenp,abires-. d from
otidecompositionshouldbethesameforallthreepositions
of the codons, and should be equal to the nucleotide equi- causeofselection,therealizednucleotidefrequenciesatequi- http
libriumfrequenciesimpliedbytheunderlyingnucleotidesub- libriumwillbelessAT-biasedthanexpectedunderthepure s://a
stitution rate matrix. In reality, the nucleotide composition mutationprocess.However,thisimpliesthat,atequilibrium, ca
differs:thethirdpositionshowsmoreextremeGCcomposi- therewillbeanetmutationpressuretowardAT,whichhasto de
m
tion,reflectingtheunderlyingmutationbias,comparedwith becompensatedforbyanetselectiondifferentialtowardGC. ic
.o
the first and second positions, which are typically closer to In order for a codon model to correctly formalize this u
p
50%GC(SingerandHickey2000). subtleinterplaybetweenmutationandselection,theparam- .co
m
Thesemodulationsacrossthethreecodingpositionshave eterresponsibleforabsorbingtheneteffectofselection(i.e., /m
been accommodated using the so-called 3(cid:2)4 formalism x) should not be a scalar, but an array of x values (i.e., a be
/a
(Goldman and Yang 1994; Pond and Muse 2005), allowing tensor) unfolding along multiple directions. In the present rtic
for different nucleotide rate matrices at the three coding work, we address the question of whether we can derive a le
/3
positions. However, this is also problematic. For instance, it model which is able to correctly tease apart mutation rates 9
/2
has the consequence that synonymous substitutions, say andselectionwithouthavingtoexplicitlymodeltheunder- /m
s
from A to C, occur at different rates at the first and third lying fitness landscape. In order to derive a codon model a
c
positions.Yet,althoughmodulationsofthemutationprocess along those lines, our strategy is to first assume a true site- 00
5
alongthesequencecannotbeexcluded,mostoftheempir- specific evolutionary process, following the mutation– /65
ically observed compositional differences between positions selection formalism. Then, we derive the mean substitution 03
5
are likely the consequence of selection, which is stronger at 0
processimpliedacrossallsitesbythismechanisticmodeland 5
the first and second than at the third position. In principle, identify the mean fixation probabilities appearing in this by g
theseselectiveeffectsshouldnotdirectlyimpactsynonymous u
mean-field (MF) process with the array of x tensor to be e
rates.Thus,althoughthemutationalprocessmightbemore estimated.Basedonthisapproach,weshowthatthesimplest st o
complex, there is no reason to model it in terms of a 3(cid:2)4 n
model that correctly teases apart mutation and selection 1
structurewhichconflatestwolevelsofmechanismsthatare requiresadifferentvalueofxforeachdistinctpairofamino 5 F
notsupposedtoplaytogether.Simulationexperimentssug- eb
gest that the 3(cid:2)4 formalism indeed leads to less accurate acids. Similar multirate models have been introduced previ- rua
ously(Delportetal.2010),althoughneverinconnectionwith ry
estimationofx(SpielmanandWilke2015). 2
the question of how to separately infer mutation rates and 0
Themutationmatrix(1(cid:2)4)ormatrices(3(cid:2)4)estimated 2
2
themeaneffectofselection.
bycodonmodelsarethusnotcorrectlyreflectingthemuta-
tion rates between nucleotides (Rodrigue et al. 2008;
Kosakovsky Pond et al. 2010). Instead, what these matrices Results
arecapturingistheresultofthecompromisebetweenmu-
Toillustratetheproblem,wefirstconductsimulationexperi-
tation and selection at the level of the realized nucleotide
mentsunderasimplemutation–selectionsubstitutionmodel
frequencies. Conceptually, it is a clear symptom that muta-
tion rates and fixation probabilities are not correctly teased assuming site-specific aminoacidpreferences.Weusethese
apartbycurrentcodonmodels. simulation experiments to explore through summary statis-
Practically,thismisconceptioncouldhaveimportantcon- tics the intricate interplay between mutation and selection.
sequencesinthecurrentinterestininvestigatingthevariation Then,weexplorehowcodonmodelswithdifferentparame-
betweenspeciesinGCcontent,anditseffectontheevolution terizations are able to infer the mutation rates and the
ofprotein-codingsequences.Animportantfactorhereisbi- strength of selection on these simulated alignments. Finally,
ased gene conversion toward GC (called gBGC), which can thesealternativemodelsareappliedtoempiricaldata.
2
MBE
ImprovedCodonModelingApproachforAccurateEstimationoftheMutationBias . doi:10.1093/molbev/msac005
SimulationExperiments Jones et al. 2017). As expected, h2N P i is always lower
e fix
Simulations of protein-coding DNA sequences were con- than 1 for simulations at equilibrium, under a time-
ducted under an origination–fixation substitution process independent fitness landscape (Spielman and Wilke 2015).
(McCandlishandStoltzfus2014) atthelevelofcodons(see Quite expectedly h2N P i decreases with the N (fig. 1D).
e fix r
Simulation Model). We assume a simple mutation process Ontheotherhand,h2N P i dependsweaklyonthemuta-
e fix
with a global parameter controlling the mutational bias to- tionalbias(k).
ward AT, denoted k¼ðrAþrTÞ=ðrCþrGÞ, where rx is Theproxyofselectionrepresentedbyh2NePfixiconcerns
the equilibrium frequency of nucleotide x. This mutational allnonsynonymousmutations,butwecanalsoconsiderthe
processissharedbyallsitesofthesequence.Withregardsto
mean-scaled fixation probability only for thesubset of non-
selection, synonymous mutations are considered neutral,
synonymous mutations from weak nucleotides (A or T)
such that the synonymous substitution rate equal to the to strong nucleotides (G or C), called h2N PWSi.
underlying mutation rate. At the protein level, selection is Interestingly, h2N PWSi increases with the strength eoffitxhe
modeledbyintroducingsite-specificaminoacidfitnesspro- e fix D
mutational bias toward AT (fig. 1E). This distortion of the o
files(i.e.,avectorof20fitnessesforeachcodingsite),which w
selective effects toward GC is stronger under an increased n
arescaledbyarelativeeffectivepopulationsizeN.AhighN lo
induces site-specific profiles having a large varriance, withr stringencyofselection,underahigherNr.Likewise,thenon- ade
someaminoacidswithahigh-scaledfitness,whereasallother s(GynCo)ntyomwoeuaskmnuutcaletoiotnidsecso(uAldT)a.lTsohisberarteiostdrieccteredasfreosmwistthrotnhge d from
havealow-scaledfitness.Conversely,alowvalueforN indu-
cesmoreevenaminoacidfitnessprofiles(i.e.,neutral)arteach strengthofthemutationalbiastowardAT(notshown).Asa http
site. Thus, ultimately, the stringency of selection increases result,theratiobetweenh2NePWfixSi andh2NePWfixSi ishigher s://a
wthiethmNurt.aAtilotongebtihaser(,kt)heantdwothpearsatrminegteenrscyofotfhseelmecotidoenl t(uNn)e, trehsapne1ctuivnedlyefroarmabuitaasttioonwaalrbdiaGsCt)o.wItairsdmAoTno(atonndolouwslyerintchraenas1- cade
r m
respectively. All simulations presented in the main manu- ingwiththemutationalbiastowardAT(fig.1F).Altogether, ic
.o
scriptareobtainedusingthesameunderlyingtreetopology fixationprobabilitiesareopposedtomutationalbias,andthe up
and branch lengths of 61 primates from Perelman et al. realized equilibrium frequencies are thus at an equilibrium .co
m
(2011), along with the same 4,980 codon sites with amino pointbetweenthesetwoopposingforces. /m
acid fitness profiles resampled from experimentally deter- be
/a
minedprofilesinBloom(2017). ParameterInferenceonSimulatedData rtic
Simulationunderthisorigination–fixationprocessalonga le
speciestreeresultsinamultiplesequencealignmentofcod- Fromanalignmentofprotein-codingDNAsequences,with- /39
ing sequences for the extant species, from which summary out knowing the specific history of substitutions, can one /2/m
statistics can then be computed. One such straightforward estimatethemutationalbias(k)andthemean-scaledfixation sa
c
summary statistic is the frequency ofthe different nucleoti- probability h2NePfixi? In other words, can we tease apart 005
des,andtheresultingnucleotidebiasAT=GCobservedinthe mutationandselection? /6
5
alignment. This observed nucleotide bias can be computed To address this question, here we consider two codon 03
5
separately foreachcoding position(first, second,and third) models for inference, differing only by their parametrization 0
5
andcomparedwiththetrueunderlyingmutationalbiask.As ofthecodonmatrixQ,whichwefirsttestagainstsimulated by
g
can be seen from figure 1, the third position of codons data(seefig.2).Botharehomogeneousalongthesequence u
e
(fig. 1C) reflects the underlying mutational bias quite faith- (i.e., not site-specific). The first is based on Muse and Gaut st o
fully,whereasthefirstandsecondpositions(fig.1AandB)are (1994)formalismandusesascalarxparameter,whereasthe n 1
impactedbythestrengthofselectionanddisplaynucleotide secondisbasedonatensorrepresentationofx. 5 F
biases that are less extreme than the one implied by the eb
mutational process. This differential effect across the three rua
ry
codingpositionsisexplainedbynucleotidemutationsatthe !AsaScalar:TheMuseandGautFormalism 2
0
thirdcodonpositionbeingmoreoftensynonymous,whereas 2
Thismodelisdefinedintermsofageneralizedtime-reversible 2
mutations at the first and second positions are more often
nucleotideratematrixRandascalarparameterx.Thematrix
changing the amino acid and are thus more often under
Risafunctionofthenucleotidefrequenciesrandthesym-
purifyingselection.
metricexchangeabilityratesq(Tavar(cid:2)e1986):
Apartfromthenucleotidebiasobservedinthealignment,
a statistic directly relevant for measuring the intrinsic effect R ¼q r : (1)
a;b a;b b
of selection is the mean-scaled fixation probability of non-
synonymousmutations,calledh2N P i.Thissummarysta- At the level of codons, the substitution rate between the
e fix
tistic h2N P i can be quantified from the substitutions source (i) and target codons (j) depends on the underlying
e fix
recorded along the simulation trajectory (see Mean-Scaled nucleotide change between the codons Mði;jÞ (e.g.,
Fixation Probability). For very long trajectories, it identifies MðAAT;AAGÞ¼TG), and whether or not the change is
withtheratioofnonsynonymousoversynonymoussubstitu- nonsynonymous.Altogether,thesubstitutionratesbetween
tionrates(d =d orx)inducedbytheunderlyingmutation– codonsQ ,formalizedbyMuseandGaut(1994)aredefined
N S i;j
selection model (Dos Reis 2015; Spielman and Wilke 2015; asfollows:
3
MBE
LatrilleandLartillot . doi:10.1093/molbev/msac005
A B C
D E F
D
o
w
n
lo
a
d
e
d
fro
m
h
ttp
s
://a
FIG.1.Simulationsof61primatestaxa,4,980codonsites,with100replicates.Solidlinesrepresentthemeanvalueoverthereplicates,andthecolored ca
d
areathe95%interquantilerange.Toprow(A–C):observedAT=GCcompositionofsimulatedalignment(first,second,andthirdcodingpositions),as e
m
afunctiontheunderlyingmutationalbiastowardAT(k),underdifferentstringenciesofselection(differentvaluesofrelativeeffectivepopulationsize ic
Nr).Bottomrow(D,E):mean-scaledfixationprobabilityofnonsynonymousmutationsalongsimulations,h2NePfixi,forallmutations(D)andforAT- .oup
to-GCmutationsonly(E),asafunctionofthemutationalbias(k),underdifferentrelativeeffectivepopulationsizes(Nr).(F)Ratioofmean-scaled .co
fixationprobabilityforAT-to-GCoverGC-to-ATmutations,asafunctionofthemutationalbiasandunderdifferentstringenciesofselection(N). m
r /m
Mutationalbiasisbalancedbyselectionintheoppositedirection,wherethiseffectincreaseswiththestringencyofselection. b
e
/a
rtic
le
/3
9
/2
/m
s
a
c
0
0
5
/6
5
0
3
5
0
5
b
y
g
u
e
s
t o
n
1
5
F
e
b
ru
FIG.2.Overallprocedureforsimulationunderasite-specificmutation–selectioncodonmodelandinferenceusingahomogeneouscodonmodels. ary
Thevalueofthemutationalbias(k)usedforsimulationscanbecomparedwiththevalueestimatedbythecodonmodels(bk)oncefittedtothe 2
simulatedalignment.Themean-scaledfixationprobabilityofnonsynonymousmutations(h2N P i)isrecordedalongthesimulationtrajectory, 02
e fix 2
andisdirectlycomparablewithxsuchasestimatedbycodonmodels.
b
8
>><Qi;j ¼0 if codons i and j are more than one mutation away; bkMG ¼ðrcA þrbTÞ=ðrcG þrcCÞ: (3)
Q ¼R if codons i and j are synonymous;
i;j Mði;jÞ Asforthemeanstrengthofselectionh2N P i,adirectes-
>>:Qi;j ¼xRMði;jÞ if codons i and j are nonsynonymous: timateisgivenbyxb. e fix
(2) Asshowninfigure3A,estimateofthemutationalbiasis
halfway between the nucleotide bias observed in the align-
The model can be fitted by maximum likelihood. Then,
ment and the true mutational bias used during the simula-
from the estimate of Rb, one can derive a nucleotide bias tion.Thus,theMGmodelcannotreliablyinferthemutational
towardATas: bias.Ontheotherhand,x isclosetotheunderlyingmean-
b
4
MBE
ImprovedCodonModelingApproachforAccurateEstimationoftheMutationBias . doi:10.1093/molbev/msac005
A B C
FIG.3.Simulationswith61primatestaxaand4,980codonsites.Estimatedversustruemutationalbias,usingacodonmodelinwhichxismodeled
asascalar(MGformalism,MG,panelA)orasatensor(MFapproach,panelB),orbyapplyingaGTRnucleotidemodeltothe4-folddegenerate
D
third-codingpositionsonly(panelC). o
w
n
lo
scaled fixation probability h2N P i computed during the parameterizedbyasetofxvaluesforallpairsofaminoacids, ad
e fix e
d
simulation (61 primates taxa, 4,980 codon sites, 100 repli- denotedxx;y.For20aminoacids,thetotalnumberofpairsof fro
cates),withaprecisionof97.2%.Thus,thefailuretocorrectly amino acids is 190, hence 380 parameters by counting in m
h
estimate the mutation process does not seem to have a both directions. However, because of the structure of the ttp
strong impact on the estimation of the overall strength of genetic code, there are 75 pairs that are one nucleotide s://a
selection,atleastinthepresentcase. away, since some amino acids are not directly accessible c
a
d
through a single nonsynonymous mutation. As a result, e
m
!AsaTensor:MFDerivation thenumberofparametersnecessarytodetermineallnon- ic
.o
Wewouldliketoderiveacodonmodelthatwouldbemore zero entries of the tensor (x ) in both directions is 150. u
x;y p
accurate than the MG model concerning the estimation of Finally, under the assumption of a reversible process, the .co
m
themutationbias,butthatwouldstillbesite-homogeneous. number of parameters can be reduced to 75 symmetric /m
However, the true process is site-specific. The link between exchangeabilities(b )and20stationaryeffects((cid:2) ): be
x;y x /a
thetwocanbeformalizedbyprojectingthesite-specificpro- x ¼(cid:2) b ; where b ¼b : (6) rtic
cessesontoagene-wiseprocess,usingwhatcanbeseenasa x;y y x;y x;y y;x le
/3
MF approximation (Goldstein and Pollock 2016). Thegene- 9
Altogether,thesubstitutionratesbetweencodonsQi;j are /2
wiseprocessobtainedbythisprocedureisexpressedinterms definedas: /m
s
of mutation rates and mean-scaled fixation probabilities. a
c
Ffiiendalwlyi,ththtehmexea-nte-sncsaoler.d fixation probabilities can be identi- 8>><Qi;j ¼0 if codons i and j are nonneighbors; 005/6
Q ¼R if codons i and j are synonymous; 5
i;j Mði;jÞ 0
8>>>><QSðip;zjÞecif¼ica0llyif,caotdeoancshisaitnedz,jtahreemunodreertlhyainngonceodmountaptiroonceawssayis;: >>:Qi;j ¼RMði;jÞxAðiÞ;AðjÞ if codons i and j are nonsynonymo(u7s); 3505 by gu
>>>>:QQðiði;;zzjjÞÞ ¼¼RRMMððii;;jjÞÞ2iNfecPoðfdizxÞoðni;sjÞiiafncdodjoanrseisyannodnyjmaoreusn;onsynonymous: gwivheTenrheibsAyMðeFiqÞumisaottihdoeenlai(ms6)fii.nttoedacbidyemnacxoidmeudmbylickoeldihoonoida,ngdivxingx;yains est on 15
(4) estimateforitsparameters,Rb; bb andb(cid:2).Then,fromtheesti- Feb
Where2N PðzÞði;jÞisthescaledfixationprobabilityofcodon mateoftheGTRnucleotidematrix(Rb),amutationbiasbkMF rua
e fix canbeestimatedaspreviously(eq.3above). ry
j against codon i, at site z. At equilibrium of the process, 2
averagingoversitesundertheequilibriumdistributiongives As shown in figure 3B, and under a variety of scenarios 02
2
(number of sites, branch lengths, tree topology) in
theMFgene-levelprocess:
Supplementary Material online, bkMF under the MF model
8
>><hQi;ji ¼0 if codons i and j aremorethanonemutationaway; provides an accurate estimate of the true mutational bias.
hQ i ¼R if codons i and j aresynonymous; Inotherwords,theMFmodelcanteaseouttheobservedAT
i;j Mði;jÞ
>>:hQ i ¼R h2NP ði;jÞi if codons i and j arenonsynonymous: =GCbiasofthealignmentandtheunderlyingmutationalbias
i;j Mði;jÞ e fix
Interestingly,inspiteofinvokingasinglemutationbiasacross
(5)
allnucleotidesites,theMFmodelpredictsdistinctnucleotide
However, because selection between codons reduces to se- frequencies at the three coding positions (Supplementary
lection between pairs of amino acids, h2N P ði;jÞi only Materialonline).Thesepredictedfrequenciesmatchthefre-
e fix
depends on the amino acids encoded by i and j (see quenciesthatareobservedonthealignment.Inotherwords,
Derivation of MF Model in Materials and Methods). Thus, the MF model is able to explain how a site-homogeneous
by identification, the inference model should be mutationalprocesscombinedwithaselectivepressureacting
5
MBE
LatrilleandLartillot . doi:10.1093/molbev/msac005
attheaminoacidlevelcanintheendproducea3(cid:2)4pattern estimated in the MF and MG models (h2N P i’x).
e fix b
ofnucleotidefrequencies. Additionally, in the MF model, h2N P i can be restricted
e fix
The mean-scaled fixation probability of nonsynonymous tomutationsfromweaknucleotides(AT)tostrong(GC),or
mutationsh2N P icanalsobecomputed.Itisnowacom- vice versa (see Mean-Scaled Fixation Probability h2N P i
e fix e fix
poundparameter,expressedasafunctionofRb; bbandb(cid:2)(see undertheMFModel).Weobservethatunderamutational
Mean-Scaled Fixation Probability h2N P i under the MF biasfavoringAT(i.e.,k>1),themeanfixationprobabilityof
e fix
Model). Under this model, h2N P i is close to the true nonsynonymousmutationsishighertowardGCthantoward
e fix
mean-scaledfixationprobabilityh2N P icomputedduring AT,h2N PWSi>h2N PSWi,asexpectedunderaAT-biased
e fix e fix e fix
the simulation, with a precision of 96.9% (61 primates taxa, mutationprocess.
4,980codonsites,100replicates).Moreover,asshowninfig- Reciprocally,foralignmentgloballybiasedtowardGC(b-
ure4,theestimatedratesx betweenpairsofaminoacidsis lactamase),theestimatedmutationbiasisstronger(toward
bx;y
congruentwiththepredictedmean-scaledfixationprobabil- GC) than the bias observed on the alignment (i.e.,
ity computed analytically as a function of the underlying bkMF < AT=GC < 1). Curiously, in b-lactamase, the MG Do
site-specific fitness profiles and the mutation matrix as in model estimates a weaker underlying mutational bias than w
n
equation(26). the observed bias (i.e., AT=GC < bkMG < 1). This effect loa
d
More analyses are shown in Supplementary Material on- couldbeduetothefirst,second,andthirdpositionshaving e
d
line,withdifferentsequencelength(498,996,2,490,4,980and compositionalbiasesindifferentdirections,whichisharderto fro
m
9,960 codon sites), different branch lengths (decreased by a disentangle (table 1). Concerning selection, we observe that
h
factor 2 and increased by a factor 2, 4, 8), and a different the fixation probability of nonsynonymous mutations is ttp
s
topology (90 mammals). These analyses suggest that the higher on average toward AT than toward GC, ://a
numberofsitesdoesnotinfluencetheestimator’saccuracy h2NePSfiWx i>h2NePWfixSi,asexpectedunderaGC-biasedmu- cad
formutationalbias(bk),norforselectionpressure(xb).Finally, tationprocess. em
for large sequence divergence (Supplementary Material on- The results obtained on empirical data are globally in ic
.o
line), saturation of sequences (multiple substitutions at the agreementwiththeobservationsgatheredfromthesimula- up
samesite)leadstolessaccurateestimation:boththeMGand tionexperiments,namelythatthepresenceofamutational .co
m
MFmodelsfailtogiveanaccurateestimatorofxb.Themu- bias results in a selection differential, taking the form of a /m
tation biasbk, on the other hand, is still correctly estimated slightly higher mean fixation probability of nonsynonymous be
/a
undertheMFmodel. mutations opposing the mutational bias. Moreover, by set- rtic
ting(cid:2)¼1andb¼x(cid:2)1inourMFmodel,weretrievethe le
EstimationonEmpiricalSequenceData nested MG model, hence, both models are directly /39
/2
The two alternative models of inference just considered, comparable. /m
namelytheMuseandGaut(MG)andtheMFcodonmod- Theempiricalfittothedatabetweenthenestedmodels, sa
c
els,werethenappliedtoempiricalprotein-codingsequence using AIC and Likelihood ratio tests (Posada and Buckley 00
5
alignments.Severalexampleswereanalyzed:thenucleopro- 2004) favor the MF model compared with the MG model /6
5
teininInfluenzaVirus(ashumanhost)assembledinBloom (table1).Ofnote,owingtoitsunreasonableassumptionthat 03
5
(2017), the b-lactamase in bacteria gathered in Bloom h2N P i is the same across all amino acid pairs, the MG 0
e fix 5
(2014), as well as orthologous AT-rich genes (such as to model is in fact very easy to improve upon (Delport et al. by
prevent the confounding effect of gBGC) in primates 2010),andthusthehigherfitofMFcomparedwithMGisnot gu
e
extracted from the OrthoMam database (Scornavacca in itself a very strong argument in favor of the use of MF. st o
etal.2019)asshownintable1. However,oursimulationssuggestthat,inspiteofthelarger n
1
For alignments globally biased toward AT (nucleoprotein estimationerrorontheindividualratesbetweenallpairsof 5 F
and AT-rich concatenate in primates), similarly to what was aminoacidsonsmalleralignments,theestimateofthemu- eb
observedinthesimulationexperimentspresentedabove,the tationbiasisalwaysreasonablyaccurate,evenonsmallalign- rua
mutational bias estimates under the two codon models are ments(SupplementaryMaterialonline). ry 2
0
greater than the observed nucleotide bias (i.e., Inanothersimulationanalysis,ithasalsobeenshownthat 2
2
1 < AT=GC < bk). This effect is, as previously, probably better fitting models could sometimes lead to less accurate
duetoselectionatthelevelofaminoacids,partiallyopposing inference (Spielman and Wilke 2015). This point was more
the mutational bias. More importantly, the mutational bias specificallymadeconcerningmodelssuchas3(cid:2)4.Weconcur
estimated by the MF model is more extreme than the MG withthisargument,whichisparticularlyrelevantinthepre-
estimate (i.e., 1 < bkMG < bkMF). These examples behave sent context. The 3(cid:2)4 model is typically better fitting than
identically to the observations made with simulated align- the1(cid:2)4model(whichisthedefaultconsideredherethrough
ments, where, compared with MG, theMF model estimates the MG model). Yet, and this is precisely one of the main
a stronger mutational bias, which was also closer to the real pointsofthepresentwork,3(cid:2)4doesnotrepresentthecor-
value. Thus, a reasonable interpretation is that MG is also rectwaytomodeltheprocessesthatarecreatingthevaria-
underestimatingtheunderlyingmutationalbiasinthepresent tion in nucleotide frequencies across the three coding
case,andthattheestimateoftheMFmodelismoreaccurate. positions and, for that reason should not be used, in spite
Concerning selection, theestimatedmean-scaledfixation ofitshigherfit.TheMFmodel,ontheotherhand,givesthe
probability of nonsynonymous mutations, is similarly correct logical solution to this problem and our simulation
6
MBE
ImprovedCodonModelingApproachforAccurateEstimationoftheMutationBias . doi:10.1093/molbev/msac005
D
o
w
n
lo
a
d
e
d
FIG.4. Trueversusestimatedvaluesofxx;y betweenpairsofaminoacidsunderourMFmodel.Thetruevaluesaregivenbyequation(26). from
Simulationson61primatestaxawith4,980codonsitesover100replicates.Verticalbarsarethe95%confidenceintervalsforthemeanvalue. h
ttp
s
Table1.MutationalBias(k)andMean-ScaledFixationProbability(h2NePfixi)EstimatedundertheMGandMFModelsonDistinctConcatenated ://a
DNAAlignmentsofOrthologousGenes. ca
d
e
m
b-Lactamase Nucleoprotein PrimatesAT-Rich ic
.o
DataSet Bloom Bloom Scornavaccaetal. u
p
Numberoftaxa 85 180 22 .c
o
Numberofsites 263 498 4,877 m
AT/GC 0.792 1.154 2.028 /m
b
AT/GCatfirstposition 0.583 1.057 1.303 e
/a
AATT//GGCCaattstehcirodnpdopsoitsiiotnion 10..177174 11..212912 22..564418 rticle
MGmutationalbias(bkMG) 0.853 1.447 2.073 /39
MMFGmxutationalbias(bkMF) 00..639302 10..714184 20..451296 /2/m
b s
MFh2NePfixi 0.336 0.116 0.525 ac
0
MFh2N PWSi 0.297 0.141 0.594 0
e fix 5
MDAFIhC2NePSfiWxi 03.471.62 01.6059.22 10,5.42877.0 /6503
Pðv2 >LRTÞ 9.2310(cid:3)13 1.2310(cid:3)31 3.9310(cid:3)296 50
df¼93 5
b
NOTE.—TheMFmodelcontains95parameterswith75aminoacidexchangeabilitiesand20aminoacidequilibriumfrequencies.However,wehavetwoconstrainsthatreduce y g
thedegreeoffreedom;thesumofall20aminoacidequilibriumfrequenciesequals1andthesumofall61codonfrequenciesequals1. u
e
s
t o
n
experimentsconfirmthatthisleadstoaccurateestimationof level of nucleotides and selection at the protein level. As a 1
5
themutationbias.Insummary,thisistheconjunctionofthe result,thenucleotidebiasobservedinthealignmentisdiffer- F
e
b
higher fit observed here on empirical data with the logical ent from the underlying mutational bias. However, current ru
a
argumentsandthesimulationexperimentspresentedabove parametriccodonmodelspredictthattheobservedandun- ry
2
that together justify the use of the MF model. Based on derlying mutationalbiasesshouldbeequal. For that reason, 0
2
these justifications, we can thus interpret the estimate of theyareinherentlymisspecifiedandareunabletoteaseapart 2
bkMF as reflecting the mutation bias, and the difference be- opposing effects of mutation and selection correctly. As
tweenh2N PSWi andh2N PWSi assuggestingthatthefixa- shown in our work, the misspecification of these models
e fix e fix
tionbiasesaredifferentinthetwodirectionsalsointhecase doesnotstronglyimpacttheestimationoftheneteffectof
ofempiricaldata. selection on nonsynonymous mutations (x). This novel re-
b
Altogether,ourMFmodelisfavoredbyempiricaldatasets, sult is important, as it is reassuring for a certain number of
and simultaneously estimates more extreme (and probably previouslypublishedanalyses,inparticularcorrelatingxwith
b
more accurate) mutational biases compared with the MG life-historytraits,inacontextwhereGCcontentalsocorre-
model. lates with life-history traits (Figuet et al. 2016; Bol(cid:2)ıvar et al.
2019).However,currentparametricmodelsdonotestimate
Discussion
themutationalprocessaccurately.
In protein-coding DNA sequences, the nucleotide composi- In this work, we sought to find the simplest parametric
tionresultsfromasubtleinterplaybetweenmutationatthe codonmodelabletocorrectlyteaseapartmutationrateson
7
MBE
LatrilleandLartillot . doi:10.1093/molbev/msac005
one hand, and net mean fixation probabilities on the other thesameinthetwodirections,eithertowardATortoward
hand,andthis,withouthavingtoexplicitlymodeltheunder- GC. Specifically, in the presence of an AT-biased mutation
lying fitness landscape. In order to derive a codon model process, the nonsynonymous GC polymorphisms are
alongthoselines,ourstrategyistofirstassumeanunderlying expectedtosegregateathigherfrequencies,comparedwith
microscopicmodelofsequenceevolution(here,amutation– nonsynonymousATpolymorphisms.
selection model based on a site-specific, time-independent These observations have some practical implications: for
fitnesslandscape).Then,wederivethegene-wisemeanfixa- instance, experiments observing a fixation (or segregation)
tionprobabilitiesbetweenallpairsofcodons,impliedbythe biastoward GC at thenonsynonymouslevel must alsorule
underlyingmicroscopicprocess.Finally,weobservethatthis outthatthisfixationbiasisnotasimpleconsequenceofthe
MFprocessshouldinfactinvokeasmanydistinctxparam- tradeoffbetweenmutationandselection.Moregenerally,our
etersastherearepairsofaminoacidsthatarenearestneigh- observations and modeling principles offer a useful prelimi-
bors in the genetic code. There are reversibility conditions, narybasistobetterunderstandhowmutationandselection
reducingthedimensionalityandallowingforaGTR-likepa- willwork togetherwithGC-biased geneconversion(gBGC), D
o
rameterizationofthistensor(95parametersforselection). and therefore will help better understand how gBGC will w
n
Inferring parameters on simulated alignments, we show impactbothnucleotidecompositionandx.Itisworthmen- lo
b a
d
that the model derived using this MF argument correctly tioningthatinourresult,wefocusedonthefixationproba- e
d
estimatestheunderlyingmutationalbiasandselectivepres- bilityfromATtoGC,h2NePWfixSi,becauseoftherelationship fro
m
sure.Inthisrespect,ourworkgivesthefirstclearexplanation togBGC.However,inpractice,thesameanalysisandmeth-
h
of how to correctly disentangle the underlying mutational odscanbeappliedtoanysubsetofnucleotidesorcodons. ttp
s
bias and the observed nucleotide frequencies. Our model OurMFparametricmodelusesgene-levelparameters(in ://a
canpredicttheaccuratenucleotidecompositionatfirst,sec- the form of a tensor) that is meant to capture the mean- ca
d
ond,andthirdcodonpositions,whereascurrentparametric scaledfixationprobabilities.Thisderivation,anditsvalidation e
m
modelsfailtopredictthem.Wearguethatparametriccodon on simulated data, shows that, even though the underlying ic
.o
modelsusingthreedifferentmutationalprocessesatthefirst, selectivelandscapeissite-specific,agene-levelapproximation u
p
second, and third coding positions (3(cid:2)4 formalism) to ac- can nonetheless accurately disentangle mutation and selec- .co
m
commodateforvariationinobservednucleotidefrequencies tion. As a result, this study demonstrates that phenomeno- /m
isnotatheoreticallysoundmodeling.Indeedthisvariationis logical models derived out of mechanistic models are more be
/a
anemergingpropertyofthebalancebetweenmutationand compact(i.e.,notsite-specific),andincertaincasesaresuffi- rtic
selectionasshowninourwork.The3(cid:2)4formalismhaspre- cienttoextracttherelevantparameters. le
/3
viously been shown to lead to inaccurate inference of x The methodology proposed here for deriving inference 9
b /2
(Spielman and Wilke 2015). Altogether, we concur in this models consists in proceeding in two steps, first assuming /m
s
direction that 3(cid:2)4 formalism is inaccurateand not mecha- anunderlyingmechanisticmodelofsequenceevolution,pa- a
c
nistically sound, and as a result should not be used to esti- rameterizedbyvariablesthatarederivedfromfirstprinciples 00
5
matex. (fitness landscape, mutation rates, ...). Subsequently, the /6
b 5
Applyingourmodeltoempiricalalignments,wealsoob- phenomenologicalinferencemodelisobtainedbymatching 03
5
servethatthereisaselectiondifferentialopposingthemuta- its parameters (here, the entries of the x tensor) with the 0
5
tional bias. This observation also points to a fundamental aggregateparametersderivedfromtheapplicationoftheMF by
g
propertyofnaturalgeneticsequences,namelythattheyare proceduretothemechanisticmodel.Altogether,webelieve u
e
not optimized but are the result of interactions between thattheapproachusedherecouldbeappliedmoregenerally: st o
evolutionary forces (Sella and Hirsh 2005). In the specific inference models can be phenomenological in practice, but n
1
case highlighted in this work, the mutational bias at the shouldnonethelessbederivedfromanunderlyingmechanis- 5 F
nucleotide-levelresultsinsuboptimalaminoacidsbeingover- ticmodel,soastocorrectlyformalizetheinterplaybetween eb
representedinthesequence,comparedwithwhatwouldbe mutation,selection,drift,andotherevolutionaryforces. rua
ry
expectedbasedon theirfitness alone. For example, under a Ourphylogeneticcodonmodelisnotthefirsttomodelx 2
0
mutational bias toward AT, AT-rich amino acids might not as a tensor. Thus, Yang et al. (1998) introduced a codon 2
2
necessarilybethefittestbutareexcessivelygeneratedbythe modelinwhichx dependsonthedistancebetweenamino
mutationalprocess,resultinginastrongerpurifyingselection acids, measured in terms of the Grantham (1974) distance.
againstAT-richaminoacids.Thiswaspointedoutpreviously Additionally, Tang and Wu (2006) leveraged x tensors in
(SingerandHickey2000),althoughneverdirectlyformalized ordertodetectpositivelyselectedgenes.Thenoveltyofthe
inaphylogeneticcodonmodel.Oneimportantconsequence present work is to formalize the articulation between the
of this tradeoff between mutation and selection is that the nucleotide composition, the mutational bias, and selection
observedhighermeanfixationprobabilitytowardGCismim- between different amino acids. Finally, this work is still pre-
icking the effect of biased gene conversion toward GC liminarysincetheMFmodelshouldbetestedagainstamore
(gBGC), although unlike gBGC, the phenomenon described diverse range of empirical data, in terms of phylogenetic
herecorrespondstoagenuineselectiveeffect.Althoughwe depth, strength of selection, and codon usage bias to assert
did not explore the consequences of this at the level of in- thevalidityofourempiricalresults.Inaddition,severalother
traspecificpolymorphism,theselectiondifferentialuncovered parametrizationofcodonmodelsaslistedinRodrigueetal.
herealsoimpliesthatthedistributionoffitnesseffectsisnot (2008)andKosakovskyPondetal.(2020)shouldbeincluded
8
MBE
ImprovedCodonModelingApproachforAccurateEstimationoftheMutationBias . doi:10.1093/molbev/msac005
inabroadercomparisonoftheaccuracyoftheestimationof SelectionattheAminoAcidLevel
theunderlyingmutationalbiasandstrengthofselectionon The substitution rate is considered null between any two
protein-codingDNAsequences. codons differing by more than one nucleotide. Otherwise,
themutationratebetweenapair ofcodons isgiven by the
Materials and Methods
mutation rate of the underlying single nucleotide change.
SimulationModel Selection is modeled at the amino acid level, that is, we as-
Weseektosimulatetheevolutionofprotein-codingsequen- sumethatallcodonsencodingforoneparticularaminoacid
cesalongaspecietree.Startingwithonesequenceattheroot areselectivelyequivalent.
of the tree, the sequences evolve independently along the To take into account the heterogeneity of selection be-
different branches of the tree by point substitutions, until tweendifferentsitesoftheprotein,weassumethateachsitez
they reach the leaves. At the end of the simulation, we get of the sequence is independently evolving under a site-
onesequenceforeachleafofthetree,meaningonesequence specific fitness landscape, characterized by a 20D frequency
vector of scaled (Wrightian) fitness parameters D
perspecies.Thesubstitutionismodeledusingtheorigination– o
fixationapproximation,thatis,substitutionratesaretheprod- wðzÞ ¼fwðazÞ; 1 (cid:4) a (cid:4) 20g. The fitness vectors wðzÞ used wnlo
uctofthemutationrateatthenucleotidelevel,andfixation in this study are extracted from Bloom (2017), which were ad
e
proTbhaebilmitiuest,abtiaosnedporonceseslsecistioanssautmtehdeahmominoogaecnideoleuvsela.cross e4x9p8ecroimdoenntsaitlleysdoefttehremninuecdleobpyrdoeteeipnminuIntafltuioennzaalsVciarnunsisntrgafinosr d from
sites.Ontheotherhand,selectionisassumedtobevarying (ashumanhost).Foreachcodonsitezofoursimulation,we h
alongthesequence.Duringthesimulation,giventhecurrent assign randomly one the 498 fitness profile (sampling with ttps
sequence,thesubstitutionratestowardallpossiblemutants replacement) experimentally determined, which altogether ://ac
(onenucleotidechange)arecomputedandthenextsubsti- determines the (Wrigthian) fitness vectors across sites. The ade
malthusian fitness (or log-fitness) of amino acid a, denoted m
triutthimon(Gevilelenstpiies 1d9ra7w7)n. randomly based on Gillespie’s algo- FaðzÞ, is scaled by the relative effective population size (Nr) ic.ou
accordingly: p.c
o
MutationalBiasattheNucleotideLevel m
FðzÞ ¼Nln(cid:4)wðzÞ(cid:5); z2f1;...;Zg; a2f1;...;20g (13) /m
The mutation rate between nucleotides is always propor- a r a b
e
tional to l. Moreover, mutations from any nucleotide to /a
anotherweaknucleotideareincreasedbythefactorkcom- At site z, the substitution rate between nonsynonymous rtic
codonsiandjisgivenbytheproductofthemutationrate le
paredwithmutationstoanotherstrongnucleotide.Themu- /3
tationratematrixisthus: andtheprobabilityoffixation: 9/2/m
FðzÞ (cid:3)FðzÞ sa
A C G T QðzÞ ¼R AðjÞ AðiÞ (14) c0
0(cid:3)lð2þkÞ l l lk 1 i;j Mði;jÞ1(cid:3)eFAðzðÞiÞ(cid:3)FAðzðÞjÞ 05/6
A 5
B C 0
B C whereAðiÞdenotestheaminoacidencodedbycodoni.At 3
B C 5
R¼ GCBBBBBB llkk (cid:3)lð1lþ2kÞ (cid:3)lð1lþ2kÞ llkk CCCCCC fptrhoðezmÞ,rwotohhteicohsftaitsthigoeinvtearnreyeb,dyfi:ostrriebaucthiosniteofz,ththeepsroeqceusesnscpeeicsifdierdawbny 05 by gue
B C s
TB@ CA 2 3 t on
lk l l (cid:3)lð2þkÞ pðizÞ ¼ZðzÞ4 Y ri½k(cid:5)5eFAðzÞðiÞ; (15) 15 F
(8) e
k2f1;2;3g b
ru
Whichhasthefollowingstationarydistribution: wherei½k(cid:5)denotesthenucleotideatpositionk2f1;2;3gof ary
2
r¼1; (9) codoni,andZðzÞ isthenormalizingconstantatsitez: 02
2
0 2 3 1(cid:3)1
()r¼(cid:2)2þk2k;2þ12k;2þ12k;2þk2k(cid:3): (10) ZðzÞ ¼@X61 4 Y rj½k(cid:5)5eFAðzÞðjÞA (16)
j¼1 k2f1;2;3g
As a result, the ratio of weak over strong nucleotide fre-
Thesubstitutionprocessisreversibleandfulfillsdetailedbal-
quenciesatstationarityisequaltok:
ance conditions at each site z and between each pair of
r þr kð2þ2kÞ(cid:3)1þkð2þ2kÞ(cid:3)1 codons(i,j):
A T ¼ ; fromequationð10Þ; (11)
rCþrG ð2þ2kÞ(cid:3)1þð2þ2kÞ(cid:3)1 pðzÞQðzÞ ¼pðzÞQðzÞ (17)
i i;j j j;i
¼k: (12)
Of note, by modeling fitness at the amino acid level, we
P
l is constrained such the expected flow ((cid:3) r R ) of assume that all codons encoding for one particular amino
a a a;a
mutationequalsto1. acid are selectively equivalent. In addition, in this modeling
9
MBE
LatrilleandLartillot . doi:10.1093/molbev/msac005
framework,thegeneticcodeisofparticularimportancesince Z
X
thenumberofcodonsencodingforaparticularaminoacid hQi;ji¼ PðzjiÞQi;j (21)
variesgreatly.Asanexample,tryptophanisencodedbyone z¼1
codon,whereasleucineisencodedbysixcodons.Intuitively, If codon i and codon j are synonymous, this equation sim-
this variation makes the mutation bias more pronounced plifies to the underlying mutation rate R . Otherwise, if
Mði;jÞ
amongcodonsencodingforthesameaminoacid,sincethere codoniandcodonjarenonsynonymous,theMFsubstitution
aremoremutationspossiblethatareselectivelyneutral(i.e., rateis:
synonymous).Ontheotherhand,themutationbiasismore
hQ i¼hR 2N P ði;jÞi; (22)
constrainediftheaminoacidisencodedbyfewcodons. i;j Mði;jÞ e fix
Mean-ScaledFixationProbability ¼RMði;jÞh2NePfixði;jÞi; (23)
ThesequenceattimetisdenotedSðtÞ andthecodonpre-
sent at site z is denoted SzðtÞ. For a given sequence, the PZ pðzÞ FAðzÞðjÞ(cid:3)FAðzÞðiÞ Do
mean-scaled fixation probability over mutations away from ¼R z¼1 i 1(cid:3)eFAðzÞðiÞ(cid:3)FAðzÞðjÞ; (24) wnlo
SðtÞ,weightedbytheirprobabilityofoccurrence,isgivenby Mði;jÞ PZ pðzÞ ade
theratio: z¼1 i d fro
m
Z
h2NePfixðtÞi¼zPzP¼¼Z11PPjj22NN(cid:4)(cid:4)SSzzððttÞÞ(cid:5)(cid:5)QlSSzzððttÞÞ!!jj; (18) ¼RMði;jÞzP¼Z1ZzP¼Zð1zÞZe(cid:3)ðzFFAÞðAðzzeÞððÞijÞÞF(cid:3)(cid:3)AðzðÞFeiÞA(cid:3)ðzFðÞiAðÞzÞðjÞ (25) https://academic
.o
owfhceordeoNnðiiaÞnidstQhðezÞsaerteodfenfionnesdynasoninymeqouuasticoond(o1n4s).nAevigehrabgoerds tAasrgaetrecsoudlto,nh2sNoleePlyfixtðhi;rojÞuighistdheepesonudrecnetaomninthoeascoidur(cxe) aanndd up.co
i;j m
overallbranchesofthetree,themean-scaledfixationprob- target amino acid (y), hence the parameter xx;y identifies /m
abilityis: withtheaveragefixationprobabilityh2N P ðx!yÞi: be
e fix /a
h2NePfixi¼ðth2NePfixðtÞidt; (19) h2N P ðx!yÞi¼zP¼Z1ZðzÞe(cid:3)FFxðyðzzÞÞ(cid:3)(cid:3)Fex(cid:3)ðzFÞyðzÞ: (26) rticle/39/2
e fix Z /m
where the integral is taken over all branches of the tree, PZðzÞeFxðzÞ sa
whereas the integrand h2NePfixðtÞi is a piece-wise function z¼1 c00
changing after every point substitution event. The mean- 5
/6
scaled fixation probability from weak (AT) to strong (GC) 50
Mean-ScaledFixationProbabilityh2N P iunderthe 3
nucleotides, denoted h2NePWfixSi, is obtained similarly by MFModel e fix 505
restrictingthesums(inthenumeratorandthedenominator) b
fromweak to strong mutations. A similar computation can TheMF model is parameterized by a GTR mutation matrix y g
u
Rðr;qÞandtheselectioncoefficientxðb;(cid:2)Þ.Asaresult,the e
bedonefromstrongtoweak. s
mean-scaled fixation probability of nonsynonymous muta- t o
n
DerivationofMFModel tionsis: 15
F
The MF codon model hQi is defined such that hQ i is the e
i;j 61 b
average rate of substitution to codon j, conditional on cur- Pp P Q ru
i i;j a
rentlybeingoncodoni,theaveragebeingtakenacrosssites. h2N P i¼ i¼1 j2NðiÞ ; (27) ry 2
Importantly,sitesdifferintheirprobabilityofbeingcurrently e fix P61 P 02
p l 2
in state i. Theaverage should therefore be weighted by this i i;j
i¼1 j2NðiÞ
probability.
Assuming an underlying site-specific mutation–selection
61
P Q P
process at equilibrium, given we know that a mutation is ½ ri½k(cid:5)(cid:5)(cid:2)AðiÞ j2NðiÞRMði;jÞ(cid:2)AðjÞbAðiÞ;AðjÞ
fromcodoni,theprobabilitythatthismutationisoccurring ¼i¼1 k2f1;2;3g ;
61
atsitezis: PQ P
½ r (cid:5)(cid:2) R
k2f1;2;3g i½k(cid:5) AðiÞ j2NðiÞ Mði;jÞ
pðzÞ i¼1
PðzjiÞ¼ i (20) (28)
Z
PpðizÞ wherei½k(cid:5) denotesthenucleotideatpositionk2f1;2;3g of
z¼1
codoni.
Thesite-averaged(MF)substitutionratefromcodonitojis Similarly, the mean-scaled fixation probability from weak
asresultgivenas: (AT) to strong (GC) nucleotides denoted h2N PWSi is
e fix
10