From b927a832a759e2ad7047ab1e23033c0a23df3c3a Mon Sep 17 00:00:00 2001 From: sdgamboa Date: Fri, 3 Jan 2025 17:39:05 -0500 Subject: [PATCH] Fix broken vigentte with BV data --- DESCRIPTION | 2 +- vignettes/articles/Figure4.pdf | Bin 26446 -> 26497 bytes vignettes/articles/Ravel_2011_16S_BV.Rmd | 699 ++++------------------- 3 files changed, 101 insertions(+), 600 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 9089c61..1d55d3b 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: MicrobiomeBenchmarkDataAnalyses Title: Benchmarking analyses with the MicrobiomeBenchmarkData package -Version: 0.99.16 +Version: 0.99.17 Authors@R: person("Samuel David", "Gamboa-Tuz", email = "Samuel.Gamboa.Tuz@gmail.com", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-6863-7943")) diff --git a/vignettes/articles/Figure4.pdf b/vignettes/articles/Figure4.pdf index 4abb018cd2987a72ef6841a394095b3a82759f2f..490c680d24427f7fa33df06476238d306304e660 100644 GIT binary patch delta 22336 zcmYgYbyQT}_Z1M37A2*myQNE{1w>G~LApzXVHlC_4rvML?i#vNx`ys<81nbf&v&ig zU#!JjZ{7Ruy=R}j_c?clh7i1l5hC8u^Kfu-aC7o;bH3qYL}ypD`yy%RZ1S2xQjm*- z>;FBA)Lua1=i!KcP{)U$i;AKAZzmp*wQ}-MZeQCAIJS z&ubCS;5sMn`>WsB(ubLghmVTCA7+Fz-?!`EWt_oJ#Gqxb_;*z=XOEgqR)lT}^w$pK zGCcKo#Y8LOcN>^WGIbY++{sp%PB$9@Wg%N7=Ttkp z(^Q3pGkPE#i%Y5i>Se7v6|TNraxdl*6CXAD1%vPLXpBV&2bK66m-`f(^OKcLiyr+xaBB?%%gNQ|$;0k#$2RB(w@CUtmA2yoO}A*6I;CI5kZ&7QMbzi9zeKS?~o- zqo2yu5E|9+6(q5DoyY3dcUyH9)TVYM)W-kbd{eB}XWCpy<^g>y-jpA zcWGlc`57_#YUzWvp(j+Z*^%dk#V@!B|6Flj-b|>rlLxF$85Hb(T2LbWd*S3*a4eK| zy&mGgBQTcQTuKRUryAs!sZqdiMUdOS%`<(0Ouku>2*ebJdhOm_Z%Y22w(f zH+z~zu&{Q1jpt+_RcEV%^+9`5?inbwuXs7i`QI!`gMnFC7Zv9Z(urnPQDx2*G|tA( z-))pUbalzjpxLQci@fG4Mqfy>r_w2Kee&h1MXS~_C9_hn)sOY^-&qj-W1{VgYfJ&Y!?i)Te6i@jm-HqmS`3n8Rd32#7CKteSLHnls+0Dh zOylK!5r^lz+k{y;xDbzIB@uvMFTULrv*eTABDyDf&^fzq6TWn}q|<&7cbvMfi~S&3 zME(@tWiLW$oF9lQ*Ky8OOiaV?Rlg?DVd)uRFh@=Fz4n7lF&j^a6YMTe=Il#7AWvl^{oxNovQf`xwhqnt5bcUxQ;n2bjJ zhHiHd1mWdjwlhU{sqll|!AapzjrPGO2SX{y`nS}z`t2@bUQ@% zpog@ElC=lNX75TyIPVqa(wt@7oB~IiD&6{X!pj3_kzw`mD7%h0g`wAvTL~&FhcL!J zuIjxSS^xH!V09R;+x6$eV&ei&t!}lZxi1EX1cH!JI;Bb{S2FEz!3(*h`5`SJvJRAF zu=EdvYfUL9(eypZ#S9x2qm*oD>+P?jGxxcyTH5zj z-%D#;F=vy+=ns}N{+;k7M4MI3tpKa^tSZt8{A+F#q2gMhqj@3p_+wDXO=9yQ@kt<@ zq%1>+ub7VhmmeaiX{S+wk6}=}(7M-xuoS03z^O?~F$uD6Qn7g%U9pe5Q4pswtXI5z z$iv;#^|K-1#EcKTMGc?O*0c1%uX>}3GgO-NoY(ae6DSr(~?rN$+b`!?_qjZ&P7yBsFm}M?=3*Zpa=DRzQNKzlp9arv<(t z^{dr}U{M>~!`^A*{g5hq@G!l=`!WF{Ve_nUl53bUbYFr1fT+n7p`7U{+C2p#DPaBs zI?BB3DUobiZN=~3XN>hd)^ z^w2|1fw;INS*&yN^YZEZu30bMDjo!!(Cik7=-y-#3uJfdJZV3He{A!~6jC!=F_Gm@ zr3@G46N7PZe>_aZKn{Q6w^*w&LRCucvw{7xbz*;Fv>tvxXWv z)sit*c# zNcTaJw_?Bv-qU!UK&MaPQj#%}25$>q z^qn%)?k5wPV#mjcpAv=6wLN)C!6=VgHlNo{sA*8>#^+bd(1svN|y^4d}H9Noc>kx3z%Zqx2{gn&ij z;DcjsF8MgAjgqa-cBzfO&ACECGI#N^L2~Drf3 z96P*bOl26wsyiszmn+rcZ1Qv7iD_4_jJ`ulzq5CqO@1911edN~zt_Zg z{8Xs&tz10rbb)jnLA`jV2etoS=(_P_Cnb!>KQuAcY2ti8b35z6e$bTvXHn-uqbW9- zoZ0{a@Q{N}0sN}`yrE07N1WXw|F=r*c2Tg3(cOK%Ue9jP4-_ACU&e-DIym*>p}pz$b)=sKSGY11|s@LPttRRc8G*KcAwZ;THV{07gPWP_#`ETWK$S)95~)bA@%E1z}kZ zjI>X0RZ@GcR=N%eOCMq%1kVpC_j@Bir@tzj6p~%aMyS7yUdR3e83>%oi6Jj;#OZRd znj|c1!8|Y&{U78kUK)|yLeOqWw*zmnb6QJ5O$iQJ5RR zh&FWY5<;SFqlQ;sKozS*Q6lrjbH$NcyJ6r%Xl7F;95ps$F6l%6hWN zlXZkIs7kZy6W5pC^aR>eUp%i*Y5R{=S=9Vsk1HY5+T*gi^!eADdYzzeAz@^HmJk%61X<%(+AMXdd2JhNm#0z<+b#afbHo1u7Co`o1;eYF3 zL}2x)5RT9)-Gw3;W2)bU`av?vH30JzguPr0&2qIsKK84wz~eN>vd00KT-E7G zOg?>R_zr-SAp*vnze#sik3hnPMt~2NZS`P|^-)v-16D|(dZY)$U1L~8$}vAVnhJ+_ zHW80iVO1Ur$Sxd~w2w^KA9ne=#PFHmK_{n{K4+^HZ&{apM;}SXirxlmTEY`VAhGW1 z@zS+5pj&wlrDP?2o|ng|+k3?GqG4c`C#>C6*kxIh z_1fZ&jgDvy;7-L*L+D<{fc-{KM=8M>Yy15IbEB z=tp@)ckhmT*|A+&)fzGd z0JxEmG6dKsf8Oz)40a*bmMtv3O%C>s?QgE?t!W+xE!PZ#ej#zj-po>_x*ux1UX4uM?Z*c^iA4f3)vSwI4UFPNsXe%VW3+U{{Om*P47nie`O|0- zX&S!E)W{Rg50^zkQLY}G9YI?8!>z~W=7!M?`B_4~pA>)zKc2oKU+9m{W&(0#hu`5s zuTpQ%8K-@+TJh1}2Rm$h$!R1!WQ7FcPqbye{&4Y}!iaUN89+!&!=nGVl?yxhX%Fq# zW%O8)hBj%XGEctR_B1z_ah8^~L{I&CyVYxenAODbhjI(2i%&0G^^N``zWf5v5K8!F=J%^(w7}Em?8$}QVfGRzXS4uUZA}O|5<-1&N$Zl;oiPD z(ns=}bQ-!lw)?EW(4(Gz0#%*hYsK>{W1&?R8nfUN!KP*2_{N({78ZsrfLUBc6I45s zT+8w?=!E{IG@^TPsMv;j1SS2KRM^)wYh$lZOEvYu;1y4H>Nr8iwI^wN2)1=r;cE7NH$v9(!>W3z&&uG#K+JXHooMe3*eKQf)Z!*_fC0hsiEjk zV3+bYk~MrCdzhUQk_y8oL5glJ^iic5DhjNxUhe}fV*Dh1n36JPYaas z*}9)j<_?r5!*ep;7~k=`xlA6F7)=kJue#&(Cl#ix>xq51elHVw(0yHoO4M;QoOG)p z_Mz-u!7{~K`J<@ro+YQFw1F&Iff53~bZA=_+W%?dh+=@rYt|&=j%#xMF_HhvObX0e z^s0;T08%rnhgQde3}72Q=%lZ%I=)Qld)dx?T%6GajInA+!)1}D6tBQmm%J206Dcxn%Bj@G?A(>PgPW#rG3%ud3phijG=tyIA*#o)rr@AU6*~d?COPKP zY&q9PLVnFt?7jEyo6dC`pu2B_5&bF#0Aa8ibRl!#_{JWhTJnt@Ja1L}L1a>g@sh*N zP-dOn@wFbkI9aqV~^GcUqlZH3xEA(Wv+BoR}= zprOtoXBT7)AqSj?4P4hneB$jeWdl~>Z$JnF$`7!D9=~*1i!`l^0du^p56rY`Yw}4F z8C$oRDsn(zpMuhChKe#R$ixqV9*(gN#cXw)eld$sevD%RWLBJc zS%L?cN_sByv4UYD@v|MgoRy>PM+Tq;K3vwy8&C^UQ`@#>lb%>>+ug~6A$}yKKR>!%d?Rvup7-g_Iv-S&iwkgq` zzw~o@P2tN6b>SCD$$<*ao!n);KYM73ztQjtQ1tnO z5eLD}8+osCgI1lGGkBiCC5wlP+juDp=-)uTp7Vu!HDuaa!>#uwj`q6;H!ikwP9o?TQ3J~-x*(4Ux? z5?+0SB=xYqc)K4j>=zD_3=PlqJPAma_8uw6;TRZbhEL2PKi0WVQ_-i0rEsy@`qs8}f8biB&( zpCGP5Wx&U?<47 z@szDEL5&(tFQJHM4pN8H=h@7%56APy9v#f>tBMVBX|7-IHX+NM_1q(y_?gtaQtIC? zm}lL;;824XG!uzZc8A`ksfQfP!lbVoFGq3dvq7kLnWi|q)q(l#z0b@(+0R3!SQCq{ zP^zsE=eW6pJAZek3aN20T<*NkP|;KtUNO+1w7v|m1~g;AKj!NXQ)r%W>_}?K!%ena z*jotVHHMgF=Kmg@_nS`gxAbc?w!$sT3Ytf0b-+N}ZPMNf8yo@N)J#G$z3#Au$$xr7 z{F0h9c{hp2lKN!xRdyn)Q}8>tkhxzTR0Ao~=?O%X{7|(Q)3z5r2cC-pe<>Acwg@CQ zJwk!T{A2YG*C>&)X#kxXM1YR>-e}3P)Agd4}R;hZ*>d7ynu_p6Yqa{BouV~ zB=sOL=|!?4su{ODzIy^1wGFX-XQHOmM^KLlrIdm93f(&}mv8nQYjUYlJ`{}55fn6tzCr`J3_uIs zGT%l)hB1(k&cqsBcG%SX(2FNv56aN%B_z2OB`?3x zuT{>}t5n_9Rh96+1^&zKuk|9VpA=vMgifUji;lBax-Pt{1}feZ2jDHP!$VEyj@X5 z{34vX(Z=4EL6OhGP{p23zP$E5ki8;&q}l-@`x~Hbq<-m7A5NF9?9*Lw?s!mRXWW}K zyL?Xc++~6=dckSB26Bo5gE$}h_`Xipkk$J`?-(FSC#B5@TSg+ZU$$M2>HvJfqt(564ddrd04>iKUHkdxh7CE#u>iuqdSqz8q#>rk z3mETslSbQ+QPZk?xwufb-_#OMor1O|1A zQT&@g`3^g|x(u0-SF<}Hb)P}row`%)5@R5hwsVO z_D4PjFngZgUrQm0)m>K>Bc`(OCF$Z6tQdDf(%}0i*y``1wYPB@UB? zrcCmREsmi)GV))x@wiCeJJso%Zat`K zgI@Tmhp*6vj@C%quot(e0jr69!03swI_77+Whl;e-Wr?-_f=*G!m-?1Ql)HWH$fWHqknoyjOgyUNUo0SAcXgGzV@Z`NvJYHX$Dni9{=E}4ceRsMu@X_ z(2N}V-*4uVdCO}euO~7+7BmLl%z#_a^zkM7E;o;5S2;5U>Xf!xB%5)CNaw2tt(knC{o3P=MT zlBim?$13R)Bk}WRkf#`7^cA?(H>ot+2W2;Mj~oW@_P{qAOl?vZn2V%i8?@hCUVJk8 zk|w9sLz+7#K@0EV2-Q%*#eJ1*#GTq=3ex-~RG=SL&Nw{#a|I|JEd;Bg3fd)XrX}y^ z$-0}A4jHoZr6Kj8)?)m`-B%pR3)R_a=9x0W#g^Ll0cv^(j0i>X!cJ6EW- zNE}4*ajz=17v)>AgNp8@k~`g1@Mle~vGprG-Xh(+Wi`JW^CrU@X?fIy^#waA>v(pOmq9@`boliM-)5V_M<0H|9?7dbyhHfb%C&IQR#yO9kfVihIhv{pXI<04d;Du2-$JUw_hDjZ#N4!OoEO6Q+~2pIkr0N3t~ z`wdq!Z;CAzWa>&=#SEI#0<$OCA?|UtIZe_*^G`IHwlG!QX>9 z&8vp{B5p-^Aa7o7SD1)}UsR_DK!YmK$7J00z#j@)lgS_Xa=!&L{_|84QRzG+G)YK* zc(Xu~d6#lU`A+t&UQV*Ys|N+cAACGUHeV16e~JTmu3dLq6n*K`eD2x^*89_PWjldA z08RJ~O+I{^i*wg<;5C+cnOm#?BR?Y$s;~P-5C4$f_qVe9VXSYo>!%RN4_NDzHMB2! ztI)?vQh@<9v#VT+LzhAzt}cg4{4GR=E;;C#9!Qc=g)A}2y+|7hE|g%Qlr)f6WQ5U`A;|>$MX0?1tPX3gz-m_9#o(5?4q+ zdLNd9E$C~Fd``3XhvfC*;(UF-|2fGL9#~VAeO+hmIIV7C{IcAJ<88?;-R>mRiJ&H^ zkv8JG`syx-^*_7bWI}#=+3jZ*NOFNL_Meui>DuBni(T{~CAl|d5mu&Tv87y=e+uz+ zI6IMkcY-y=Rt9KN0E=>j(X?R{EPyQG*G|}t8&fe-YTzk87d*)O!`)rA#xelqMZAn3 zrWeG%e?lvUoY`YcT73UN+EMZCVj!%P|Q8?DU4+~&aD<{aJPb;fU) zs;fU2`kRt}`JNjTC229PyTd92NpCDoPuDJ}Mv|1x*`e#Vxc9`DEMUkpsaT+N<9x*^G}FU^$<0yykIZQ0 z6dfHX+B`!KV@F81WH)vS7*t_JsSSKrMBK-|d_p;#ZYWibvGi-Ny~eCr#>EQ>&a(I-)fc)$xfC@u6@Yle)_*zB|t6MZ*dn* zW`$+hvttCN#`eNOXOp$`hv^ha zTPAo=OGS-BwV|Avba9_UOXQP{GVfci>1rV(WOMtE#P!h`2av4g)x`wTIJuQMa)A+A z_r#2QkcpMQ(;E}zGCwt2c#DN_OVEzs4;p@AD2hQ@sC@mQjUKRUGsr-)6;c`Gj?vuh7iY7nKADkV5(zb$LxV_|tw^eQiZ^P7AiuF)NHZ zDzM-`V{RrxN$hSU+}t3TEo}4NiWlcVbhK&$Z*7Zaz34y{r;E}q0FI>gIfx)7d1Zm! zFjy!sn*?w2=f}VHKTv=tW0oziavv{i18+dE`YnzST*kedr1HW-!Tll)O{Ez-WCCx) zC}ybl5`s)-#Ph-?@cyo_tyCUc)@K_*A9{Z1_@BL_u!9BfKm#``mACK;`BY)mlaRgr zjv!0RR`2MDw#Yzr+&VN^R{C5&9sQnB_4%0XWZBvcq{wEsme8S`9;b?bIwbQaskdRYREARcp|5_d!h?}8hJBo6+QNXhXMuxdD<$z^P&^?|HhEnx zVm`?WzuWLsn8p#|yi3z~wk8=>@QJ*obl2e08PB4P{B~(*542ou}hmvy5PNZC70i(hZT@QC;?075EB27 z%vg~p9HoM;zIPnQau+OzV3vvy*Qm>Yc-rG2^N8?-6QC8 zBpxv%u(|zq$(bm1llLxr7MJ_OD_-ABe8l<%L=(?r4c6B7G-!*ZJ{MAjYlo3qv_-u? zNebrn2vnW+s&a+wkYCOUKe7q5N2yG_^OXsVGaqWxMh zp>o#5KhcHY#Tqc3evq-i_@7cP;!i&iaXKYy0aWKK8K@D9KoK?Lht|h?@JdX4gsY?? z0&{-v9092qw%YWL(-T0hIiCvZ#>_N{Ifq3QH}Hu_XHtEeEiSq59Ax1 z34a83;rad~Tk)n{?9qE#n0qvVOb^O=zo)^!{}P~@U;nn0JW|cgN*w^6YWbl!QyZ~R z^4jObB4y;YYh~rXJtiWw+wuuClpQ2<`Mp6b-^uKyRe*G60Tzik$--~p6)n@gEEfr0 z*%3VB=AWZaMuOb1@hX*(Y6!4i&o;^xQJQT?8A(a`4KIC*b1>YLzQ6~%vV=BMaeDq@ zoeMiL{5)u)4Ir9T*ZU!YQG+i4KmR?0fOqc~FETQr2(=F7#T#y6 zL0`w|$Bq%uB7zmNj8m()Mwkgn>iruoL_2WxOGvlM&&|?8Z7A+#DRsk_xC1`uX^^b{ zPb4eP4xsJF&=Vl`@sA#N96jmwx{9h238@m!pWp_(kacLkVP?*KEuOkYCHZRdTPWvUYUXD^VfsEjm!yJX4Vz6$#Kzv^xZRX^!=D|> znC+5PjN}VfkkDi7%Ngzpz*jXI^h(~-(>5#`8~F8@5Di!3{??gKoA7Y#PxfNwaj3Ah z&UcWo7i-i!+BPa4NI@1%VeOD~(wlu>{E0WrSlHFCUka7J5r?fFXlnjfdSrr>sn^*pvTNty?5`ru{Qolx{mAMF@7#K{YH%ElUy= z=^FB3(kvtL6rGm9jqVk^B8T7;!Ash6cT*mLiUi4vG8wgK?qE5H&T)^WqbcVf{NX}6 zvG+S1B(Y-L-@#uGO8vYuAhIW=y!LHNICA70(cyQn-uC+ev=)}&F{Rii17R08u8|0s z6bC6DyNA&G=6m>*Si>G*-af!*&R^+6@9zNXac=Q`$LwK$JsxuBHq#=EKXBy*{AsoF zG)20BE5%_AKsM9i1E?1`!yX2nHk#%iAtzr97L5m~Unb-v5u>4^-W8mD2Z#R-??>G5 zofa%Dx*#MCK_0!KBHaA^TtZf0d-X>f>rVZXGY14-l?B4@qHG$Nfq#(Gn@*1!#lK-6 zX6?0%hGOHJHOLJ_FW5R9coHE zu~^Y>iWS98I8^v~d3p}Qb!|F){u&Al;lg$5GsX3K{g0db1o{Wz z^ynvzEu>s<=#w?La?Q4X&r4JI`dR+xAiYl2%3UV6C;}@GAb)iS^E@P*5BqVpvwRV zxycz~s1=jaBI(Wj8cTBA^5v*O0Xuut50_%WHS$%`>GZ3UAD?bNQW=6RXq-3sRb>BX z#6_~ChW3T6XPYTt4%<|Fn3dtiBlnN}h(6lc%@LB@xxw|Nj$C_5O%0Nzhd9h+rsnd< ze8lU2>OJ}7vJ`)Vdie)!hRJ&#u6%Z9>I#EuEG2b`)iMp5!q)-)R*ls%x(lwVuCyZXICB5!cQ6oIWGpGD=yvII4EV2Bqr95e}xw4{@vSpy_-) zCpa3!E5HXNCNN#2C9F5lRk88Fid;@$rmCS;Dj`Zf?Xab`*m{b7H$ZPN$M52;LNRCi z&ssv|@Wji9Ev!88!xSqtg-69tcw6xXP?O4DeQdwOjIu2eXbH>>l6u)59rIk>LkSlc zx5kgFHG(G$G~M^RK3zZ#z4>Lkv8pMgyWun;YdcBhG^GaJpER-zB|mmml7>h73*KUu zg3htGRW82I6itjAV<7_2Oe^&MD}`stA>V2}iaE)PnY(LY4n}F zT`1rsXCg^cu^2Uv$+O9se z6B#uy7^Zfh%3q1?X=zNV?Z*_N#?cZlJsH$>Wc?r zB8X`|glo+`r6U)#$pO!*=rbW#VZtjxs0zwffjum7XedDYwxLy}aH#9&iIX>?cPr85 z#hVRaTz@JG4d`K{$w9T|G}!XYc%p|iM;w`8&nwzc2Fpz z@$2@Wa)_v69qcyOyQE$v!hS_>3{Otr$CaDo>}5!m4y;`f6|`i9JnTGzf0ukiy6B``O`4jhwX>eF;xt}Aw4GfrnZ*&!*&TrS9TghO> ziQ$9k$G`O^ilUK1&@7aM@MGuEnxRvM=v6&EVyIkWnT0lG%GDb(dMa#2I%*|BEX;Tj zhD0r=4uuK@j(_bukoyZ5s6Qv_kZ->it%hAebFpqFd!)|i%C0S51 zGpA=4|9O>pxH7?)QiOzdD8E)^>qE;oUC!HY+uOv2PA4PT+jAIlr7#CRf9`?vQQw!i zH@{kIEho%+m!C@sw{2PRIr``5w%YANHC|cFE+vdX6TH|#mlAB^oV+rJ$o?fP z6hqS$kewg5VAA|@ReeTKm6{s6UldhedDo} zZn_ydG0aRkrJL8)C3e<}U+wQ$y_*G!`SkN3V;A;XRI-v&Bhr`1kH$oNO~d6hUruN5H5@2U z{*4tTuzwR^aoS<{$Hd}+o_NCQ;C%vbLyoZyyjR(g8CU*(?BMdhRwYPr>WQNRG zy{P=LwDpCL^UEU+m&>Ir#Q$XefvOm6EYD%Ttx@4Rp=HIj+2J)oo%p%10oP_j&~~lP}5?(L%m1l`C-k9hxbySuAu4?hd(w8!}Dt==5OrfAUm>*$rnPw8gxXJC_r{jhJR^*;FE+o`3Dkh@xZt1 zkz~?ru3Yy97HPUBe6m@R1G1WEzu<}8_|~p$rf+2y&sJvOkf>-N&9_eyjl+F7Mk#`q zBWpB7dwYV!%2*J8cne;J@wbmE(=JO&t{8Kbn@Lcy-PHXWNB?p6EA}6PS`A#2bloOC zjaaT|o~c?Hfy~!M`#-xNEeI+Tc8=r6-On`(W@gu{?Wc2A{Eh$A>`q6Uld;HQzFjSh zRd^<)i40{Uc)9qZtjv8SnzG+yPRjmyg{3SXrgL%tE^nl=2QP-s1$%;))7pC@0DE5n z*c+-erB-P*MwLPY%|U6%uzBM_Y)%b-{h3-?R-`r#!jG$HQ%KZ%;fbrhwBxR&9NO^} zz4er+oBXhntUtX51ltiOt*dBWxmPi%&5AdtQo(S&=~^{6F3&q_>=|fZ`A$6MVR+c_#(g0Cv9GhBSk^UYO=Q>%}a5dItVhB;;%Xs>F_ZH-#EI$p zTAZultcghoC7m=Va%br{&2Vh9IF>P+kC4rFQ4uEv>~A{DgSp=l=Ad*wF*6bk?o!Xc z1;I_f%oR5KP~#596A`N{MU_-2TdC*-^!db0K5?DRuY=klRT1}%>ZB$App4(tI%x?% zP^_r(0~}4q?nHySnqd-oCT2|?`xm6qQuu6x+ODHHuT1C#cE+H9*CvQDAKkPr0#OD3 zsECPgZA1eZ(tD>(DfUf;|IU$rR(F_;!IGDkk4O5&5O1>Lm)}9LPd+1f1)-! z(i4S!b#z@vXV7*l$%~?oFboLIa;WRr^ZBgQIUqgploxqLdc=D3pNe5xATKM5^gYe& zyfh~pve+fNTEYov@Rb8K{t4-!(N}Mf_P%1>l#L${8R}Saog7cC45#(dgm$5WnYmWd zS^5l+IQ+pB<=9Oe2^Yk>H?2qGGGjOy?Mw`v9`MR?%W0`9bep~mlTEaD+VObQC}$(+L4I1V)F z6E}KX9_t~sD;2i?VFH)k(vokx+I2FWw2{108m9cq$i`5oeD8N>QlX?Kli7|FeQN?< zO!UYTI|V`z>CpwsC)!*LTWs(Hk&nFmomxT9nrJxedZiu8W zl0(A`_Q+N9QIzrGq@`!h%80R5y|oU{_2T3p|9F)Mr{`&$L10)ytBggJr@%8W(}NC8 zTr*w|5@C}R`HBlx=i>H1JUFkSpBhA@qFo;)-N#cATR+*TWZ&susVvNc+9OoJ8NMw? z8o5DHTZ@aHw%)Oy5P z9J$doha7H;A#p(`2_U_jjeVxLSmjHJdd8)FFRq=8#c@TpfHKma#A{H@Hl6r!%71Z> zT!tR&G5ZGhI(KrJOq#cotDPXwyvqe+FMn!Cl(P4UwUi|^4Qqi|Nco*Ij>4$JNKN|b zQ<#pE{vd)Qe;bAggoY4IvhH5Xb{s{_9Q+jYK&-y9mc42}2A4-&FK2fIZUK+J>(4mM zeM)-A=vk2s^Q~0yD!6_Vqc3|h;q?Y@L0bCaOEPaet7y^QJ)>59Vm~gm@t<%57T28T zKmrzM+a<~xpoQr4?LAkp3i_rQC&xI5bWj-!KIIVOZma+Z*`F4Ex!I^xH*-ku?S6`l zqgmA0oKg>pUo-Xkvu(>UStAd`_u$CNSQ8W8CyGFPPdezUR`|*Y+l_jDK0czkBluXS z7ePwatVw1?Ez(44a0OzJ;Ui7BLQ)7-DTdV3B6}kos)CTJI?0*Aga6 z#c7j*F9$r)!3RMb$G}C$QNIwAC6E6L13$i}ixGy=%e)33CrJ~DzrC7+*X6b5E$uwA zppG<~FDEtg9Gm;MgS5?Yw9E%(Uoj#6By8q;W)^yzOJ>_C%`?WY-=$y7{T`>d*xtU} zY8Og_hzgJy^7O0C3wWIO_@^KFydtHJ!8xnc_43$p-Oi)zUO@sGRG67X+kPYhrT-_| z+}J0yQHACav1o}eu00MmP@z{-VVbkE8DB_A+~Sq=UGOFubGteXI&A1EoDLZ~U)$%! zGc#J#9B^DosO9=gD+ojMLOzIv+T9*Oqy|-=@08;u8{70^ZeO6_13^PS9IY%2LvrYp zpwGj9ib?}LZYTan7$Qw2rKE#6IfdVv7{}(p@xYmu-#>bIi|E^XsH?x=%KeJk<}$c1 zgJymF`j(6*wZMrbK_zjm#lN+zc5AgTAb=nfPM79~)2?)`c_2*RYHC0S;Vf1mj@M9o zj*F`(eO!dlF_XGOQ}nvldX<7J3zN+rNcis&>WdVO)#q&36-(_osp~KWmC`aZ`?T%N|~~bN^<>&Fo z$w;0WP6Y%MGiSXt)CqeT2^q(TQl{HAF?nqsr}liW{uBOR7l3C{x{8FEDgz#CL}_S9 z5m0nv9V|rk3vP;QjHNJ6MxP;NfJOW>6XyVBT3$=szRU<&V)6^t!er6;<}HfQ0?lbP zO7}}=4gGB!-@`x0URy7+tDY6lgX9-W<)3OON7ks4_O~x0S_xc3TmgfVo&2#E7U)hP za}u(v_3QHd@;UrO_0x`%2WgAMJ4FjLMqSbE*Xgl}781RnXYpkrt!m;| zv(=s-Ol+FU$XbK~blKCvGviRv>gc5zJ(Fnkf!+-Zb*KLL@Zccg+6gxY8GNQL@OvfW z65SgJuXAm~@9TRn_P%V%wWCi)BGChcR3*v~_kv*>1&XAvxmuffZ_v$DkcSPIR7{ll zotWlFWiyU?^6YzJC6*&qmLr&b+o&_!DAiGXW#0LUfjsHf%(e!OJT`3kqj+R*dtml! zme1fn>z8FGUc;%K+QXw_Ew1Q{=j8dJEcDSr4x3)dX`#vd4-64$DFb=0tv96)(iq4f z(Zf`%pK^t$j{9^=4)=fh$km127iAe<>w(Av6C%L75?q|% zwES@S&JGp2A?w|-u|01|y|+F~bo^xy22{5m*6B8x(?;MTIkF>k@v)_G>wvnc zj7d;N7h$$UgM2n)Pn0HnPt@khS&U22OPKA*`~KkRJDP`YUD+8h(cRK#YZm2hb+d+C zTc<&zoRK5=z}4~DwgY@!gzLt(46jWGWZ=R6ZtEp}dSv2uqadEsHZaewvZ{dp!CA6k;JhKl@NnsZ*slV_+@-DqP5 zEMIz2ZQX*9Es?$Cn-ty~Vmag#^+^t>F6N&}qHew-)%Fw|?p(om9!HW;jJ@7F?8&6r z6HiZu%gZ3u8MIm^?{yv{@(UuQaxFO%16b8c4dDtMvkh1G23Ir92WQ$)&53ZVA)4D? zGXlRARY;{2b)L(i1-IJg1lc3YNZR;iIF&T^u?W7k%xWR4w(NdJ$~ldC>FM~J;#;kV zjd@pGe{~r$0p7c?_+dj@nO8)Ha#%L@7-19;LJfU9OjzUs&{zOej3GePr(+dKv%S8W z5XEfkS#?;lF>>Ak*j2k$j`i+q5@oYxrJ_Ip{=}-eu*0i-;j#&w?{NpOP$digs~P93 zdHxF3raI;qfFiE@m$qBAw>P+dWGFEV^*#x_Vlh)S7TdJqPdv}& zYvNfLLlKV*&%ktDZ%w8&X;-dv(#`h7Dq>fbb3Tu!&%ca;-QYQ1x;>5iY;Eu?40MP% zXVXVk-%#Ev^WPe~kPn@JryH!M4DEk|U7DyM!_S9&_v7d+1Kw#c&p$-rKrA|fb{v;c zyZL@T<-HMT9ZflO^3eo(86I2y1jWa?wErg08K35u0LILA7) zixfB{@*}OSjG#8`Wb+y0GJSfj^8Z$f!ikOPRjTj`DVk#t1%h2wJ@t8T8svH!e@bBd zDVI6UNh{0*X!f@STMp%Uj^LFM`!A@th7~&huR^}xx%cc%c*NzBwU~F_jo#G{tRCU> zgqW!euR@O)95K*L!3wLh0q+Amvsso)A-in30t?``IPkfv#eb0HEPLHSzD=wXM4&#UVpAOHS( z07-HwhD?&k-ej&IrHj;`8j57E!eRVzHEi=?93rcnuT1`CW#8w^C8}D2f;tHt@@pR2 zCcR#=`M;nqtNHj}2i&3yi<5U}fZ#?jEqRpDVVwcbY#7}xQslxpX)-6||JpdyaH!t* zfmbS$v1S{4vXp%vd&w?ii!ez>vV@tjgcwIzlQpt$AxoBskbPgWX3${lB1?oAOZXkt z_szfmyg1j)eVyw%FU~yAb3f1h`P}gw7AB?bq=Hw0-wI2Bxv4|}h-JOCenW6EKe5%l zgQgF9s?YfQw_r3>&^8-oTJk%U{)5wwlP+iZL2*W4>1` z@n)`q-+}O&`8f)B`}4lWoWN3E(2MFp zM@0PwU;)+6qV!ABKKD9rHg$eW*!PU+l0wU+5orbBYx$k`dGj0Y%xV%kRAOeV`kQB_ zAr8GnPIjA3Q%X}3fMCC(*ih0oFF6k<>uVRy-J9!GfZB1%6%s5c;&;7%QHg)15m;7w zf?$oL=6;MCqPYW6DJxmItYv@6k1EUE5Vbuqcrzc@n($+ob3ci`rv^j4IXpyfsNAOx zqvl=+dLS7_=~n=2hL~r8rN-E1=qx2-IPtR9>hOJ0)O|D%r`-Lhd}ztBG;4|Pd^}C$ zq~RdD7s5raUceupQT;gPY)1^)#biPgC%9&-o1UPzd_Tc11O&OX97ncA9WVCkXHV5^ z_TMk0uP0;y@aHx`u@rVPVYx=JWF;D<^MOqymG}Onx>#VQf-ZlO z8!mdD9lF^pq~MsY;y18madU7+K0JimU_Iy*MSrdaVK!mp!?S_-RC7&mNeC_?UUji| zwOl@EvRiidM)rK*pu}t1l^b|yDt-!q(E`!RVH7%gw*!*dG`c>HFekg4^1dkvs=}Eh zrYAm$+6!;P%e~nWciMMWmk*GTp+(giIX0m4(=YIW`r35^vGa1W9aaseK19q-i-;H* zn%*Q>tzZ{yG|Dp3b9CBJSV9qowyX-ItMLvNr2Ntnu zte;{pjZk~w_h~4fnoQ(cb-kiRoW}TAKW@nzcNwwtuN{fH?$z#GHWJ;Eycxkd+!_jf za<^&2PsU0<@&(H$JKfx<|H}CjS+=l}t_j0qGI$>pTE$rI#bIZ^TfiTt97WXHJg0Xu z$(%rD7gz#y6~Mj{H~>p6(EW#SLl6gJCvEI{&QV@Q6d{}uH{SVD`p8)$4<)zVx^tC< zz|nJJ;m=+A-q~ecZ?`O5N&^3jtM+;)QY@vHvZ}$iBZkwi#1JS`VEwD=E^J^pc`x9c4592RUxyUOGwJMPWM1sz>ZrJ6 zrM7Ly)^?R?aq}$o)Ef|J9$GAP>HgjKV0}9(Lgc!$0fvE>iv=f=vJX=bV!vt{e0Jj{ z>Q=y0Q0f!9M%f498KKS21|dL5!7JU5tJ-qAKPwK3t=;ZAD5|nnX~gc~fh~JeMSiu| zW&BvBj};ACDG!<|EITTJAWtUc2URb*g$SNb5YD2iK`pAcd9O*Nu+%pjYzUX?4)l4? z`8|A2zSvIoQC8#x9rS?7ps)9iJ=JC35uV5=WEnUcWT(XKfI>LBjyvds_>cAJ-5*`; zBO}$4dt(PqO%c29q4e+)KjP06N~L9IQ{d^Q zD?ZcJJK7TIyUovl3Pr5oJE%425SEq;LL8r#z5@6ReBW78W}bj;vW)>Yn)&}gI2#-j z8Ouk1NUJ?N(GJTEH{p+W$Mv)Fe*Ib$b{*ZxXu|lM-kp>sB!*^+tZ;Y{WVXS!Gl;*y zP+<#mVi0`#Lf=mba$LZ8=T=yGdC9QgV3eeMjcl+2pQpGVy;S)cG$jH+yCIYNN#E48 z`562ypDNW#r|CSLe9_!{$btjCiTydI^TaA)+I?uT#+7`IDsI*g=5PsS62*&XwiiDJ z#@3~toRr77)=G>gA5R6YP(_kF!Vqk24fuMVIlYr!_kUS=7!@DpxYQV$Bv!)E5!~{0 z+cbY*;=>?=jS=w26IHME8l<4I@@!xpmtosY%8`dU4wd>Ootq8sRg5`p!)HR7LO1YM zC%S}?12hUnvb-Llx8e~rGEC3_$NK~+kLR;x@ETDNWJHH(lzDPsp!&Kbyqq=)m13!0 znCHExo$|m{(kh0dnRg5+>}d^2{|RAZf2$~|A? zAE+{lKbpRteGuki+Y|1;&bWB2#PPUDV&{C>Ev;xjzp56KBuhoQqe&P+oQZ*m@b`NYJ?dv77wtJ(&ou>>;`dS`8#WX$cb>j%ivU?VdS;V#bs=8{CNF zMj|h8CL}})zLBHrc*)BJ&6|#ra{cx*6#bMzB8lDB%u)4nADD>^u)Q$GM}IN7?;7OK zTJ5^3`X6V$f-oQRdKD?Ewe}e(UaPs|D9jM6u3t)4dQ<5|#r>iuMq0u2Ax6{r+!eRO zk0fn}zljKR#6=3;Nilk*T#7h)_I!TV7d&|EVwR?>q>pVjbJ%ZbpZiPE_ZsIC~hw)$VZ@;lbnZtPoBs*^+1tM8KIdP!?t1^Z>nH~7|D7|n} zh`=?CFh*G$u$Ymq)ChR z={(bl3lWN4)hDR-cD4XRvJtuC7_2a@fI8ayx03i#X09BuBE#&unfY_3X(Nv}*-Ff_ z%e#F>f?xj|!2EG#SXTND!Wf?s`Z9!+&kI4dn)dCP#dgHpeY5x3lA(<_JJ~|U_tg|U z`m`PA;poq<^kZ{z)yNB>(4;pl_2hB8kHM2_*1g;Qb6pvMZFXP=Sci!n_ow^ax`(4% zUrA#POBvBWA<%!QQj(z-D~35>H#ju9xTn5Ajs(<(O-DHS_HL=+RXa%Ca$7c=zbdCok$vV?0nJn56o7t!oi0G<A9kgM<6U}?s^SC*2nwZ<|wx(mfxu3L=k;g?EZ7YVA88`g`} zo=~()RWRQ?z58W~Z3C&>u^^9H6L-W0?x)`MmYI3Jv~<^<-`#p^x*j)9%W@B1nNvhq zjk2_)AgKXgGsB%{a~y1?<8(lyzEp^e?>4O3rg)tH6o>j0&M0*`0r_N^b;{Z~rakcc zes^BOEt#m?iusK@6^!E367psh z4wA)aMAitENa;o z1Uc91xmGZJkNW)C$=d-Z`?b#ERX|EWTvD>3UMnYS2*Vt!c$M5};-ELde5=lk)W(qf zQHZ>w@!DMYkOTKDS?h)?7NH~e5hH1do6A2f*z?Ps@O;`O+nPFVC5~FrAbKNXX(c*R~U7^HU~|OYSysajJCz zh1~q%fiWa)_kX_j#tjU^=YuA^X59t5MZbt*I5hAu-c$2V*WKDz+=dTufy$&`V(e#7lbc)u7$W z$||AI55m<%VG@_=C~j$|uzB!b*{Klz#v|XKSKj_de|1xmz2KjRM z9LYkByDTp63^@9SUkJoIQwpl{4Pqnj-3C|h6V03WpuCxf7Grqxh8w+Jfr}RTceQJo zw0DP}S0+5rLs5~T5AL{7+g1Yhs7{u)A0f8f$qDx7n^V$Ti^+U(c)r$8_g0Qq8B z)HQf9yLOETe9PE-EB_Y7j1xwTEBz9GmG!T})+jywrECXBA|)KN5+VeG=ken-ZSR{d zfT}mtV*xv*g3>rm$wLK+KzT29Ve8|lyIxB0ot1g|epAlm_LYvSwauYW7Xi}Tu?O1t zuQU6DkWcR(g7xy>;$mGa5Li1c5vST&n$=kR#H>`+2TL{5p3La_qIs}QrG#N#Xyt)= z%gaw&_vpAfw-@xf(Hdc|g2$_yGrU722wyY^$v7D|#Z z3m-nXA9NcAUl>z3x5N0AlFqH_U?o~%)%EpwvXTD?2j4niy5YpG&AZtt^y~^dv;6nO zTl0x`=WA8ERMWcD`UGPCzn}eSfnoXDv3#r1o%5P^FYj?jc12gJZfZbdIC%OKAlbgA znWZL#jsudxKrhPYSWegd=bM3}Ii2bJmZrLUQyR^pKKg0Qk2tOImPf!*I(4n0hrIx2*nYRqBZX+Bzb+u6|Iu!ZtH7zg$i45&JY zIjpT4jg@KZ#Vx?u8lV)WvpR)6tZX>?zS0f+3Qo5m*@=H=IfvMv&9(QR!a{tprmOk5 zAH+RLR5~vGNvxX^8+I-?9_wLVScWQhZT5;vMk=>5t#!9FxsF!Ok18}hE-2EJuf=}J z%C0xlfKlrja3FU_6PueNZx(5>Fy8)?!=~;S`x8ZC}Z*4 z{k8PLXg5V^Qu^^hmQ(_)TEoVqTXA{*cA{ zwA5bhPv6tXTag*s^<@paw<^rZ{nVQIf;DGU?bS=<&{t;9h0Ju%weAT!U9+H>3Z}j4 zfAv&qP#PdIO$za^m;btT-7~D{&aPyhWfcl?Ub(n-h6e57b&n8vCoMcxVO9HU?|N;6 z$UA&XWcY$zXC=?C75VP<@A;116b`uI%)R^D#D;^8sbAn(@$rO5uafc#uOoKrv0@Bjm-7s`s0g(pj2I-EWyK6|17KVnA8oKMf=<|Ko zk3X1M%$oV!Irqdqd!K{P$o=8r``_Qvyk%o!<6-0G^lZM^!=LaZP z-^zI^G8$dF53yq>thgF^13EFA{*38<9LFZKRaQ|k0)NhrY=vUwfYXiJoAuj^&RZ`& zA=2JFcdydxF&@*G*JkM>=Nqlpn`6M)(dmM{HNQ~2y9o1V@pHgayX8#x3b}AylM;PdAQLkAiQ$cWN}?}`wDy_@cWP> z-1ZvXTzF&zIJ`RAmucQa$jjWWC#P@uI_m=53HX8Y%ymo}#*JT&oZCwdZtjmf&R+4Y z^hUy4_e6|@dk{q*bdZgn+9QuMPZ>&&>nCup_Bzz?)Y@AWp7=zV<3R6Gkrkk@N^p4Qf1Q;%Nax>&AFmy$Iv z#6v^_+>vM-*TI*ZTC;HcX6=2Jn~9Mv9^vb$1<@A*v+&c4uGZsKkCQTQofggh8@ISq z%aaZ7vQ5NgCRgd%d(u~?%WM;q!XqoWzY&IW=FC^1ivnNvAKvbT)l*`n+ zlmdXRc#oTIQj5Z^W+CeI+cu~1w%(u(ToA60R6}Nup=>sssxMVvt~fLXeM3l)|K{U& zPtAS^f6s&YGgzuSKYdDC>srgy;_<%!lOHmY#u)aEC2xT_7m6o=xMg#_opnfrim)hv z3a|6whQ=KO4bVW zmVvX@Ip=FPXYlk6XDwoH?&08~{!wAl$DR8;?@IJH(MnquejWFjI{J47``YL*4Fl)r zL_y63E0XYRSI*fQE$3&AvqK9CrV*N>lkkUpx)u)ZP4=8N6K;wK!4GqC`y{J*bWN_K zMFv-MYD4klis9NN&u{STk$Lzh zb1Uy#cj4FTB&UDeT?X9}4LoA;J_Fvu=ibyRI8KKa_BH`#)F<9$D?*>R>9k5}nBsG~ zTDI41gKXCXziW?0bcfEO{7J`aC94f^c=5id+*Ft3hfR-rHs`&bpagF+leyGF8XHe$ zs%qc191Co8&4+UXDuThkf6DY?+K%)NFY%~~pRu#)EL>;Ilm5|+pbbX^h@S!9Yq7`- zW!`g_QRnIIsef&5u4sI=MeScSQL|PhH60HpDv@>kLa>Y`&cG%n(C?w;Q^}jfKi|bC zqo<|LW)@=aVuxPSmiR@*^Iq2%Lu;*&h^%VgpKU561;uY{YB>iN;9u-?R-DxLcLmv7 z;9qQ|Ix!Z_esb<<8VM^44V3+kY`lJn;F$TFKoERWf-o?WB7T?Pm(JY0~Yk5Nf~-zfFfSI1GKI8!25_j z_p1>D7apjr$v#mK}4eVI2bO0%M_LXL)HO!h0%F+mw zMF6oYvtHZ3dsd@43bqYZFP%vne{)b{Uxxl%(=&VC`(0C_sZKz9LGZo2)#*C z54(;zJln4Ncii7EIxAQGL4-%md$y|{SorDh);7b6-9Fp+=f&%D=Y_*lr`#*zbX6(P z3D+lnl(;wqiQ3} zk)p3q=i<^zuiC#`W^5hA>Eu?JmO?1pG0onwQf{MAr(vzh+!*LJvx9)Nkv0LsII3d~ zn246z&@%t2)ke0gm5}9T&6J(B)aSQSYKl8<+}*G|;^|qVv$x^foj{P*0j{r zk!KJk-lTMS?z+N2h>j)4(Q_eOzhAf0S@7VouA-Xa&mG<8vc^ZT)3bc8qeC^n2f)-c z%WR-Se)kQ6t2=<6G%|C+Kq*R(b=St!hx)U;j1+~Ntj_e_w7`=@V$jYP z!WP4+ZRSt;;VgxKNmTR#_L=1m0hbTiF$b4A5gcGqnRUo(S2Oy!{I&UZ3JmVGGEwb_ zl)DhU+_!v%->%`&`$g89<7L1`Z+kw*^4rWBg3tZ9$8SJ?Z#8(w;EPW9FlmX{2;J;A z1&rY3TKIiT@SSeV`>!JiC141wl-Wyv(?Om=N#TZ6Aoof$9Lz=v#(7|ruKk#LjcKgE zd)_sw$W14VKP^_zw8OG_rBdwUVj&2v?`RCtabhR$dS^QL?71TE0uKqVj-{cfyvr~* zJ<>5XsJT2_51tKakyo&h3i6t{7QjF6KJG_U{843br<1bf2hcKAnOeDP1xIa}RJ#Jz z?FdPQIM9oHVTB?bab`aYoZb^Ei8Gx84+RFhlxEx9v81ePA*(q zhXl1sV$!q2Q-6yK+T(;#oZK*$FgQi`1E*`!^-Lu!l3)&J``{_qO2u=UN^->%0_FC4 zXu=umwe7k}cyT{!Rnj5T&J5(?0ts5Ao_CXG`j&zFxbK_UgI-CX=P@sG;C%@-D zJe1Hcbbcnet~U?OhXzxB!>d+P&b$*L@egLT_2=3%Ju$SjCOL*v6VyDZi!X6D*4!f> zXWQ-%=3rGLsjI6F>fZLVc}mBhw{J$tBK{?ma{^*uAsbZP&_$APq`^(Q`45=5G_Li} zdoSA`lCSwc)P(>RU4_0HP#z5F`s=Qyh#=}=Ny%mfbV8-x^xHO-0(2I(YZ zBkDyAG^H05Tn<07&EDnS#@U?(h#sMnKGp%$VX-ql7%=+Y-~MARdNshnV~0w-juo2u zEnsccvEP3Odh)&D{r={`9h`D$O!c(aKOrXVRfOR$>nZ4l`{$p0DAc4IlzQWTf~rcv zhzO1cOX94mpJGR{#5U#&N+&;q1OMN-bIzdiYpQ?3$LljNPKDZXNJU;q_!(jij7=w3a)s@IY{2Gg5Ic==&Qj`5#b<|l0l;bz0$mGk@irDqj_Onr3hlq$U zA6DaZ*HX4iwRHCK})@MPwbZ%bzL31AD!xp^ipu_ zFy=ZT&dyI(k`1=#HwNxNU9G>6j7&O0AyJl%@a3q%!78N{RV?NnN)0Fo)N&G>wuksAM`tAD@pKwX6?6kqXNv!6nkmyJec z+(|I&2hSg^o}MFzY#W7hO^g1&u;W9_ur|i=yvX3?V(K7|~ z0+)44so@wweVH(HWOC=rfbuN82Q>fC!(V5*OkJ$}{_#+9kI+-MD#6oO*ev&;zDm4w z&Cnpdd-XB@JFOD8*EPGsV@Cc5Y*U1)Nx6tZWkU)c9JPcU$ z@3vqy%X({SN`ZCbr18bD-Y;}0^iKtNKUyzGYt|}az=EUM@qlJ!M}1ZJ?aAmXUV@K2 z37$UPOH17rkCFx`#5>ZUnfd&F<4vcPlw8CO;@k!33a2rDOD(Vk*S)kIZ7f^Jytgc$ z?g9wuBhbl6oPg6Q4y@0oTL7H&d4tT!7kj~7?;-ThC@<-1b2?kMp_fIJjB~pO1yz>b zsZo{R`MesFqUmO9In-_9k1TQ*R~*mY#l83@;?L<+FHD0n^RtKKsj~WM?uFq2cusy z{S@135n6WvJKJi!Ea3nBoN)vvh zH05^JB8V@}G^cA;EADz0oz=3c9FxLm-iP{|jm;t6;*<{}#z;cr18DGtDboXmUs)Ee zW99K+&!UC`VTq%$4HD->5YL$#Zatmq{fjP`>Y46>mcouE-p}Q!u*UVs-rnX`Q*H>dgQa3e5bb+WD!Y&cRzg`-n38qUs^5s^9j~1D{DUVZgX6n zFjx%Ib{$qW^k|u&Z@99G-?UP&>*ykga2iLC7=lsE8BH`tom26x4n7H({?Qm5*;^e- z!De{k<+i9jxSKD*>ym8+F=0lvEHo8o3usr@M{zL$ds0MrmAK7c=$ev-N>jMfcnozQ zQlYUlb>`nSP8}|SrI2ikLQ`HjhH*Nrjw_@%FIAE;o?fYQ(`}mES=SIuOe(G!DLVQ$ z(>%})WkL^wav;MVgUk?M#_-anI#o zq>ZiB-ky9?98(c59tm&Xa5lg-gl!@;=dG{R(~nm=eWKM%E~~vmANe`nvaTX7_9Gdp zoNsCW^}4IbC-0!h$|fN+Qu;1k z)Yqsg{zEY(f^JT_vkY@IKeU_#at7)*pGsZe1lNhLRSS=he! zvyys?ynsoq-srIc&Z)Vf(HbwJM!1~HpxgvZuJLyk~_fG?2m`2b-cB))ruYc(O$@c?|}(%0ThkXLZ=Z#pT}k#zxNbN8Xos!vo1;A4t9ZFBj=>{g zGG%q&Mh^hk-`}9A{gapOz3;McRo42B0wA8HuYf%8SZE{-dQFi2kl)*}asPv?-bkMI z0-10|39_FQP-<0m@gG&2!gq^(_1^v^GTp^3Qnfl#2DwRIWbEM-?GQ=<>XJ z8PMf!j;W=}rU;{QpBCic=kdH+uf07MzFx%_zCF*oT~j=l6TYrd2NcFMJ%H*&`-S&i z%7De$(&&-y?JRI=Xwx~>D5_nX@3n>Jv;O5L8gVG=;lj+8|x84@%-XL9CuMf6#Z#v6vPEJWDgpqM?)yGaRNAaIMK23+Jnf8fDw) z8wb3l2hF1 zb1h7mtfn>V=p_q3!Tmv9G56NxG{_ND`@Q?^j+0#razXZDL_^T3i<*8Cbu6xA^(f&s z%(svT>a&y97OoIU$m%e*OzU^tMTzJUF%-o$-A@py?U?onwE~;QqM83{5`F6gnrfHj zAhwXRY>3&(pNTC5ReoGDFwDC|tUcxc3*g2+4l8{AiV_8n5>K~wJT*C^(i#SGUX38s zsu$Ju$K)Em9)40yAX(QuUM~7rql(wS(T28MK*YF>I;Z0)>^%-Qlww-~zH+{1q4|#^ z00!EbwAIVWdHZC87IZ(s%k<|8R#?|OT#m(ut=D>j`+ac%So;S}0|ubp(bSG>HtuR~ zokP#1@m%vKFEfq}AA!|xo)BfILHw+LY%%MBMsCaECXCNE%CI4MI=c$S%^Ph`1r`U` zr&GH&%DYM=v-{9>d(t+jYe@|b8$yM>zaKVD=xndN8LCgrPb2!Kgk!#a^+;3Z=#9S2 zn#p5_a3I=`2s#*EW1}!WR0D%*_f~os2sI=g_3z}f=|U!pc!*69?H`1u%jamMJR%K( z>G`bs<1OXyGZfaB>O;cI+oqj+LVx{Nf$tD#?E=0h?kiP7us6ZQCxzipHTmN@@$pmfMjJ+!pPc|X$fd!4tN#p^CY5a`=Sq}n2!WWO2MQ|m^AB*w zCNGOofRs+ix%z5o)*jvcWF?pIkC5+I-_^DaWZ~+2EdMLP2zKiDFjq1@ zCuq%vgX9QXY2@43y*oZb2RfT+QTv8uyxnu?%X1qp9UzNOUJZ~xpS;F#v-^+E*(y*8 ze|svUpd$-mKHBtf=ObhxI&Cfv%nt5IEz{bqy>iztli@DIL%AFdeT}=D>lnQ}VKsZO z;xDMv0I+2i;VjH(<1zDkmQcjXD~qw@UMeZpxkN8f!)VZ|Q7!2YF`^aA-%CvA|Cs?v z0U^;l9_Y@{-XT@u*B_4A?CT;ZR2Cz%=6Uq&fu;Um;cT=+GT|0_J;9kE!!9f0sWH7} zryLfG1MuEe43JR;DJ1UQqR%Zcyrw`Wy%b__!-h*?jDj&HW+C2L5VFm|en-Z7-49@R z+%ddK@9vPed>K{A&`Z6i>}B`nBp<8Z{Zy$zUJ)34aD zEjRZn!g*sACG;&RMB!ksVJ-EygV|i0rBdsL7`@@bAuhe z*deV_Bh{-Xft{}_J1h=;Ks!0;|7T-iFLZ`8ri;GbVV_xo5yB#w5W1hmuS9|;IkRKc zg%XXViwjPU`Zi=lY46Ewr%2bR3?*$4Q`9nnlo?pdIl}3nH3BQsRT2|9QW=*GdRZp}{%{u!#$bzi@mcIAT#_4D#WYAnycyAY>qwtexWIgNv zxn!7h;7|a2zrP);rZ}N)yDn{o%fY2$W4TFvEZTg{N6Nc%cji>z`s*xLv|EdV^sWMX z7IDwV39^tE@|fov0dyBUopeOsO=mvp0iy3~1v5K*NayjO^1N$c%W@}!q&O2c&X&mA zAw?yGou8xn8h&FvPwG11L1Iw&N7CyQANQ%(NvbKTFOR1LlIaBPOg}ILJ6C6R#lL1| zdhsW%LdIFp%Y2A%x}LzE$n(?ZA6RTD2Eq13lx2A-Yd=H(;b)~2PN0SY?Pr4oZ;NA% zIJ1;LeCNcWojS6wd;EhJ6*aOjev}@{>;$p<+WdiCSY{d#cP& z!t2aIt)gV|l^SKD6L$<*+6`gLn)oIySi3;2Gdb|B; ze--Dlf0!WUH67{zqoJ1&zmpdsz3t!ZH|e!D{{eIl;BVL9ErpY->#8D|@|I?u<8;%9qcA!{*G$j_fV06y( z+9rF~5t;H$$B>cu=eTy7yxiY-0BwEg3!oEp z*%3b^M<@u(M`0aXL(O`>P@|bjls>o($ZDtkSbgOs@q9CgU_j!#Q#9mO6#dgbm2fJE zQ)@nWd|Nok@8Q$sxaaGNBuD65#=#_`YrjQMDhsA0iwin~+Sc>Kq2-a7lVIZ5WG)qa z>v7k%+l2kW&BLo|Tv0{V0aC+(qqz&|-=>Ek0FL2k5u`>0HOr>RwUMfI0a1?yc51cK z#5O9dC-1I5{9j$_z;brg97}EdW|}uOSvsx=raSmab+{Dq%Inu|tf##QNTXJmRr97{ z$k=jOGdPV%`#Jtkn5&;_=JNIskpS78F~Ghb%_QkDl0Cn>!@pA$FpP5#k6E~0CB749 zz!opGzlSTPs}>p8eF_JA@v2tZjhn<>1~C1c>OtzxCYPiprNntz{Z8pL0*uO6yFTlT zbR}8L3!2H+rBFc~e#a6>oYU!Wto`jv(@;oK5O}tBpN1E|Rrk1TFMsv~+;iWOZmGz* zQ96Ndq@Q2Z8#dCG(BG|>X98s}gt6-^Ro!syuW$S56EsI51oj2WEs-Rf7yHh)qg<9!6^7VPp zckfKQfApZw|B?Yl`S7-?hClxQ@X7{|GFOR=t^ubkAEf~&VN3oNz z0lSJ@gKw-C_lT0HDC_hsyf0ZX-wg3Ef@rDBsR2_X8iY(-0cc$vI>f`-$v93R>iNsk zv`Q9d=%7oMC1`tD5yvGy-CSk7*Fs`Ch^W+AsM8}m7d)Gc1PNr2<*$EWwEa<=dMi!6 z7WQ8wV@;acExao~Z;x$%p|3l=kStkNpO8uF1>KTtA$ajJ4%AyW*CzS{zXDHZSutvD z6m;Ke^7K-~=YZ_yY?=0SC>HMHkSh&)=HV7R>{FIOnkO%A$kupg~)ac6c<&FlP>xDlYmTrWe-H?Ur=%7tqGJBo5~JEn1pfbL1roXal04~GtQo9 zv7C3O4t50hD?D=hLqbbMo!3= zM!J)vjSH;ra;I#1d6&}@2NuWBEzk17fC?myZvVd=0PLLD@~Rh*b6W;tKf^MHP}Om0 zFL%+iaakHsdVV}Ol|1kICS-_M>K93~Pwt#>HjsZ|s$NTb6Bpv2D%Wc`5(j3@8>ZK~ z*OwH<8X~awU#O8PrV{-V_B6^A;%s?8-D1SOE0)8{Gjk1RVeMV~(7m6EF7A8(P%+$Q z2p~XxE*m4<*|ucp4(oWck6#8NU1L~G3fg;Pw%oA`glQh z3oLUE{bqk5+kg$AYxkUpSOb)j7Oe(ZU}ex;SdqvdHS7~gpZhC6 zxr@g#JkvTpmhh_LIb|8RxJsn`fAW_Hva0->*Bi%6r~@;Qcgm?^rj$`@m~fDCujn?W ztiSqZOBBQQos4XI2$j{U+yvO5Hz@anL1yoi#}RHuFv_tqHlZ!RN!2k(U(} z&7-{hk;DVvdy)$0SFOT_@i3ha5+zSP)9U||Ji!>&4bBxjJd`s2#v`mdfigp`mfZWB z5BNP9g+1)IDK~Qb`5M;nr__WA<%fN0jZT4w37h>LnRf^E#Exj0y9*mhMu*<-v<|!3 zxi&zPI*3v;SvPeCb!(c&LZ~OlBPkkqa7S4T8T$1#xamqcpntm;mPSyCCX}idJAoEf z>he2e6Fel^MzKtQkS4E0rkCyk8AHB@RW_$=$`w=MBFf0vfr=;ysp9>{UnGB-CPw8r zb}e|o9wmyoerh8Rxw{=Z6@UZcQxriLSXPD99k+lfz>5-nh>s{9oXfV69}lbPmS|M6 z1{D;q?RWhaq1Jw_-6_7JbEV2*zM@ytB6d2tX0L-aG)Y7?NQwd4O;F_3+nAm97;FZG z#yb-6ALSa(sU_fc;%R`f3ieorsYCsI#nV1oGHfJvYi#=h+VA=GIE{(m4ZQuqV?5S0 zTs>il4HEL){~_GVqAnsnDbAD$9T7-52APSC-&XXr#1AzQm)sG8{CE1i@qwY7$~&Bt z0LpO)_CViPC)4)uyDsGJ9-E4!3CR9&%prF$&}kqTS`1zZbt9oK4Y|lvU$Ga0ls2I< zYz_t|!#0*y2?scui!TSIIYB7#7fEDY5F6BOY#AbyOryImF%|p5c|5mriSx?VB02Ci z)D#1<7}u4)wuO~U>l=fT>;z!?UsdGaDjgH%Kir~}?Ka7oPms`n#=!YaZ8%4$k%JYI zym6tSxXSGA--%PQEw%-}bh)wpk_oQF2BYcku?Qa&u_Om0^giO^6ww!|ChxXk1O(qd zx;=-_QXuH#p_x6{53%y>xmc{a)5y0<>k($Fn*c+Eb_alT=C}E^?P^|4#rt$Ue6o4g zlRY@|*xq!3yWud(rG~7yfWuT{eUU)JP~0>W;%+Tym;m-J>yia+pboXq_?lhSw5PJo zhdU%zuO;2DsEg+ZZPDer&3zaPi8p0&ey6>I1GNHYXsVV+Jo3DFsK|Q2JRW-2QwImW z89#;s4vQfuPLT~f%xRr5C7nKuQ5+}w#rbJ7r|VT`UDP$~h9V6`t+M~iSBU!%7B_XMmvY-`Dy+wU&Sqc1J#E4QQkGR=FOQyzfh3JF>5cOxC( zAju6q^czo>Eo4^}R-M#6_UNZqU8T=)UY5w#)3;ILhsOUQ;{xWgQxpa}2W{!BfGt7-@eF0S~b&p!YZRos7fuylLyep}n z;=@qlZ^QT7Myv`A5^k~~=Z&j>C65~>w2AK0B@yB1%9q?p)8r&YebH{ec%(BgtUHt7 zi3iwAo`!&ovH!^GkN!F+wEriF5^o@&Oa1NRfe_hWA11m$l|LZ=h;CH^U*LFDroLH2 zk!0Y637TtB`ZbtCSvqc-H{nSEgwAamOwU>f=YqOno72kOcg-K>wIFfaZTji5PIV5H zYkm4LT~JRp2E+OWBf*l(kTf@At`C&=o&m4yyFeXy`ReLGDAi5yul5jWM^f+g7BG+Y zJ9yblmix7E0JB64hAkhe;t*Ps!a3#787tb})$aDbTdf!&w}sbD$hx8)1@hbAWJ7(5 z(3%n5B-!xw5%oeHEK-mCR|dae`lr54I+qg(INj!ar6Tp|7OIB}C)LBrv2g%D>_d4N zM3YY^#?j%#4<{q=XiF-Yn9Qjc2bDs-WRW>V?2_${wM z*Zb+9-BhO$#6aPlcL8_qHA*#b^XdG1tKCpV!cxJ#x$BqmRoelR!dN&caYZIG5&k_> zS@$nIVxwDRN7>T&Q~gEf{l(O{zS1Z*miIOO?DO_V3u!_NyNZs-bHK&MKXzB;oGD5$ zI)hw9TxX$0g(}-V?rtG95cf-eDYkXM^?LEQKIT|AATcf6;;I(0JiQ{vK=KXN zdkpXHya_`4B8r>Ws&v1c=l3oUjH>5ppT4YlyH9hA04^uwY%f7l^Y(n>ngD#Pz1=HX zxLS`Fz8(QNTw8B7PSzOTM8JTP6~~on*FsPDGDP`dCvnyGaof$h_6zhWv5fkgkqmM~GBu0vqvlxZ?81$8-a%TkFi1S-<)daF>$V!qJch4< z6l!lzoAU))GutqW79tPrI9Bv)8Hwm#BO%&asXxbPN~Xb+Ns1s%+8B|XD}H99kSG|= z7eMD5gv5cVHq+wFWbxDb)6!o*JUXgB8*6?=rx(nbPPR@7T<1G8`LUW{|JX#*avw%+ z^WSKpayav+*&veqOruvc>z^Jg`aljmrxQk;6i3xco<372c`PkvKq@Jv0$cjTv~3hu z?-<=NVKMu!!lkGBKq%L`)~6!?Zp5-V;XeDmT>#NQ$QVvZ=kui)!>EMm?4!LXMa*Va zJSmET#?{CVHh=M+O>MZ5xe0;a(&muY+RwUWy}(YVGo~jHK>;pLJd~ zdW+d`dCl$@@aQRxuLkGv7awEUB_8}Z6|9S#(cEnI7vP|V5XWXPkJ?6ebxM?cV$D!9HKwC#*0qM zAK=f*C!|MJ__H^XFN|fa@7KxeAX_?dUu6T)Sb@3sDM5i1&_PQO+3&*sW4lTa?B%4) z-suOAq7@B)sgW=I=!guBw9^w6$3bdKw{i0hGe{N#`O5{?g!O<8*hsIi2sx0xdbmGFYXNh)>|q)=#=Ng;*nt`8xMV4t611-MmgBSGhq_=v z(O;4eVC za0jkBl(z$-Rj~d@SSc!Kh|B>4Dqgh5=WVvvt^}u?A8Y#;LfjVRr9CsG>_|+s9BqQI zp>PE@4Nt?bMYR6(C6q_xACrxy1r~K`3J;I$dX^F9hmC}4=S%8+IZG27JT6L)+#jKS z)d^;sH^>qla`Qr5eRg z=(672+Y#IJ-*tq)VEK{cv#d{pv9!9wcb1k`syldAoI|54vnaj#rrQ9=yFHNMQ*hj5 z6z`ZtG`hM77^+D?gflk-SBbqh15xZfJX2RyPqQ)Z9z0Y}Hg5i@wkRXzzrEdm-vx#7njDA>LA8VGPM8L^Clz@*!&*NNBUUo-Q>@3Eb$B4OzCQ6c zlh$D*B7`(mZGU#hU^?+Mji#9E9h<1f`?hY3hd4~HzRc_l5W1dy>2>@B@~OLmxtWjO z++fe3laJd@zFpb$tY}N@+oUr2(U)i^F8`ULio{JN?RSjQ%HA1)eD?h^(RX!|2*RF*j?! zI=iVI4Z{lyn?*zGK$}2E0s}d%Yf}yaq)oi1pT)li=}+~b&vCSWC!&s>Fz>PmXHObE zBFYQeq?$G!ZuTO55oTy1clFahr3?+Qthmx_nZ1^eL67|wIJz)lEcubbgc6;7Eoavb zm&={65Y`r7dS_=YM^hI1{^M24hsP4rEKm(39_yH%0+X<}GIM&-Y{NgFUp6&cRh=HJecn+b=kdAt0Vmt1)cf z^_NLXR;@KKVPULAQ9JExRTbv+>ZzteYyEcE7xpVz5neqOjicZu$-ez z8nl$B7S|slNJ23=td{6&rX?+hUe|%X`1fS|eUt{h#hYmYf#Ao=oF9~55@d@0$mssU zL4d7k?9UO|{u$tq4BhR3FC3EP$3Cj6kYTrynR#X&ONw$!7@IB3^`zgfTea2D#$ox@ zc$D4vaiq>w#{<*M3w{Q%qa|~#?#VaiCH0k`{)@Q6L^*xR<-OXn#eq`U$K&VQ?=1h)88X{j_@UhO9e&(0thlNhFzo~OzGUeDsHdEl)%asl-YrFFd zGQJ9)0Os)Uk7R3aSBO&7!3(v$-X8ClVk_0Uc^3y?7=Lpw2So-vM3>Zq6DUEA~0Qh zgDw;@4Qbz|*gy-g5+5=|XOtg7`~-^f{#?=8>C%aGiy548ADh3#AQE}SCh#ll;_MW^ z#>08E;#mT`B^1K=h#WriiQ6rC^!4XvWLM1I;BRLdS=}h~LwtiHmmwQ5B3#a(L}H#F zw(E86HJ=8a=s91m8WRV0h#sRhRG2XK5PufQE!UCtv@ysO4c4L0%8J6DGkFh)YTbmH z=*`A{v4kAIc2^aP8!ygz`>Aedj#!r>)U|#%(JbITjn+_bP;b*wEz(%+-AtgqCc%x@ zxfiZi4&l0CF;m6fn#x)bf#l>1mTE@kR&_#f$f$!N#yRW=5v;l)4YVtoWWg>iDZ0V9=@^?U3#f6uC;{xh?+x0db33i*$F3(?nn=f-cIV<2 zPl&&!0ytT0d5y{Zeqqm3+f!?(Xi}`S(FfRJ0tu6n4@3lK{|$0F$$2hb?ZQ0Nt_qI= zJsG!S0aHQ2o1MkJ3M+bC2G*-aR|^)R7D_gJa=dXDV^QK7s=r3vEh-jCcT8k6yq+e& zZxJUJ;A2m!&r;_bS+qx*xoi&> z{fAV#Q(hvT->SRlhSFkZ$g!>vh8bF}l2sMUdSir_bw!hX$wzkiBg7NJb+ez z&qQjZm7&a1_d1`X)}{5nM9Yef_Du?@OHs!2ghIb}uW`hr;XvzWjOHr8s{N)J{s`E( zVP$0whrt~@TXJH*;@e06#R`QpX8s~`g_e;?J{0HV6(jy5&qxiZeCSL>iCz#|+U269 z1Z`HmWup&g#R=CI(MUhmsivbVd4n`AB45#CG#o}EM_+l3yn#*ZreM)*#bmMl>mq}X^N;< zEVPVLG~V2o_hKy|C>E#HilHc;l>BhLCi$)S?!2ZPg?vjl1-<8G!QG4kl;2$YVZ{MU z`4&yeUb#(~o&nHIi#Uo;!-0Zs;20HNzCHR0dU>kA#Z9{Z4og(jf{RwLbv zcpvRph!JldaLIVzqg7^pW=9=N^T0r9bq=-nEu}ft(r4Fy$gNc}G8fmf0I|dyVAfSn zP1A%srW6CrAQn(1=?PI-j^I9Dmli;G6J)Ex#2%fd_qx(d?^PEadt%^xSy8LZtcd4Q z(@@9yz}Apf_qAdppS0OQ$8^6NJ2=Fs)L87y4$gx)Zt)^Qu_G`C(}gTd zStDXFUO73cUtrl;^VmAgj<)dy!An7COn``6)*nKJ%T9pDi>&(t8}6@=phZPX^&*M% zESNNarQ7p^vFONnRHa(^?Y$~yFuX5#P&d1L=$aNXBe~wPId&oM*(E5hS1&z`czqP) zWN2`7>KJ4&s&)_eiTLwRx+I@|)X-+#p#-;+2#Zqa@rM9(ssQyV3 z`9581SGUw(+?vqIHBepwqM~6|tzUB^5JPgxx7iWER!BR#4;M*t_<~< z7utk$vJ%WpQlW|V-T9kgTRJ}yIH9ER$aq4+e7ivBeI*RY=MAdb2C`^_$Pkpjzrce&}UjU#JZnN{+flm%)RTk6b|Lnd2MxubkrpcRlqPX_Z6iu_0)CLHRN7cVo?YqfjE_Rf(9b9f_gfs1jxxuS<*7i$2%y2i?5sVj>s*u z0ZoFTCFeiC?2mM=SWGrwi@Sw%>=5DlaKU=@>>u3i?v~*NgdjV_63oCk)xQg6f2CIP zk*k}kx z*tGComId82LQa3R7aP4b^>#HYMNh}3nE1yH>U)fFPBKyR^?z0+_|c5phHc#7B3}1qX$rWVjYywQKadkXC)M>nZUji>zdeK9 zAaaqtG5Wz7La3Pni>*e4)RX)4oe{eOjZdR3Utc@WbsjaJG@-IVFD6wBx{kaGTXg8C zt04!2HSLSM5|fXOP?7CqajTTvZj_pB$y0Y*;Rg6CcN2uI24rFvTHxcd#TCy5IIg(- z;dx`xW#C-l0kinZw+HPA=ZzEs1)<| zX-hIno@65#H=_R=Eb*qte6gJ*;Q1w*5IlS+4_%V#sPmt+K+5& z1TY9)FQ~KRZ(!KK+V@?8L`UJ{scT;Uo0Vi181_8jxt2hAoyg(!XjV@d+OKeW&ED^divtzpWW~NkWs-rlPABaq+vQm(J9$H}^`} zq@m5o;Ac7KpB)WlCi3>`qT*OcZpF5(N8G^J+t&O4shq4yGFusc6CcICi#86h7ao1D z88H^K>zGQHUaU)5Tp&Qq!y8Kzaj%h8JKh583T`%|Kq$I60Q^BnhZZnd@fCFu4?Y%q zA7}a5ayceA4u||bG~FWFaLO92SWmvKk~IoI`xqtQofC`y`>H}O0tPUl`wFJJNq&?K zqj;?;Aqe@5FJ3dUKXy|6Gbj-a`Q^FaM9yIEKl4EI#o`(U9Ky^0TZo#aF=H3GtZNzp zIC%vs7dHnR+b||Zp_xB1m`i?2gh!|uh-&Nj%08Ey;c!0qD7BF`5)}8}RO#syFaEK? zxA8#-%B~>PdF*Bl-Ch{QHwy*b78nN)W4gcY-+wQl%W41J&%lb;l@;A7zE~1*M|TXN zbP70ou&nE4y5%;}m@c8``T)&)Y}gwCa)NmNI&aluY9)@x{34yE4l%o1sBZ3X0 zOBOY$2TL($?w$ikgl9xgasAR8yZSWUwI7q<3Qe9xU9q;yHz$WR{(i$SBXDokX~n{U z=pv>3)Twz?Fxqfoe@6$lnd8ya-(Dn1q?$&)Vqo zp9Byc-3Xp|*Pm|+tX+o1Aek@1VnO0!M~I6DFE(-7_en`|%%vLdzx}aO4w3x1slu+s z)xKSxx}Tc%+cYnHu#FMho~wt8v*5hql2!Do9*}Fjea*ywaF0(uP36D!j_XiEl-Cke z|F4bnj%s4-+PH!eMT(Rwf+3)E0qGEmAktgt5K5>bBE5tTc_|{zAjQy$QbeRi5J7tH zARrJxinItwH*o3ifcLKTuFo%jP4>)MCpj}`pE+kgzsE3DpBi$_IJLawP`1|iF1E)e zTVOH&{Zoi6y2NsH&ArG`pyw67Sz#ovW&tAiZvh(H<23>)KIM2$5GA(dCN`7)&9tK% zoh+PQe=(uM>oo^ZhPrFzL9Xq^xm)$aQ4io>1}Fx8VJ)Lml^*Jqpn_hB2S~EP$t5d!6Z}WK%ImanduyK_462)(#8UY% z*P;qhq_%?-g0(+|OB*v&VBKlaBW2=hgymRVjy{xDgb-ol+a;JQ$7(4=uGA+W$|q2T#MPlc4ljB5XFqLMpkAx{-;r>CIiMYx@+ zd-D>B4rp3mAj-J{V?ajpDV^8n^_=QWcI;+(BKRk21L>~MARHO6@D8TmtE3|wx-e0z z0Dwj+n`xwkX?&=;fyxjTTj{Y zxV^YDB8yl!$$IFn`cOeQK*jbhqq^sPj7p?(Pz+oC zsQYyUN}C6*C5~=gR@4>Yl=EoR3O~t3x`0q`v@{%ExTO~lq$X#M2#Vc}()DBdu1FR5 zsip#R#V5<#R7xo1I!uua-3JYZ!RWdJ0ljM$4;0O9oZw!p1t)zq^jv+KHiZ$o&^ zGrpyF^l!esrQ9)SdNsS!s+2^(I@&b#w*=+fFC!Ieyh0_ z%mG_?y*+JY%3#rG2p1<_l|pczrfkR11C-yH=fhSX-hVSN0N}#&;r73gMyxrUXu#k zT(;dC?vP5Ju7%i8wXkuKMp6<*nFh@MT9oyEr4sBnOoV+r`Vp0l{{On^RK1DU779Jq&pIxDI_Tk3dp5$Q^0(;<^{zImNTWhzVFXA z8^KC*Yh=p>KK>RI>@C?DU6Ljct5S(q?YwU|um5Vh=iiGuHZdusc%~zZRF7=VXsLM_ z7+5M9>B?j#YWsTqLm&6QY9PzJ*M_dL{0}-joiYM~X;m8Tq%#p;`m^vo!K%=8e1g=7 zvSUk@FFWZm*(bC_)3(pGsZFKIXE|$3n-sy=fU6>yC=g6X@2F#WN}|0$K)7imLoiUpBqd%!Yueao6RL!{X5V;;ru z)XSNZzmy?IWJZaJ)vc#xGpl`k?qvvk{s$Zxm~6kmJ~`wpEwa4}XmXexJ!h%7ZfH0a1o-ChpAise)MgkI^XeWuJXj|7nHFfR zd+x*K7;MMw=E66Qx%57r(Cc3hlq@6Y9OQ+L+14LUM^;#*qKvY|`Sbjmohv5-Ik)Iz zD}b`$^RJeCE_BN$M-y=_(Ju;yZZB$tnA#XeIM)T*QN0alZ-4Vc1$g*pWCm0Kz4_Ey zFlQmdiho5~;z-Qmz#1wb*yvn|>OP9t=5gF|>v^O&t>J6>BETw#cX|(FHmB%rQHsDM?e!xttqIvrEi}oz z=dLk=R+x=6Ns>PxG#uF6Ynr^pN_}%3Ka-U)a*bKI|7W(D-;n_mpt|WkA8cmowSHCV zcGNOc<6{b#ZDh-B?$L?iXZq<$-DmSv1B~RNgi%2w+od(Vepx|kUm%CV?K z<7{BwtN3-_n>@AGVQRU{RU zBevn4HO`U%&{wKIv%p2|!a8cu1r~5$=^8g17^)FNCT<5wLH=PNI_gFDq&?6+sF1PrbDWwX1C4fn2rdhc+@_-#;f)J3?7?(% zU-ZdgF021mXZWBiBRc6?QJp)2$xoxD`qBt#UmNORFKJ*7WNN0zv0A{HZ)0_Rw{ax| zKosu?e3JWZ#nZr-NUTiTljX9X!)3)3HsK#IGTQM=?7T&Zvq%mI(2pAQ0Jc|Dd~&w} zxHsjeI&S=J^^c}bF6L>B-H}0#v?N7W`#K!w zwNDQx&I9blk$Y3m!}IR?WElwi)vcO8p3!zh3K7=9VOcHO;*(m0iJJQvd~T{m1@R^2 zl9v;~*3{qv`>9+RLUjo+8lt(3;-)M);5yZ6HKJn-T2v&y^88Eh6jK_qD7+rIJSpw{ zQ4uKVNdnj9&!1nK)xn}#F;y`IHGkg&bb!>ZWTb$Y*fl*NwLH?&!>u9Pm_?)oQ^Z>H zT&&u|cAvv18LnSll>j!t(u%8BnHy=%#Y`XAzBN#$`CIrPHS55jgDpQlbNNXoou1;{ z^@*Qi5*wBruHZ^@wcu=ptBxyWKGv}Ldl#&8!U+I#?e5=Q$$>g|5ZWyZ`I$T7n|0g* zO%*`rP>kAe!p9^$Z6G6Fs*8bMk;f~j@-C9McDqdY37y1|+l@*M>X`0Di;#=PD7uAVe&l*w?gDzy} z7ncMTc~YVtFVR#Ga>KR?e12+%*pK7hj#2BQ$40h0aSp7wY~XClurj~UO#D~j*1(@|t_ajsoSX0_T(#jxr6#r?(Io#?Dtg0`=GcAZZM1zH|oIE z@2sxw7(VG_sGTMg9q2VirAF0{iAF>vt~Sw(zaX7Uw5j;jm=?dlM+6Kvvg8Flbq0Fz z=wgP^DNn};Ko578lmpBf?9^(U_iv~~B0lGF(?X5jno*9y$pwPJ7-b0Rn;=Zl4}7u` z8|Cn=l}y}$BVM~WucDBOW!qfbtU>*y)w_t^P6QcSVYjCL*KT;G>0{0Ax+ORdR_byl zTz`Jtj+KhkqZtyDC4)@-u;B^$&F>%j<=?qhlTlLX3{!$oyR81QXh|7u*t^r-;ls&Q zKeVHH=5=5k=;44W$FvSfj^P37%8 z#jDtyz|EKZ)>`-d&`V&7WFvtLFktW1);aHnQ|&~RL;b9PMnkxRQap<1x4G!8LlRP3 zDN9$YH6^QC4HHr+hDux0RtAbR!MSdXg_<_Y)50%CDn;V$&WL?%o1f24C7Xd`F-rRGJZkTrRdmpeGm3~+*cfKlg0l6}1 zmt@_maqVAYo3bp8Ecnk;xVa1!i_=5~Y%*w*XMrM&#_|24_r}6&oeCT>K70M(QubGI zNe71@-JOD37PB90*L#OxKqsSaVLO^QW}jt4!U|Kfv+HPh@5a@7K{$W}vs~ukVd}fy zu1s+P>~8|*S-+nrvjL8Nn&HkIcoNi|9zP{uwslBr8>lFh%ySR%E>?E2&UeGWO6k>( z2^|mmha!bCK9sg`6p*mT47)u1&C}=bqJvQT=mxhgoHxqKgW1e(GOh^B`Qu^mm$*Xk zDEUBb3HRmG9QEgn6S^daDkD@PUUL+Q##?19Jr!R$`yA z{beQ@9QX!jp({Yt=4Kd8T{&g#*`70bJZ^SXt^>h-z zsYljyP8U!|$E7DrXZwrB@+$jfmZy7hEV}pkDhrj>fZ5ab{qwd$&vox(rNZZiSp}a) zs1Gl$ot}PgU@4*dEx^AG4CTk>44iqlnQJh5{nI2<$!=UQ*+*N2Zn<)M z9cc;s1rb5gsVZ-yB$aGU=65Vc9-cbP(8{lRM%kP&MOji2H*|^HyI=fyCLnWuJx={A z6O_k^kO$J&Q5JvO?_QV)@mt=x(A#*&OuDSR!G?1(`%)tZNumPAEun(7Mj1-851ibt za)a6;)JuCeY-17z=|jC0_CT9l)$%il?bP@&HV-jDN+yZ*y^N;PvFV@ROuz4b`Gzg6 z+ss_See?+?Y#eSer|VkyDRhb-b)fhf_B8%7dyG#`l z@16EWOmD7ixYx8$WcAzI-18N)oI`l5*rx#Q(_xj)4U3$EI-2b|I8g=LG(N+Y%UGMt z1&d1|AY}c9`gbVcOJDRLTEL4rfibcvfCs$sp+a`x!TAiZidvEvZvKML|1Q5^e|ju&}VG#2+zH z(Lb*fm5@j@dMt?&k`VsC*Nce?Oa2iP6Bhp?CN2#9BPJmt@<$9RB=l!YT<|Q${=Yw7 bOj7XAh9$+IiNf~ER|O@=IXINGmC64LZH^Hi diff --git a/vignettes/articles/Ravel_2011_16S_BV.Rmd b/vignettes/articles/Ravel_2011_16S_BV.Rmd index ce0ec87..a936252 100644 --- a/vignettes/articles/Ravel_2011_16S_BV.Rmd +++ b/vignettes/articles/Ravel_2011_16S_BV.Rmd @@ -23,17 +23,13 @@ library(ggpubr) library(tidySummarizedExperiment) ``` -# Introduction - In this vignette, several differential abundance (DA) methods will be compared using the *Ravel\_2011\_16S\_BV* dataset. Lactobacillus is expected to -be enriched in healthy vagina (HV) samples and other taxa, such as Gardnerella -and Prevotella, are expected to be more abundant or enriched in -bacterial vaginosis (BV) samples. - -# Data +be enriched in healthy vagina (HV) samples. BV-associated taxa, +such as Gardnerella and Prevotella, are expected to be more abundant or +enriched in bacterial vaginosis (BV) samples. -## Import, summarize by genus, and filter +## Data + Select equal number of samples per ethnicity group. This was based on the minimum number of samples in an ethnicity. @@ -45,12 +41,20 @@ conditions_col <- 'study_condition' conditions <- c(condB = 'healthy', condA = 'bacterial_vaginosis') tse <- getBenchmarkData(dat_name, dryrun = FALSE)[[1]] - -## Select equal number of samples per ethnicity group col_data <- tse |> colData() |> - as.data.frame() |> + as.data.frame() |> dplyr::filter(study_condition %in% conditions) +col_data |> + summarise( + .by = c("ethnicity", "nugent_score_category", "study_condition"), + range = paste0(min(nugent_score), "-", max(nugent_score)), + n = n() + ) |> + arrange(ethnicity, study_condition, n) +``` + +```{r} row_names_list <- col_data |> {\(y) split(y, factor(y$ethnicity))}() |> {\(y) map(y, ~split(.x, .x$study_condition))}() |> @@ -63,6 +67,7 @@ set.seed(4567) select_samples <- row_names_list |> {\(y) map(y, ~ sample(.x, min_n, replace = FALSE))}() |> unlist(use.names = FALSE) +# select_samples <- unique(unlist(row_names_list)) tse_subset <- tse[, select_samples] ## Summarize by genus @@ -81,7 +86,19 @@ colData(tse_genus)$study_condition <- tse_genus ``` -## Get prior info +```{r} +tse_genus |> + colData() |> + as_tibble() |> + summarise( + .by = c("ethnicity", "nugent_score_category", "study_condition"), + range = paste0(min(nugent_score), "-", max(nugent_score)), + n = n() + ) |> + arrange(ethnicity, study_condition, n) +``` + +## Prior info - biological annotations ```{r prior info} prior_info <- tse_genus |> @@ -98,30 +115,25 @@ prior_info <- tse_genus |> head(prior_info) ``` -## Convert to phyloseq +## DA analysis + +Convert to phyloseq ```{r convert to phyloseq} -ps <- makePhyloseqFromTreeSummarizedExperiment(tse_genus) -sample_data(ps)[[conditions_col]] <- - factor(sample_data(ps)[[conditions_col]], levels = conditions) +ps <- convertToPhyloseq(tse_genus) +sample_data(ps)[[conditions_col]] <- factor( + sample_data(ps)[[conditions_col]], levels = conditions +) ps ``` -# Benchdamic workflow - -## Set DA methods - -```{r} -## Normalization methods supported in benchdamic +```{r normalization, weights, DA methods} norm_methods <- set_norm_list() -# norm_methods <- norm_methods[names(norm_methods) != "norm_CSS"] ps <- runNormalizations(norm_methods, ps, verbose = FALSE) zw <- weights_ZINB(ps, design = conditions_col) DA_methods <- set_DA_methods_list(conditions_col, conditions) -## The following chunk of code was written for compatibility with -## a more recent version of Seurat implemented in benchdamic for (i in seq_along(DA_methods)) { if (grepl("Seurat", names(DA_methods)[i])) { names(DA_methods[[i]]$contrast) <- NULL @@ -129,16 +141,12 @@ for (i in seq_along(DA_methods)) { next } } -# These methods throw an error, so they must be removed -# DA_methods <- DA_methods[!names(DA_methods) == 'DA_ALDEx2.1'] -# DA_methods <- DA_methods[!names(DA_methods) == 'DA_corncob.1'] -# DA_methods <- DA_methods[!names(DA_methods) == 'DA_edgeR.1'] names(DA_methods) ``` -## Run DA methods +Run DA analysis: -```{r, warning=FALSE, message=FALSE} +```{r run DA, warning=FALSE, message=FALSE} tim <- system.time({ DA_output <- vector("list", length(DA_methods)) for (i in seq_along(DA_output)) { @@ -156,53 +164,43 @@ tim <- system.time({ tim ``` -# Enrichment - -Get direction - -```{r get direction (effect size)} -direction <- get_direction_cols(DA_output, conditions_col, conditions) -head(direction) -``` - - - -```{r} -hist(abs(DA_output$lefse.CLR$statInfo$LDA_scores)) -``` - +## Enrichment -```{r} -hist(abs(DA_output$lefse.TSS$statInfo$LDA_scores)) -``` +We need a threshold for DA for lefse-CLR +(It can't be the same as when using lefse-TSS): - -```{r} +```{r median lefser} c( lefse.TSS = median(DA_output$lefse.TSS$statInfo$abs_score), lefse.CLR = median(DA_output$lefse.CLR$statInfo$abs_score) ) ``` - Create some variables for selecting and ranking differentially abundant features: -```{r} +```{r variables for enrichment funs} +direction <- get_direction_cols(DA_output, conditions_col, conditions) + adjThr<- rep(0.1, length(DA_output)) names(adjThr) <- names(DA_output) esThr <- rep(0, length(DA_output)) names(esThr) <- names(DA_output) + +## Lefse needs an effect size threshold because the +## method requires both LDA and p-value. +## Otherwise, it would only be a wilcox-test esThr[grep("lefse.TSS", names(esThr))] <- 2 esThr[grep("lefse.CLR", names(esThr))] <- 0.15 +## Lefse will be based on LDA rather than p-value slotV <- ifelse(grepl("lefse", names(DA_output)), "statInfo", "pValMat") colNameV <- ifelse(grepl("lefse", names(DA_output)), "LDA_scores", "adjP") typeV <- ifelse(grepl("lefse", names(DA_output)), "logfc", "pvalue") ``` -Create enrichment. Threshold values is based on adjusted p-values +Create enrichment: ```{r enrichment} enrichment <- createEnrichment( @@ -218,37 +216,20 @@ enrichment <- createEnrichment( alternative = "greater", verbose = FALSE ) - -# enrichment <- createEnrichment( -# object = DA_output, -# priorKnowledge = prior_info, -# enrichmentCol = "taxon_annotation", -# namesCol = "taxon_name", -# slot = "pValMat", colName = "adjP", type = "pvalue", -# direction = direction, -# threshold_pvalue = 0.1, -# threshold_logfc = 0, -# top = NULL, -# alternative = "greater", -# verbose = FALSE -# ) - ``` Create enrichment summary: -```{r} +```{r enrichment summary} enrichmentSummary <- purrr::map(enrichment, ~ { .x$summaries |> purrr::map(function(x) { - pos <- which(colnames(x) != "pvalue") x |> tibble::rownames_to_column(var = "direction") |> tidyr::pivot_longer( names_to = "annotation", values_to = "n", cols = 2 ) - }) |> dplyr::bind_rows() |> dplyr::relocate(pvalue) @@ -273,7 +254,7 @@ enrichmentSummary <- purrr::map(enrichment, ~ { Create enrichment plot: -```{r} +```{r enrichment plot} enPlot <- enrichmentSummary |> left_join(get_meth_class(), by = "method") |> mutate( @@ -319,43 +300,11 @@ enPlot <- enrichmentSummary |> ) ``` - - - - - - - -```{r enrichment plot, warning=FALSE, message=FALSE} -# Create enrichment plot -## Not a plot. This is a data.frame that should be used as input for -## the plot_enrichment2 function. -# enrich_plot <- plot_enrichment( -# enrichment = enrichment, -# enrichment_col = "taxon_annotation", -# levels_to_plot = c("hv-associated", "bv-associated"), -# conditions = c(condB = "HV", condA = "BV") -# ) - -## The actual plot. -# enrich_plot2 <- plot_enrichment_2( -# enrich_plot, -# dir = c(up = 'BV', down = 'HV') -# ) + -# theme( -# axis.title = element_text(size = 17), -# axis.text = element_text(size = 15), -# legend.text = element_text(size = 13), -# strip.text = element_text(size = 17) -# ) -# enrich_plot2 -``` - -# Plot putative true positves and true negatives ratio +## Plot putative true positves and true negatives ratio Create 'positives' object. No thresholds were added. -```{r putative positives} +```{r putative positives data} positives <- map(1:length(DA_output), .f = function(i) { positives <- createPositives( object = DA_output[i], @@ -372,38 +321,9 @@ positives <- map(1:length(DA_output), .f = function(i) { FP = list(c("DOWN Abundant", "bv-associated"), c("UP Abundant", "hv-associated")) ) |> left_join(get_meth_class(), by = 'method') -}) |> bind_rows() - - - - - - - - -# positives <- createPositives( -# object = DA_output, -# priorKnowledge = prior_info, -# enrichmentCol = "taxon_annotation", namesCol = "taxon_name", -# slot = "pValMat", colName = "rawP", type = "pvalue", -# direction = direction, -# threshold_pvalue = 1, -# threshold_logfc = 0, -# top = seq.int(from = 0, to = 20, by = 2), -# alternative = "greater", -# verbose = FALSE, -# TP = list(c("DOWN Abundant", "hv-associated"), c("UP Abundant", "bv-associated")), -# FP = list(c("DOWN Abundant", "bv-associated"), c("UP Abundant", "hv-associated")) -# ) |> -# left_join(get_meth_class(), by = 'method') |> -# relocate(method_class) -``` - -Create putative positives plot - -```{r plot positves, fig.width=15, fig.height=10} -positives <- positives |> - mutate( +}) |> + bind_rows() |> + mutate( base_method = case_when( grepl("lefse", base_method) ~ sub("lefse", "LEfSe", base_method), grepl("wilcox", base_method) ~ sub("wilcox", "Wilcox", base_method), @@ -415,11 +335,15 @@ positives <- positives |> TRUE ~ method ) ) +``` + +Create "positives" plot: + +```{r plot positves, fig.width=15, fig.height=10} vec <- positives$color names(vec) <- positives$base_method posPlot <- positives |> - # mutate(diff = jitter(TP - FP, amount = 1.5, factor = 2)) |> mutate(diff = TP - FP) |> ggplot(aes(top, diff)) + geom_line( @@ -441,75 +365,37 @@ posPlot <- positives |> scale_color_manual(values = vec, name = "Base method") + theme_minimal() + theme(legend.position = "bottom") -# plots <- plot_positives(positives) |> -# map( ~ { -# .x + -# theme( -# axis.title = element_text(size = 17), -# axis.text = element_text(size = 15), -# legend.text = element_text(size = 13), -# strip.text = element_text(size = 17) -# ) -# }) -# k <- grid.arrange(grobs = plots, ncol = 3) ``` -```{r, fig.width=15, fig.height=15, warning=FALSE, eval=FALSE} -# ePlot <- ggarrange( -# enrich_plot2, k, ncol = 1, labels = c("a)", "b)"), heights = c(8, 10) -# ) +Combine enrichment and "positives" plot: -``` - - -```{r, fig.height=9, fig.width=10} +```{r combine enrichment and pos plots, fig.height=9, fig.width=10} pp <- ggarrange( plotlist = list(enPlot, posPlot), ncol = 1, heights = c(1.5, 1) ) pp ``` - - - -```{r} -# ggsave( -# filename = "Figure2.pdf", plot = ePlot, -# dpi = 300, height = 15, width = 15 -# ) +```{r save plots, eval=FALSE, echo=FALSE} ggsave( filename = "Figure2.pdf", plot = pp, dpi = 300, height = 9, width = 11 ) ``` -# Perform DA with lefse, Wilcox, and ZINQ-Cauchy manually +# Perform DA with lefse and Wilcox ```{r manual transformations} tssFun <- function(x) { (x) / sum(x) * 1e6 } + clrFun <- function(x) { log(x / exp(mean(log(x)))) } - -## Relative abundance (TSS - total sum scaling) assay(tse_genus, "TSS") <- apply(assay(tse_genus, "counts") + 1, 2, tssFun) assay(tse_genus, "CLR") <- apply(assay(tse_genus, "counts") + 1, 2, clrFun) -## No need for pseudocount in the next line -assay(tse_genus, "TSS + CLR") <- apply(assay(tse_genus, "TSS"), 2, clrFun) - -## CLR transform -# assay(tse_genus, 'CLR') <- apply(assay(tse_genus), 2, function(x) { -# log((x + 1) / exp(mean(log(x + 1)))) -# }) - -## Relative abundance + CLR transform -# assay(tse_genus, 'TSS + CLR') <- apply(assay(tse_genus, 'TSS'), 2, function(x) { -# # x / exp(mean(log(x))) -# log(x / exp(mean(log(x)))) -# }) data <- tse_genus |> as_tibble() |> @@ -532,8 +418,9 @@ calcWilcox <- function(dat, val_col, log = FALSE) { ## Separate components taxa <- split(dat, factor(dat$taxon_name)) taxa_names <- names(taxa) - taxa_annotations <- - dplyr::distinct(dplyr::select(data, dplyr::starts_with('taxon'))) + taxa_annotations <- data |> + dplyr::select(tidyselect::starts_with('taxon')) |> + dplyr::distinct() ## Perform Wilcoxon test pvalues <- vector('double', length(taxa)) @@ -589,11 +476,10 @@ calcWilcox <- function(dat, val_col, log = FALSE) { Perform statistical test: ```{r run calcWilcox, warning=FALSE} -wilcox <- list( +wilcoxRes <- list( wilcox_counts = calcWilcox(data, 'counts'), wilcox_relab = calcWilcox(data, 'TSS'), - wilcox_clr = calcWilcox(data, 'CLR', log = TRUE), - wilcox_relab_clr = calcWilcox(data, 'TSS + CLR', log = TRUE) + wilcox_clr = calcWilcox(data, 'CLR', log = TRUE) ) |> bind_rows(.id = 'method') ``` @@ -601,7 +487,7 @@ wilcox <- list( Filter DA taxa ```{r filter taxa} -wilcox_DA <- wilcox |> +wilcox_DA <- wilcoxRes |> dplyr::filter(adjP <= 0.1, abs(logFC) > 0) |> mutate(DA = ifelse(logFC > 0, "OA", "UA")) ``` @@ -639,7 +525,7 @@ incorrect_taxa_wilcox_clr Let's plot their values for each matrix ```{r get abundance values} -transformations <- c('counts', 'TSS', 'CLR', 'TSS + CLR') +transformations <- c('counts', 'TSS', 'CLR') l1 <- vector('list', length(transformations)) names(l1) <- transformations for (i in seq_along(transformations)) { @@ -655,7 +541,10 @@ wilcox_raw <- bind_rows(l1, .id = 'transformation') |> {\(y) pivot_longer( y, cols = 3:ncol(y), values_to = 'value', names_to = 'sample' )}() |> - left_join(data[,c('sample', 'study_condition')], by = 'sample') + left_join( + data[,c('sample', 'study_condition')], by = 'sample', + relationship = "many-to-many" + ) head(wilcox_raw) ``` @@ -666,29 +555,15 @@ Box plot of incorrect values: l <- wilcox_raw |> mutate(taxon_name = sub('genus:', '', taxon_name)) |> {\(y) split(y, y$transformation)}() -l$`TSS + CLR` <- NULL l$counts$value <- log(l$counts$value + 1) -l$TSS$value <- log(l$TSS$value + 1) +# l$TSS$value <- log(l$TSS$value + 1) # this (pseudocount) should not be necessary since a pseudocount was added above +l$TSS$value <- log(l$TSS$value) # this should not be necessary since a pseudocount was added above wilcox_raw <- reduce(l, bind_rows) - - - wilcox_genus_plot <- wilcox_raw |> - # mutate(taxon_name = sub('genus:', '', taxon_name)) |> - # mutate( - # value = case_when( - # transformation %in% c("counts", "TSS") ~ log(value + 1), - # TRUE ~ value - # ) - # ) |> - # mutate(value = log(value + 1)) |> - # filter(transformation != "TSS + CLR") |> mutate(transformation = factor( transformation, levels = c('counts', 'TSS', 'CLR'), labels = c('log(counts + 1)', 'log(TSS + 1)', 'CLR') - # transformation, levels = c('counts', 'TSS', 'CLR', 'TSS + CLR' ), - # labels = c('log(counts + 1)', 'log(TSS + 1)', 'CLR', 'TSS + CLR') )) |> mutate(study_condition = factor( study_condition, levels = c('healthy', 'bacterial_vaginosis'), @@ -713,7 +588,7 @@ wilcox_genus_plot <- wilcox_raw |> wilcox_genus_plot ``` -```{r} +```{r, eval=FALSE, echo=FALSE} ggsave( file = "Figure3.pdf", plot = wilcox_genus_plot, dpi = 300, width = 9, height = 4 @@ -721,118 +596,21 @@ ggsave( ``` - -```{r, message=FALSE} -stats <- data |> - mutate(taxon_name = sub('genus:', '', taxon_name)) |> - filter(taxon_name %in% c('Actinomyces', 'Corynebacterium')) |> - group_by(study_condition, taxon_name) |> - summarise( - mean_counts = mean(counts), - sd_counts = sd(counts), - median_counts = median(counts), - - mean_TSS = mean(TSS), - sd_TSS = sd(TSS), - median_TSS = median(TSS), - - mean_CLR = mean(CLR), - sd_CLR = sd(CLR), - median_CLR = median(CLR), - - mean_TSS_CLR = mean(`TSS + CLR`), - sd_TSS_CLR = sd(`TSS + CLR`), - median_TSS_CLR = median(`TSS + CLR`) - ) |> - ungroup() |> - arrange(taxon_name) |> - modify_if(.p = is.numeric, .f = ~ round(.x, 2)) |> - select(-starts_with("median")) -stats -``` - - -```{r} -types_names <- c("counts$", "TSS$", "[^(TSS)]_CLR$", "TSS_CLR$") -new_stats <- select(stats, taxon_name, study_condition) -for (i in seq_along(types_names)) { - pos <- grep(types_names[i], colnames(stats), value = TRUE) - mean_vals <- stats[,grep("mean", pos, value = TRUE), drop = TRUE] - sd_vals <- stats[,grep("sd", pos, value = TRUE), drop = TRUE] - new_col_name <- sub("mean_", "", pos[1]) - new_col <- paste0(mean_vals, "\u00B1", sd_vals) - new_stats[[new_col_name]] <- new_col -} -new_stats <- new_stats |> - rename( - Taxon = taxon_name, Condition = study_condition, - Counts = counts, `TSS+CLR` = TSS_CLR - ) |> - mutate( - Condition = case_when( - Condition == "healthy" ~ "HV", - Condition == "bacterial_vaginosis" ~ "BV" - ) - ) - -new_stats <- new_stats |> - pivot_longer( - names_to = "Data type", values_to = "Value", cols = Counts:last_col() - ) |> - pivot_wider( - names_from = "Condition", values_from = "Value" - ) - # filter( - # `Data type` != "TSS+CLR" - # ) -DT::datatable( - data = new_stats, - rownames = FALSE, - extensions = "Buttons", - options = list( - dom = 'Bfrtip', - buttons = c('copy', 'csv', 'excel', 'pdf', 'print') - ) -) -``` - - - -```{r plot fold change} -wilcox |> - mutate( - sig = ifelse(adjP <= 0.1, '*', '') - ) |> - mutate(sig2 = paste0(round(logFC, 2), ' ', sig)) |> - mutate(taxon_name = sub('genus:', '', taxon_name)) |> - mutate(taxon_name = as.factor(taxon_name)) |> - filter(taxon_name %in% c('Actinomyces', 'Corynebacterium')) |> - ggplot(aes(taxon_name, logFC)) + - geom_col(aes(fill = method), position = position_dodge(width = 0.9)) + - geom_text( - aes(label = sig2, group = method), - position = position_dodge(width = 0.9), vjust = -0.5 - ) + - labs( - title = 'LogFC of taxa identified as significant (adjP <= 0.1) by CLR', - subtitle = 'logFC is indicated on top of bars. * means significant' - ) -``` - ## Lefse Define a function for running Lefse: ```{r lefse function} -calcLefse <- function(dat, assay) { - res <- lefser2( - dat, kruskal.threshold = 0.05, wilcox.threshold = 0.05, - lda.threshold = 0, groupCol = 'study_condition', assay = assay +calcLefse <- function(se, assay) { + res <- lefser::lefser( + se, kruskal.threshold = 0.05, wilcox.threshold = 0.05, + lda.threshold = 0, classCol = 'study_condition', assay = assay ) - adj_pvalues <- p.adjust(res$kw_pvalues) + return(res) + # adj_pvalues <- p.adjust(res$kw_pvalues) - dplyr::mutate(res, rawP = kw_pvalues, adjP = adj_pvalues) + # dplyr::mutate(res, rawP = kw_pvalues, adjP = adj_pvalues) # res <- lefser2( # dat, kruskal.threshold = 0.05, wilcox.threshold = 0.05, @@ -858,23 +636,23 @@ taxa_annotations <- lefse <- list( lefse_counts = calcLefse(tse_genus, 'counts'), lefse_relab = calcLefse(tse_genus, 'TSS'), - lefse_clr = calcLefse(tse_genus, 'CLR'), - lefse_relab_clr = calcLefse(tse_genus, 'TSS + CLR') + lefse_clr = calcLefse(tse_genus, 'CLR') ) |> bind_rows(.id = 'method') |> mutate( DA = ifelse(scores > 0, 'OA', 'UA') ) |> - rename(taxon_name = 'Names') |> + rename(taxon_name = 'features') |> left_join(taxa_annotations, by = 'taxon_name') head(lefse) ``` ```{r} -lefse_DA <- lefse |> - dplyr::filter(adjP <= 0.1, abs(scores) > 0) |> - mutate(DA = ifelse(scores > 0, "OA", "UA")) +lefse_DA <- lefse +# lefse_DA <- lefse |> +# dplyr::filter(adjP <= 0.1, abs(scores) > 0) |> +# mutate(DA = ifelse(scores > 0, "OA", "UA")) ``` Plot lefse results: @@ -907,178 +685,7 @@ incorrect_taxa_lefse_clr <- lefse_DA |> incorrect_taxa_lefse_clr ## the same as in wilcox. ``` -## ZINQ - -```{r} -calcZINQ <- function(dat, val_col, y_Cord = 'D', log = FALSE) { - taxa <- split(dat, dat$taxon_name) - taxa_names <- names(taxa) - - taxa_annotations <- - dplyr::distinct(dplyr::select(dat, dplyr::starts_with('taxon'))) - - pvalues <- vector('double', length(taxa)) - names(pvalues) <- taxa_names - form <- paste0(val_col, ' ~ study_condition') - for (i in seq_along(pvalues)) { - df <- taxa[[i]] - - res <- tryCatch( - error = function(e) NULL, { - ZINQ::ZINQ_tests( - formula.logistic = as.formula(form), - formula.quantile = as.formula(form), - C = 'study_condition', y_CorD = y_Cord, data = df - ) - } - ) - - if (is.null(res)) { - pvalues[i] <- NA - } else { - pvalues[i] <- ZINQ::ZINQ_combination(res, method = 'Cauchy') - } - - } - - adj_pvalues <- p.adjust(pvalues, method = 'fdr') - - log_fold_change <- vector('double', length(taxa)) - for (i in seq_along(log_fold_change)) { - df <- taxa[[i]] - healthy <- df |> - dplyr::filter(study_condition == 'healthy') |> - {\(y) y[[val_col]]}() - bv <- df |> - dplyr::filter(study_condition == 'bacterial_vaginosis') |> - {\(y) y[[val_col]]}() - - if (log) { # If log, revert with exp - healthy <- mean(exp(healthy)) - bv <- mean(exp(bv)) - } else{ - healthy <- mean(healthy) - bv <- mean(bv) - } - - if (bv >= healthy) { # control is healthy, condition of interest is bv - log_fold_change[i] <- log2(bv / healthy) - } else if (bv < healthy) { - log_fold_change[i] <- -log2(healthy / bv) - } - } - - ## Combine results and annotations - output <- data.frame( - taxon_name = taxa_names, - rawP = pvalues, - adjP = adj_pvalues, - logFC = log_fold_change - ) - - return(output) - # dplyr::left_join(output, taxa_annotations, by = 'taxon_name') - -} - -``` - -Run ZINQ - -```{r, warning=FALSE} -zinq <- list( - zinq_counts = calcZINQ(data, 'counts', y_Cord = 'D'), - zinq_relab = calcZINQ(data, 'TSS', y_Cord = 'C'), - zinq_clr = calcZINQ(data, 'CLR', y_Cord = 'C'), - zinq_relab_clr = calcZINQ(data, 'TSS + CLR', y_Cord = 'C') -) |> - bind_rows(.id = 'method') |> - mutate( - DA = ifelse(logFC > 0, 'OA', 'UA') - ) |> - left_join(taxa_annotations, by = 'taxon_name') - -zinq_DA <- zinq |> - dplyr::filter(adjP <= 0.1, abs(logFC) > 0) |> - mutate(DA = ifelse(logFC > 0, "OA", "UA")) -``` - - -Plot ZINQ results - -```{r} -zinq_plot <- zinq_DA |> - dplyr::filter(taxon_annotation != 'Unannotated') |> - count(method, taxon_annotation, DA) |> - mutate(n = ifelse(DA == 'UA', -n, n)) |> - mutate(method = sub('lefse_', '', method)) |> - ggplot(aes(method, n)) + - geom_col(aes(fill = taxon_annotation), position = 'dodge') + - geom_hline(yintercept = 0) + - labs( - title = 'ZINQ test', - y = 'Number of DA taxa', x = 'Transformation method' - ) - # scale_y_continuous(limits = c(-3, 13), breaks = seq(-3, 13, 2)) -zinq_plot -``` - -```{r} -incorrect_taxa_lefse_clr <- zinq_DA |> - dplyr::filter( - method %in% c('zinq_clr', 'zinq_relab_clr'), DA == 'UA', - taxon_annotation == 'bv-associated' - ) |> - pull(taxon_name) |> - unique() -incorrect_taxa_lefse_clr ## the same as in wilcox. -``` - -# ANCOM-BC, MetagenomeSeq, and DESEQ2 - -ANCOM-BC -```{r ancombc} -# ancombc <- as.data.frame(DA_output$ancombc.none$statInfo) -# ancombc$taxon_name <- rownames(ancombc) -# ancombc <- left_join(ancombc, taxa_annotations, by = "taxon_name") |> -# relocate(taxon_name, taxon_annotation) -# ancombc |> -# filter(q_val <= 0.1, lfc < 0, taxon_annotation == 'bv-associated') |> -# pull(taxon_name) - -``` - -MetagenomeSeq - -```{r, eval=FALSE, echo=FALSE} -# metagenomeseq <- -# as.data.frame(DA_output$metagenomeSeq.CSS.fitFeatureModel$statInfo) -# metagenomeseq$taxon_name <- rownames(metagenomeseq) -# metagenomeseq <- left_join(metagenomeseq, taxa_annotations, by = "taxon_name") |> -# relocate(taxon_name, taxon_annotation) -# metagenomeseq |> -# filter( -# adjPvalues <= 0.1, logFC < 0, -# taxon_annotation == 'bv-associated' -# ) |> -# pull(taxon_name) -``` - -DESEQ2 - -```{r} -deseq <- as.data.frame(DA_output$DESeq2.poscounts$statInfo) -deseq$taxon_name <- rownames(deseq) -deseq <- left_join(deseq, taxa_annotations, by = "taxon_name") |> - relocate(taxon_name, taxon_annotation) -deseq |> - filter( - padj <= 0.1, log2FoldChange < 0, - taxon_annotation == 'bv-associated' - ) |> - pull(taxon_name) -``` # Plots of BV-associated genera @@ -1226,112 +833,6 @@ sample_sizes <- filter(data, taxon_name == 'genus:Lactobacillus') |> data_with_lact <- left_join(data, sample_sizes, by = 'sample') ``` - -Relative abundance vs CLR all taxa - -```{r} -data_with_lact |> - ggplot(aes(log(TSS + 1), CLR)) + - geom_point( - aes(color = study_condition, size = lact_tss), - alpha = 0.3, position = 'jitter' - ) + - labs( - title = 'Relative abundace vs CLR per genus', - x = 'log(TSS + 1)' - ) + - theme_bw() -``` - -Relative abundance vs CLR BV-associated - -```{r} -data_with_lact |> - filter(taxon_annotation == 'bv-associated') |> - ggplot(aes(log(TSS + 1), CLR)) + - geom_point( - aes(color = study_condition, size = lact_tss), - alpha = 0.3, position = 'jitter' - ) + - labs( - title = 'Relative abundace vs CLR', - subtitle = 'BV-associated genera only', - x = 'log(TSS + 1)' - ) + - scale_size(name = 'Lactobacillus Rel. Ab.') + - theme_bw() -``` - -```{r, eval=FALSE, include=FALSE, echo=FALSE} -# Plotting CLR vs Relab of Lactobacillus, Prevotella, Actinomyces, and -# Corynebacterium - -plot_1 <- data_with_lact |> - filter(taxon_name == 'genus:Actinomyces') |> - ggplot(aes(log(TSS + 1), CLR)) + - geom_point( - aes(color = study_condition, size = lact_tss), - alpha = 0.3, position = 'jitter' - ) + - labs( - # title = 'Relative abundace vs CLR', - title = 'Actinomyces (BV-associated)', - x = 'log(TSS + 1)' - ) + - scale_size(name = 'Lactobacillus Rel. Ab.') + - theme_bw() - -plot_2 <- data_with_lact |> - filter(taxon_name == 'genus:Corynebacterium') |> - ggplot(aes(log(TSS + 1), CLR)) + - geom_point( - aes(color = study_condition, size = lact_tss), - alpha = 0.3, position = 'jitter' - ) + - labs( - # title = 'Relative abundace vs CLR', - title = 'Corynebacterium (BV-associated)', - x = 'Relative abundance' - ) + - scale_size(name = 'Lactobacillus Rel. Ab.') + - theme_bw() - -plot_3 <- data_with_lact |> - filter(taxon_name == 'genus:Prevotella') |> - ggplot(aes(log(TSS + 1), CLR)) + - geom_point( - aes(color = study_condition, size = lact_tss), - alpha = 0.3, position = 'jitter' - ) + - labs( - # title = 'Relative abundace vs CLR', - title = 'Prevotella (BV-associated)', - x = 'Relative abundance' - ) + - scale_size(name = 'Lactobacillus Rel. Ab.') + - theme_bw() - -plot_4 <- data_with_lact |> - filter(taxon_name == 'genus:Lactobacillus') |> - ggplot(aes(log(TSS + 1), CLR)) + - geom_point( - aes(color = study_condition, size = lact_tss), - alpha = 0.3, position = 'jitter' - ) + - labs( - # title = 'Relative abundace vs CLR', - title = 'Lactobacillus (HV)', - x = 'Relative abundance' - ) + - scale_size(name = 'Lactobacillus Rel. Ab.') + - theme_bw() - -plots <- ggpubr::ggarrange( - plot_4, plot_3, plot_1, plot_2, align = 'hv', - common.legend = TRUE, legend = 'bottom' -) -``` - Plotting log(CLR) vs log(Relab) of Lactobacillus, Prevotella, Actinomyces, and Corynebacterium. @@ -1344,7 +845,7 @@ plot_1b <- data_with_lact |> labels = c('BV', 'HV') ) ) |> - ggplot(aes(log(TSS + 1), CLR)) + + ggplot(aes(log(TSS), CLR)) + geom_point( aes(color = study_condition, size = lact_tss), alpha = 0.3 @@ -1367,7 +868,7 @@ plot_2b <- data_with_lact |> labels = c('BV', 'HV') ) ) |> - ggplot(aes(log(TSS + 1), CLR)) + + ggplot(aes(log(TSS), CLR)) + geom_point( aes(color = study_condition, size = lact_tss), alpha = 0.3 @@ -1390,7 +891,7 @@ plot_3b <- data_with_lact |> labels = c('BV', 'HV') ) ) |> - ggplot(aes(log(TSS + 1), CLR)) + + ggplot(aes(log(TSS), CLR)) + geom_point( aes(color = study_condition, size = lact_tss), alpha = 0.3 @@ -1413,7 +914,7 @@ plot_4b <- data_with_lact |> labels = c('BV', 'HV') ) ) |> - ggplot(aes(log(TSS + 1), CLR)) + + ggplot(aes(log(TSS), CLR)) + geom_point( aes(color = study_condition, size = lact_tss), alpha = 0.3