From 3f95a1febd2e5a4d6c67c9890596436cabcebf6a Mon Sep 17 00:00:00 2001 From: Ryan Gutenkunst Date: Sat, 28 Sep 2024 15:10:42 -0700 Subject: [PATCH] Add preliminary figure describing Vaquita hs fitting --- analysis2_manuscript.tex | 17 +++-- figures/PhoSin/Vaquita2Epoch_1R22/hs_fig.py | 70 ++++++++++++++++++ figures/PhoSin/Vaquita2Epoch_1R22/hs_hist.pdf | Bin 0 -> 21494 bytes 3 files changed, 81 insertions(+), 6 deletions(-) create mode 100644 figures/PhoSin/Vaquita2Epoch_1R22/hs_fig.py create mode 100644 figures/PhoSin/Vaquita2Epoch_1R22/hs_hist.pdf diff --git a/analysis2_manuscript.tex b/analysis2_manuscript.tex index a16d118..4c0e319 100644 --- a/analysis2_manuscript.tex +++ b/analysis2_manuscript.tex @@ -498,15 +498,19 @@ \section*{Estimation of the Distribution of Fitness effects} \begin{figure} \centering \includegraphics[width=\textwidth]{figures/PhoSin/Vaquita2Epoch_1R22/PhoSin_Vaquita2Epoch_1R22_Gamma_R22_Phocoena_sinus.mPhoSin1.pri.110_exons_DFE_plot.pdf} +\includegraphics{figures/PhoSin/Vaquita2Epoch_1R22/hs_hist.pdf} \caption{ \label{fig:vaquita-dfe} A comparison of methods for inferring the distribution of fitness effects (DFE) from genetic data. Simulations were performed using a model of vaquita porpoise demography \citep{robinson2022critically} with a Gamma distributed DFE - acting on exons parameterized by a mean selection coefficient and a shape parameter. Estimates of the + acting on exons parameterized by a mean selection coefficient, a shape parameter, and a relationship between selection and domiance. Estimates of the parameters of the DFE are shown in the left ($\lvert E(s) \rvert $) and right hand panels (shape) respectively. This vaquita model has a single extant population, labelled pop 0, and parameter estimates from each of the three different methods are shown: 1) GRAPES \cite{galtier2016adaptive}, 2) polyDFE \citep{tataru2020polydfe}, - and 3) dadi \citep{gutenkunst2009inferring}.} + and 3) dadi \citep{gutenkunst2009inferring,kim2017inference}. + C) Shown is the simulated distribution of $2 h s$, where $h$ is the dominance coefficient, along with the simulated distribution of $s$ and the median inferred DFEs. The distribution of $2 h s$ is multimodal because the simulated relationship between $h$ and $s$ is not continuous. + % RNG: For C, need to extract inferred values from pipeline. + } \end{figure} In Figures \ref{fig:vaquita-dfe.constant} and \ref{fig:vaquita-dfe} we show a comparison of estimates @@ -514,10 +518,11 @@ \section*{Estimation of the Distribution of Fitness effects} vaquita porpoise genome. In the constant size model (fig \ref{fig:vaquita-dfe.constant}) we see that all methods perform uniformly worse that in the human genome simulations, which consistent underestimation of the mean selection coefficient and overestimation of the shape parameter. - In context of our model of vaquita porpoise demography \citep{robinson2022critically}, which models a population - decline towards the present, we see that methods to infer the DFE perform no better-- again each method - signficiantly underestimates the mean selection coefficient and overestimates the shape parameter. - %ADK: add why we think this is the case-- is this a bug? or a feature of the genome? +A critical distinction is that the previous human simulations employed a DFE that assumed all selected mutations had additive effects (no dominance, $h = 1/2$), whereas the vaquita simulations use a DFE in which more deleterious mutations are more recessive ($h < 1/2$). +For all three software packages, the inferred mean and shape of the distribution of selection coefficients is consistent with the mean and shape of the distribution of $2 h s$, the product of the dominance coefficient and the selection coefficient (Fig.~\ref{fig:vaquita-dfe}C). +This is likely because deleterious alleles are typically at low frequency and thus heterozygous, where their selective effect is $h s$, and all these tools assume additivity ($h = 1/2$). +Strongly deleterious alleles are often observed to be recessive in many species \citep{}, so these results suggest caution in interpreting DFE inferences. +%RNG: Could use Kirk's help with the best references here \section*{Detection of Selective Sweeps} \label{sweeps} diff --git a/figures/PhoSin/Vaquita2Epoch_1R22/hs_fig.py b/figures/PhoSin/Vaquita2Epoch_1R22/hs_fig.py new file mode 100644 index 0000000..ae7a1ce --- /dev/null +++ b/figures/PhoSin/Vaquita2Epoch_1R22/hs_fig.py @@ -0,0 +1,70 @@ +import numpy as np +import scipy.stats.distributions as ssd + +# From https://github.com/popsim-consortium/stdpopsim/blob/main/stdpopsim/catalog/PhoSin/dfes.py +gamma_shape = 0.131 +gamma_mean = -0.0257 +dominance_coeff_list=[0.0, 0.01, 0.1, 0.4] +dominance_coeff_breaks=[-0.1, -0.01, -0.001] + +# Construct the corresponding gamma distribution +gamma_scale = -gamma_mean / gamma_shape +dfe = ssd.gamma(gamma_shape, scale=gamma_scale) +assert(dfe.mean() == abs(gamma_mean)) + +# Simulate draws from the DFE +N = 100000 +s_array = -dfe.rvs(N) + +# Sanity check, using method of moments to recover simulated shape and scale. +# (Fixing location parameter to zero.) +# mean = shape * scale, var = shape * scale**2 +# So scale = var/mean, and shape = mean/scale +fit_s_scale = s_array.var()/s_array.mean() +fit_s_shape = s_array.mean()/fit_s_scale +# Or just use fit_s_shape, _, fit_s_scale = ssd.gamma.fit(s_array, method='MM', floc=0) +print('Sanity check: Inferred s shape and mean: {0:.4f}, {1:.5f}'.format(fit_s_shape, fit_s_shape*fit_s_scale)) + +# Convert to 2*hs +breaks = [-np.inf] + list(dominance_coeff_breaks) + [0] +hs_arrays = [] +for h, s_start, s_end in zip(dominance_coeff_list, breaks[:-1], breaks[1:]): + hs_arrays.append(2*h*s_array[np.logical_and(s_start < s_array, s_array < s_end)]) +hs_array = np.concatenate(hs_arrays) + +# Calculate mean hs +mean_hs = hs_array.mean() +print('Simulated mean 2hs: {0:.6f}'.format(mean_hs)) + +#fit_hs_shape, _, fit_hs_scale = ssd.gamma.fit(hs_array, method='MM', floc=0) +#dfe_hs = ssd.gamma(fit_hs_shape, scale=-fit_hs_scale) +# +#print('Inferred 2hs shape and mean: {0:.4f}, {1:.5f}'.format(fit_hs_shape, fit_hs_shape*fit_hs_scale)) + +# Results from inferences +dfe_dadi= ssd.gamma(0.33, scale=0.002) +dfe_polydfe = ssd.gamma(0.33, scale=0.002) +dfe_grapes = ssd.gamma(0.35, scale=0.0035) + +# Plotting +import matplotlib.pyplot as plt +plt.rc('font', size=10) + +fig = plt.figure(2132, figsize=(3.5,2.5)) +fig.clear() +ax = fig.add_subplot(1,1,1) +# Histogram of 2*hs, scaled by 10^3 +ax.hist(hs_array*1e3, bins=101, density=True, label='Simulated 2hs distribution') +# DFE inferred for 2hs +xx = np.linspace(-2.2, 0, 1001) +ax.plot(xx, dfe.pdf(-xx/1e3)/1e3, '-k', label='Original s DFE') +ax.plot(xx, dfe_dadi.pdf(-xx/1e3)/1e3, '-', label='Inferred DFE (GRAPES)') +ax.plot(xx, dfe_polydfe.pdf(-xx/1e3)/1e3, '-', label='Inferred DFE (polyDFE)') +ax.plot(xx, dfe_grapes.pdf(-xx/1e3)/1e3, '-', label='Inferred DFE (dadi)') +ax.set_xlabel(r'$2hs\,\,(\times 10^{-3})$') +ax.set_ylabel('Probability density') +ax.legend(loc='upper left') +ax.set_ylim(0, 2) +ax.set_xlim(-2.2,0) +fig.tight_layout(pad=0) +fig.savefig('hs_hist.pdf') diff --git a/figures/PhoSin/Vaquita2Epoch_1R22/hs_hist.pdf b/figures/PhoSin/Vaquita2Epoch_1R22/hs_hist.pdf new file mode 100644 index 0000000000000000000000000000000000000000..0cae2df8eb027bd45baef8101177d344c0016acd GIT binary patch literal 21494 zcmeIa1yqz>@Hk9}2qGPV!~!CcyKKQy(jh4zB_K#k!;*^9A>EAU!C)v56|_vd!HLKbLZZfJ9p;cP?wS86X1sva}>V@71a|9 zfFK};YuAZIML}RKPe*eQSlR?_V(VZD0;`)?n!A9YKm`quxHz%7y&0yW&`%BI9PH5` z7`6gfSMAz$b5k@(@bIUs2U=beZGtuj!4C`6P0(m_XL}F=^OG2?X=P$&ZEp!eet+xi zV5(`323-M$m5~ELG50`&z>0PN2vXml(%+woAj6-)5%?Ve4BQy#oy`I04#5X&n7cT* zI-3IcVBpsPfmO`StWBgGJb)1)z&|Jg3K9^4fDDPj(g0FGOBWF8u%@KFy#w$|08{N} zLSSD1pq#R~y(QWTB=CJSIcr-0R}fgv79fC(xv7JhIR<1Gw6nR19kECHN=8pCQ%y|J z<_2}uS-N8}wHy!Q<)Ceo;apG)+8I6dv$9p$X^+MR&+q2mVmctwqIu+4KVhaoA)~B7 zTK{b~pMQ(E>!do)&Z3UP!Rm*N35QJYy8X_uIqQS<9k;9v!rdD+uYmu*x!&KP-+z3C zSNshxzwmBj=>@JH2e)tAdyj6`i0GEyJqeFl*-&cHE;PJ?a8P=0SD6mSc&%gZ=`8gSYu0+aMx^;!v6n>2 zayXE8Nmh0y?P!I;XUS5tG3)A_#tgX4XU#ve#$6N6)z5v}YpHQxwHX+UE3H>Z&vTVO zOAfDH8GY{Zhhxg@&Cgqf&UfvRjhfd(+nKjn56nac>FzKn~QngQ<#DSRaPNR4(Oy1WN zvmQ+OV12As?-2qCEeBfq&3ULI^I!!U7cLAD;3RiZTe1_6rI0G>r$p6l0mWFrXHQvP z^}`Gkq298+uUvSDqih`B4W{(^00zZUb*>-ygsPnFsZVv2?PP67T}Cs6k=?_yA*+4> zK;uM!r0&F@XH}JeCmvXEnK8n(_wx9&);X@Aape1dP&26bxnlbSP@E#BKd#{Q07H9Y z;=J4`?D&h)i9MwE8ia5MbZ?S#_)l|`W!okB$=l?l52k2h=gIh)qp_Ky>;%8{YMw&= znkOBO;K|guP=WJz-Op>daWcI+`ND8-SH|=+!HGo4UextmX=(P8iT)RS-RzEMm+!Vx zhfB*NXr803&R^ip54h>;)dh2mAZbk?3G}(xZxDRWexiO^-hiX8838?^gB)@Gtmnd% zs->(f{tZ{B8(Z}aoJDBDj`vu7qpX)fsM2ywa1$ z-$_x#Jf|5n2+ngpp(3Jb;iL%CqL^2o$}J#Ox#7_+oDkInBIia<2Hy00#qpAe(=C)x zi=$5hL4HC@Ai&g*w8PNf@}7Y+qj7>T{%{BmuY8A`)PaGC>nWU;o%j5z4;zGWhY5Vu zvzFusVn?-tkDugDm%4*r9BVnnKp((!I`I+>Q$z!8b|}*M1b*rF(v$!L@RHSQ1db!g zg15mG^BRc(yT7Bk=XyLE{!xF-sABHqYw1j18#A1wnvQ(NoX8GCDrw@XEB8&vobd9( z)=tL3UT6+{k%&Kjmky?3aN!=e=kd}Q$yz^Umi{(6Lt)J~L?ZaZ5+L<~_~h{qAGJsE zyx8bH+i{9s(m#+p)^aT_qe4%9rWCh%PT5s=QTfK<3|uBz=6D#{X-Fk+R9xhr%5X^S zDPps0U8w5^lL#$k{>UqtS&j*!xT6IAj#*K@By3js167*dT(9dsEcr>+X1AntFP18f zstYJf8tG*?C5YmWc7Y6R`mD^3~x2v;iz+gIL<8Cc2T#1|#2PYt7$KkCQ7 zQ@fiMwVDrfQW}G6pxnXfPy&#{|4P)PGqcrqmhzl2Sd+c6RoX3;6z_8WAC_WRUo>a+ zoffWyjb_vyE*<62%-D~ft|k1_pLW7Gwf*u0nTm$wi@Jzzp6{zQ2~FayAk5}|3a2iY5p3n|db}`3;=vs* z4JiPOM7(sCeVyaF*oEy=M9%2rY9v2fiVJe8+@{=&=A5cGCg*n;V+|VyhT~#w=UB7; zfr%{p>W$)$ZmxXhieuvu-=~xx$@)sOW+^(h5zbYx?Lkm9OaI_Fo-cQY=JQX~!&zCo zPgL`mO~?QkDBf8(C7M}ftn=ev#P=!U@#eQ8+y3HW2tAP`h#Y(1&Y;d9;(uLngCJVi zum?-fKuAE|oURKK17U8h*-8k|b*HW17A*Q2-A+iC)TJZY^oV zxOy-D+;?`FVoUcgHdqzC2M-eAoHt{XDH5}!ts4K1$H7kT;r!Rj8x+hqgkxpoa-SroJznphh6&@ zPs##pU($6=;e~jHQx#^SaQNjHH>w`FpWMG9A0%mMz=z!)sYTOf zvIV)+^&c(B%N07%?2PMUcm^l2#^cj$N#rv_B~sk5u87KunLe1kVd|kdTa$f&p^}xy zC%-2oKNhBO@}uD!wC5wt+3@pWwse-p3_W0F#k`>4+#5!QGiSlo8r=AW-?z;)a;Xtb zM0_HwekW#W$?-eXW;RMTjqQ>E->HnS$@kc)HVPl)s!6y>vYNdh*BrcY`n<~fy5sNb zQPIb>F)O1@IxXfnno-Y4F_Gf(Z6|fHz1FHd)wd;Qzuix5+~;n%duPxMw@OmY=Mtqx zh1b$tPvFSaq`RLE+`O-i)>m@Q#Tte>--~0CfiaIPrOovu^$cB54AZ}QtHx{wMVI!8 zM!ZRH-_mkKp=J-_$=n zWhZy~);-@VmsK}H?%tgudW#z8|A3~DJfEK3JWaWm8ufZ+V~$J-3MJ?XI!o|yg}jDk zZSPjswx?KqiSVwwS*6o@UhAz7P8&^Pg`HjDqHCdxwYxdD?K#Otb8p=F(lr7$n<7j( z84KC&NLz-wo-N&6-%yyvCCwyD!KvHm_Lwzl?A#c8{D^P0%5F)Dfm2$E5wE#Fw|e@f z7S(CwjjxYJB{xgmE60ud+2>x#+zegksw!Iz(c&L-eC_mYYl>-(FL&;?d#A@+rfxZQ zT$84`{q2?;u^e#&#>|aeZ=&>Oa(njjCgy`3xaYlg+yx)stz3EP{o3Vi^0%xjP+_T- z3pbSn87d5>Lk)P9G(go9VNR2pHn%nn=gXp;E@+Z@t-t&H;Zw-mC`*&~#>()ZmazY2=@;+=*iqTPU;M{#@h6U8!+MPbF)Lk0E*OH1gtH?AE9A zRF&SSpvuN&{L=$FyeO6$T*_FxdRBSO`Q?{SGNgUk6xBq}dGfGSD|NLl=vNC{^Isvp z(!?}sRZ&r#P?r&#L{(}awbGL`WuSyJU>zW$JV+7waNSK1yshWTF`jEh<+2*8NHr@w z<6D>4{)GSDRr*jC*`>jeI!Wqi=d3ETajjY!+~+4~&pO9nHjU{SFb!Edf2YXlng|0~ zqIYz)4aIU7japh}pmcR6k=-TDo)_%x=aF|j!cw!`XvWzYGvq$gDt=vGg32e(hsTrc zfyZy&Pzy0kEZ~Y?rB0-4dfs6bqa&nebWS7Qc9Lx-k%0ZdeTU2v_q4lkrTEA8mrn^j zs6$Sidbr?lJb78uFDuy<~IEy``PsZNUHM7>U}A3^N9}J;cl0z*l6#x zEPnTzNhFVCT2?>p0Y?g>&T2=;9Hk~gLUGR@ZaPNPL!4ZX2la$2ug4pcld6@M-D?G|vqtvyseYS_o?l3o2i_K#URkJkX=wWDMk zp7KvTkHaOxZW+Us0~kk^qeX(fYl^d+_*LATMYamw)PXWoj;J5p>l^WM;xNtUCbD34WHdsQ%X^A%wl z%=byx;)-lR+^EJazO8edO(HbZfwkiqwF~dt&6@DeElRVoOB8`7_bwI|Hm0GW?Axcq z?_6^`dy;fuIOln=ZUFAsgCs9CEhL@PBoX-C^#?=mBTOEiFQK71cXwCYtn#9kCf{A& z1JiX9H$gnTR1V=vM~bJ2)UZg-;q00%?yZp8Z3@343xsIN2 zV);r(&xd|NZbrM3k4ieHJ&h3#p7INEqDtd@y+Mf{Y@?EMdi02Bf^S{psT5hM#bmxq z*FT;A`jR#@M>C>VV_8H=3h3i&MPUIWrJCioP2>Bk?4ZfaRx z?iiVMzoJS|AjUi?zUm)CG#a05Xfx*VfrAkDJ!s;r*>nk=v1$CWj;K-jgoo#wYf2x& z!eUAYr6zgl_w{cEY_&*b-OpbjJT`u($)bNDsNDcLQ+2;Oo@D-5@Y5R?KtX(Ea!i93 zxocm0La#_a@Nog7Y^zV!*SfxEfAtKR2~J47o<|aVqD@)yadJ#{SG6BR&*S1SOLWF$ zY@dC_V2(77W=WaC7yIEl+Q92{LpA)FrycH;lIpa@wNThH;}czIwS>7=MohhIzDcb_ zDCXs|lF4S7X4ZI1X^zmO0&U*z?eWbvZDt{gS^U9hoqK>f3*2DtTKK-u?CIp%C2qB* zPj)g!>5<{jbdoCePF2XN-m3y-&I_O59KzKbSqR}=UK<&e6_RJ`H>cP%eYfI|&QZJe zD(dFr^Lit4$g(1MV)MJtZwu`?dQREtUQR*IZD_Bp%xtiI<96|WSFT9@VYAa*Fzd~? zjY09Rlt!!1a0kKTp`G&^`bOS2J$F-TqP$OP9z+{>^JTx#WxnZvKg!b{yh@xkti=$& zWY`jANJRJStL36+)Hh`j4uRhZ!avi4?hQwmfgDb&BFEV+fs*8``i1D8!O-o zSNFDz$dJR$I{uU62Gvr}5*b5=cx|2HNrm{bJe3-%Z!|oq*ph^1dbRG&9rTgutd>qT z)f+jCmhG=^yfM$)pz@kN6B{0`oA2ZHu6ePQpx1i?H(2>=%ZIHm-&QAY)!u-UUq05l z;SiPS{Y5vHuwYdF>h8u89hCYO-PeNE&qR-7SBMlGzoglO=AcCh_`KI$q9Y6@6q{BQ zclSD%Zco*GW06qnA5l~2uOyM}9oYf?YS5RerSk3vmAaG|Vfi zv@wIbJx`!Y;L0A?4I-Hs4L6`$fec$tiAoz^(sW8sGnpEk3Ar~kzuUKU zrpWrr!#JaSZ3w(Jxf8h#&^&=uwhcNtd>^VP5q^Hg6iN=&@PtA5eaG9ygr^n^2` z!MPXH$39kq?H+Pk78);bNt5d=h*YxErawn2I{1Bf!h5lND*OWxb)d5c{};-_CUcfY zs)HuH_Lmu?6UpG8WvvD?pfA&JJq#kJQxhJ;!N293MZbBVSABeZ!Z9GjLQ8TU-sdsH zR2$x(f=_c@`>8=}t4i*j7_XOp+r+BIyl-i`Z!vD8;>62OK-hJOaW6&ee`|T1QNt#$ znhV&#+Lq=PZe=k(PN`qEh!m-J7{5=?MIf}iTJSjYps$G~9O=b(#Drk%**`e;(14D$ zii4$0T+9!1^3poGGTJ;c=GRTMT{TVYUHH_l*;+ff0!C^%YiAdfz=n@JMDKs300=DUVv4yeKtWL;FlHpo4?ZXi*m$H(92Ly1 zEv?X)8Pffc1+F51OA>2SNqb9MbAVJ}O|-e4HYTHD;(^5sKnwx-`JMuU=izRLgg`+5 z7bN2U-9BIfLLfLW$8dlJgkVA-Bn(N60)Q0+?mwV{2oMYf!6pO%C+0n-J`5>v_=E{! z>i1YR1`-H~ zK!Okmpc5nvBnSuU3IPB>fIQ}jc@F~;a0w&`w8!*`5<&sb!!pbhiwsPV81uvcfk6tB zkO1Rh+9CkN5BmV>0^?%yFgQR6A#6ejqOebtAOJ5=jue0b2p*0L2gXG}0m7g#?cg8; z@Du`wg?VBV6jKk0B?9J&86P`o7y<$?0E8vRjEMl+V2SbrILs5vKbVo=m>&=fi7;Se zP{TZdgq;`6aEFN)0icL&iJ|Ym834e;E&wrsX@rT!_Nr=7xp*M2ZlF(Bn;+Ah!8M`0IV1uI3l5dF|d3A>{rO2 z(hr{a@dS7QfDcZLePZ%ICG1=uzQ@%0Ct>&kGpm^AziJO=OZ~Cu0sG|lJ5*qIJiM*^ zcXx#T*hVm$?sxf*68H=8AKq}P8FE7AMqVw&mXOb8Xg2lNeY+4oy4PZhp<5uq4K^#hq77o$i(0Q^yhi*BKiW<4RZs$I_E<0B!~ZRWhybYLS9XJE!;yjdXP z8+|Dl)@C`~?LSlwcXXf$w>eS9zJ5XlFRuA%OR5LSN!)ys7*qbYmlH4NCzPf}WR5;I zl^nKCLY%UD<-1lAb)gk8ne&NZnOTak{5;y;T>oT>uBc7ANN-N>{4Lc9cd^qYc`u6v zoySMTgnM498|i-wcwtq!$2kga5HMm1q(+@M0@U2PLD10?90Ai& z+C5zG&zzxhNzwi(grwmc&2BSJVsX2t+8s%5orkXpO4XN2P*i z+T@Ig(tKC_N7A<#%BL^9lFLoFpmhaAC2;91`jQ~Srh;H2@xZNIKb^u&-7ggGJkPHP zyOlR>wp6~4;VzOAogHHwSI{4o8Tc$#Iq4#eeg^Ut7i1-rv5qtf%x(JG=4uCTs6tk7 z7W&KQ1*|DnFr*G~S#;Sx=(OY3h7t=FhTr+7C8tlehI|3_iJJ+4>Z`dGpUd8T$(a;Q zN9EMYFCs#Fjq2o?;7XPR@$%i9OV78`tV`}kJ4v#^+wZ@%Z0wx(IKW9zGdOny{fJsm^g72ArnwQ`$!ieT26I6r;}lX;CB8l}dqw z*~C2c6N&p0lqcp>3cb=MZ-0Y42%{L8V)?w3EEJT)yG6OB-DqAs2GP9Pslk`)X1pZf z`p%#oHcK376r!?RIb~2?r>;Nw9#>KyswM{ccE;VneLWQanWIamUvO<0&%r;I8PQ zR3#!6`1rtT7xfEK!QAK z;&=K0(&YC|50Ke^q>os%5@j2RM0zt_N*?FH=S$_#qTlP=S9(&={?w*DU>P0sazb!t zM|{Mlm8n_)pTQiOK~mDu_? z_}Pxoj25bh+a?A{ks@nzrbUi<`C(s9kyVc#_iMp$6+^0owzT$-dT+X?@2r*?%nb9k z^lrA*G%>R1yjkGx9D5cpt`YZs5VcX3II!bb85ZdEPX5L@1xME+(B=ChFJ5VS?P)gLB9w_82m5RFG53PFirBrQj_KaUZ&n;k7s&!?v~2|5>JGN_n{Qt zARb|gM<{FIrP3`9HHZ=CC-=vq5qCQqqV+#%<+9Ww?Pm6LE6x!o%JBw{cRpd^nDw~6 zdj?O8&bl>RS&sB9zEji4e4kewSMXaOxSc8shay=vLBo}{NbpIhf|o{HY{u@xQUYzd zEFu?Em9y;m+IjD{$6y8>RHkPiq!OH6q}nX-qpFxTh0@7-F z@6GZhiL^@6abU%xPavv8azyUe>mT!c-T2abHhhB7P5WYA3ySokc^t~mtW4ybZR-sA zmq{_dm}J4GFg`z?B-E(cPC1M0nMMEUo%ABE3L1`G=>*a8X@-NlZ@+DS+1)1$KYxq; zh}k(pvw}hXLO6-VD%s-;0PHLt$WPw9GoJY+_k4ZwS8SF@PbG%p1P>-|yLx6>F&B7%V9U)Ss@NNTowH;qbiu65 zU26|M&E<^QNcC+!l6%I(mOH;+RrYsTJ#&sWST!e&a?>VcBbD zs?NpMVV6}m)v1T-+8&dTZ`jI8bGgI}f4JKA;-pQYEv?7J9RlW@H|LZ{wZ|tOsIHcX zj`we}=#QdVj-dS!DkBvB7b>G{jEXY>>xqEP%xv(d!dV@$zFZAG<&RCLo;RN5(GPr& zKZHjm|A|t1IW&Qa1$q;_0KHiQ75OAQp5mxhYj4nI-Uc(dG9|*@BVr`&UafH>)z&r7Hn?mV>hkbK#FmZIuZ7-ymj>!t0&GeyGu#&I0e>prjVC0N|_(om#; z_}?t^<}C?Y=jY))4{*+e=;eWBY5bboumL+U+>Cw?18 zY2Wbp5yU$}1%(2qt3U3;2sKnJun+fm28u&2pPrCt%FnMxj&0wrk4;gH{>0inM8%Zi zk;2aY(4uof-wb3irK#NB^tFyDlg z`BAhMn!F@RcVKh%cMrD7jba+Y~1C z*?xM-zZFy@tpR%BMHYfy&O^pf2{Wa`E|>_KMmEa|98X_F=ZO)XLw#!{4sbbRzO1dl z8Cu8|`2iNmmC1eMt_LmeN*D49XRO#P=$(c}X!rfQ@y{;UJ&Pw6wir0ZOIbH%Mx`Ru z)1askP9G2bz?|pC=3Ok>l_R`hrt+egJ!O1IZ6alh(Xr<0D>K!EPBe4xR(j;nAmknQY?Vs?wO-&`IWQAKOdEeK~?Xg3r7gm>_#Ai=?rf zNx!PWb&WMy&G9ZolQDwJea$q=$0`0gZ@!g>kOzz+V2a|c4 z9SK&WxW;_x6zTbc(PyNAsAEKBv2lsXx+lp50$6y(wqjUa1ozI|${ugf^NQ^rKS-Wu z)4a5$@qquDlwru3zATDQVp%qndgYHlqj+;tT%HiJ1lAsdW#O87D?;lYyiZ`IQYCR^ zB-PJtt+lzT+D`-Fg&@E)D!jYlG>pn5OxS59@G(l;J0 zICCk?k0*YF2%W9{w6J%o;{Ax5pfumw)C2F&J676B#(@J$8@e7>1I3dIQrteqN_o9g z5n=0PCeDS5@bOgLc>2NKcjrQn#;3jR%`K)fvgf-7=|@oTi1Wc;={B+&O7=3q`JhSD z(O2PJQuLR~LX7adlBPWTe2U1})}u!=&+OU5Iu7=;qw#9qY|_PfyKv{TYqrdTS;QsUMOV6Kmp}NSkGHL1`-FQhl{>jTgp?)>~ zIV)q&fx0iR`5R|q4#_M1?#mQ!$(Tz0c>GKr+VK0Ym_uJ4>;%&qta@)=>*V*LEuMeU zLW=VspQ@6ei~iL|B$pM=<@zwi)|@pjd$(-coXd;bj5sh-h;wrY#}%e01Cuz3-i8aZ z5qI`@BrfSm>~;;iwM+_hD>%}=KOXtmQM5h9QAwx+PvXww>C;}HIL6+CbxK20KM}ds zn}Hi!Ghe-*mx?>#UV=kvWY$(vl=C2QWl_}ctEH^7T=>f{Az{{7{Xvx9^x!>a{TTzn z@@YHv@;aTjA6ed|gAKke(J3tlJpFxk;76TK{>CID8{O$d0O&TWg?JngAL5PT&M#z3 z`(;7vLy}&nrAV%+wLUb^aWX--Mb+@ZPTiKL|JbnuHEkT$1#vFmlf;A;4dm{)wbK`_ z&kkQJ&B-w}Ql{k|W(O#L_1jWJlf(qLnF638j}| zO2>`qnYxIN!!beIeFQv5sQrIq^f(fpL~?C@I%QtMUc7-sM`E_jmW1Otb(W_%j~Mm{ zeHM!Ni+NBM`??kA)uB$pWJbh3e22x!d=;5_m(F$Cmt(yVYHokX4iI_Q!cIKBk5)+NR3ldzI_EW{)vz?!aJde$1%h9 z>wpr3^w^#hRq||@#1e^b$8*i~S<2;^8Mk>~V+T(Pg>mw(7}nEFer_2C2h$_lj}M%J zMEa`eOBq+E-IGV$I!S>cISfwmlgtoC^!%p0)qkfUOz0h zFcfer^@khn87Lg02|*$*cg;Zsk$4r(&AjX`?4Xlp!;B@(V6}rdG+&5s9s%nS+8_${ zXVHv^_+cP6J_hMX_Rc;ZSaJIdMAhAEGYuLO-L8H*;`A_~5y_R~-dPhQYL0gAsjs2C zxBOm2svxKvXrC1@#Y*NUov-0eARMN)u`C`tsGIWrhTpo)lzaq`N2p)|0ImK=6xnDI zN_@a{a6>hVeTyAmyw3v4h1_&79kbPF9?}K*H&FS#eB~Fq;1^uY@rGcvB4iego??6)GI!R{vNq7&E!hmDS)zq{MQ8 zUE?|QG$(~`CO#$g44ShyLvnD-g5w@?Z1xj>_6~hdz4>N1u|gW2WQDqQ62@J4hG|F% znd-<*Re$ew8*Jy49U5Odn7RB4s6W?uv>N2xNp_l~C50-CjpU-r@q<$5sJljCXT;k+ zOH;?hlDX8}CIVQtTJ{*hK$kX_;N*15FTs$4L zgo`J|U~vR!N9Z+B*k2qLRMflefZ6-j0NkZDFEBPR6NZUC+Pf|Nj3_qAN@+T4E&fX4 zGxJS$ltXv@oOG~ig7unx|MRl){TvcJw*z%8@)C~ni@gSe_T_Lz@OYb0M3PNa)Z3;8#~7fcZ~#pryFpq z`r9E(04V%3*Q=pO+aM$pd&IXuKCZgelW}Yd(fhYc7Va?}qqyhM^(e=bfK#(B&P-H$-;Z_eCNJ*NQlq*qh@P9E#;dnBrM=Uy=Sz=k8AfFv2x%{fl%9}_OE zY`mOo`i?_Otwj7RdriF}io8Hw@R*vFG6^g_S7|C?@rk);jR&gaB|%1X*Bd=MxD5Ss z(V(3CWg68TuuI&71bTDQu42tMs$I|gXp-q(I=}6!Epg7GVy@|AvU5;CI3l0D zhzDn^r&W38MQtU%ST2`U-*2N7$-N^BasA;gnFz;--O$D3TDxi=4B|hj_?N{B?3~tW6!8yN+%l4VtuL^$F9Rv3 z|4QT_dyXI=w7g*aNuP00?U}=90uGJ{}E9>4JGGxz=6z# z*zg`~i*QzS5nb=`_@FvAPF9Z8cwJyFDD>h@Jy1=~$kUpLUWveWWktP#eF_s@!XZ7; zVYRV*fj4A$X9b>$KJ{BI`pnF{UEJGULD&HvcPHcd`ZBLg_2^N8=KaWrVL6bfua&km%vZCo zu(`%w@PEn;emzsq)OCNPVKO`H$?aTR0iq+YIYO_20SC=LVuR2?*#nXhaMkn%`Pfga zj)N3pcAHKWj8o$~p*7+~^tIf2f&;axmeJdl39l#xXWh*MpAvN!QLbnNDLiP@W*uHl zftJ-EERS_*+;5%DoptK=yU=z&|77QOb$Qub7gH3oE)9E@Of$6o=>kLIi;oX2uG~tr zZjSN@Naq>h>(nLCIP*qm<<4f7C*PyvA&*lM?+`O-htl4X{CrP7!u#uKEnU#S3ZtT9 zc&K`ekDKo+1rpQ^Ubjmu>#*|zR12@p$i9~J<%ms2wb!H9-){An-`gJx7i<#dUc92O zXF4c;W+lFZ7rN-ByN>=c(&>K^zM8!|KTJXiH;;(zAS`)>3d&Bor0r279e3)vKqu=L z(#$q`qKiAC)J;(wR7~+*Z%I#O-M-446P}Ko#CxB8!H$NIbh@I=sLq$>BArd_CT@Z9 zg74x%+2^V!qoXWLD?S}ZOwJJ+74&a5aK1~5&iIli8v1zF11$32nkMlC&mN4Pct&{q zRG0)>fBIhdT@bHM<6KOURoh#bfZtn!;9ESScDl0?S|CMJbwrliXRmeUF{fmr-nfy6 zBo@VFyWw1~7uw$m&|B7wn!hApHY^iMk?bxuc*;RVr#x)JZb^C}1fwlEQ#ewGl#2K#F9MS7|Lf7v8COlg8snv%WhZ{c(k$ie2>W?`Q*}_IK6h zZ*D2YqFMmW?S%B|Ht~T(?ig3rO||^%!5GC?!X*?b zaN=Cvq}i7^8xHnR%Rc{nzv+4SfHrfBK~QBpYs#8N+ymA*O>wUj{`C6SqrKPR^NE#w zE*LHtGy!L?$)uGxU7Ghh5>NXc$JG{M99cKO+dojLZ_Nhbpdc|>_oqgQ9 z@>+?%(+99t{5$JMyD>|tT)OlhH)hOA_{V)2uxGt+?PS$ z-Isss`mp^Uj46M1TLwPn@sHaw6r<7tA0GL~NBMWRW#V7lmNCu@E{-Oqz&A_)x#o!Z zpop!x1^Ro|8S51c2d+IF&7G|s%s_xABG}8^*#QK0um^my(e6MFZRKokP7LU@U<+$E zbHJMw@HGbf9RV+65Lgx zuo(ypbZS8i2ENQ<1p-?Gr9d~fATVHcvIm9)zyKVqf$hl|<6{j5pmYVyPQa1G6NF`S zJ6Bt@wWF;k;B^SN6aM3VjoIXW{;COfqJB6j{L2VGP32GS*T2mhfDzX18i?}1R{O&Y z{@>zu4Tt{hcKu@$zw(dU><{st-%yfPrH3p0?D1-=V)D7WN5 zP|WEKD3b@m0Sr%0YU)|bJsw? zNgrEg2810P=FEX<#^i{>79b#O0{a##H2?;{{;&ih4Gwcw7>R@}vjzdW88&AN0yr0& zvjf4g!w`f2+4SMx$FK*YJ`O86fPhOVY;}yZ0%R99?+5}sA%Gn5PvQ=1%L6dN3cy3D zhqZrUiVh_o5W|3N?F{hOVa^4Fz~V>zU4Q}+P1qV}5D?FCn8Qd$AmrowFF;rVms8j` zu7Gp=Va^j+IhY*rp>GRm#RM-TOKWLuFuT zNuUwf)C44e83QZ{{a;oC#OlA*VE@Z%5Pm>J0~Ik(aq*u3#{l}B@&3gILzv$RG4%Pt zqH+R!KLwH5S=+n1fB-M4|I!3A6+apr%J0MPxPg_d&9Jdphm^sDbUC;JVNO_FzlUEU ze{L)de5?!;d3GKM4Rf@0K-*eh1G&NYh4|sTAS*Q5(M1>xw)^K5zk{V zz_-tSZsllZ0lH>lY6Fb&y$7HswtH&_dl}3(;?BzmLm|M07YZr_Sc%~RFg}PNCj`RD z^%KoJoXstWF##sTKzz&hKR}^E0L#Y$^iLTkIt}v&vj1KtfI?vw(=TO63@`msCJ02S z{aPl3SuVfS6GHyppO6s1y#J{O`~bAQU&~;a?cu-mfOwqW+d={5?zeis;Q@?*(bxXl z9}wB{+qei6pmY9GPY?kpp1+hKFrh8KmLYz}7nmhriDLf#rL!Okvkm-G2FGlI|1Cp7 zfv}kWmO%tDVK%>(3H}a?5ajoILIS_zCItJ94(cU%+-IQacjhITfwwl#OgexFm*+6y>B0;fd{ p2M5f?jg^*)_7)Bx?8c3ij4o&sXEb)JgaT^?_{=AVtePD0{{e>eiyr_0 literal 0 HcmV?d00001