From 282f1ae290775ceab0c6b257f0d2dc7485fa13e0 Mon Sep 17 00:00:00 2001 From: Greg Caporaso Date: Fri, 12 Apr 2019 12:06:59 -0700 Subject: [PATCH] added dynamic tree visualization (#314) --- book/applications/biological-diversity.md | 98 ++++++++++++++-------- book/applications/images/pd_calc_tree.png | Bin 17889 -> 0 bytes 2 files changed, 63 insertions(+), 35 deletions(-) delete mode 100644 book/applications/images/pd_calc_tree.png diff --git a/book/applications/biological-diversity.md b/book/applications/biological-diversity.md index b7ff615..567d43f 100644 --- a/book/applications/biological-diversity.md +++ b/book/applications/biological-diversity.md @@ -144,7 +144,34 @@ We could compute this in python as follows: Imagine that we have the same table, but some additional information about the OTUs in the table. Specifically, we've computed the following phylogenetic tree. And, for the sake of illustration, imagine that we've also assigned taxonomy to each of the OTUs and found that our samples contain representatives from the archaea, bacteria, and eukaryotes (their labels begin with `A`, `B`, and `E`, respectively). - +First, let's define a phylogenetic tree using the Newick format (which is described [here](http://evolution.genetics.washington.edu/phylip/newicktree.html), and more formally defined [here](http://evolution.genetics.washington.edu/phylip/newick_doc.html)). We'll then load that up using [scikit-bio](http://scikit-bio.org)'s [TreeNode](http://scikit-bio.org/generated/skbio.core.tree.TreeNode.html#skbio.core.tree.TreeNode) object, and visualize it with [ete3](http://etetoolkit.org). + +```python +>>> import ete3 +>>> ts = ete3.TreeStyle() +>>> ts.show_leaf_name = True +>>> ts.scale = 250 +>>> ts.branch_vertical_margin = 15 +>>> ts.show_branch_length = True +``` + +```python +>>> from io import StringIO +>>> newick_tree = StringIO('((B1:0.2,B2:0.3):0.3,((B3:0.5,B4:0.3):0.2,B5:0.9):0.3,' +... '((A1:0.2,A2:0.3):0.3,' +... ' (E1:0.3,E2:0.4):0.7):0.55);') +... +>>> from skbio.tree import TreeNode +... +>>> tree = TreeNode.read(newick_tree) +>>> tree = tree.root_at_midpoint() +``` + +```python +>>> t = ete3.Tree.from_skbio(tree, map_attributes=["value"]) +>>> t.render("%%inline", tree_style=ts) + +``` Pairing this with the table we defined above (displayed again in the cell below), given what you now know about these OTUs, which would you consider the most diverse? Are you happy with the $\alpha$ diversity conclusion that you obtained when computing the number of observed OTUs in each sample? @@ -166,18 +193,6 @@ Phylogenetic Diversity (PD) is a metric that was developed by Dan Faith in the e PD is relatively simple to calculate. It is computed simply as the sum of the branch length in a phylogenetic tree that is "covered" or represented in a given sample. Let's look at an example to see how this works. -First, let's define a phylogenetic tree using the Newick format (which is described [here](http://evolution.genetics.washington.edu/phylip/newicktree.html), and more formally defined [here](http://evolution.genetics.washington.edu/phylip/newick_doc.html)). We'll then load that up using [scikit-bio](http://scikit-bio.org)'s [TreeNode](http://scikit-bio.org/generated/skbio.core.tree.TreeNode.html#skbio.core.tree.TreeNode) object. - -```python ->>> from io import StringIO ->>> newick_tree = StringIO('(((B1:0.2,B2:0.3):0.3,((B3:0.5,B4:0.3):0.2,B5:0.9):0.3):0.35,(((A1:0.2,A2:0.3):0.3,(E1:0.3,E2:0.4):0.7):0.2):0.05)root;') -... ->>> from skbio.tree import TreeNode -... ->>> tree = TreeNode.read(newick_tree) ->>> tree = tree.root_at_midpoint() -``` - I'll now define a couple of functions that we'll use to compute PD. ```python @@ -213,29 +228,29 @@ And then apply those to compute the PD of our three samples. For each computatio ```python >>> pd_A = phylogenetic_diversity(tree, table2, 'A', verbose=True) >>> print(pd_A) -B1 0.2 0.3 0.25 -B2 0.3 -B3 0.5 0.2 0.3 -2.05 +B1 0.2 0.3 0.2250000000000001 +B2 0.3 +B3 0.5 0.2 0.3 +2.0250000000000004 ``` ```python >>> pd_B = phylogenetic_diversity(tree, table2, 'B', verbose=True) >>> print(pd_B) -B1 0.2 0.3 0.25 -B2 0.3 -B3 0.5 0.2 0.3 -B4 0.3 -2.35 +B1 0.2 0.3 0.2250000000000001 +B2 0.3 +B3 0.5 0.2 0.3 +B4 0.3 +2.325 ``` ```python >>> pd_C = phylogenetic_diversity(tree, table2, 'C', verbose=True) >>> print(pd_C) -B1 0.2 0.3 0.25 -A1 0.2 0.3 0.2 0.05 0.1 -E2 0.4 0.7 -2.7 +B1 0.2 0.3 0.2250000000000001 +A1 0.2 0.3 0.32499999999999996 +E2 0.4 0.7 +2.65 ``` How does this result compare to what we observed above with the Observed OTUs metric? Based on your knowledge of biology, which do you think is a better representation of the relative diversities of these samples? @@ -596,10 +611,7 @@ First, let's look at the analysis presented in panels E and F. Instead of genera >>> ax.set_xticklabels(['same body habitat', 'different body habitat']) >>> ax.set_ylabel('Unweighted UniFrac Distance') >>> _ = ax.set_ylim(0.0, 1.0) -/Users/caporaso/miniconda3/envs/iab/lib/python3.4/site-packages/matplotlib/__init__.py:872: UserWarning: axes.color_cycle is deprecated and replaced with axes.prop_cycle; please use the latter. - warnings.warn(self.msg_depr % (key, alt_key)) -/Users/caporaso/miniconda3/envs/iab/lib/python3.4/site-packages/matplotlib/__init__.py:892: UserWarning: axes.color_cycle is deprecated and replaced with axes.prop_cycle; please use the latter. - warnings.warn(self.msg_depr % (key, alt_key)) +
``` ```python @@ -610,7 +622,7 @@ test statistic name R sample size 6 number of groups 3 test statistic 1 -p-value 0.054 +p-value 0.065 number of permutations 999 Name: ANOSIM results, dtype: object ``` @@ -630,8 +642,7 @@ If we run through these same steps, but base our analysis on a different metadat >>> ax.set_xticklabels(['same person', 'different person']) >>> ax.set_ylabel('Unweighted UniFrac Distance') >>> _ = ax.set_ylim(0.0, 1.0) -/Users/caporaso/miniconda3/envs/iab/lib/python3.4/site-packages/matplotlib/__init__.py:892: UserWarning: axes.color_cycle is deprecated and replaced with axes.prop_cycle; please use the latter. - warnings.warn(self.msg_depr % (key, alt_key)) +
``` ```python @@ -641,7 +652,7 @@ test statistic name R sample size 6 number of groups 2 test statistic -0.333333 -p-value 0.889 +p-value 0.869 number of permutations 999 Name: ANOSIM results, dtype: object ``` @@ -659,6 +670,7 @@ Next, let's look at a hierarchical clustering analysis, similar to that presente >>> lm = average(human_microbiome_dm.condensed_form()) >>> d = dendrogram(lm, labels=human_microbiome_dm.ids, orientation='right', ... link_color_func=lambda x: 'black') +
``` Again, we can see how the data really only becomes interpretable in the context of metadata: @@ -667,12 +679,14 @@ Again, we can see how the data really only becomes interpretable in the context >>> labels = [human_microbiome_sample_md['body site'][sid] for sid in sample_ids] >>> d = dendrogram(lm, labels=labels, orientation='right', ... link_color_func=lambda x: 'black') +
``` ```python >>> labels = [human_microbiome_sample_md['individual'][sid] for sid in sample_ids] >>> d = dendrogram(lm, labels=labels, orientation='right', ... link_color_func=lambda x: 'black') +
``` ### Ordination @@ -762,6 +776,7 @@ F 0.737954545455 ```python >>> from pylab import scatter >>> ord_plot = scatter(a1_values, a2_values, s=40) +
``` And again, let's look at how including metadata helps us to interpret our results. @@ -772,6 +787,7 @@ First, we'll color the points by the body habitat that they're derived from: >>> colors = {'tongue': 'red', 'gut':'yellow', 'skin':'blue'} >>> c = [colors[human_microbiome_sample_md['body site'][e]] for e in human_microbiome_dm.ids] >>> ord_plot = scatter(a1_values, a2_values, s=40, c=c) +
``` And next we'll color the samples by the person that they're derived from. Notice that this plot and the one above are identical except for coloring. Think about how the colors (and therefore the sample metadata) help you to interpret these plots. @@ -780,6 +796,7 @@ And next we'll color the samples by the person that they're derived from. Notice >>> person_colors = {'subject 1': 'red', 'subject 2':'yellow'} >>> person_c = [person_colors[human_microbiome_sample_md['individual'][e]] for e in human_microbiome_dm.ids] >>> ord_plot = scatter(a1_values, a2_values, s=40, c=person_c) +
``` #### Determining the most important axes in polar ordination @@ -822,14 +839,17 @@ So why do we care about axes being uncorrelated? And why do we care about explai ```python >>> ord_plot = scatter(data[0][4], data[1][4], s=40, c=c) +
``` ```python >>> ord_plot = scatter(data[0][4], data[13][4], s=40, c=c) +
``` ```python >>> ord_plot = scatter(data[0][4], data[14][4], s=40, c=c) +
``` #### Interpreting ordination plots @@ -848,10 +868,12 @@ One thing that you may have notices as you computed the polar ordination above i ```python >>> ord_plot = scatter(a1_values, a2_values, s=40, c=c) +
``` ```python >>> ord_plot = scatter(alt_a1_values, a2_values, s=40, c=c) +
``` Some other important features: @@ -950,6 +972,7 @@ Next we'll load our distance matrix. This is similar to ``human_microbiome_dm_da ```python >>> from iab.data import lauber_soil_unweighted_unifrac_dm >>> _ = lauber_soil_unweighted_unifrac_dm.plot(cmap='Greens') +
``` Does this visualization help you to interpret the results? Probably not. Generally we'll need to apply some approaches that will help us with interpretation. Let's use ordination here. We'll run Principal Coordinates Analysis on our ``DistanceMatrix`` object. This gives us a matrix of coordinate values for each sample, which we can then plot. We can use ``scikit-bio``'s implementation of PCoA as follows: @@ -964,6 +987,7 @@ What does the following ordination plot tell you about the relationship between ```python >>> _ = lauber_soil_unweighted_unifrac_pc.plot(lauber_soil_sample_md, 'Latitude', cmap='Greens', title="Samples colored by Latitude", axis_labels=('PC1', 'PC2', 'PC3')) +
``` If the answer to the above question is that there doesn't seem to be much association, you're on the right track. We can quantify this, for example, by testing for correlation between pH and value on PC 1. @@ -982,6 +1006,7 @@ In the next plot, we'll color the points by the pH of the soil sample they repre ```python >>> _ = lauber_soil_unweighted_unifrac_pc.plot(lauber_soil_sample_md, 'pH', cmap='Greens', title="Samples colored by pH", axis_labels=('PC1', 'PC2', 'PC3')) +
``` ```python @@ -1010,6 +1035,7 @@ Imagine you ran three different beta diversity metrics on your BIOM table: unwei >>> _ = lauber_soil_unweighted_unifrac_pc.plot(lauber_soil_sample_md, 'pH', cmap='Greens', ... title="Unweighted UniFrac, samples colored by pH", ... axis_labels=('PC1', 'PC2', 'PC3')) +
``` ```python @@ -1020,6 +1046,7 @@ Imagine you ran three different beta diversity metrics on your BIOM table: unwei >>> _ = lauber_soil_bray_curtis_pcoa.plot(lauber_soil_sample_md, 'pH', cmap='Greens', ... title="Bray-Curtis, samples colored by pH", ... axis_labels=('PC1', 'PC2', 'PC3')) +
``` ```python @@ -1030,8 +1057,9 @@ Imagine you ran three different beta diversity metrics on your BIOM table: unwei >>> _ = lauber_soil_weighted_unifrac_pcoa.plot(lauber_soil_sample_md, 'pH', cmap='Greens', ... title="Weighted UniFrac, samples colored by pH", ... axis_labels=('PC1', 'PC2', 'PC3')) -/Users/caporaso/miniconda3/envs/iab/lib/python3.4/site-packages/skbio/stats/ordination/_principal_coordinate_analysis.py:102: RuntimeWarning: The result contains negative eigenvalues. Please compare their magnitude with the magnitude of some of the largest positive eigenvalues. If the negative ones are smaller, it's probably safe to ignore them, but if they are large in magnitude, the results won't be useful. See the Notes section for more details. The smallest eigenvalue is -0.010291669756329344 and the largest is 3.8374200744108204. +/Users/gregcaporaso/miniconda3/envs/iab/lib/python3.5/site-packages/skbio/stats/ordination/_principal_coordinate_analysis.py:102: RuntimeWarning: The result contains negative eigenvalues. Please compare their magnitude with the magnitude of some of the largest positive eigenvalues. If the negative ones are smaller, it's probably safe to ignore them, but if they are large in magnitude, the results won't be useful. See the Notes section for more details. The smallest eigenvalue is -0.010291669756329357 and the largest is 3.8374200744108204. RuntimeWarning +
``` Specifically, what we want to ask when comparing these results is **given a pair of ordination plots, is their shape (in two or three dimensions) the same?** The reason we care is that we want to know, **given a pair of ordination plots, would we derive the same biological conclusions regardless of which plot we look at?** diff --git a/book/applications/images/pd_calc_tree.png b/book/applications/images/pd_calc_tree.png deleted file mode 100644 index dbefdfea115e4382273cf5f804691b9f9d35c2a8..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 17889 zcmb8X2{@GB8$UcLkt{9NY%OFkp|WJ(zABWZED|@OMpJ()4e*gEquJ?K`d1lUbpZnbR=U&b^kDlmgU1Oj-MF)ey7}T%-bqfYN zC_%9Ci=}Q_EBTOIgonXmJ$de2y5%i?=;k{aNwn_s7Q{&=X<71yizN z!Y%9xRf(sM{M1aa>V5zD!vWLdD-w-``0AQMn?v{1vhTJ2{OjR?1AcAP&mLaSeDv%@ zB6XDO#a8AtQR!^>t<1p1Vy}dl)V5UN$P-)1ms_gLh`qa!RMvtrQ9)7pwCb355?xL$ zavc>6Yy6~XDGZ@Js0UtZ-sRe|t-GndJiR@*Y{xk*?3DOBlTlwktUz}-{CC!u@b~5p z>|OpXXvwEf-kYCNByUeeiEtX#le_P9I2v;9!n4ZGcN@Je9%{SiN(G6#~-*T^| z{-efqch2c#M1tnYwF*-tQi>R1W;Y1$=nW^B@GAJfxR|znKbu{Gw+vBc-eueJ*IYWG zTI^yVdQeRADWFjQn;~{xBq@5$+^Qup(Yw|>aC;)6q?Eqb@knRi=Wx`513Two&QIYV z0DoQE$W87P{Y@-cGhx~4U_;Gbb5&{2gQpBFgU*RtUUt=haeT+x$WXbtLt@#OSW@9U z?l$N2h~$@gdh-$cL9r6TH6!T`WdHbK%cl$>VD=vn6YdE4?W8ZA)Kdw{v*8;&j83#(FNVp57{!;5}| zU%V(u4^i`q3DJD$_>#MY4dtJd&jI)_^hXky^*?>4r(4AtfW!?bUQ){$vD^5R5qrL}N!k%zI(wJ+_msnn!UHEc!@R3?_mwxotk2n+eJ|@? z)(;9x18;pYUPr6w4#8*EMC>?>Q0i$&mKLq$g2Mt9p4joAAK{DQcO<6 z;itr+9MJAUF`Y3U4(I_Ptx6g)gVh^KGeoxj(X-x$JK=Mg<)6CDvk->M_f>z{GS9cp z)=OgTDs~ES{TSH+HIE@z=eR0DmO$hd=it@fRiX~eoYdf_^{RB$-VdKT)k3qUj$oF@ zWpWHD;-$Ly;>3Yzd)~U;4aQHxtA^h9Haj_IJDr#*!&N9>iEQyuosG9m3 zeU39MN${qswJ-N1Uzu)KQ8wP2|1SR{9xwsnyyUgy3l8Y2g1+SlU*5)c-}XXg`i_8M zdr|t5e5-qzU*g6uB9QfrV|psN(d_YtIAExqhD2)?8a6IXDJgREJ9bw34ykx!6vypefo&@^ral8XDNvrO{$>N2?@zBOM{?+c~8^ ztiTK^8ec`3xvdXll{h6pnM5 z$v6M-ZrlFHaQwQ6YJfbdBFIn*5l-E*st}GNE@P_v-{7`__O=@A`H#l5zFUb5w={aH zQsC{!)TS=*zLz1&T$&JmohAdZLqrDtFtWWv-|a@XPN(#yg&%xT5?iSC>@H+T1LhceGDaE7R>AF=-!4@a0 z7-D=@U=^>_)N(&cl7{w?IV$^}^WmHtOr+y8^I&fo1S)k{h@G|Nk@-ag?*-`})}_f_ z4JQFhgds_y)n@KN(cFxFr!Y-Jk|ZLTC;nN1xjhM^|5>!Sh`EI-yw%_Xq^2mae7QBI zCzBy@*bZ3s#ZY$mXUQwK3zMwB#9lO?3KlIywF3y-xl7x>1G z)!@deWbS`^=d#&LmYkd1yz;bU;DZz|#1baFp~ZY~cXDcaEKW-Ch^VIX6`)Z3E3e6~ z&#MZ&FFHjaI-0)#1}NM^Qh7JfaEHe9A--Rr&I$!nS7h!Uqe`Ge_% zE6IkPZM8zG_E)y>WEDv^34bPU$-uJqQj;+=gAx0nq}%JmJ6Je1{ZV!x z!{mks?)d8x_B_to93{}oy~8K&A$2ivlGh?RS+7VV=J59sMoq{TNHG&vaVGa$nlEcw zmPu#hjeIl=?Gf za}8(@LY+|x0;-ir@2@+6M?xqm%(xHWyzktr!UW|aHme8I-lOYi5(7vZ_tqz)Z@&sjE0 z_bc+Ci@$1;`Nwn@73)ii#t9)c_Wm5`lEpUOAi3c;Z{?ECz3CP&Mu87(7aXtkeh*d>uV*}~ClWdq5@OyyshsqN;M3k=0#;3aO`Dv&7?9BMgD4KvGv&VsM>FxEywFg@Mj{ zWbjNK$Egqto$7@y#=)elhZo@^<`IR;6Hhlpox7vBy1SRVIu7IF2on*!W8YRgkmS{0 z$TcUtt*y>arDZpcd0)SUl{|?NMM-#i7MLJ}3R3I7NI>B*HB{&_$KjtpTrx$d_GlL! z3u&HH9xNJjV=XRX1mDlR^?EyrHV$cKLvKdu4Cp5o2{?_|hXuYby+u!#d&{_c%X}Z5b zg@Jdm;RZoN^qQUl1PfUrkamwQJskkfJ)QsS11n(1P0p_!8ylVP4LHFNsVp6z`Xe`X z^QiDZIOoCo~Zu&9U4z1KyLHKVJM-1{M@} z09iO&=oq&zsb$Z^%JCJA?dWopBFelj!XqX_RQW}Sw(!+26DQfxs!8t8ox&%l=by%! zkNFoCA?dMNo;&OHxZVqNW@$^&m2;aPJ6=p{!EaAX|5kODWpBsR{wlzl?+DP~6v=IS z>;?t>Ql{IL;eEDi19J{yM|FG6m^3;xYGmOqPFSNrzZ@&0Py=D&`O$%Nm(7q90_{C> zohQORl>}}G3s2nBdEF@P5T;zo8K3P81RPz{*ccCIXSU>_KqP-`*VLD|wU8%I*T8oP zemE>vcxz?Hm{(T`ND3u2X`oo%7PaxtA*kY&<=a1y1p^Bo%@B+%StwH<^hrd7*qv^Au&bGk!+;RQ8m;Q1cxC+U-c&@dYP+ha&@|`<0 zPon%dn2nI{oGiN)v72~*kRJQ0n#)5)p=K!XmT3gMogT}^j$Uf(M9##Uf4c-P(K17u zOIC9~uTN>rS4Odbao#**h&Re8)1j;av8D6FmxfU_RQ5Xfy?*=mNzaq|JQ%RTl|uF* zkCajNVCAdX$i1Bb>$=^fpf%_E`}A}t_F?x2egA$uf>wo&@Yo@dp800N8Fl>}S~!_Q z+u>7%c7KbS`QW)sdM0%(74IMPw&lZn39r&j0EkZh!zB0!W*bT`#RGjKGZ#d6}KVB zuEAqMFcsV!OR4OhxfzdCylY7X$+%a#APAX3MavROcUlj0u;vxgvN$R9Q0KUi3bj%m zEekB;kq_gxc>-Wr<;=k@lr&z#r$;=@jtt z)ios#{79}oWhZT%{aW7j+EHatX%CvD_=?>PFGZ;88+H<1-v>6(f$M#c1 zjuEpK(ZwOS-ax(#suk&s81c?lmWD98 zrpqZi8w2;u#)6!3@*x#Ll7#{0V1HCbM^PM9`_P+9NO6zQ9UZSoo>Hs~hke!aNwNw- z{N-9uQ;l3+T`>qPSh8$i3ta2DprDl3)aU#IE<@OMss>|I%WD$U`n2;cf;d+5ht^}l z#DM5ZZcB4=A?X5q{iT{Z9leRI*hm#1S0AmsUuQa_My>TvA_nT(Mkmy&Pd!kf20Zzw zgUC2K%orE=e$QXQCan7MV30yb7@dk5Ftb!{p*9S8e*M}tkD1Y)Wm~(`IMI0DD z=m*F^k_5k3`gGp1O&wEnwaM-2QR7WI@FQ8Qw(jl&$G@#`?n902;FC2UClR_Mlu@bT zc53GJHeIZv5!c+SKKoW;{C@T41$YShQAn@i9h~3LgV_M56Zp<1Nes3Qp~_x4DK=;- zTW%brjb~2p*$7|ntuJXpaM?A)@g8sLRo^e*e@By=+)XkhR2Z-5j zd!Yz`9>khZ=w`}HrpG&GrWl{=6_U#!UqwHv$x&R~Uy|JRPfpRgXxOgE2fY~ifDt9Dd@t2eJ?LWo)QkSTh@esoVI z4_2S4W5kkFhqpBUdUI<))Q;-DXOF)_Md=GmquHCR=Cw`|y1k4oW|httRj3=2EiVpr znLGQ>2wPt7TXHyW+ZH+LuYwyn!ZKD=N~^wHin@c;R&jGVLO$F5di=Blm31w=NEk32 zETGp)uA;H*VOk2&#A}CFxxXh=cse@%PS#P*YjUfRBlw#ZeMjCkQr5O~-)KB*2qgUT zvYdx^mxUZHI5N>WY|DIjZR@wa3Ixa;CSHY+I`eoON#t6_>gQ4EBputgeaBwidbgE4NFjg<$?v~} z?jP8jfg%qu7b>8qdg`Fvwji#?Ju*?MY`1H9_xUr(Sp0>79i`*|-vyYUxUYXdux+(m z9+d-euep%bVtNC^3)=jzq##oQ7+5X0lcYtxm#=e|{M_sk=P#-y_bFe5+< ztB1lm@ZkF{;7^Q)0>M849~5vG3TDnk`gOMSt1ZuJ#n>{3CCOG+bFc%EH}~bA%|^Ws zs80^hMsfwZde7|?*03CeqPuqOYX|4oEY`!HZgm|#Ga9NH(178^RwUhku+n6U+YvoF zhj-=nwPh~-uxXgq$#7sn2yMxM&1V$?D3#7|^Y`xRaj0k@IxtI(yGq+E^B!e^rS|B} zPveubw1yUC3#MB(+Oi4Hah8oi2tC-g#UFnUMONSo@A|}!N!S0rfB;zw@ZbmYjo7~w ze3i9fEF*ZjrxdjJ&I7m23wC!IBZt)YHu;JdY^AVOjKV2DT+}Z+R&{kaVoe!r9SvAd>-gO9J$82Rpa9PC z(KL3H1Pg9dhnzZ-#sn ziik4$b!{M5gBZp_P0MT+?Zg$9?dQ5oxY@K36~}UjiaIq*0o6L@TGB5kpre%DRQ4sJ zXxh*o<;KpVKg%fE+<*Ec{It%`SI%9YM{OBt@~k`C<0GtUlnw+pM-f64w?nFkyj8W; z-dXuZ{9ycP9fe*iYyaiR&BOd_%B-2$qWhJr&j?NGaa9NU=xUBe1}iD2!>Fi(sZQ%0 zqI>vvffl`0c?n2mhcrDkv82E{)wNT@R~}3BzQ_!RQQ_{$GwF`l^AwBt0;^iah;wgN zSvyD##A{Fk5~#oQLw`-jEcP#M{?;;De?8rsDnhB3)9LkRJKHhWpfp(WZr0!6mC}Jj zSk!qenJRS3iz~d-(P_a7so~|qLD>mHWu!xB;LXf-t)6=_N~DvO|(w><8wuH zV8M^TU^*&lrF2^6{Q@>NaA!`SLjZvx)?pU7>1mn8nIRd#4n-;fwbJevpmE3G(72UI zDoelLUjha0203sC@Ea+nc~*IUNRF%Ipk+RUJ!geF^$GYDqO-_GP$Fl8)0jUZD)2Ku z9n^UL@Z1u~LNa>!)78qlb;Z*<$KcV}Shz_%6AiX9(Hr{iW`d1<72ln8i7$NCFYxIkPm{h%T^gtQ)qp2!WdFky1C00w%+(+xv(# zj1CIgr3PykmypKlkdC_9Lp^P7$hgklC`sz>7(D$kWe#5um-P*jL~5mU-f6q&foN8N zkSvU?=-}5#s^C>GQpw&1ra4-XY(ogeImtdL^|GlN?wpG1UmB2UpEfjqE5%96?4u&@ ztdd=f)rq%EjFw&KmhDWKJa5A_PY+;e?ta$F{t9N`4tb`!L3V`T1!NCE0}G~lsB{24 zo9vIGBbIxsF3!C7z&AXs5$y(b214O?NhfV;Bp~0kwZ9E2WQYjFX3hjw%udyBo1oO`MEq3~3ls;RQPe zhXKVyT&ShPfRZWYpBpy(Nl`VQ%spoagVp0EzbZ!4?x@o;|8$Mj;###h$=`8a&&g3i zJLwuoQq^ot&*(h7mff?I*VaRNk1Q%qiZ3nCKL&4dKgchC9Nsslc2c~hCK3RK=Oy>Z zu80ghdPioxZEdNg6`oCM#;T;lnC`M9aW5OOd-Rb~^qe1AuDSZV_W0ME7#V^|1@$@N z3;Hm<#&LoDVM+(as@>wHmgZdAnHaBWeoxyv+OLrg19slK`_1sROS?tzSsj={W=6Jc zt!c~^F7hL%cx$k>CRz4(vY7n$azWB`=xR1>>4>N@P_&P8BM~nvJi&ho+m1vXVGgFE z3X5hy*D;(2ngj*vLv*Ksu>ktkL`8q`eiSyj6Qp5T3;#4JOuE}v$1fdc@?LG-gm)}5eXTM&w>UWf zFE~Q5n}1vPg6E6!N<}5eK^(%~=Bb^s?$LhTy`{O;8-5k=h5B*ddv|0H)&p!5K(&=y z@i_S0o?VwC>yKHfuT2ZDsW+13zJG?;N!1%8#*J;xS&;lo1AHVl?G>OkYty5%oNFLY z(%AE`C`gI4(!vyf*$jUeJJDMwSAVt1b8Q5NB^sy0{4G@9zD9bgWgI(CMg8+W)tlml z(Ne;%d{kgX(Aet}a9FgAD|*296>hL$=3i7(Q?R{nVaNZ$`hRO7y93fwg5dwRH8w3U zq0{zuXNWk+Z0wg{=^3E)=Y*64b6DV*qXG}M5B`^}?1SGs_x&8yf|Y_J;?xu$5DHOp zPHh=*XG}cTp&QI^3|IlceprkgY3R^j^NF|4Ko;Tyy5F&cRm)bC0#%-4u#AuJH^3jX z^w8%HXh?DhT7ptCXNjc?BNamZt^kG+LVFn%?vsD--JnK>IJ=4;Gz<`!d0F6Pzo+XH zoz{GM9m8&^LipcKRFXF7#czLTIsoiun;u-I`baV0q|-tG;#Y0bzj<=Awq6ISR!;Z# zO&>wH#)PYJ063@tjQtZkAm#sw9d(kN56(p1?R=K)o6r$bg6=%7hQjz6Tw{gsOocO% z`lf)G&nlSv-eJVPi&*z`+*DS`4EV2m%a4vwn(+mZFV>?*rv{g%Qc7nfI^W2f4%twsqy7`znVLC96*!>jxDN8`Exdc2 zbgioDdT~;K;GoQ3#dg4a)Y56u)*y=vDgs$#fFRD{Snti-Z=IJPq2A{?h*OqPf7`Qx z5JdU9^7O!{L|)PxiDZdxfzFwB{Cz96GtiPI=9Fq2`N~EoI&p!Rubo?H0~Uxhb1EI0 z&q(U!6(=61={BNov8XR290&wmwG1vH%Ke*Hp_ADwbm-2Dw&0F`kWY~=dN9v@c@dh9 zn#Td2tHA4R(Rm2S?&lO0&c`-F2kFS$R@FTm0Xfk>AqUEq{uF?0fGPf`urF|k+J1fl z6am=%-wykstG->yh&5}PHK%B7yUTF?1>VWk7_l(@N>|{hQKeC@7*S(Ojs*@YIjb{& zUQBXToURYIWLTSYdNAx$=hle-H$WKs9@EhAa zhp^|BJ{k5;bE?$EN_CrVmc5w145e632l+K<>PmZ=mU{Zl18kdrxmE4)ed`cQ{Njup z;5iCHwA+XhbaMUP;SoR6;62uu>4CepWhWJG#*nt2=dg8Nlf_L4*YCx4a#3Qgv+F~8 zDx6V#2ab(k7EG1q)cZu$%pvo^U?_GJL)&b07;Ci0a5H^H>kUoD@$n-sDE)!ebJC*M z;W~|-2QER9EYYhRE0uQap#~P%S%&W@7hXXW{kQ`TChvX7Ga-n2b zjvHR}?7{bB#i|-D^IT8&Lna$ejAg7qVcTDn)MNPgTu0=1z7&uKM0|(CdsFM9kxX>V z!Ohopf4;jT*`Y}4Q{-+Co!{8g*JrVc-3{b?TOt8MQX6?vl-#WG^A}PaA7sWOk3rR+ z0C5Wh&;gU-s_xvB7p*-`7zhY}hh_oe3njdAi^c%36cdGlRU~j&5R#x^M~g0Vm!1{B z*8MBCU`q7cK&QN{)~KltrJYMS>R4;aw7J_%HTmXehi>ggvPFaq9mQEa%d-l3!8WWa zW*vf{XA*IBwLe2~2Rz>3dJ*5c`TBa9+)>gbUbM>K=It7?%AP;PBI}P^zKqzksat+9 zbW_FEr}0O=P4h3!RF(74V06`bxM{?vb&Gn8Rxd)k#j++qq09J8M=7|?rHo!FclOLx zqI}(VK2Yd{oEZ2em{=NxgL^OX%?m+g=EC-`NTs!z z&s>Cv4`m_Nr1=HMz4YhUK@ zTF4bu-;XM(26-DmcLO#v}#J^FRQh@exP&oW=0h(ri^l6zj|v&3e=3bnNKh8$~WJox5y(ZQU4xnz$M z0JzHOE4xbSnNz$h2_ag2#;Y&GDRrT`!gNz=4S{;gK{$Q=B>D3N8#7DOYZGKYD24S| zr1>awU9zLTp5%Z3t5#B*t+fsR8IWl<1%x1m_SDWs#)=_F{Vk7|B8D&6!l}VRe8hEn z;U@UrOykb2qH@H}nK><0qxjCu4B*5B4!-%3Vxo7!nyAE4F?o!CaP!{Nxj4WYOq&Kv zHZ<#uk?eMUzQbuEB5C>tMgI1Ior?P5RJhv;XHlp0_1J}}l|-h}ucy9H)4`%?u(*6i zPtypy3Y($Ng+n|OQbrCzY?Ubg_h6Iv@ER3v7G84Y>1(9BBVwz2>r+{52?#y)wEyOr zFr|MtTl-tCAKS5tgjiWag#c2=9x-UIJ6tUtMoCMX5=w;%vbqasoE!lhd~}9JhCDcG zrKG&_s(789P*B}ilHxpXKz%^tecVH(;Ao$74)4BC9+8{+lkT{d3d#m``D$|XIQoBsv%STJNhZFdSQVjmpKz=n3DPdGpwEA^-TEDmN{|y+ToT;#wblA zhmz4`2CD@&lTf}|X>$dJA_LY^hq`R_6TT;BvSBa{x|I|8mkH3xZ^e(oz=iV>@nxmWj1Y(WT8TXfeDMRgLqx_ zwRhzPfGrq*nD!qn3kl-%6SYK}ADG^n9KE>|?p}9`!dX~1T+pbD=D{(>0yU%cq~;zT z)6t-t0o`| z%?yMFp8Mk4e~wiK?Je&GY?F!07vf7Qp;L6N?TDalnid+M1i5!-&lWr?=}3i^r>3q< zk^A3d_^BzP*ZrdoG&5z5YJW7MdrRlmSd<&js@Kq$DJ%6-r|Rlro}l)AagSb~z1QXa zM_WrPMF9Gs;RLyXSbSxoM=o(@y2WeF#NNkqRbuT!mYJO=!@; z1=I5BTO7&+_!!zoGwDM-NWgykPX?)&7B%Z%q%$hJ#?LsVW?4b?FeVd#X!R98KR(Rz zis>z^fEn^e?u3O5eUXp%55v7)J8Ec62l#P)rkgPvT<0eyT@uvgnTk47UuFlnX=D9m zXh-p4CgXw#l&@V0JrD<6Fg@GRU#9SWRIM)YmvsSlp|wvi)YM>ytd_rd5L{j+ zYp=e*ckS86Utu6TaH6Jyri2TwxPf;tD?P+DWlcsKa9=vlcrTo|oI0Mme%tJscf^o1 z8|}Vwf)jK)y=D%*nJ6_5GG4Dz!*L`8ja8zfOi@M{Gj$3tJ}H3RB#8I8YLThZLEEOI3P)FX<2zd2LSkCx5uF7GS{?c zQ2XJwnIt`8OI>9|7z_3YOd~ReYUMk?Z($(4u)5K*%p@bN>qqB$>rg7 zsvqWc2Be}q=^(Yh{PRO%fATiQ3|N)4y98YsSlR1H9a*@PC@r|J6AIc6Kv)G3`}6b> z<$seLQ6Jzp*ul3D98+ucGEu=jr-ybFY9^!|JDE+~fz%RxizIZFP4>EKS#Z zUzwjUE32j-UH>U`$+j3+Gk`=5R_se@3tCs4g)9aVEsO#N z80@EGJQ}k{#4zoVFne24Wxx}{Rw6aA9$zRx_!R4XT=HsBI{*-c&b7E9X-D#~+>Fw@uUDk$Cp137&d|z}M0y*sGnx@9r4{7c`8BFHjCPz*TXOMqp~Cj(6=m-)P;|=9&Hsl% z@<#LS&yg56$txnKH!eN(E6N)9OI^SOAtZc~BGjJ}`O_^^F68cpB^pa!$*>|YZoL&xIg>|r`^HNQqO47u^Atx-uG*w`CKLM{_% znV-5ba>i^_7#l9IRsMx$_~h(bX89#$NJOBO0(O9Rh!fzbDI@+=178gZ@@EpbW*=p_ zzputoGx_uQ8_n%UO%2Mv)_S|YA8rGGrf>@ogNAN@3%@kZ`nhSz)?lxwoXL91CS{%e&3t4mK6*mGBB(e&Kg(TjzMyB5 z3qAIx|8We(9XEB6&hV;mdKlI25td6*XBN?#<|7IGC{!vTKS)){=W)j%tGA= z!8!q6H7&h^VC@h#T9sqMZ24gX+7F=7Rd;W}d`uT1*-@khQ2=86N2*XN^Y1^3fI+*D z?uNz4$RX+fHH~f`(nfimU%;uz-kn2~VG0EqUsUaHi%pu{;gM^o1$jc2M6%-Qq&~ZN z*!d$x2m?ww!XL(_!*a$n;pM-^;E|NnAUgPNaQITkZGLLXv_3U4G1K1I&u*t;Xf&|? zKZaP**=*-yUD-<#0)?Fz*%lbCPFOUSn?b1!NGhEEwU1zmeMEj<$sbZQ)YmT?$Tl{O zXtj86R@}0&an?k%q$N8^GU$^{C2fa32!X)A2RNO9i?2F zg3I;*qWe3tyZMIV0m=%0YUURR>FGPFfm{*%6MaLEI>o&tC|i?_svbmn_1xU-7BCVn zFCFggBfrrHD2MW7Gj)zd4Gqj4dD3W|Dsi%Fy~{A4ly*6>?kOi_@)at6OLy%g@teyR z)RRL$EWV#~<=&ABx++O8L)RxxGQo{HF>Hi<(p)gFNyov87xRL{xXX_HtZ z4%Vs7*ztsJu9dKTcL;+L;(y)dKQ0oi6Y`dj69sjY64Z}Uf~x<2jcSUDNAK+*_ZC7r z()I{x5sJhz#j;I6t;Fya&wXvWCU?CFx35QBnZs<#KH7d3FpS1zH!^0KEo@5Yat)Ex z?-;$Jc*G`a!KfD!t=T7w=aa}6kF-zM@F^O}f4L+ATRArx0mo?9MzaN1_GmI~Q z@JDZHX$pzAX%KEr2-NOT)wI7NOByQpy@|EymPw<5a(ac#>}=8`nQg)7dRCp`MPbz- z@5d(S2WO%{ehYbspGVQ?P6nY^z9YYV)hD-Yw6Q;5YF<{{!`eWqMD}cKY{W*;8#R(v zwfEMwnc~sq-p46NZv6m0hzPfdwl{9%v+kwX(1=Nf^fKrDa{;HQ8w1Ah0PP&U1vS)-YK((+xTXBU?V z$66sZYuEphd0JaEEX>7k#BPEp)BtOQypqCgqtHUZ^C>R+&T0820%vR^ukmP^Ua2s9 zTG`6S8S1t-;&azcCRXQe0nj8A(*5)8ByH%fnW3+BcOM=m%F^@J9{6}u`YraI$rqI| ztt#&^o~|Ojzp0P*=Z&JH@U5%W*>Y>=JqiexiDoxcF)53TEDMKe0zm3QaqsO|3%_ei zo3q$g=NTna= z8ONV|g-dFwiE*hBys)8fv#gcaT^r)^KJcc%Xo}MKsc@{2{w+#|^OGMOaorI~$k=@x zuM1w4sR$g+_o%H0#s?u4Xb(UGvO*2I{ipK)ln6=fE(DQES}50m{&AGQdTHS$_3&~i zP(Z2Gu_jk;d_nk8$tSQ+fxJIuE20E#SN;ftvQ^OmJCc1h{ZBho4Yn#TEicj)@zLp) z8%=(lf7UYJWBAQn5)V!W0&dnRq9!GPw4hDtYBuF$KpdP5SU@KOQEeWm(%w+B+c$b$ z-aGRTWN(bU;1`ge+oTUyi0Mzdn3 z&25Xok}0Iy^ZWQM_jYgdzp1AKx5$Zwni^b=Or%p?OBx_7LaV`D*3kIAcz`3c#ZH%% zl;cu<&2S5O8yVq>;oO-b{O*g=U3-nXSOuIk`k^gHOgJ3^bvQnyJRsJ)2LoGxetcfW#Ef5H(@%U>R~OUEFrd9F zozS2M85s!DOZ=8=Y_-sVr3Jwb`fIKAdmtH8ubW!=>PZ{VMjBpTK%R1<$o%J;oJaM* zGr~Q;*7EAh(*FAF-mwAFwhW8%6jCs=`fP`-N8daf?&iOPUp-SZs2SrdO7xeoTkdGT zS6RDuN}%f?wCRGSncm`QtTjx|&gkO}Yh8L~r?T@g!(3OAUgQu1q%bI6K+lIUj{tu6 zpA@-6@TDZif45}*r#EDCP-NME{(sSltcduBHo#tvCiF)+@W}6)Tr#9Un_S&B)X1yQ zo!|kJXG3)Zd>O;ts#S4|3f(HY_{O)d!w`{P-Mk`mOIy%{>DgRTXq1g(r0U4gj*whK zwuT&?^6)iciBUqdR;kFvBHcX3v{a_o1AGk~(`z_j8HnYs^k%nWcdCd+Zqj6~bT#~hvV-9_5UKz2`0ci!< z5?IR5G7W5&6c+!yg~ZLLpXdR-*YS8OBA}wVwIUumG@(jfDL|5hc3Y9VVJrva>aw$$ zOlo}_>!?gdUX6T@oewC|X=w7dlb@4N;wq##)VmwvFly(jBWYCeQh#(VkR9LI`>s)m6WLJuUVVUdaN!jrv zPt^&x?~wZfj8Ieo;0v9+0JEal7x)IbODJ4WPDeh0zwsp>zE7&J!=SB}-3U#>24}Rr z4GYwV7YOs?W*>6o8B1`imVZVkXOzxL;m+C@F9l^WPHj#blWpNxV#(X-VLz&Mk>{E> zdbD}{!3i9gK!Z-M!(wD`NkY1D>o;N?Zwk(JYXnw>=)I*xqU!n#*T%tt0%XEX9JMu3 zZVJKx<30KtKcZBnYNQtnc#DyPL3MxDU_E#%@`aPe)m|oINLBvX2c0NEu&Z^y|Erc! z-}Ckz){Ht?nm6)rv!$Tgzem_}gr|7aU3*YfzSEa=@@JKwLFG)+$fBim!xaIrF#w01 zl)Vy=CeSIQT2M#9CHMxp;E9{tP-0`aH;)X1ZvWkKI)Zkz_XKwdX?rzN!W3C9=VtqF zjd~k=p~9eoy997pgQ03nmD=Ge9B|#^pF_(79a2T=89WDX9Pl_H}G^@@7HH>8O!K z*C$gXx9!W)P7I_;S^yXH0Ow<;p`2{k%-Utb?>F4jOytB1&Rf|w0?p8R{O3Z{}=e}$yVRJul?ynoPHy8Zw-rMb-7$?;Qv5{MXwwYKv42{U^+I!-=6TIY% ziGu*}Kjk1-ssTwdCMgaMh(En$8IM2WKb%-T)qJBi)C`{LnI74ssrhnhP`DcQQ zh1xqHN&GgJF~KUy^{d$P)x}}*9mXQSC%*TQ12M9ok7j|O_9q)Wr2qT3UMjT zrhgq16F&Kov5$a%#2&%co5)Z)T{DK88KqX=C-5@nMA#KAl^BPFD+)hkV;EIE*OoTU zdN5*G+dAL0ZyByGAcX{*hCKoTNFSZ7qrj=h9R&Da1vZ`cHqY_%Gq80CebZ;jw-Rh| zyLc2OJjlE1>FNI3$ekKfH<%0YKjauG~OJBHq-xUTsOAY&~{;ZZg3nkxt&bIt&p!%(0=A0xR+dVU$nJ>o- v`i;@