From de319d0497b215270f4ac0ecb8536d6e636acb26 Mon Sep 17 00:00:00 2001 From: sdgamboa Date: Mon, 23 Sep 2024 12:19:39 -0400 Subject: [PATCH] Lefse DA are defined by both effect size and p-value (raw) thresholds --- R/new_DA_methods.R | 7 +- R/set_DA_methods_list.R | 11 ++- vignettes/articles/Figure1.pdf | Bin 17034 -> 16844 bytes .../articles/HMP_2012_16S_gingival_V35.Rmd | 91 ++++++++++++++---- 4 files changed, 85 insertions(+), 24 deletions(-) diff --git a/R/new_DA_methods.R b/R/new_DA_methods.R index 1771f47..f7ea010 100644 --- a/R/new_DA_methods.R +++ b/R/new_DA_methods.R @@ -257,9 +257,12 @@ DA_lefse <- function( statInfo <- statInfo |> dplyr::mutate(abs_score = abs(.data$scores)) |> dplyr::arrange(abs_score) + ## These raw p-values and adjusted p-values are artificial. + ## I'm only leaving them here for now to help in the ordering of the DA features statInfo$rawP <- seq(0.04, 0, length.out = nrow(statInfo)) - statInfo$adj_pval <- stats::p.adjust(statInfo$rawP, method = "fdr") - rownames(statInfo) <- statInfo[["Names"]] + statInfo$adjP <- seq(0.09, 0, length.out = nrow(statInfo)) + # statInfo$adj_pval <- stats::p.adjust(statInfo$rawP, method = "fdr") + rownames(statInfo) <- statInfo[["features"]] ## check names colnames(statInfo) <- c("Taxa", "LDA_scores", "abs_score", "rawP", "adjP") pValMat <- statInfo[, c("rawP", "adjP")] diff --git a/R/set_DA_methods_list.R b/R/set_DA_methods_list.R index 75b7096..d4c6683 100644 --- a/R/set_DA_methods_list.R +++ b/R/set_DA_methods_list.R @@ -190,13 +190,16 @@ set_DA_methods_list <- function(conditions_col, conditions) { # wilcox.threshold = 1, # lda.threshold = 0 # ), + + ## I really need p-values to help define DA abundant features + ## So I'm using a p-value threshold only for the lefse method lefse.12 = list( method = 'DA_lefse', conditions = conditions, norm = 'CLR', groupCol = conditions_col, - kruskal.threshold = 1, - wilcox.threshold = 1, + kruskal.threshold = 0.05, + wilcox.threshold = 0.05, lda.threshold = 0 ), lefse.13 = list( @@ -204,8 +207,8 @@ set_DA_methods_list <- function(conditions_col, conditions) { conditions = conditions, norm = 'TSS', groupCol = conditions_col, - kruskal.threshold = 1, - wilcox.threshold = 1, + kruskal.threshold = 0.05, + wilcox.threshold = 0.05, lda.threshold = 0 ) ) diff --git a/vignettes/articles/Figure1.pdf b/vignettes/articles/Figure1.pdf index 404413c904dce32e36109973996d202d34819614..4d943ffbb526b470c248e50d4643a3e9a72b0266 100644 GIT binary patch delta 13147 zcmZX41yCDp*Ea4F+#$HTOOfDKoZ{|IaScvzhvIGpic2Z(#idAbD_W$*p}?O$&-?!K z&Age(?C#m?=zZnv-e(dU2|F1Mn=FjP&BM*l!66hdVP;=ZhA4vrulNaEH%Uj_TDXZDr0EFOSwnBw~JQ@4RND=~=KO zej~Wz^qVD*y1sP$%6;AZb2eOv{p!5F^jUSxnQXlRscg5&iKzYINAOjzAQe7wma~C? z@W`2_V5>`3i+ylz)Aq}cy*u#hXk}$BjCon{{G|7J!`?BQe~V2-Iz%z$j-R z3O5!4cRajYl~*66bSYFfX9}L#{R46s;7JlNC=Q4qMhPQE3Ui|;40GdWLy&^R03qfy z_xQ!B4)9#g@i1HEfoMolEm8J}V67~}8cq7cPqlpBZt;*zF0{0`{ z-F>lZevtB%-+Pfa*rJxjDt~Gt*}6R$&XEtZpY7GG4I^Cy!)Gs+tT~qBo{|G<@FzDjyB!K9 z)*t`;wCp&o?{p~5_-ID9^I1F6z?%&t!I}|;`wW3wMN|sn%m^`Gi&-JudpmT?z4ulL z|MaaA`BFY`s0(k1>jHuLHXj%j#Ap-~%4mcMhF1ev?4Vlc`PPzti@;mtb`sLXgzx`o z3Do0f0%=D}B|8wK=u#j_eKZH^c^Nq=KdMWqm{_Yka*%bdwEG3VVSE?Z4(b#Tbs#I4 zwOr#=fh$qNFU~-=@rF`q`L3&co{{WEw1^g+R%;50&z`BR> z&oy_ZFgKEVWM`MzDjTOowcQ{ylrluX##8v-_f#%}_R=EWs@e-IYH!X^AUKhS8m0F@ zSK+>)?9IKlA1&80#&~rse2SmLC%sg)-@R<`WSX}d3^t-0fw$JG(NCLP{Wh)8M^{-F z>Y{#roAI?$lE^AdQGVx39pTRZ`&!c@y>eN@ls2pCfGNu^ zFCkvtBOPaCZ48B9p2U1TVKO79bTZ?Y;6Hj8P;5g^D9>;-$(WWe&^Ol9iRpKGkCN?G zdVv%CH?TO-q>2;1(X-}f{OJ5re?T@VP zI`_qk(w@Ya2JcpuNsW)ib}$ZhWtk6M%GsJVZ6cVgrHIQXQ8dQGxRkO1AJa+PjhL0o z^_i7pX(Ior?y#)VN!YI${1Q^6$q}Z3fmm)^m3$R+H$%ONF0z!H$ z@D|MT6Cca5FwYPy8;mMRv`kI2L>xy@Y<^S=7gz`cMQ4-#{E}zIk>M3N?b=vjw9J3f z4r3_XH@j7An>YTZa+lukd7f-_Qz#-NQE2_Jy16eG-PRcOoP7AJaoGE)8S#>-4&TJ;oh`AQ( zJ1G3!UnEftVI%mE&rDY+84#P zJP&EC_tS&&;LpE9-WuCs2Qtg=-s%OuSw^bI=WOFtmIqrtp@E!Cr7+gs0W6--ELM7G zMu%?+go7t<|l*S-LZYTzF#tJciLVE|ry&JqGh*;@?Rw=NIRStvdG}XMl z=@T)^5CkUrR86(VdSh9hG8)cpVD%rPq%Jv%dnhPrvWzZ5T^hjv@{o%|OPs{@Z z_A!EX!o-ZN>3l-&T9JbUNguh!l!J~~VO{dH{LR|DH^oAqc+T}*p|Ilwomy*Q+B`*f z@?X9M+T`y(Tk5;37+&I7kZ{bB%bR$50PPBmj)q6gy>8~IiPB}NYXIOuEJqh~Ar-S? zBiAvL<8?h`-h?$FP#Aoegr+?+;+kY;40Kq0(;PjJW?dqstuXjIQuRvR*Chmh2cg z_6AjT97g5(ci28u55WtnfCO& zu4z{`8yUqn-ELoDO6eKhoqR9|3JfSRIeh3Kj#Ti%PEJCgPRkCEZEfH&6{y3VU)*^p z+C@9cGpVOC>3R9HY}&KpnQ;00yz-7J*#8yWb9AfBNZE5gGnFf$q&oB7n|TawxB%Sd zZv&CN)|~0htKV`@4f>-&e>v+KZDlfYT<+yLi9zvN8I!+1lN< zCRlIk7*uB=s;4a#Al%k;6B+yX^c3B0tmIYw=l61TEm-@EW%ZF;G_mY)=f^wn>FwGG_Wk1q29J>o8L3Dg z6d}j1U4%mQU*2beP|S3O1x7Bkzn-F81sGR9t2S45KJrDlK zND5M1M$_$~*Lo0#ynY5}Bq1)(;25Fo-`zh8c2AqHvFEoZ6-op?`u-YIE}lu;Ol)YC zbdS1|c=3FFzMX2{GBJ_+1Qe1)Mvjg1_3d5T_|Eg9!46-$zbRRb{Sw*f)NAkd_BkJa z`T40Z_{DQ7ex?DFya63fS;)C@4gD-`*+o*iA$c#QLesWsW1u&$$QwXHQtw8GO%t5@NqVYwHxhP?rE1y7O5j1^u6GntlvA!w?L#l zi@J@{Lrje7Y1;%hXpC{}p01Hi`?RmEiu!md-JY|b1@&N)_i6=mqFO3r5sJ?3d_qS3 z%qbAW{QA2_^<9{R#7l?Z>nGnzzqzb~=LSIN{kRKCgW2=>)j8)~uM_2#^YvNbd2sg^ z12S>_F9=B*{Kcn^sp_1)!LRq1_jta10uq}+ouW$twnksTZk@#-g`UAy#i%KyNBY<0 z9klfZGThFo%t|z)f#3Z_6D@tKyGog^sMV7=-lFwZ4`!a>oIrKxNoGiOh)c8D`zQmfT^(Wxad(TDo3-mZnTl?ie1X(D5Vs#-_7A>dP zqv6X-_T#z3Lw8+#*uTEPOJ2kDrSOTss=-t+e{#zGsO>&VSG(-%^9lcQ(NQFYq1*5e zf1)Rist6wun<`>oi>xQ>7yfYID$M~z?4t(^jwSd+rNus-j^Gp>9pK z(U0bsS95Mr(r$f0i+$@t6`H%BQd|Ika09lvGJw6ef&r+7z%4~-6MwRXGZjcd*e2}H zTyc&ta>Iowx+5O@BEIBt6p{yQ!uJ2n`HK-ooCIbFy{I0!BSRFpAWqFM0@pzX$JggV zOShR42}8{AM%k3l##dIq;fL06A1`_zBYcQ2oQ&I<+cFId8{D})g4-4{hQtcWXEU%E z1x%-rf|-m}L!@`^@Zd4)&mkW8%darN0+^r@D-tydfEn=XB?tz)BVxyci~)TL2*Y<> z-&c#;NB2|9GezI+&;xEe4OUVN&Mcs79C!{Rmc(L7s=Td1X|C}dRweZtLWmD48>LXw zMiF(Fp&eB*C~X=LXa)QbNsXgxCI${bY4XiiMc6?l+&Lnb{9t0nnhWfK%Wcp=rl#)b zCk#jnHKeWs-j7K#IYa=7Hz#?0qSJ~63n#6@F88MaM6GKm8*5(_5@^nQkU9>k;L2)3 zIf$fP(b;8zF{X5YNl@}EFg3b-B(XwKlNnUrVCD`K^&u-) ziAzVNf-X}+%?P8mr-_tnDW;JcWc0&uA|eeU8hh86y8CuxfI%UpN&Fsygh@-L^;HwB zDB?VT$$?}LY!DMWEUqMb%rzYf`Zt6)$cEf;1^iV*1s-**EpIfyI0|6h#}j8hq^6n4 zSFEH@$H=VQU&dlhHUh5?rp~dygaDjHyb)&Yj!=N7V>rxlCq4@L^CPZcy_sW}V= zvtdiKOhhKyCJ8O5`=d@qQvnMNV4O!$pQk1&-I{ko>HEw1us$0anbM2%f(-XXQQL0&cwSJ&va%%-~2I09?k#=_{J z%;_*T@-rs{{iExnfv0GNK7NI$_^M_4|4PBx9 zBDtwqVk~mm#|h`x$7L&3zVNC|x;(3Vf$4ESs}Ce{lfv%#qKTS$VG>_>?>%|AO(tLR zJAo*YxQEC<6n%?}(V9JIyt2{Zc}!4G8O%|Mptk;kl!+m$fIYIdJf0M9je`}DT0-E)9_6U!^tCX!|N2GjcjL#0wbCIVY zuH87di5*+%yLDJ3wds)g-$#S=k6`oJgp>MvfmPZ*2eiq7J}*V`eBurc5YeZ&?YL_6 z+9)xDK=kSnL^+WK?C?)w4G1%QnXon45E`VXHZ_3{H!2bso$rkiOjV7|_)}*GPILvC zwlU6WE2xmz;Xs!vmL!^YwzL9e!pr`e%I6Feb!Jj&!e>5kZ`kJ_DCbvo4Jv>v|e(!gkxPj_4Ry%p4@MVp7_pP}k2hIinXtJ=bC%~^?@obQS&cWrb|mNZ&we3M zOQ_ADjNhD+8>Ovzb5*vJ+H0{+SCvLCSmFfnkF&D^E%>!~m{c@4fv)g8>y%A8Wz-2s zG$Sf}ltr$F{n)o$O}pT`CgWrdcP7&?-ZJRbig2}{i^R%k&_rHBE%aSS@;p+4p(Q*E=C?<8M?P0y zflPFA^nIW4(`n8au=7tPl@oL`(c=|m&t(EFe4)?3PwmQd{S;>>%$MPPfQtAMRvte% zY>8@71dO~04?ULycPf>iq?a_IPN@+ZcxjnjZfXBlRrtq0s$(?@ap+fcp`LIb(VfXV z8=}nj(3k^Fv%s54o*9YB@0i#OXYV8KsD$gIS?s?nbX&i@+(_pL)Yy!H`W@5f)9|@f zj3gaI%nr5kh&tt2zO&9m&FWb%bop0w03D z`h+V_JyaP{5k0H91f?%+(cv4kViR_4b|Uu3Unl<7D>!x5n$GVRUjz4tH+P@Aau>2W zryU;-I+H1P%KN9q`{}vW6FhZ((nF*Y@q6lb%Q;iu8Yulni?l{3Eu{0bQZxatMrJjo z_)p~VNz7++j*I@9@-hz^ITw3EU{9LO@J8AFurQ3EdEc|CVr*?MVz2xkx0($HsYLZ- z%T3J+Zb5>Tyry}}zo)!#-sEXBOr*uXZ{5M-LLo{1+yr1@NehTh7ns+NgEAy_+{!O3 zv*ch{Cn=KxiiD;L{wLR}i5JcHGo*dP`&O7`!v^q098vs$iB-}*&sL;)i@MPLr6riE zG;_={r@GMe%&6= z_WTW2tuPh-B|fVoOg>{8e~76aAOWNOL3?pUM&D6#UVlw%vr%*0ybXy92$rr3g(EC? ze`y@9UV20WF5YK6cI-yddBEP_Yqy#`PGv!|USLjk2d(-=1*pD1gb1at~4VJK#TvC6cVr``JpycNFUyPh7^*<4;YWV$lwRG zaQsG)w$!HG<9GYmn6@O#Cju??VV4SGU7T1kyYSbzEDyWyZ%=714m=Oyt_=QGT114= z5oV81o7S={!mA=K*fwGMZw~dJE_083f;HW|6*uF-&C26)4oyfx|B}8D$)^2foD^aa zeUgUoUYh?sp$9p?H|1kKzkeob>}tW)N7MK-^}WN$ZTPQohoH$N1&6;CIKfNw`b!5$ z(-L>hIc^lcqz0k|rj`C_)!&<#Od|r&J9CVoVyJn`X^c|sUcPmu7q?9^Qaxa%i=n@z zYY7$ljy$agfB#I~3=6bWgr-;X;xj_p*96U>J7J`LHlF`qs{GT^SRMOPw%+sF8}E^! zq`3lrJ|GS3lU*(eN?NP+gbF+$WquWGfC{09*;-f0agzrmCF3U~hqfMA6m_y5*d3R9 zB&C)f*fCLXAiVVo)-IXbl3<8v%L!SP`59S==N_5cU#?2hOPpI#^?BFWS#@bC@(XCD zxrQw0yjq$pKAY%1ov>e;?+;^Yvy)YSCEfPCLAr0)4VWQS@M=R+P{DEpB06H0Wk415 z*E~#Gk%lVBmL#2)NXIq};QfA(LFYC|7KXpN!A(eRU202`5&a)Yg+nED`Mq($bq=0N zM9yk{k)+Ma1!mh(2VO@zDNXtQoBJDi(|;Qqln>jBN@V_oWzlw4sgxy+DLJ*UZ{Iy{ z@B>)a!m8T#2-W|9Mwa}*QwvLcHAs5--2k4)G>Sju*a3`b z1`W=t2g4C&VnK*`{OE@SK`n!1q$NEW?DG zejviL8@m40EMnw?%9q*+&K(qqPnM*W`mpeyclneFr`7YHr^fg3Jidyi-i zft1_Xuu@sXHc^Dz83)~k2?JrS?U|gzP3{TD83!jH0bc4F-wvOKysp;QfLp;V5rNd) z=?4m5Y6?8g;zSkEYQDM$i6G!3dZ~s3s>W!;C{t+#{%Qs|OJpFJCEk^U7d>_wxlSX? zhuDCccF0zwkoE?4LOdUhL>P>s=#a_f9r!KYr>xut07LFP#gb1xpQH1ZihuA6tLhY& zuiWFOxYi}7DVnn~2Q@zuty|)Yj}u2^>_(x8r9VZDlXn2jISliH>lsn(nsrt?RkB8htNCn@Hz}S>?C-K>Z zt6KuN$T9Md9ZvMdv;k05={vFdS^Gb=yu?8SQ316b%f);LJBju}e3B~mf!uFXR#tq* z`S$0`nTgg!<2(B^V7v=x326BaQm`r{x8KTt1FMS9wfTf`8zCq7`2wzMo%@d`+*>1& zD2#LoSBw1C0)ieVEQ{VCDM?qc$79GJqLCy2WI>=JQJA&F;R-HUMhydWlqc9lT>`H0 zdFws_b%PKqK-b6CaRU-3pzVxm}^qeNPG+>arRu!|JhH{nJ^3+OYF6(5@s zf~EXp_!oBA#G02D(DG3)oug(;VT;)?q`nV^!2)sKHL@`QXtXLy_DeoJK5&1qoyXPa)Oz_Sd{Rxa>d3X>M zESTI(v$E4_)rDzTmL&7eO61N)COz*;ak~U^wCm{$Y*sFkJm@5Eutbs2L;>5Hqxe-Q zp!Ubm#rNZs%K( zJ!%>mnJPx+a?d%-Dj*J7WaCaY=^q7T;(QJ4WR!o(t^t`T3H0%mrj)1cpaBq0p<^;H zF5|P>&2W~bCve|iPoROBO*096(bXF!(4@wiKmbd9Y7~846bm+1;^GN3xewNe>ztHR zO=0$qiD57N6KEBfaOcMy44nn`aS&gfBT0ZW4t{_LK!(a#S7F0ipW2J+e+-p9Tj^7u z31$+6C0j=VN1A9tr!0w7WY!3gNX#?A1S)Gwu%d%(tRLDS5jafnQZ2IDf4-Be|ywg6EXM8c9fL3e|!H&adgi9&*I%dTy)x};AVBY zt|i7a3dWmi619+oM?R55;En#hCD*(qrn9vcRab!Bx5}ACwuv&i&|_!r3GtmTS0v46 z(Xr7p85f)9aX-E!m%ZyDQY$sNujr#%z>eWi;tCu*_W6;E6`36j;ep4x8|<0N-Tp1o zJ0p4CMSg!fT(ZySDzOF8h0seyS8T^mFn&gCC4(CfRDv}0DnGUdL&S}bfK{WvK26zzKU+=Nr?Xzc1Z^| z<1Hf#-O2IL&&v@Ww~jlUp5W2E&NVA}N3P@oa~&9Qp5&?u=i3hd{zY)%8u+j3Wrt08 zBa4G6@Y+t237sm~#S=(?8P*N;DB2|~)DXFH0CzMOaw&`u(MJ=G(Q4f*aiA^Kuj&Yr zFP^D@2`&RngOc-p#8fLNn7s}AfLMx{D4BO9!aO@ojWOUht-gDrDIk(kQ zKNzp0gZrC0Sbv}H1@=uEuCqdecz=EC6B}QGKF(aCD<4j0Rd->m{9b+wKt)Yl`~3P_ zGuo|ErJjyiqUBI^NqCsoGS&sgP;;`XYn8f}x;_@<17@d#0(xh!{@&T8DBJdk6($5w zD${_EGd-Pmc4y|lBZ`X1C%L@rx$189<^uwTpn{)Of-fF3J<{2^I5C<&Z~-6q=f}clO)zfZ8Mj2pPylO=T5p~}I?`%|LYmrDlHE<5f8J0m*rq#NfA0vEr;bxWjiT!oko$(SVdJ4v_;pdRq*il{vxu7vj%0bS zN`$2WC{Z!P$he>%3PSVb2?~2w z>0|uGC5|=X`B(<1OkS=qAudFlLiPjk=ru_>m^4WD`aO9hey4>zwZr*TQFw{1cWs@uQ>ZS$?Bk`o+ho0}YKw3tmRjO2fZ^Y;)Lmn60VBRw+$Lxoe z&iEYhnMR4S#K=!^3*60@zFDNVa|xr&);G7o>%bu))d7^ZD6ey^4ph zC*d}#N1*{$#&Yq<(a>*|43A@)He?~cA8t8qp|quUq(CbGmQt_|T-^6T7ZLhsb9V^N2nmK4_v9sj25s>dL%>0}lDy z<5a! z`F+vPKEiMP31Fq0!)gRYA#XSI%%BHXX-+zZN%fWte ziUCAql$0jumC{)g=-k(me`n>2VS0%?Bg(OLn%8#l(5b^e2QAReQiP#z1nBV~Rug(C zej)~u#}-9Jq!hUEDg)L~7{V$`H8fn4bz0BH7;CRPJ`hB8CV?$9`qf82AT%}r29#{{ z?RF~0C91DJ_q_{iS~2JQ@MA1-a?N~oJM+uDr!zk0vWidYZ%XeQa(!Nbszso~K8cnz zEs5*CIeN@DjlPg<%g1=lgr@f*BHD$*w-wybVe$#nu9fh=?drqg_o8GK4)jRhZXdlx zjiH%1Xx9U$M8-1Jd~DZYDNDsxrFMrfXvm{R9difCRm&B!pP;CZfzA;k>AYBLL_mu{ zDuo0ZKTBq~4m>zus7N5(*fJFSZX=Y7wh0&6HrA~&nEE&jef83y{?4K_L#3$nE~ahqi<_aj_5m_^S&3WU#? zjTwl+=j21Ut6uzU7O{ceDnev3#GMb%*+$ICm)YM8uo!mhA8Ag84oPhv=fCVir=(!t zOoDG2-53SdG7=a$^(tCo3wlU1C9&xEZfmn_ADhf@L>3Mj|01n2XH^1=@nRTi<}r1i z!gnK~B0}x;pZ>*o<(e2g3%CaKzU7S=BtJ5GojLJBIe!r4&}30U!^!L~?GTjcd~b!1 z2OYx}q1O5q7h_PA%dY#2C8L9d-$Wz^mKSbER=~{P{GLzTOg{Qqc~Ue)jSu1^4+H)k zYhq$?1~1JiT@!axyHB7O*o*q(olhmqDMJO9jonsq^|RjB3_Z-s%VL3Ane&0XnY&Y} z`aGnx!=5L0^#1XiL#8=^E+|H>=ZBJ<SDy013ss4MhfqXx{l5~jWU(}J@ zC*on==c|K_6HC#b#`v!fc)=Rk>soHNF2XnYzdekM(3*&-B}~gD!yiBE%J*sqXJ*A| z-2Pm=T4W~#$F~W1t-WeL=u>jYeeat&tif~o_2h8sxm~(>__>4Fy1k;HUYCAWs!^k6 zkcw4eVsyBQ@zH5_vGCJRG(?cp<0`MrOYD9*Rah-5t)SeJ+k+ro{-{D+il9k>;C7um zgKMXlt7-Uf%;BG;7hR@cg`$n6U9;ML`v=1uG07EKa8KP=X;a;{A8icS`^pct0)7}B z=3ZK9qY=44L4Lnl&y0R{)dF&g_Fm||-bfViQmu+*+QKd<;CJFJDwD7=zkf;M(6Hrd z|C{_5T2*e27}1~9S&8PaytRZI&Z=0CpVn9Q>G$&Op|+}td7*~iMVr3ftE#@J7)zOY zMR9_^4*71}6&>%5H)z+sjEnEr5@3s&R3{$WjB>o!{H#l+wFJ3KnxeRm*b;m*&)ptK zCt|bz=5l!e?8|H#4E}R^Z^jEAx_LQYI5a(%5MQ$`HTR(3jy^N(dAOU(%9Z%@ROft` zyGE~dogW|9&T`(aSsT2Nzr^hsyP2Cxu1mLw*8hJ zkF1MFT-1C}tk06CVFSl8OVjrxz;XuJx5(bW+0GmNE$a}&H!M8ME z)`FZwlm+V4SXGV(2#UWgi$sx~%k-g}TDr2Ae=%uqnr1fHxUJut@6?rh<8Qqou7LC! zsIf+prr-&OC~}M40W5QYk?@dT$Lm6FVj#gg2p0uS^>sPZw7ky;^XdxeKOpWR=vXuB z14&l^A=Fde_)iV-iPLN>VFF^G4zTq0ksA>=O%ssZLD5!b>5;t#>oXgv{W$<-;v}w- zx@jLamhy6#Gcpp_iJ;X^hLF_{Wh#G-PZSOPIb1%e2!sDh6HD3c}_WQW8xT9f2#j0nlG>_?t8M0ivElk9}0j zeJEW_SN?G<)kl-XL;8WO-d1E^9oR2RK|>QKy5BNp&}+R&9sQxBK{{$|7n#n#!(D;N zeE~CN7?ZVnWj_@==HF1KYvthjA7;=t9l^#?I*=18ee;yF(il`a0U>8FZGn(; zwXJ@^z) zd$*AP6M&uc> zsp%hSlsCUt2eV-?nm9)r?NkUKIhCg2UhWLYP1n@FTtB%6W4_+p-C?^@eYL$gb^WCa z&iAz+QqXt947%BqQhWH>PC2|NSyrc-{*x z$ciTn=m6jSdhs*iSv}o(cJaCU#%YgbQ7PvU^t5;?ICGUaXYl1o%|~Jb?7O}(vv#I- z&*R(*{jXTh>z@fCvMt-`jmpXb5%+2USjhh1)Z?vo@S5G&>wWvhOW=6%pWOTrw>`Zx z^h+ezYQ5eI$^gq?lJ^(GU!w0`FR)6yZ%_Or+3QS$P1voY1&Lzvlcht9s}{TmidkF) zzLUA0O^Zop zp_~-mw~uqFM1RH@ZZIN>akerYp^w{ooh*y`im+`y8|q5c8WCd z3j1ID(9%%Y5QI?FC%KW3#`IB;WZvZb&Rw$!D!fk7O+DY!^*@msxS^}jJ5{V`cSiB& zfpH)iIfN- bxA6bK1^D>VL_l(A{Cxaq^z^c7a%le#YrXL% delta 13339 zcmZX)1ymeC*DXrW!3pj#xI=IY?(VL^-CYKEcXubadvJm~!QB#qI|1JC{ojA@TkFnR zJ#)IIYFC|oq`RxnWbmh`&`$|Guq>=BENrZag4q8n>25-?vG64GeWC9r1z0 z1yrDDsQ<-zrqJ9O;miMRj`{&f&JMdYT%;B(m!?C}*7EXa<*(~aMZ{R}_j50!MuvYb zsBD6hlT-1IU4rKu=Y1X^pR40TP?xbF5?_8<4~hFQIo61iAWGNvxyQTDgF&G2#;8-e z$gg-a&2j&juAhHR&^P)({*PDNy_-FqUSH0rFFfuaj+RjGtPr~T`W{c2yNn<1hwESV zqsLbe!a(b8&00falJ(}nX7>L1Vc)+A z%`UF@*S!}mwF7i|>ZnOC!9|@(f5TyYnrj3mSu9>nrhNfert79v*A|c@Zl}}5k?(yc zRakfa9@}YF2B9k{k-A))F;4Wh;xt4AX8^5;pe01n)T{?Q=vkNGS*-3f3zrCKMuT|n zD73{Y7zWJ_!sHkA#OJxNV&Hi+m^G-1PUFVkF(+j<5rmnGNO~0*9El%;% zGDoPIN3GONnoQdi2l1-PCQMpJI?EcgVXJeNvE(1pAYkv2#7)?S{me#BGD|gcaqZOm z7v-F)!q^!>F*_I?)#M3uxGO?vdjx=@ zC)gmgV)e|tFeuD-LWC|k9 zir)DK0bGJJTI2pyThL$!C_*04#1Oz246=%jzyMj^gpsb1S&vE$lST+^<2ktjkb8jl zZeq4^qMDG4YWdW|f1R}OZaOBCg-Y!z0h*(lu1ybeX31QFaxPxKquac0mfo{|vIQxx z?M;OnbslAQtK0D>AG15qP95@0=@J#&pXJhDMss>W`SDL(hgfB|RV7$lV=9qU=_)=k zn1EWW3tWG*Nk(1zQvUv9bj#9Fxg=na!=o;7(~?;pz+}upYc_H(=hL+6&-Axu8v76R zzAd~+H={7SgDl#fJ$zAzZjtoxHoAQ0k%M&Po}+Z+)Gm66r(TdHaXRVHA-b4LM+p5@ zM@aRN3%p1UlQ4UEbW4l3deKz2?(kYCPz#jJ1#M#^9>K{4m|sTjU)vy+(gnBNG;Q)7 zBU;vOMdP==@@axkiLWWr&UJUV1O8tc{@Aneb8CDSs;uscXX^E4`;LlVH(h;0Dq(k0%DyC_faOswi$uYB; zaO85PNuPOdj6&e$SvA4GwQ{yeppOQzY_E`=e$1-TndYlUQKl0O+Mtu0W9Y&)D6Z&} zzpt$h1XYxDPHE+~GgG24;mzmeOxJ1r(Bz#DNH5(rej+VmY=OLtZDV>ZVrJORQ+DRA zc4y3Lp?|Eus+;BJd*^9}zKK*0sMP_H@t<|^=4Y^;_ve?0N)>EOKrtXCkCe|;kD~xL zNj&2s5P$a~(cx|j{@XtK644>L{J;3mzvEa3pds)I&uPfzpc6V&JW@JXfkF}$eB%*E zao1=XTuNb{E>KS;`v;MfaT+uFLrFHJ==M8*BQMu*khudH^@rc zxnIiq99dW%ktVo|lK+7Cm|K4qdRyT>xbq?!)c zqii06YbK0aJPG(R9>gM%iDsRG*I=I_Zpk^WQuJj!2CbL@V3L9dl+Z>Sp<&bfmJ&+; z-_t3Y(W(#$3>4F~3_Yc=tq`cp)`%~`L^gmTW-U!Zr83Zx(wxdl+$`LUS=yK7CpRX7 z=)~yBP}4OOHu-wY#-DnVE`b4?fAQqy%$Ssy{PtU6qlN>MbRtb_9Jx9>#LBDfuN5go zOl~qK@P7}?o$0<6rS7cBk8+d;g@riZ&OcnRpRlrIB7;12k;jl2b`_7e60rmHsV)B)bJ%zQ!rIM!0!uj!z{dmY8nwr< zDw%^=IxgmqUh?@hwQ2-mrfrCBI9f5a{F|lTi$%zzT9BncD(O%;y4X!6ZlM4cxrQGd zR+b+fZNrmV5L&(r+Eh9Ee|tGn3ZyzN`O*nTC{hBN5b2EqYV8Ke)*eNbl-og(QfOiX zV#G;Sd&RTZpW)W7C;&x_MGAG{b07JTS=Ju%;)JY^;R)Vcr}}(;UMS41d6xM_JvF~v zo+BB3rBQp-wOlV>W{R@v-;o_b1zj@JYfd*qbMe87RV>wvHN`WdX& zK@EOxCcpW;zm6ugmT&bF`Tb-0^dGD@$E<|l6qV??dcXGwCtto>EK5TpINQXcVSYkJRFhH=JbwRHe z2|mF8e~cSAt0LE(Yh7w!wDPB4j|wx=okz2!#{{+xk&bDwL$9L;lj_c)*~(7e=dZf! z9l$jHJ<;2BrDGo%DSpVG@=rrvcLHZfUUR1|W8LBU(-R!XaE!26hEeg!hTW5Yw?BWc zA`dzxf?Q_Y_ra;rezH~Y2xB}GX;fy&Wx-JoGtO*qj@RO5$HXB~8Oiaw%s;|48SwHQaV#ALag$cC;O+_8~I z*{A~v!di{8(fzu35@G=X{@#m2-l>UInl9%D+ZwC^=uo^l1Hp~g`(GSOwAj-ys}+S2 z3d)1bw8}^*w!w)fqa?jg&qr?#xpLGU+yrK?{l|ZTSUPI=a4Q8i7{}bQ;+vEs7B3-& zK~iY*9UV9kp##2OE2X5;`-}BNVY5{apHB~`$_+s;)|!k16YKdTTkroaLA_OPvXm%> zuVxQN@nomc)6+281&=o~R2D%!5ef(Y{#ZRN^YG%h-t`$A1PrupE12p%{>ZOe&;L|D z@F`6kL)y0;^sw8@X_xB$esm;p*b-*f0IG<2bIxGR{{xE*cwS@(9k|HRPmZUbl5^0p4|rET0|qc+Nv0xw~cM{yHS$@@;{-U!UGy_pf@Z)?TZ#vtPTy<3Xgb z7Y>5{&wt$D?TU^(xdfm=4u9=^=^a0NSX$`t`FH!a zAKv6+aO>u}=p9lDvu5G>et-M|B6)By=o9FRQ4?|sK)!Eg-n(?<-&$jRQ|5sBUZNx1 z;Qp4Fb!N4+lkv=lSNjs44|?+a8KGE@bk+dRZ^K>F@CUIbdF8jzpm|}SPnO2+Qm6=S zNqbc`0yEP4?e!fOmxHxcFLs;YR-a$HA4u2p_%XtwV2{^0!9tT0cOd>$^ z?)|=gezmEv{dRT#upV8ri&Kqr;iwr~cb~<+IN?A-j2bT}DrWLZGmAFh!N3o*gB$4K=(H^clhD3e7aLe$8acvVN)q93`py5cf#TDCNfj7 zNKe9{4X2zJ5}$2@udt>b($Zp_R#m(XI<&?W@~f=Ey0;=G9qys z2e8}nHQu#(~#6htPOIf$#d4{?;;ZK+GZU*cc)=e99>EGkehz8zN6; zG2CD3$N_59or-1MK3AdZ@Fr3bRUwv@K^m=*Ri%HGhZSxXzva7@I-(P4NEKV60^Ije4Z`hhpr;7(GL7Jr0sA3@vN%1kmml218u7dr|-cMzP;KHoj!$H1_IQ3}nK2O!xk}%b81t74jX$ z`&~GKS4laM87!kRE!0}CZs9gg++Ru2(n&xRgPZO||H<40;;S;mE_{M-ooLV^kiaau+`$Ik8-$oLJb?yRhI(gyqWXwnypcs7bc7BYMdhr;7Kbeu&%k7I=e z5`tztp`l~gq5;S-e2e{-aK0F}kbx-tfRh7dHMWsZWKc0|@Jql$gL?QUkM#D+ zfCXK)FhJHMDhgMN@hbGLB*HEObmhH|^Bxw_+^;3UA=x68UryuoX7Hp3^C-l=vo*Y? zD8^G8*(;(Mw628LA1eW)TXUkU-Zc|0XOZadqdi0qP4RvQnOj6G!O1BdbC?VbqqNhD zEhdXxCI%W|hwa2c-&4d}ngd0C8jc;^7^WmjXn}5>TgP# zG3-2c-r-Dh-nokgt2iAidIGRt`#xyh>Nv&pgQ<~xO;pjwbMFUjNdjsG+T=ij9KjGq zrA_+R7@w;V5c~j&_m8J-h11W9Y_^NQ9^zCXtGaN6Ym-r+AZvo~tNl?zIs15_oD))b zZ5FZXyBl`MI>l$QFGxzohj_ubD!pP~1_H9S#J)fSjgf?#S;9yu;(!s@qedaHthmOI zct$r-k)|wm=3TvOnU>%1eimrGDX0*dD zV1FWLlGg*DL%bO6;mdhA#fDTXGqSG^bFgH6x>wc-@K$I6V1}FU*t|+)yvg=S2a@5E zfz0kMoc{FOX6Xb%%`=#AAQ=5?o>4Oa2!zW(b3qCJkwSGM-@tTM@~T272)e~hF^wre zy}^_w-d!(7SL&@6KrTGy^O|@7yjAV;&0et`5*Th@>GSqV=BRP@(wfdbj)-?IT$O77 zrFON?!c6SIrk$5Vh%1Z0Bu#Bo*&)J^2SxF&^ye|X+7QlhDba=%rWfJQqH{vZjC<^4 z&xSfwe%&QFCJ(!S_vKN2K_N;Q@z`0$=%tQgh2LIsjG4a9rTJar>%-#=*Fxfc-?mTC zrj6T+_zO75#gjhn;k5C4fJ^=fe;ew-ewRSdqq6=&G*2j@oH0L z@W8pDNMeik=`+Ih%=&RvP_all(Rc95ZP#}aLC-@g*$}&SoexXXm{C>3$t<9Q^rwg&Hl)+V^2tTi=J6f zH_>)Bxom>vEg$?BidqNI&vwjNIXi})kesbGOKZ%uO=(;j$bK)s6~XR%R4=a5w!D~z zPZJN%V?jTp|rO1uz3>HoElim}1O19&V2 z%=Tn*%oVy9V*9gp2kE#B%_&9TftgqRnZf{MgGp-r&Umd!l+-6GAcmzCzm1B^6 zVo@RMOxI|) k$!ZGr6&DyGwafDnpUksNhi%ZC`c_DE_cI&#c#pcwNX_4r2d*_;pnsaYx?r#wm}YS+TA3Q-o` zCZ_PLduuZE@K&nVQo{PfY+$#?UjZsrH@B(VKxn4==pa_RaAH6o(m-u#lIV0ghk3QY zTCav2QFEKDbfdHsr~Hqb{SRkgyHJER^1^@dECo4Sp=EFcjsS1boq1({eB_Rrkme|Q zFJrf5n^iGwc*@+zr0%p2Wvr=?dqO@v#6*XDw(FC^UE)H6=jpjypKH~^f94!gT(a8w zw6d-9i>0&|nMdD-_(rFfO9ERC&XNczorhoUYxiupgyPvCeBNnn;*t!p1!a6w(!_Vr z!Jc9)!-z`7gJs@qDdy!|n3$?l2E3Hk)~OtDcX%pmI>*tVsKQp{2BW)A`dVXiuP99( zv(O82`gD@^!LURE zH*6eVAQdiW^B;DO#0ybl^E|aBF6F?{By4^;BL6e9C`zO&Fnt1-8SH3e49r=Ce=+pV z;Sjh=O*bTgvRG2aQ&`UYwizBEo{&)z~hH1JSFkgv~IeZQsYm z^h*-<+4m6gJQq!PxEzA05Gglwl@E)CepqzJrj*KfWHmr4!>yI{g9ULfKCpR~Pe7K| zv0qewu>=SUmZ8eL1mZyAeI$aUByw@Q=7TZW1d;e=4}&W?FwZ0N*ZXq5<0e zg9o=!(*u{oWwbW*MC2UM71(;g$&=IxLz)^mjdRQBiJ0V^=yoX9egMuIH>@97V=qny zG$m^>y>=hPY9p$!k~QNs!X;H92{himzm2f~cKEJJ;6mAraD)T#fa?GlV+B!Gk=D-$ z1BBreA@<~^Azx9o#Yrab*~+GlSh=|jKO{y-k!o`S-~zVtW=$oZdt|flk65}@TXIwH zS2c(}lKO=PZ}WY`x>!Sq{_yiRp?7xvZjzd8`BI*I8<>^xD(yJQ`!sJbty= zh0MEjQ`?X`!){$tzEV*|4M(O#az5MC{S360F=Ze|HL=S7GzjmfP7wmy{Dp>SZARw(4!g??1? z8avV{292haHyFh8nVVI7Hmu#FLeNGQ96y302b<#CrBKrG2EG-~+U6Zwuo7*>#!_SL zCP1aI&fc@|R07ew7fDkT!&IMWf~zN zz(=C}#(h$pSDe&!Da9Fm(A5ba9nEzYnTQOK58UXNUZq;Vk~i40>O7?OVxNHuNfyOKVgXh`@Vz z-nC*P1PPi{i8n-dPD!?ZI%i{D2N?cALDNiZ&?I;fHo;9h3EvK3(D;L?#=o0P2Zllb z7Jf-SMl@CNZ=k_x5}6}7HZ_*R$Ip%l+Uj28LYL<8fp%Z*JON_X3ZxkE{DIBIBT9`P zMYi#tKgr+%4*yj(uzRBDSaKH+=ucKy=I{<28~(^x9$+CPhI6Gzg?-x(AXqR)50~Ty z#+O@F6YiM+1XTbQ#(?vdZ~G9(6L6a(s=P-(Soa#l%sBUO!~H;tQU|<%A-AiEpOH4#>}DOMxC5xPTzsY!ymsssoJRtT75vpvN^X z!%9RFcHxX{J*Tl`V5>1&Bi)GxtmfnzuaQ9-jBlWc=|Z28usIleyXgQMCz?e}9l#)C z%_3-aIzT|C2B1iuzT^!Ic(+o96_fx5Uuz}{8g9g(DzZZkR^V7bD6{+f1608U78uFC zsg;5^og#ur)WNT&UUW52c9C_jW+hN*|Ljuct7%U1*3fCOImEGWznGv5X z0OxI@ir1yYBr1eTI%ed2b)+eR9BD7fW$kvx!g9oIjUCiiLtN#q$F209>8(8rfwwm0L3qGa)B(gVN6II01Z)Y;ZwL*;u^YvcoC11-wWPfiRT=0?$ zu@O0Sl)|DBsub>3bYqjP8PB#4cXSP9KJ^F$Z%kZyt9f?Bk4jB|(l*jBECxrV7WCF)JJQMX$x(#lxyD9*Xt}(N{QM&w% zO~2vOw&CaMiT~sRvuc4R(Fdl>%4;bfPH`@!Rb_=MQho2Cb&;!} z{O>r76`jK_=Suf9<+3tc-xZCkb{J3TRg|qwc=&p#L`iSQ0CC!_zni9EA0hGE&%NG` zm|C#;=IjW)EZ+!4Sq-r$6HFb(wZAK_X^ym>L~Iugpn=Y3OVIBO4(`5l1{)wZspJo+ z_4CbWmJ}_Hu<8oaV|`@}!cm8xxXXl9)cheQOJ1z_^?Y84tqa}N6c9qp0We#wP;3u$ zCnT7-{QT>Lw{mFmAn_s#y$NHA2w(QDQqC;nxmH@dGS&>2v?D!=7_Yxcp-%HRMv!SY z(9HAshzZn!{l~J52ev0j1{i4HpfXl4=w2#7lwT%h2(X{RKba@Jj;*kDPbbs&FEYC84p9x8J$k>cUd~TBcy60gjsS|_+ zVx#_F$FY&XF=J9?8%UptnzpO3%{!;K77VS$*s^Zn&Ezp9kTZ3{&dn;WWh9Hf}p%$@Fu`A=JL66Zi2}4VM(V(eMwvdwOtD7qT z+tX}ODiBKBcffpC8Fu(RFH{7?%0q&1wU$VMkZRO+eCex9Lu**1G;EBtAXW$sRHhMb zrb!{CK%k7k4*eD)QHX0?q+~UeoWaN`+}v}X-V3vK{zAD%MyAkA5Oi`rzVn&<-`{DY z62>BYXon&txjKnLG51K3$TLO+N!@My9*p*#}(r)ql}a4Vx@T(etm^9=VdE2rc;@| z$V{eJELsft=;oD09-xG%m_$}Rar4Y<+Ig@H_Z)+^57LG>RPz{=a*ne&rp#gC#HCtZ zKz`Nghcfk5_UK*_do*tgWDeKh8U@c}30gr;j>4>vq5D zt=tonp8Y<7W@GAwi0t9GNB?QSwf-L>)GWkp)dDxoQ|;jLVo6pNBA&uElzMG-`-!Xe zo`7d8+wApDpPyq3chy5AaX4P+UN3Zt3iZ?L2Zg$z$83HK9$=LZZ10rhK_IezaB6t5 zGZ`v_v#eSTJ~}SY{=g5WmFS>rpZl7Tx9lxdKLtW(DzSt#+Ipdq z`|oIgTxF)lx>(k7BXp)R(~#sn@b=0cW-|~Fx9rlbOROJmX+thf*@7nUT^~~8dM6g( zBG+XdPlr+hKIz}6AcnUi4#;bD$Rm8Y!!P+GNh;YKQL?rHR`WTV;U zUN=ldP(_|dmO#W}Ri0D^{)ARfdk)Zm%feE^&@=gFo}utf0{p1rO@db9v1h3yE!LLw z^&(S+e}Ta72h1~7VYs!wimtwO4ZD-Dnlt-Lu9?xbRi4Xdtvb= zL#Lws0tEDD;>b;1ovUM@@t4Y_n^8fEevsG0oX-I)scDd?epSX^^G`)0R9S9xLk8E~ zzu>q{P%^GUO)0=TgJ0f*EWC?RS@a2Dk^(6;Qh=8BWZUdzZ`D8pt&yv@EtHeL* z3py$}&)RWA{%2>u^;zggb%0kR1fCmPmTmLrcX%Dhm-NRL>0*0&KcS2PP--4WCC!o) zht;Gy5JTJ^D@hK*r2*AZ;@1p#xoqfumaa#%;X5wt15NXAxOfPg-PDwyH)U$wLNvt@ z3{GMN5+HMxA6NkSiyG-V)nr#Aao=BtFq5p>ibT4fScbomNafYeF=HZ6L@F#?g(Qw2 z@lb4yp`zslBX{Pzm7=nMRNXj9mP^5xo9m4AbV}%?r+rMwtIo2nc{qxFsl04+!k;b} z-h1Mpe#D+R@B*B^Ey~2PDx-TBZ61U7MT${PqoPt9o9nmAm#Yf?$Qc#US6luc%||W_ z9hD$hc{}0@{PaThVGVkj2AUat8cV!9ie?eb@TAzWB4{ZgK7X-*XolIxEW_vMv7Q%yAk zP=dzjYDCCt?&wUF$lE`37UQzB9*OL6E2|z1v)BRe%1MyyL@C@hz|uxqG2-=6vb7|o zZk<5Y3n-NL+$`AP`GF%_4RxhmK{%H+aQ5Rhimsp0fGP^%zv39rkuJ-6!+op7dv{qq zIGu0nB&1~pN9jjl`c;ice+yhh!?;0jS^+}f8x4|#7|%DTpdCH;BPF7RdMbIK=)ec# z?L`=VA;L1>bY#QqaDlB!sRkENcgXF9CIh-pM=eX`#2&-Qw_XF%f^P*0k-eye38uae zr2V2Ey$Ok&+%6Cq<3|}R!xtX-*|$ec$g+L`nL)`ah^+;R>S<1YK#z;k!|_NI zdG#QqyT1ea>|LP@q~Ap9u#ZC&P2*1$Nr8zCI3e8E$X9 zPszAtMK>wC?>;@1jOU?XbN%0H4#Sh+BOH(#v4%pGK(+Y*+*|KUN4_PqCrQ-L-!t!> zFFTImsuKj$Gl)Jp(NhF7Gl+pLamTJ#Ak9O;rU56DHCO%E{beGMeyR^It)qIX&!g?O zr+(^773sBL+qmc*;y|i8G1-=rhk9zb$jVOx)kCKOo7dXP>ys<0DUf;UXoxI_rn(m) z+qqe9*|_SRj`8bJ7M-H;K*j}!dfrab4#+;0`cY}#Q9m`fXF^lMD=Zn|;~PZtH&xvL zHU|x)kt(RFI$KC{%pwGxlj~3*fxxkQsF~gZw5Z{wewLlOk}7N>-I|73&S)GM$>O&8 zF=H39fV^;CW8J1n+ctf7hE60Rvh@|3m%dN1$M*(n!>GaM>G98bS6$WCX3szKhx6*+ zCKm=<+y4J<*I&cy-Wq*wB z@KLGS8lb)?Q9BB|);72>@a@_nr6=j-D;F+0X0A~QM(fxho$%<+HQJkmHRhZQ7J~cS z;7sOqrRFA<0h_4y-T9|vNcsiE7!#+*llT=W5vPvewxP{lIZNgvu;%wwS zsx{qyWRZh4KxUsAr$`zogS$yMBliBsJ|09Z;jinwlHU<5$Yhnbn@ZD?h8+DW^16#-uZy5~s z3h&5{y0I>+rSTWY3fa8IWu~H^8RrKnYdh-sIA%q8)JX}Dt!;t|{{_^(ANabm(WI|2 zp&WXs(76;FnrRp)#wH|hck5Z9DQj<#>eAUf?e$~G<*bs*oBf1R`>MxNlG6S=YkV&% z^?x+4V(=33P<0eglToPx%haW;Tq!fx-bKP`-`*@ZA}Sc{6YNqsKoj4E(>^o?kAy4n zZg-Xf(K;=0qz+@zS85$5p+=T99`s9%inEZ=QZltvGWb)Mv-*7C$e7ukp&42w@b~|2 zuloL(ODa!E8xfI-I)i(a6y_n3Aaxj3>P{#_ELEkyi!XMZ4o0J%0;`Xb=_$gJE}A6a z_pUE^y^CdykGVNA8;$~le=1>D+m9Wagh6U%%E7kP-%G2X6j+?hgK__Zq;y!XQdGr8 z7fq+3(lAa*@>$;O3QFZ`PJ9Y_`kERx*Z`3P(-xj!1ISvSOk1O=YcGniaVQaOt!jMyR7pHODibNgK^MJ6TL~H>K4$CA zvrHIP6msB{8uKL9P{Y=P#2JTAZE+>Pg{M;BkfN#0#~lZgUbquTr6NN9ZJ8Zq`WYn_uP9gM-v z<^$IQ=*IAou^zxzq7lq|vHj-zzJJqQvh-{B4YV2la_l8tT{k#*<@r)-$cpUv*{r&D zWqYIBE6<~OjraeSe;L}p2Qe?`ynrj?Udy+8*E)Ee?FIc$3o(ow4lPIjxS*){>sF7w z|9c>L>&8-gJ1h9Ip``sAnsCOgz*mredHmm6`C3$4c=)>qVwZc~Scg&C5wjr3E^VG= z;P2NxYKnNnn)Q2(G_j4pNo;d##a;nad_4T)QX_=+z3`l{NNrO?1(lvjO4pwY%M6GM zex~_aPL`^1++8aua5U`4Y2FkD%KmZnM|_d_8#Q;lw#oS{+cHrk{`L>c04jvar9>8R z#Y+I3I_ii%%p^+i_fCb*bD!M;u@+MMcsWSd|Ix@oh)dYx%82BriR{@eFrpIF?RKyG z5#_n2N3mt?Qq$|O4j*hXRD6o5n&TA|dG9^T%bP0~#^%4YFBM!@d$}$ zzpQL&xkts9Tb50Fa_964CmF!7gS@zL%Lhf{%9sifl$Y>))|m+c8I`7iN}k%&dR?iF zd4AGs@cg`M#S&J;H=0eBhj0wu#dx5?6>u40UX+In9Pi@e+%^QveZ!Rl3ZUtgR9O^ zaE3mM!NJ}>&&6c${R9AFeiPgxT$^qhWTQtGVCUc_W%>A!+LI#uNBF|V!T#m{Y2)Nz l`+wWGlNp40vDvu)pYEKTY{|aD5(r#u90-(@;))Uo{~sR=c?$pl diff --git a/vignettes/articles/HMP_2012_16S_gingival_V35.Rmd b/vignettes/articles/HMP_2012_16S_gingival_V35.Rmd index 0ba40d7..fd80987 100644 --- a/vignettes/articles/HMP_2012_16S_gingival_V35.Rmd +++ b/vignettes/articles/HMP_2012_16S_gingival_V35.Rmd @@ -200,6 +200,40 @@ direction <- get_direction_cols(DA_output, conditions_col, conditions) ## Enrichment (adjP <= 0.1) +Lefse is the strictest method, defining DA with p-values and effect +size thresholds (based on linear discriminant analysis). CLR affects +the magnitud of the effect size, so let's define a new one. + +```{r} +DA_output$lefse.TSS$statInfo$abs_score |> hist() +``` +```{r} +DA_output$lefse.CLR$statInfo$abs_score |> hist() +``` +```{r} +c( + lefse.TSS = median(DA_output$lefse.TSS$statInfo$abs_score), + lefse.CLR = median(DA_output$lefse.CLR$statInfo$abs_score) +) +``` + + +Create variables of thresholds: + +```{r} +adjThr<- rep(0.1, length(DA_output)) +names(adjThr) <- names(DA_output) + +esThr <- rep(0, length(DA_output)) +names(esThr) <- names(DA_output) +esThr[grep("lefse.TSS", names(esThr))] <- 2 +esThr[grep("lefse.CLR", names(esThr))] <- 0.06 + +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") +``` + Run enrichment: ```{r enrichment} @@ -208,10 +242,11 @@ enrichment <- createEnrichment( priorKnowledge = prior_info, enrichmentCol = "taxon_annotation", namesCol = "new_names", - slot = "pValMat", colName = "adjP", type = "pvalue", + # slot = "pValMat", colName = "adjP", type = "pvalue", + slot = slotV, colName = colNameV, type = typeV, direction = direction, - threshold_pvalue = 0.1, - threshold_logfc = 0, + threshold_pvalue = adjThr, + threshold_logfc = esThr, top = NULL, # No top feature selected alternative = "greater", verbose = FALSE @@ -293,21 +328,41 @@ enPlot <- enrichmentSummary |> ## Calculate TP - FP ratio (no threshold) ```{r} -positives <- createPositives( - object = DA_output, - priorKnowledge = prior_info, - enrichmentCol = "taxon_annotation", namesCol = "new_names", - slot = "pValMat", colName = "rawP", type = "pvalue", - direction = direction, - threshold_pvalue = 1, - threshold_logfc = 0, - top = seq.int(from = 0, to = 50, by = 5), - alternative = "greater", - verbose = FALSE, - TP = list(c("DOWN Abundant", "anaerobic"), c("UP Abundant", "aerobic")), - FP = list(c("DOWN Abundant", "aerobic"), c("UP Abundant", "anaerobic")) -) |> - left_join(get_meth_class(), by = 'method') +positives <- map(1:length(DA_output), .f = function(i) { + positives <- createPositives( + object = DA_output[i], + priorKnowledge = prior_info, + enrichmentCol = "taxon_annotation", namesCol = "new_names", + slot = slotV[i], colName = colNameV[i], type = typeV[i], + direction = direction[i], + threshold_pvalue = 1, + threshold_logfc = 0, + top = seq.int(from = 0, to = 50, by = 5), + alternative = "greater", + verbose = FALSE, + TP = list(c("DOWN Abundant", "anaerobic"), c("UP Abundant", "aerobic")), + FP = list(c("DOWN Abundant", "aerobic"), c("UP Abundant", "anaerobic")) + ) |> + left_join(get_meth_class(), by = 'method') +}) |> bind_rows() +# names(positives) <- names(DA_output) + + +# positives <- createPositives( +# object = DA_output, +# priorKnowledge = prior_info, +# enrichmentCol = "taxon_annotation", namesCol = "new_names", +# slot = slotV, colName = colNameV, type = typeV, +# direction = direction, +# threshold_pvalue = 1, +# threshold_logfc = 0, +# top = seq.int(from = 0, to = 50, by = 5), +# alternative = "greater", +# verbose = FALSE, +# TP = list(c("DOWN Abundant", "anaerobic"), c("UP Abundant", "aerobic")), +# FP = list(c("DOWN Abundant", "aerobic"), c("UP Abundant", "anaerobic")) +# ) |> +# left_join(get_meth_class(), by = 'method') ``` ## Plot TP - FP