From 3edb69376da291a9a864ba8b0e62ce3d44eec3fc Mon Sep 17 00:00:00 2001 From: kevinlkx Date: Thu, 11 Jul 2024 11:28:26 -0500 Subject: [PATCH] update tests and documentations --- DESCRIPTION | 2 +- R/ctwas_finemapping.R | 2 + R/ctwas_merge_regions.R | 5 + R/ctwas_summarize_finemap_res.R | 5 +- README.md | 17 +- .../LDL_example.ctwas_parameters.RDS | Bin .../LDL_example.ctwas_sumstats_noLD_res.RDS | Bin .../LDL_example.ctwas_sumstats_res.RDS | Bin man/merge_finemap_regions.Rd | 5 + .../LDL_example.annotated_finemap_res.RDS | Bin 0 -> 91210 bytes tests/testthat/test-ctwas_est_parameters.R | 3 +- tests/testthat/test-ctwas_finemapping.R | 6 +- tests/testthat/test-ctwas_get_region_cor.R | 3 +- .../testthat/test-ctwas_preprocess_regions.R | 2 - tests/testthat/test-ctwas_region_data.R | 3 +- tests/testthat/test-ctwas_screen_regions.R | 5 +- .../test-ctwas_summarize_finemap_res.R | 38 +-- .../test-ctwas_summarize_parameters.R | 6 +- tests/testthat/test-ctwas_sumstats.R | 2 +- tests/testthat/test-ctwas_sumstats_noLD.R | 2 +- vignettes/ctwas_main_function.Rmd | 100 ++------ vignettes/ctwas_modules.Rmd | 220 ++++++------------ vignettes/minimal_ctwas_tutorial.Rmd | 123 ++++------ vignettes/postprocessing_ctwas_results.Rmd | 89 +++---- vignettes/preparing_ctwas_input_data.Rmd | 204 ++++------------ vignettes/summarizing_ctwas_results.Rmd | 117 +++++----- 26 files changed, 346 insertions(+), 613 deletions(-) rename {tests/testthat => inst/extdata/sample_data}/LDL_example.ctwas_parameters.RDS (100%) rename {tests/testthat => inst/extdata/sample_data}/LDL_example.ctwas_sumstats_noLD_res.RDS (100%) rename {tests/testthat => inst/extdata/sample_data}/LDL_example.ctwas_sumstats_res.RDS (100%) create mode 100644 tests/testthat/LDL_example.annotated_finemap_res.RDS diff --git a/DESCRIPTION b/DESCRIPTION index 6383c2e4..0406da0f 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -4,7 +4,7 @@ Type: Package Title: Adjusting for Genetic Confounders in Transcriptome-Wide Association Studies Improves Discovery of Risk Genes of Complex Traits Date: 2024-06-14 -Version: 0.3.5 +Version: 0.3.6 Authors@R: c(person("Siming","Zhao",role="aut",email="siming.zhao06@gmail.com"), person("Wesley","Crouse",role="aut"), person("Sheng","Qian",role="aut",email="shengqian@uchicago.edu"), diff --git a/R/ctwas_finemapping.R b/R/ctwas_finemapping.R index ac2db69a..8d3cf20a 100644 --- a/R/ctwas_finemapping.R +++ b/R/ctwas_finemapping.R @@ -111,6 +111,8 @@ finemap_region <- function(region_data, rm(res) if (!use_LD) { + # set L = 1 in no-LD version + L <- 1 # use an identity matrix as R in no-LD version R <- diag(length(z)) # do not include cs_index in no-LD version diff --git a/R/ctwas_merge_regions.R b/R/ctwas_merge_regions.R index 5e0d2308..66008d9c 100644 --- a/R/ctwas_merge_regions.R +++ b/R/ctwas_merge_regions.R @@ -45,6 +45,11 @@ #' #' @param cor_dir a string, the directory to store correlation (R) matrices #' +#' @param LD_format file format for LD matrix. If "custom", use a user defined +#' \code{LD_loader_fun()} function to load LD matrix. +#' +#' @param LD_loader_fun a user defined function to load LD matrix when \code{LD_format = "custom"}. +#' #' @param ncore The number of cores used to parallelize computation over regions #' #' @param verbose TRUE/FALSE. If TRUE, print detail messages diff --git a/R/ctwas_summarize_finemap_res.R b/R/ctwas_summarize_finemap_res.R index 54ed342f..5dbe1e61 100644 --- a/R/ctwas_summarize_finemap_res.R +++ b/R/ctwas_summarize_finemap_res.R @@ -70,7 +70,7 @@ anno_finemap_res <- function(finemap_res, } # extract gene ids - finemap_gene_res <- finemap_res[finemap_res$type!="SNP",] + finemap_gene_res <- finemap_res[finemap_res$group!="SNP",] if (is.null(finemap_gene_res$gene_id)) { finemap_gene_res$gene_id <- sapply(strsplit(finemap_gene_res$id, split = "[|]"), "[[", 1) @@ -115,7 +115,7 @@ anno_finemap_res <- function(finemap_res, # add SNP positions loginfo("add SNP positions from snp_info") - finemap_snp_res <- finemap_res[finemap_res$type=="SNP",] + finemap_snp_res <- finemap_res[finemap_res$group=="SNP",] snp_idx <- match(finemap_snp_res$id, snp_info$id) finemap_snp_res$chrom <- snp_info$chrom[snp_idx] finemap_snp_res$chrom <- parse_number(as.character(finemap_snp_res$chrom)) @@ -142,6 +142,7 @@ anno_finemap_res <- function(finemap_res, #' #' @export get_gene_annot_from_ens_db <- function(ens_db, gene_ids) { + gene_ids <- unique(na.omit(gene_ids)) if (any(grep("[.]", gene_ids))) { gene_ids_trimmed <- sapply(strsplit(gene_ids, split = "[.]"), "[[", 1) diff --git a/README.md b/README.md index 6a1c0061..8dd6d5ff 100644 --- a/README.md +++ b/README.md @@ -26,29 +26,30 @@ Running a cTWAS analysis involves four main steps: 1. Preparing the input data. -2. Computing the gene z-scores. +2. Computing associations of molecular traits with the phenotype (Z-scores). 3. Estimating the model parameters. -4. Fine-mapping. +4. Fine-mapping causal molecular traits The outputs of cTWAS are posterior inclusion probabilities (PIPs) for all variants and molecular traits. To learn more about the ctwas R package, we recommend starting with this introductory tutorial: -[A minimal tutorial of how to run cTWAS without LD](https://xinhe-lab.github.io/multigroup_ctwas/articles/simple_ctwas_tutorial.html) +[A minimal tutorial of how to run cTWAS without LD](https://xinhe-lab.github.io/singlegroup_ctwas +/articles/minimal_ctwas_tutorial.html) To run the full cTWAS, a few more tutorials including: -- [Preparing input data](https://xinhe-lab.github.io/multigroup_ctwas/articles/preparing_ctwas_input_data.html) +- [Preparing input data](https://xinhe-lab.github.io/singlegroup_ctwas/articles/preparing_ctwas_input_data.html) -- [Running cTWAS main function](https://xinhe-lab.github.io/multigroup_ctwas/articles/ctwas_main_function.html) +- [Running cTWAS main function](https://xinhe-lab.github.io/singlegroup_ctwas/articles/ctwas_main_function.html) -- [Running cTWAS by modules](https://xinhe-lab.github.io/multigroup_ctwas/articles/ctwas_modules.html) +- [Running cTWAS by modules](https://xinhe-lab.github.io/singlegroup_ctwas/articles/ctwas_modules.html) -- [Summarizing and visualizing cTWAS results](https://xinhe-lab.github.io/multigroup_ctwas/articles/summarizing_ctwas_results.html) +- [Summarizing and visualizing cTWAS results](https://xinhe-lab.github.io/singlegroup_ctwas/articles/summarizing_ctwas_results.html) -- [Post-processing cTWAS results](https://xinhe-lab.github.io/multigroup_ctwas/articles/postprocessing_ctwas_results.html) +- [Post-processing cTWAS results](https://xinhe-lab.github.io/singlegroup_ctwas/articles/postprocessing_ctwas_results.html) You can [browse source code](https://github.com/xinhe-lab/ctwas/tree/single_group) and [report a bug](https://github.com/xinhe-lab/ctwas/issues) here. diff --git a/tests/testthat/LDL_example.ctwas_parameters.RDS b/inst/extdata/sample_data/LDL_example.ctwas_parameters.RDS similarity index 100% rename from tests/testthat/LDL_example.ctwas_parameters.RDS rename to inst/extdata/sample_data/LDL_example.ctwas_parameters.RDS diff --git a/tests/testthat/LDL_example.ctwas_sumstats_noLD_res.RDS b/inst/extdata/sample_data/LDL_example.ctwas_sumstats_noLD_res.RDS similarity index 100% rename from tests/testthat/LDL_example.ctwas_sumstats_noLD_res.RDS rename to inst/extdata/sample_data/LDL_example.ctwas_sumstats_noLD_res.RDS diff --git a/tests/testthat/LDL_example.ctwas_sumstats_res.RDS b/inst/extdata/sample_data/LDL_example.ctwas_sumstats_res.RDS similarity index 100% rename from tests/testthat/LDL_example.ctwas_sumstats_res.RDS rename to inst/extdata/sample_data/LDL_example.ctwas_sumstats_res.RDS diff --git a/man/merge_finemap_regions.Rd b/man/merge_finemap_regions.Rd index 08635049..114c7583 100644 --- a/man/merge_finemap_regions.Rd +++ b/man/merge_finemap_regions.Rd @@ -78,6 +78,11 @@ genotype data in genetic studies.} \item{cor_dir}{a string, the directory to store correlation (R) matrices} +\item{LD_format}{file format for LD matrix. If "custom", use a user defined +\code{LD_loader_fun()} function to load LD matrix.} + +\item{LD_loader_fun}{a user defined function to load LD matrix when \code{LD_format = "custom"}.} + \item{ncore}{The number of cores used to parallelize computation over regions} \item{verbose}{TRUE/FALSE. If TRUE, print detail messages} diff --git a/tests/testthat/LDL_example.annotated_finemap_res.RDS b/tests/testthat/LDL_example.annotated_finemap_res.RDS new file mode 100644 index 0000000000000000000000000000000000000000..6a71af1835d9641d352e477464f262e8549bcd81 GIT binary patch literal 91210 zcmd>lhf|YJ)NTrdPLPh$lnw$B1f+wAD2S+ZkrJtb$WQ4tp$H0S1VoxrL_r`FK?G?D zprMmk=t80>L|Q@%CGGOv`%io`Hx zrf4;FHsyTi#0tOYvMZLpr2FS(I7F#dHB(+*{<5_N@Vd?;WG?3g#c+X3kfH0>FAFPa z7tJ#x9a_KXO>HyB!^igq;>Y&{7nU0qJ_Q&#xw^VKyE?bw#e{`(H5KC}57^RgBy%h6 zS3AV`w_LrpY@9B+UAlYe(p@zV?)&TmLWbE+1vHUtC3pD0>n98IkiAgg=4-oTE_l#_ z`|(nM_6fFhIbPx-(e?p6p)E*aDUaUK6tj9#>-4b`l)h$D5;JpP;q`yRuhuI8UJq7F z9vJMo0_ysdW;wByEKQ{HQ6Xw!HsveD9<(#=C8p5WhI@Zevcr)$ZOa`^zRP{97w zMPyqBpopv+qfLrQ-Tc1QpxHieOD2g#8wH3Hw&O<$$ ziU6`U%z<4$5gn2~;yV;e5b2lYjVNSI9>9gG%akmEm)}|L#U=pEkx3@|yi5OjZ;_W% z#n!s>O@a4PW8;L}z*O0IiY#ff*O$E9FBC+9|GR|r_c4?bg!59x|>45p&gEFzaTCh>+K z_3%PmR<9ah4#Xzby6~ zd;u5-jeZl;&<8IpCHKBhod?h|Hce#sc1>PER>}XOFZ1`An0byRnRD;Q#uX-c;;e+vN%CGx-!*oNa8qg7 znFO;sz|I!>6lu4JdzIoYGMN-o-Wj;=ZhqLjM@)~6(}$4BGm2PBOA(S7xayODm&C+> zfdYeT%&-I_&^&;li|aI6V%i|*O~S(&?~b}y*$hE>cF-cdbDvTH8u%5uni)h zmO0^xrj39dYL_r!C|qiZ9vexPCK2%k;6zjJfx}3lLpFoIjxC2}ZyUqLmr--z1RlKigT`xJe54TRZVYFSFi3reP;wP!4-M?D2 zy?c`j_eloY4!sPF3E$4%mNaSmuk|$E%L~*Ie*hg_90^4CNLbit>(0dJRuH9rQKK)O zB^f9PQcJGOOj=*GC)-?1B;ln50OnHmJ|k`n*x^y0h`b^RvYM_rNv%abq0F6sGBUtW zGCCj@8(77`P;23(Hjd6B+eQzYv~)c{rMe!#SN{z}K1t0+?z_i`Ho1?0GjfM|1C+K46e)K;JlzYnzaeIbBY`%Y58!t~X4^bj` zP6khI+?*svvSEWgz6`Apk<|h5$;Nt%Dfc3P8;bnf19F3# zU~TyTI!w(AMLgGhiBfCPzDgF`ycV6t*qi!T$X}rTUBM=$)j)bP!1xE(^y4R#))qj# z0c1K*0c$Wr56b`~R4b<5O|U8jNe|DmpCuX;PZ&S_75eiWvW&)Y?xY;{GP-*P-4lyu z^G&yl%W}5I0I<&d7;~M3yFSFw+c7gkxJ1Ct;eebv@gAf0wV%Uo)PG7##@U{xDpTpuOZ-@wqWLPu?i-Phpc}k zmZwQS9!eXkfky($-~tuiaQg}jTyvTd>K4B~4qoHl7hLU2qCgg-X80d6fsJrVKO-Yi z6fa*F{Sx?}37XQx1CxV`%tSjToArsmKb!);N(%*aWrX@8C@BZrb0Od^&K+r5Vd5S? zSNVDcWU*_@8voeA)qh=sl(V8!?Y(@oO)`Ae+$w|~%rHmp@hJ?wYF>8F& z+YkBWMmB(qoZd1ZV#c7)B!JjA8$R9ss!=j67T>ccl)A=buEpapEpmQ~=OC+!^!w#g zYk@N+w=yD2m~u4PHMV1mIml>?#gBWXz)`x)&Dpgyi1Ue(jDGqOxhMVB=hd?6oQ#js zJA(|sc9B)bQ@9>U#c(YM85xRv0>5S00Z0mmCB3Ya*fF;8Cg~b3_)xm}R@DaR=+a(f zhIBs?`bYA}3gg)#*yyK^*C>~<#0OWr0RO3@+oZoIkB5&UcVGA(TtBq#$jr6~2uBW{Xunq~R+e;=Rr(fhm7F#-vij6}DBU;R~pqj-y-&Ap_?2q`b49{Fji z7;J;;Y|;2@)yFeE3}li<0Ws3pxI;wI6ZY}dJxv!Pcyc5Boxv#;0{ntg6a+4N>UOOT z_m0ppzHQO#LR86|gy+3{fDr)m0?bh9HRB4z-I%wI2r?774&J|V0bHK(=%R(gF?%@U z*emzW`7B_94E4m2#SQymj*=`sh+u*mS+rTl(dkO^Pd4V6kc)+o5 z`#F$(UcW)D33cIg#zD550NvrUk$!`{1=q0*cT*dSdr!`d;bO(fCx2^~w>F_X~?_;R_NZKUiejP8-qp>Bbgs9WIur#}HPnLi=SrzY6aGJm3}7Z?V{ z$qmskFDt0g<3WX#m-k`bhIvTA2ja-cvKFBo<3M#n^i{F_>sLi)N8DRZr^3cxj!kdV z2>6#W@a{irt|gWv4?`kH!csy(tFG|UnPnBMN_f(6DT?b5Zzw>`35Tri@w(uZv99C3 zYu$WkIf>+Z54<2fueP;T@d~7jPYiw(+Ugn^G%U%fB?7@U@ zQ$@c#!LFf!aIYcS$1GeR{a{etu@R#6ux8`Y|!QSW|R$X!>O z2?rmF`1hM!h-~L-;v0>2f=(k|FSG)CxK@>dVISo5_zz`vPsX#PM&o~cnB_l@c)bKl zw5d;LT?{=va^&|6;wgZmwS1e*D#rol*&Tf=2w;NZiqfc^6J$r%2yKkL6l7o|ks$Yq zYe}W}-|E-D$c@*Lg%olvvYPw?7WS2i^mAYr^spq z4B2{(;*T|`L?Brg8-ZOH8~vHn-1|=#MrRE{5n_oY-Ft<*3SjDqJOQeV2yk06kFN*y zIDPZWXRvwpWn`p`qxAmuE0<|W>U2{9vdTqd)a}oT5LcbSQajj!Bqv#^3XEPXl!al= zP~Nln(+)DFK>Jw-k4($BBSf|hsr7~bU>k+8=+fP$joTml4j;{o%y#o^OKby{l=LC{ z|J^k!|32*p^x-7CqJgU>;^Q+_xJgI?PCtdIL)-O(z`br0)1z*4ooC<9=uN;~1(=B& z3;d~?mTE$M(MpRHUTtDAg!|;{!k}~}RTOE1EcL}>+i&ARUcSAudU{&?hhh?mK%WE? z?B|pG=jr#Al1qQvf){?@08;Kl2$8o!vA(w)Nxt{T+?eY=CiffP$I&h5_wSl#(>V@h z_SDd+?Ym>Zd;Y{waL3a*!Nn(ZF|zW8+K#cAVyR_$d_3kg4KRn2kgECQpahFV{DLiSB(&l_*g(R~b zb5{J(V)<)(a&G4OAhSip!vmE0aH*UJfO?K}hQ&d+&rqY$t@SrXLL~hb6L}IeqPeLR zKWBWB9t>b(9u%Z*{?&p~mhX+0wt-j$of;&yXlUWy&%|=(PYTHNZxWPjzGp#8kgQ`{ zxe%#({C#m_rnyYkv5;t!wSrGYU1jkfvGRoI4MSwv&v0PaCyYN{qVI6HTX520BNT0J zoUHYL0%sbH_pn8lH1|wN#-GaiSw+TtxB?Y^vZlkmXA0Y*943F%_#<+v6-NV0zfluQ zo8?`X;0daz8!s&18K58uDh|#;1`c}V11935A*FE^+OaFT)6ZI`GIlGq;{}($J8pN4 ztRTt~k4WC-Qs^luj-KMcAF2>{(?y@ezeTu}qFz?%cNvmN`IHGCJaO5*ug9bJ5W`D- z8E6{M(mgSNU+Wx05=}g&DRa~3q|?h1(XSNH-qn&xcpZ~hzDz^tPTsPEi}_1)`?2sc zCJQhbS~?o9HCiIT8b85Fzv6|f-ok-B|JtV>40{N!jCk;@42`Uf$JPN~{jCSDs)ps= z29)#}jm`1n3flbRA@GD)kN%jgkKnkGIn$$?W3emIy_=bQb*l%Opg{?mM5!G`hE0N7 zN#Ip=N(Udnpe61jh&UWCzJFn1dkn||Z|Cf3trD&D8zN{6+d%&F%!Aj-=6^<@F?d(* zN&lW31bl*+4`SF=Y-K3W()3{EJ>~>hcDprY_j#kijy4fiK9Hb*l-<_!BWi*ee?36k zWt0IRe{&R68xZpo7$Zn6f{r^>EzvbfE#0CIv_Oe}JHb1J<{Z>wb1s4fZ-{& zba@Gb*u+X0>GAk{_$fdDJBCJ-z{`dyO;Liq2}|Z&3=MXqB!X3qE)7H@+>v`}cVR@H zhmxmxcM7gSCPP!WnfCAax1SOCsjmn)$g~vsWCa&BD{TR*EqvIQpu`WSJV}eaEXlce zMUpR~5at0Fq82xaP6qbzO4{_ z26-wg_qB^;ftZToTqa5{9>!ip1DIDueN0t`*TMuQ5x1S)0v}HxSt6@V#FAJ3R1j7i zl1gBIc0rbX*0&%&_Jo@}Hcb{3rZ=&60(W~VVS}?W{0VCJkC6B~$c19ibX8d5ecUA2 z`_fO`3J<_jS?BQ5&w)dP(2~NY(5S=jFn%be-{=a1j||&9MQ5pfD0nV28FoXqeC|zJ zDSuykUgl4*8$YhdNr0?eqnB-7qe5MNV*qF6Li)NC6M41^&x+u3&XzD)0;DrMt-{fQh8+2qweA%2P;*#D+rb6UY zeow&w7-EI|F{(|Hhu%IzqR2WQ8BuPew2KdPFq)Irf`K;E1e=qjAxTao)m64(ArXXl zAE(vC^3I6WqCyGRlRMD}C@_rE3VD_HBCwQCJSC3rMOv(y8>)LsXhm}xG8x;e1uIQ4 zMVj%_QZ^-|`PW6Id70|_hS`mzf+jQGaju?g)DxsBPFJM#wf~gj!!%d%afM|mB}i@Y zpZttElK7-{Mk0cnwwqF&?8RFOHWBOkRA^bamja!MmqPBHztW&73SF%hYVj7cmE>#! z#G}62odnNn_<>9%YJpp~Dsaex^Yzo)QAo=9SHQ4`UKx=dHo<56eH8;<-oLms5E!@D`RQ#3K7QHz(>)k)=J#yH}L&gvUu^|=R&Qcr+ zqDG*7^@xWP(Ok3kPg4(bw*lQsBiOLq-58|aux?DQrDWIxFgE^yu8Tvs8`n^M^3~kV zdOxAb`wyba3;d`eCB;AE&!t-6hnOAUgw$@L+>%-?rxrfFv`?R=Jy?J=NaypDO9of06{d zU)*GeUvF+7`G3`0+&IT`Qsjyz$hxTh3~^#6TCW`Q@WffJX^#D?zml=$d?-kVmPr)v z5iI)h;4-46oyeoj|N0BVJTz$8`)I}Hqj1EQQ?-%!zKrBBQ8V6sjfCvh@PKov=rRdF z<3n};Gaj(*SXkfE2X&{Tkg~jiy>S}Cpcr*w&K>7BiXS}Q00|B`=~(Z~h_oc}SQ&EV zB{^C&0*lDuA>dE7jOY%ocmQIwwlvYJr|PD+^_HdlcF|xHi5dPy_2v^e=d{iFQ#Nkh zYd@=lBuOY+$H zf_WS=Lpdv%CsRaMD=ZVbfb<#I&i z-r(2|vE}`ccM9lPHE3b3<_y23Hda3W{3N>c9}@pu3O4@y6LR2rqJ@gua&p)khI_sx z)@pd-Hb01Uqu~oUfn&RLt8iBvywn(TdK9tQ+;laI4>|B#=Z9k#@6Pwg%JL>PdZ*NA z;HExo{P!z(q9`RAN|Gl#Y#t+Eet(~d-u7W@n8}gr_CZ*iO{FGQmM-0LQvhcE#o#?F zG*s|E6*T=f0?1ZN%5NPQhui&0x&OjvSwgAsmd7bMg2&?&Y?!TRN|fDmV3KobN;mrt)q5EOw*1Il$KT-cHeoq?;NY5?60UdO_6nED^Fp2&?K4oSTgty= zrAgSn0>Pgr%FFY%%qrf2eVEk!NKg(KT?b6^yuJb05wViH~hssxAwJ2ad{ah|n`b=89C#v7z+(^jYcAV`JW z+gVG5T_osGV}g_Fv7(6$tYINhH^9@;{cadDjziVmGlcs!c%Z}EDuAr*m>+X`@bG=T zN@}NM+Wi6}{us4LQF7A{>vH;y0ynyAvY@RgS3>^q(cudUz+sOqv3CDbECS#1f#(af zfK07=T{jVOc|_v0ww(c&k@~#&{*8YPucSD-o^Jpnt{)1x5a+N&TWH`qstfR+XkFM@ zg7?4oYn`a%Bs*>7{TD28s#Ba=cH<7p0`V3MD|;IaAgXttpx%trj@N$mgu)Y#I`Ds< z98*){KYWh=JbM@xA6nFGDwP&B?8K3fU&Fyv>t9p#zfg$M5gv6ds^y;qC*f4ukqjyGi&*HPdDe zt`}9J=lNiA{#syj)$oVAi|YLYnlz@BTcE!zkgd_bmgC*4y3Aal0-v*Sal}M0V@1ELeynC{|vDvynhdaXYBL}_oAnSSkwd~#{fQIaC zA*$m?;n8UKygL%kdjz^*SMV6?0H<5 zJHnl&mtHXCZccq1PJY^1gKt}mTs7uFUgd*FEagva7Dx8zBe4}bhw$e`P0!Cqac-1qXRxV zsc6&G*08=qhkD4CM=j^h^P5={-vLX;t(x>6rUySr)^R4GDvHzS()y5(++>n@A>L_8 z`?=8|@LU>`SpoQugC2i$)BX5rVNf5t@MSgKQ76RS^e465ah|3-d+2f}XpD(V zrfpX1u_T&lCBUR_H-CbskCV^To7S}{#5d|vKb^&?eCx13di&2O3lgUao`eaVp!n#Q zRy5uFmfp*Fn+RfPuJe%w?tBNZ@96p@z8OnDu=m|cI6l|-T^uBL*`pWNDZXSL%7@|>Z>g2HJb;i|ZY7d!K0yeS2Tzg_&qD^1C;+C`z)tCdZl9Nb4|`4W zSD~w2Kj&aiAILQ&JpnMDmvHZE7>Ln^9_S`t9oB?1Zw)Va^X+TxiO>bNxUXooiR@qf zF)kXxod)WPYVN^jhNnEj05esnQ`EpOY5_?u2Yi^ zB+J=b?~%l<*lGHKPtS;oC1zLxa;Q}Xig4+<9R-=$1C1gMZlYg7U>jRk#%=i?V5Hk{^xkMB6K95*$WE0!dS_@hu`Kv4H)Q*`2-{qOiIw z?q=|}ttDwAfNAs@7ypQv^m3QY5zm6I#s6)Yrnn-igGYrD!RFzl#h%zKxJM%b5?4paPn;5&Z6wI;2R{4?YvA>DH&+HOs3!pT z(J(QpBiuN|9I)dA=OVa21jDq&#(`cM3>?c010Z>$nHT^@Obf(xZtq*^#^K`msJDK} z6;AQ)O-C0|yKWy5q|CQ~FfnV|PvoAu{`41lbVg=qd3R>FkF;2Pe}d~8biZal${8{f zU%3z{M9ayn1X9FQfJqwTsj-dik9K_p)*Qu(UcEcD&2DVxB3Jyjq|swt@SKJC%x6=+ z!CyjjiNQ<0B`gUz|#90zO;A^Ne87DDsR7?K(h?onPI?0rxpC6t|ZR!&gs=D zQLNQlZTHS_AGwAPNM_HoiEYzStG7xp?R=M(q((T#7|o7HUR6JX7_yb90(3C^hQg?iI*qeOn^y9fbPHlBC`$BYCZL$6&zWYzOYYLH*_PO{{WB-< z18DO+VA=1}*p_p9AnZYt*JmYc+2;uh*Hgw0B74pEHv)HcyRxdM+sFG89v#&pd>Mh2 z{T3_d+dN*!)}n9RGfK7K`mKAQ-qVPo=x(GNmNTn}ViyiLh%KeH$1UQWzqCK5jy zr>@ayFtY~K?Nm967Kc`&(Ayvg4H zyV7&;Qqf1tWwp0bBHzLcZ3L1d zD+0iWt}|wp%2;}^Frn^d3=sResdVBt!Mjy)$3fHig(bOQ`1g&8Lqh2+KlywpIHH0E zFhkQ)BHz%`nHSc5}!)1PuS?e&N$K^>P ze*E%TaVWYKZ1(Pqx#7A%xx~6DNn%zWi_+(qe5_x*8!kQCqMtlr{*c*IIDwHQTX_s0 zwx*bRbFN;dCqKahd&<`#9g|mE4jZru!R)3ux%Wx_=M&>b|A(5zdC^RO3sWOAfD=3Z zLbw_ONK(0hFR?8>YL{iPcW6PQUvQf`%{)bD8zU!9{CwMg_oBN_l7Q9gwLXd>)g1Na*=0v@{>Yx%xQ>yl61n+jV^2 zrzg8~4A}E!_5>s%=rk~P&<%eW4o<%M_{+Of0}s3Oh=_TPUZdbwhVG)nyz6%Jd5zbG^XD#ngbF8kgb?JW;fqH z`qPR1%ugAnL}-jI3glU&0>CdF#2^vx@3RDA-n|5DUH+S5LCm-7Ww-{=+J`|2uI>+G z>+f2!)u+-U9b;8Tz5o6?Qobm|-q?0^2Oz`63)@3~VY$?HtH{Ofmsa_}Fjq3g#f)Ry z`QRL>eoMdw?>R}yuyD;U_O)fFb0PD6rod_q?;nTan{R=nx3+fWsp_87d#AlG=I6j1oN++8cREK& zP(0kL_H_EDi@V@{!xoqjmM?>5hOaTR*7io4%xZ7J7dr)Z-2bX2xa#p9dbDzkhDP2!EIIaH$d*R~)R>Kfpb_g%8iMF_hIQ)EUh z*AY4vj7Mip#}E4%-|ARHrQRHeHN$Ix&Se*MxEbBS;-vbUaBRHZN2s%gC%+FTiEqNM zkR6;4_o|!0;%g@$agDgO-ihL7N{J`Y4qzj6L(B3+McIUA;na#=h55t@fT>d$ z_v)uOJCeeL2e{jr0@}XGG%Z|Sq4b}tgQnbf30zj4DDw8)9ct2cDid7&aCmx>b+-_; zvetz*`aV#{=HKy5)J10<$u9Vb=q6$`u-hpxq47|b^nC&fGIcdSk~3b-Ow70mJCh#Z6|(cA z;{0(ivkBKs!~hc99pMY;JLsOqWX#R)K!-^OEEFw(o*!!(uMb7fudgL?k#V#575L`d zq~)_!bU#vi)g*)>d}2?J#_b;7yODC7dNBGo?mLOa&&Z z$JU_wBXa(D|J@|6oe=HRUB5hB{b!b^MW(9G4f?xRFk z@J|TiM+_jQW@t@8FI2DA5&!7bV1sTZv;&rwE~(*?)dS;%1{Oh>;Jl?JJ7dt{-aLF?EiTpZ{Z?;8V5q4y^J3RuX?MoWh*ABn-l)a^vHfV zYxN>8AiOU=vWR_uKE=W&@~C2ARXk&JN;6*6M=1z%Ppo{ZL!4kge-rk!W4D6YVeiwk zqVJ1o0~)q$7&rB-3i|8{rfg2i;8sc?JGY{b5klrEZV2~;Tb==X=uR=JmqAmOW%>6G z!BLO-S9hi3k!CBhz=TOtV1(Nq|ETLc-(5!A5q;fpFZc_r8`!e{F#=)M6!5Lo4>06c zzVvlXj^Gy@Qttoojc{o^MR4!7g+X`kmZRc!rzx68vIMh%>%gsga+R0v)-AiX*MnXF zBJpvHaOM=9ALmWso1R*=Pmir{qk3~MHnbJ@h9XgMb5?*Vu;-5qlXs4*Eir9}qwLBa z-qnsIFrG?JkH-tYXEmoI0UYnW<(O@#b-mYjFA&RQe%T*Le!t6#PSwM zKOjsXU3XvY$9JDIAHMtW3=#z>GTsS^$_tWLKEh6!X@vQOx+k?F#z>hvBk;f_@R@$SAj z^&oxp-ok&sMIjYL0Sqz`{p)Obb@~-Eio8e2~hjVFIIZZlRZ$`oh~_YoY;Kn@fEl3Oo}?e|=ceAjka=4gbJqm;87iW9gK+J~ywI(aNG_B|RHX`zEhnL{+y%Sbi z^$$5vr44KxYwnL8+X9*2d4b@?DqgUai-?tA1IK!pKYTw5NgFN+C2TL&< zfHcOEV@V3{39?!q>s-FPG`c`|`KvS5qkMJgRkF5zVCkB2j`$ zNzw~s&5c#}B@q~QF2pNXkMO}Mv(^KY1%Axk?8%ZbvP_RSa>&-G+tVRj-^R=wWn0QV zyaVife=SU17grud@V&slx0b*nvhH0&WEsG{Jotfe!J9lgjbARfX=^S4Og|J3Bzf~Q z@4XiwqrYri2RQVfyZI8nZVJbjW%~+Ir<*0$8a2O@m@mXCx?k8IF+k8ZA@?uy?hxaa z$5T$mhg*X^M0}6AORgm;PtyTV>fx1=^q~tHb0;5+n^6ZHPN<7fXIDJ2R&~kzf*FzD zRutR@>MsJjTZ?;3trqrb>R*6mKHQEKxETsC3T_pl5`&=xghtf(MCM=Ooydo+_}_&$ z=edmT?s6>EIpZ)w`;M!*NcQAgAL;rL%kFyn^3+KRu;(TMK(}A5fQ?U%hP~N6jdBL>#5+vU*i4xI+8ROQs#%?f<V(Vm^HeW-AEAu)@RgfScXlLpm7y2$w)sfDS8dM9__ciaG!%zuP{aYKrpB7?b=*td{h7a>;mJg{_VM=I+=_yC01J zQ@=7bX>|EQ%7$E1+|gAL9mv|4#~dYec;dsPNR{+p7rCF>Hb*5@BIxQ34Z%!YHOs0iD=TrYmZRDIF*KlRI zQh<9C_UqS$C&N+A@?pD|z$QBcA6C2@vQ4(t$Ha7{uY>$mh`Nya7rN5B#kEW3-kDn5 ztpw143vpfe1-Z{FGM4{L6|<|fvo*zXP16nqmTSI$a5hb5i3s$H?qnb4D{?A7^CU397wzU6i#v{D!R_zniM4mZ{~$kTX)gR8`?5non|w z@I3Q5^EMS){EqLo^%IlPB5T(F)()zEHs8sfg)=tw%Jm+$Y}VKJcRVk%CS>`U=;|W~z2wT&3B6!LeD^Jvwg?6gU2`g?AsP_lX8A6JZNm z>~mTP|NGj+{%_aXFSg+$id(h}b-|djx3U-2Hbwk5Ls>Utm&Fy=JAYqF`l@|sjUh`g zUwiLl{eSeZ$*#Tab}8Y1n*C_+F1)I$U*Xv)3*xSN=PJ7Ut2>=Prf~|J%v^}=^jJ(v zjqjFV)K(n@1Q^#WOq}yEe)swymbSWcq_vJl<_i0_^LK(%#cNfD_FMt0Y!0zo_5aX8 zV9ABLN zxQ?lHz{$D2-Ic1Z{od8~>8a$%9#7d{c^n*^awn5_uSo45=)fUC&n5NBr*{t~2Ak2m z#iORj_?=HJA|_-!&)8c{?N{)29P&QA$7TQhhf`7fW#`DpJ2qeXZ&?4E09#!X^UU9W z!Qrg`smgk*CwQe%Hi*|xg3_zPQ4x-LPL~TFx^{n8w5;0Rs4pth z7=tMmISovN236TkW{;wNK2vngKbX0}+05Z*Z`1AMf7!bGd8C!l2Up4^kh`C-ck5l- z*O!Cx;%!{?EMFXWl~ko%vQ~~buROT(G&|lxUEeRS&gEo&?ePrfA(&rk*np6zH9^rt zw-_ss^PyabcPGMrX}QeDVuGru>KCd}+F;8zD-ZcD#Mw)ASyAF1vhwe!jgfoDa2vG! zb}TuM$B_`HaQDhY*YkImO*9`hzviu7QaX8g(FUEHUwr%T=O z=yzW2NS`iNF5#d*5*+mKtN$;^iTfh|2GBEC=Vo%+>Pm6&e=(7UWqxr+mD?{-R0ZqJ zt>_7o^0N((s`A&Y+;P%nmFeH^4OfaPMx8aJim~sNB`P+SlrJ&tEYy!O?d&ZAbO1kZkA|v3)>3{ZT&&4*?%dTe(@1rq@61z{}>*LE+(EDxxqn<2iNq@ z>!dVLv(lR+QenzZ<%c#d)jSEjk}EH)eo^(g-cey*@&Lp=Zn%#2;Ol$aT1KXF_*c&x zr|3{($k%eov+XtVAOUMr2QGE@*n}@>eca1An-Mum{S+;~?o5;3f6cZRR6c(Y^PI9~ ziCs3BaP5@5(qa2dwKARK`kx_h6*pIctX;`coL$Vb*YAvoLOG&?FBB@S+Z}4CzkPJ_ z`;!j4D92NK!wUwA6!2F%QpAeqaTW;wlsm(~@0ABn+IPuJ^;L#QvupWV`5Olw$Mh&7 z4?TDyq7D30<#hkm1r0YERSE@L(nWJdMZ4yXExmihFv$bvOlZc)QXMazx%{XE3D@~|T1t69C&F4`^p;c1m1p@EcKFk- zHovDT`W4oy-=F@WcWQ{I^cS@`=i*19W12yyJWD`cR<)-gdat+7Rm>@~K2tZuaT*VS zZ)}c)Cj&k&*lQEJZwE3juRhtyeE;cg>Yz}U?0Ks7Ly@O$x7?w+cl)dd7v$rf-HxfU zqcO(o&-4{Wp1QngEIiK{p(xOK7RmkiN?NsVZT7CXwHT*%evFHb<;$g-A5_NJu4s%R z1|uvmN_hXsdPUKmYu$b!;Oq})$Fn)28vosV)2TR=sVJRB#9vyrcNX6l!2R7a@qCu} z_-@c8Wrfo(mvf!kokXurE8LWVUwECeRQSWU>n1n*o>S)M#dNg$vy5W2Z-n2YOmKG@2c$^9Q^*;6-HfpQngGe zug~Z&`*joDbtTL7pvI}{@`6#$&#R&UYbx#h?$az{qbI9>+a>tZVHk$h#w%F#e6QF@ z+pazCxQq4U4})0&1oT=WQd53b_y+J+Xih-OfJo)oPFKMD}D1_lg58<>~{c9?zUb!^@E(S!jL?hEOF_P&Q*Q) zoc@=m9BD)Ll9vQc^27(V(r6${-Lp^Yvfe20hJCI)%zLenmvynd{b?KMo9={hYVI;8 zK3By5EW(d{L8SJ2|I-4MkxR0=N;;CgcH$<#+W!pkoby{tKPgrCfZK*=_D1ZD>uL@R zlq+*C{-^up8!@TbRXOt_nxfIJ=WSnVsirM&HuL2hIi)V|T+m-Ge(>++HBr=jPWgg{ zUhIO9_wE-loSCQQRQxLL?B#r8a9o9Ia`Cm1Vc3rhF2=%zA1;sEbg><0uP_k=GdgOvENIwpE?$@+8(yl8IzC(qQop9-6As=kYGrdL!c(oyyr5usE?Ns%@#?=z;7}Ldnv8<=nt#rEv6^-# ze`AtiwnnFUK{tIb{O1(AEz^|fmFn_1T}Aa$)SIf%Q};Ga=H9$G+C|<~(X;o;Bc}Av z3Y7J5J9++IN;z9Ubiw(r{ngrFRbl;~!nXx)?3#3?ei-{Iuh*g4s3<)3QSn&7)^N`@ zk#r38P=#^HUxyobl=JwY9`tm`_0kc*3RKT_4IOv6PcU-rhKib~pWm$7dz_YCpIgEE z2$AY5Td~@z#+PX;^#8l!D z9tuXDUll!V{@D6CjKS{JyVl2GCk6nSZgBuwU6#%dkAv3Q#luwe4!$s$m03Z0AWD?#*{_%m<0w zAsfD%r))ndD$LlabYMP@bYFO`*qjobTqOJH8SJc|jBCxtm;>L*K(PD&VehQln()Ft zE-|{BQDcl6NRALe8EkZoMhPhyA>APzo74ztgaIQ3#h_DSfOJdzD3Ov91jTrK{&?QS z^9Ig8=Q`JY&i(z|pOd$*O(AX1Gh!y-2cO>i62K?)znm`VCKYx3nTZrDnbk#SSk zmgH5V^+-wI8!$glaZzb*KKoH_#8{3akkHqzH@zZD64z@YFg%xo%bfR4SzuH*6(nEv zoaG?7xctL15qRonc*=5YOkNz;NUCpY;)M&d`_c}*r@GazU1RIUj}<%yJ~|gBNIO(~ zLjdt6hK5)eqS5I_iroqOaGEj&Y$9d46I=KhH;p-p{H2$-)=0CeT~}u&H5y>p!uO{$ zM@%Db;)VXRBqK%}98?fB*2ywZsQN37iLd$-Gs$=Mzt>kj2xF|I?gJ*W9c=`oE&2MY z@mp2A`8PL4b1)7tuuc_K0IWV_=Ck-5CST@o5t!B#mq~y5X!Avl+`MC}dkFiTQgd62 z$p8>Rkh|cZ7P}XAs%MtWO#6>~csNoDPZRCDoy$@jqvSZs5YP2^N<3#js`r~vVuKm0 zu2HkNY1@nY29?_TIt0-rMNCeL0dcN0ZVi}n4$A+9V|m5U$pfS?4nP#m4_fgH`1Y%7 zy)3vEa^Vs^bkOpaQV|p|>({#_C=b_uz`HCvJ%hi;NZO=;swWU?bZlo*lt=ZyxpWBL z1$`XW*$89`L0aI zbeZYE8CNExL^YFc(3(2Tm-zPM8=W{uLu9h99A6tjaG}yP8W8wO^I(&&f_Q@M<(YS6 z{G3|k(&$iN%4Id<3R{^g8K!%{b)ETd=kn5-Vfs4tB2GO7kl?c;NmbUw2c6Z&q;ij^ zoW1y#?e{rEYD$-_>}wA}6Yj*7w2f&Pd90G_KEh<%5Zv(Oscs#@Sfl=#3KMlxx|gFt zQc^BVSLnkIz|`_h$OhVhK!D0$vya@|5^%%$(XK0H7HR`HGK8ii#>4|;=Ey%BZV}D} z`u=ROKVpBoiJohnly@<$C4_}P`|U_1JIj4>YmM5VPPv>c3W>^mul>d1dnvXCybDk2 z&%UeKtu4Tb6h|Lel>UcNSx%VhC-s}!{#l0qMswh$z2kO7_Q<1iI{pYp@J91LKqm1t z3u4EEzNe#DfUbL_?q3^yQC$cB9Zc1l( zk&=a2T)ak!#YuF0GdBnfj0-hd1bwJz1MlQxLSn3b6}o!#TUkmaazJw4wS#CA#0jiguxeEp1q zg_rtHxwG$}n#T6TaLwxyw==|dnZVLJAIN6ECC2{!W1RZ;%4d2uh*{q>!R%XUeHKV3 zflseASBGlNS&f~%@}t=>+~T>KB?)Tt zzS8@1)ed=?)N$=x-728pA>lu&^qx01_+<;TFynrX0;nN-e9i0hk#pJ$Qo-(30;qyq zhd=}|7ZNxH$^LD)INdGoJTV%YwO>H(q8W7cl#kIC25?E7F(}C_(_YJ=vpFr(U%+K) zJUL78!(AFep70rpsl6S1%%=-`EnMhJV`mxaYA1wn)4?X2wzaYmc(QJ@(a-15u?7{s zPu~h4#tHwLP){|(9|#!dZwbVBPR86+cy|(pCjs}|VHNh?5 zmFIUuyM>VpF7Xa0c%yMuCmG}21sX@5td9=^?e2T9^7z*&^K}8&KsZ-|?)qHt9to9#0yN^}W+X)fgP3J3ARkE`Ys@k8P$Cp0B#d`^}tS{`9 zOJbdb?dhZj9J#>-k?zfFhU zDm_(#51dgV#4S!G_rNNfCV__XsHO6H?TW56mUKb#`w{Q+ za*TB>%A|~Y{%sIn0Bgi?R;3K|``NGOdU?bX(tF=( zWU!ml6vPJ`_f2i%Wqf*soFvXa_#(Mk3*^?w7*W)9Tmn-#KcW5*vZ2s<2ce&AJ5wUxS&q?Y&&bS*TAL`G?~T7!tMQhPgD;Ab?VD$(<>76sql$@fJw0=tq0M3^qjDW+h1SA8GbBf!EGeWI zzlRg431SJdnw{=-mFp;S`^_TVlT>e&SN_xwWM`@4}n~R7#RE zh}S%tmZ!WXQMnn7B2j&Q0>_!Ynir|_&S5#P-W$Kyo)D!TU|8hl{mKtktB?(w0Utf6 z7G=rr*jM@FSfg2*0v5tRG{_*Tmep4GHp@251UMjWQ&X@bN`ztJ1Eyik(h--h`zr56 zfiBnl6r%c_x^Zfom%Pm-h0tbOmI64uQT&wIX*(f;tCZ}q<~&>EuGY5tlJm};U z8q6aIKE#zjnw1Wkn^nz`sF1f-EV55zL&okkW;2~k5zItwSoay`Y6J_`-*Bsn~7~B%Lc5{n68$&96 z#zM8cfk9-L5N639ua}!wrBBK=zcGE1uVr#JmUJLkYN7F7lf+ogUjaN&QC3V>=hRx= zlfO%J($=m*EKT7i*;82~65Y+?2OGl^T+vIOYoKhg^P?DIHp{qEU#W?)2aZX(IejL^ zuM%X!e5(Mv?B#gHnKt^weiX4@-Ya)q5?gxZ|CI8$UD)Y$dO@KwJyLWUk8-fNrHBs-?(*uBSrcB&;Is zk3#Qr`4mWwy8>%ouPY`wQwC=G^c{mBL<}gC@JVCq4t21L$xzQu>~O;)xOP08a+u~= zmiu{Gs9#_5Hw}fPE0X8ld!?~yO>)cWzAmJf6`T++j%))LyIFU8%Juz`C;AAD;&M+p zx$J3HHy1aU0dyO@t~VJ@cZAWNPW)wwVKeC7vSM?8P~KVOAT=k!^py(cp0a5{cb&K& z4$%YC?1;89;gm?|%rEY#x6#(v)<%eG(0xjs3N}*GXMM4m9gNzM!l8)XY3aXerU;Jf z8adS)-7PJb94T}YnOLn_9=ou|-_eAL8dvuzPw=!h9w`bL`6af*tjgOZjnv3f9B z`9-p)32QNYs(fr?qDnh5@JaH5Tf5U)Aql(0&v-{y%bku0&35+(W!RY^7 z#w_foD?O`l$i<>Ux8gvd4Gvi~wam60P~&w4)9q5uA0sN9Fu>hh);43#GctN^3(N)j z^cnM!vb9lO&Lf|P^2E@w75>_8QK0J2M>_Clf<(J~CCCRQ>lCItA)Xs>*du=9mO~4? zGBBXIw0MaG~k5FOOlzeIEkQ>jLz;*P!?m?3P%%}K|Fz)64d0Ce6< zwV(i|wvVWT{!Q_TNh3(aogEBjc~(G9 zqBQse({cgnE-OkQP$>7sRA%0p?V+duxB7cZ?y%QD^8yxt*B|rg<6CQk_H;>Cou)Ic zljp*6;@BmtF2R=wJw#G)!DD1jW)rXh&E}*o)c&U=%Xem~K4653^Qx10d*++){7rAf zLG^D&f2`YyVrJQ+sWx~zZ35KQ27!t^MNo&wO-8=bm#~=lp&HejYF?*x{zaYq-O?ZZ zY_BLk~(np8vHwXJKaz7w&osJ@U@xsPTA$x$vdhE3bCQB8xR zmBD4)a+9wD>uw71!zySoWkFJ$7Seu#%}(Gx+pi6+ca&{nWtwvBqA8`fSEr_G@mkLH zWx@B`Va9%HzcpQ{GbFQ&{hzTM`DFCbGn6t1Xfr z{nFgFQB$eKsH(d(5BW*5sBOrSzqMrt;xP--ijPr;3OIvl6bF6&rhU1a-|DSL*`VeO zg5#ZhNmW1(Qi-8h6SGT+QVq8E>I&Z1+7iu{jA;MZwJc)+sWdrk=YWWYP70p!>nrlQ z%)#9Ad3*L;a*nQh^ypMql0bqaw!;4O(jn85UO@jB1fz3`(J(4rsyphYwcJ&&D@&(D zSafXrm|Wo%Uy`$83-Oz=Iu9xsbXnxSv{z(?dQB5uaQE#+XCOfxXC?{^O%}+VQ_Z|e z5D@+b98nogLVbHq4_Xg&%D8~s63$~1^W?xGF|g0CgH6YhhP6im{B;>G$O}ZBF(MY6 zsK93TsiU6Kx7hZW61aO3*FALu{AMr3dHi#`jKE6x><6|ZC3S0A!?F06V48OXe9J1L zk%r@ssvA~H!Y>90#=A-`s|_WM1V@rVAzm_ulpU>8d$Sm>0H@jeQ4?mecKk z9gE$YQ(MeJtkIP$hJrDd3h@E*PG6J|`-rD9u*8@zUgmhmj-|`zX!>TWhFmphOu}r+ z)IVjG?`!*Q0(SF93S9p3k^FP7LY=Q>f$dN=mM%z9r(=>p?~)r6oIk%3zS0f5o%HqVq>D#pKWD1yuWNGEn3`B zXhiBS&gp&xdi*-!PPtfUA=#kwSh>AilD3DHMbk19z*N6Ic6G-6u|l0~`1?D}-yFLV z`6Vls0L#L+ihHsofluwGH+%qL|FMo3ZZjJ(-p|YvTTYvbeWAh1OQo0DP#;cyuUS!z z;Z~tYL|_o}m)}A(Z`>)Os=jR{AKt%$B4xkIGVg~u7=^bAR}_<)TP(1jFo-yg@u@wS zW4lP;A3KvH4RYm#mBu;BxxxxanXY*K%n$7>-TJ}v3fVJ6GMn19A_lyg_A@;k;2=M^8zi*?V%C-4TOMAYt=j)fWt`+%oJ3 z-v%-;;R9Sc1{Rp1A-9hv!auXF-ji5tZZlHdhaBYZUxM(M5GE*gA@n;l10Cf;dX`n^qLWya{w-a z)St_lcI_vnG5GF1$4EO0(&`gq4R9wJcOed@LaOGfUik>StdX0|}wz%5VKKIo{!M=@3KK*S30gtf)utTudB{ zip{d+A|o&9%6XdQ!+z~3b3%jX1hGk*o8@z#2WGW%(u-4b?rH>Lay}AbzWw65FeSMs z_S~K=Ayo@VpU*s1szcNm__OFc*Y}3Rf=Rg8GZUzbWXxV1FAD?{_h;toLizIiEdu#i zg9?P9cl^h}176+!u@t!&k-Io1!KtZ$NDrOd2tHS_%J;Icc&XQ+2#@TP0Cd3=nLX$PnUeX7#j;u>Lvs&apZd)u-75`}%HbLWm%p)IsJYogXmDZW#QQoo#%xR_mM_N; z&Wb8j0DsFGQBtKKcR8J}9ZKz-4=lW;#oj4bgcq!*3c_pMu{Quh6BAg%TcH~d>&2G8 zHP9o}jyB*?Zmt7=sN5`Ho8i@47$alw2g+|Gm|IXsd}^^uZJi)*HJ1^=aUwlNL52s! zYuev7U{MHIWc^URIY;_RACFLI2GwVkh5d)>Vtn;civu5<7#98{ZMX@J`wW6!S(~8? zY5D4lRP=@p!oLdTzBtu9$VyFwQ;(=$oQ@$Dho`tqG>kjfQouPEOe^ne<8Mb>J7$ik zucRAD-^o%ll9khwg1;fz&x#@E1#*rEtW3p+f9MS>!d!1r|1@>UDNT4x!=Xm$5A0fI zyv;h;t3-u*t+8R@0~A9v1@zxrtxK_EtX)hM6(UjC*X>Lv&xxHE@W8NJrp-NF0Da6R}vq(vT#N$LJwXGIfN1R@N zR9|3Lsg7rLU0aD_puTvnZD{H1O$Kz!U|}}>gBQnYD?l=2IUya)_nL+2ERmc@^_$3v zfN+jD)rJ&nLOCm~@(*>Vbe*Pa4R{z#oxOS%R>R~0es(`CAwL#7Sb!QN z{oz>{uV)QyVY(UC>`Sy6&4c;KKT6G-;7OVmr}dZv!oub1pa#8->g~36nq)quZU@a+ zW~Q*ip!eo#)JqWNqQ{t3Fvh;g$l~gOV$ew#_kYa>(zF_R=OoDt2Q(7~f>V-W*2$1J zkP0%fe*4TZUNsccRGNxQ3@j)9%Cb7oM{JHEo*V#F`icpg`Xf5O#B*H%`XbKhkVILn zGMZ{IPQ~xQ23(Cb6EK;SAY$p#6E-`oyOPD&=;$5W6N9pf`YMuS0)2dxbJ7F%_Cp!jH!;Pkkzn zRiXE^3H-BlaoHdBpmJxe9*Qak<;hRb+GKaD_M&1^OsdXONWrw2c3B}H1WfWCm#G#n z!ieC|I?rPB>e%=Va4ihe<1p!0u;?F|O=&e=Zu`<}g-CJc@k#(gzX0!dA#3V$Y{Lye1JZ*>ps(ep-A~K zsIr@++LhMk`P(sU%n6TwiElhZ|Ldn3Cxhmye}?g~#hzvS_rBsfzC@h&x?@F+JJVHv zS+9tFzEz9l{9%~7-bsm^XFmd=w92ZAZ$6BnWpH-Tno%aBp`UIh6CRk!#K%W)aR;#JFk3v$02Ejf!a* z7ywD4GxQ(&bkHXFP*#L01?2k9kT17?E?YX{yEG zJ}LZ;@ZmP(Gp!@jyl&gGE9>w6j+OP#kMv%=n|oR>2xV4J=qTd+?V zI>Z{J{Qg;k57|xUG*a4;WW^+IPJ5?NUB2j0SwP-AexcRiyJXJoq;TcXY<{uQ8zXZY z4H3FRk=S+1UTJWe7DmWf4fjzYaz9JNkMiOXEwjN<67q^2wW!9XvrTy&f}V=S7H}I{ zV}3U1iS@`aGa7cDO+7UAII1pkja_sH^O_C(DDL@$ zdF6r()V$Ej3vH5s7rR3Mlq2ZJWmzhz=UmapIpvMkEdGe<)VIKe(`{jzvSZm zk<}1jS7ncxNLSEeI?GE`Mto8mVSg1>P|D=UXg53Va=g!Xt_>I{d&W{j_Tq~S7R1T3 zQB!xMTb^zv*{Y~Q9i}|BVgtzFsgE-SPBB7QeWVfP9P#I$iqP-4u{GrSeuOLP?@a#x z0!$kfdb?qxKXdvdCFEvp*5}qf<|J)2vll>>LaqYByu50Q`t*R@(z|&3xFbxlRs5Qn zIoXK|eoJkV1p@d6Zi-;MgqW-DvanejPTCtp?(>f~rpb+2(~==cZhr}) zi*=iW99YMSyvmAvhmun?I6Lfo2?$dY=!r(PBhvvU@08xte&!GG$vbZ?m@UG*)`X0k zG~#P4Yk(auG0JDL5mw~i`CzSLx4x@uzFvg+(*B@i>eJJrVeRtO>!oRL;Gr>k<&|~0 z83Zs6dyArN$0|d z*Alau+MkmeJ;&8(X77jf>i;;sHS&dp{(7pbNoSoB9+4uNT`^i@3(x8MnJ<-M{T&9G zek0K1WQl8}jb#Z;jpt3E15P%EHJQzsMze47aMTFAbr1co@8zMgS6@^NWBp#2usv#aeyEezct|f> zFx30MD$+aOUkY_(B%~ssUBQJ{+U4P@_f!;9>^(vo5Y2VlL0a?`t2@(d4k9%72=wmo zyY$`c284YN)di#DgaNxK^APh+s~5HC9hFt}x{ZY{>W!Ma)jMB7>(UAbrfC!9kSiif z#4@U1eVNAYPs-!9K{=(Ui=ScPiA*Bq=O1|{tI8kB8iw2PoV!#l7X1<3oqjl7JB-ub zw@`Je72noxs6DK#+}5%0hh;K-Uaw7JY+C!rsw}>5WUdtAO{93=YIYZPk13DcsPU(9Q~$K1$R~!{>Nl}+ zi%^<0^tOMm`<(GOt8NX*am2UeI&(wE`QxikV3z6a`$DG53Wc86((!-whZmgq#0Op_ z3sPI-zX!jq)D0jSZP>5JysCnh23|LQkWTNLp>9$SXQ6jSbf?f5noN2=opZW`f@Eie zzghVqWmJrm%v_|hxyE0jfzb5Hs*%O`Voo)AHBJYU|%&+@gfoWrTsxE(DXLR z$?^q+Osxz)bWK0GEGZneK93u%BU0#z51qOwgXkr!VHPNIY3A>HzrDq}$JtCsAmT*J^*Q(Rhf+a4 zvnJl0A<`N)?zMOf>4p2b%WN-crBTKAzbiC<2ho{DBBreL1$!8lnVhoS4bi2%k_Wj8 zkOtEXyNtJj=ZXSp>9@kF(TEKNhacj<^@WKxOO=sK;xA!*u)62=^Lqt`z_eLj71EAHjKAkwu+0Z_6-sg|^EEzhr)1+& z5b&CQCi}{6Lq4>#ib#9Cz`Mn`r1jDyPVMo+1h;jAhSws5V}oz*5y~04_-kj{$LK9p99~p+EsA?KJPdwT;Zjzh#Igdh`3A zAME|gWlHUJ02m5c*Y=SLOZbu%D^#AE22tBU)&{J_MES;potzd~Kved{<#1}e5)~7* zE8?#ul9kr7F8vM9td3v~+2B7Qg{tUUQh-0D$yYc)$y5H1*2$BfE$P1DEA!E$=|vgO zt`zdzrF;Pb*KJ)tjZ57OLt(VX=l#(DhjHguCU>GG-#b@yT%$dBQQMakHA4}?n!K8A4GI4=`(MMSzxj&azy3{ZMy zN}&nxZo5r^>T^MLU!`0-X)k+U*JftnaFynCVp*#^@0!7-BOmt+I^9uIvt*h)hrOV% z6(Uq4y;N>s?mK!5sCgTti~AWsg!_9#Y3CSGrFQ0LmLI(_rUNLkq6w35<<@-_vWqEBZ1c{Flc&NKKs&w$&4JH$GJzh zxt1HhuRGvGWv}@<=>uRN&SsS3tVcMO9i65(n^&vQz)J1eW?hlcD8-%5OIA#-Qh3*M>2gk@E>ZdR|+X91JcA)YHyeo-s$!v=65Q zPRn$8-`#W3Ap98;!xPmAgYaf^Y7MKmC?A=~6)6jE0JDQyquCd^Vq6M^jc+zQ+eo+c znRX3Q#5!tAceHDrOE`Rb&1G(>UHf(kP4m&z*leaOQc%J{^Svwxz5GlD!i+bkjG^zh zrGMK^cTims)~K@x$az6Ml<3fhW1oZb+Evms`ldby`Pw6n;{qmGnH?Yr=LZef4nfpW znrpLu8eHuZs>il+?-R2WB(`kN%!^rIY#w?taHU5Vd_j3v{5Ra~O71G>cSH6yHm$@B z*#FX?PT8j4yr>nkCG+)m`azt#b0W6n=StF{G|O~aCZre`#%S>aw*IA8t;AwFDiFpR zuhh!4wl#P_KkAqwH5ZUAdHx6&r?mjqy0G?6l;X@x(!RZ!@spXwF=_l?F`inQx%n4L!6^J?O;*;C<5dm-?}08W zf`ondFk?GYt5%_Z(c6;U#`!Z0X0S6besD%$J8QR@E^;%N1jW67oWZDo3pm|W*hDR< zL-oXVisnGvac+4PBbB~gute}vFa;pR{g?THHx_gf17UZ>?nauc^%te8b=JCbuN18@ zY<1Mbui-UiEpBDzz`1BMGqo+5ASdRoggA(Dcw>g0A$t7Reowg-9y-j_o zpCHZe=5&tC`g+^8bCze`-`HN_!0(GGV0Fh+I+a>vg@4YChywBy4lJo5&3fLi$BAEM%Q2c*Ms^bU1;{g z!4VE5p4Yhs~WtuzrlzIJ1w!pqw*IaMsCA6pz!yZ!d_z#>m zMeOk&1)b8GSkI7oTAfXPD;*p*&D7mXPq!1JUYSvU{4O;`2PRc{a=(aFPuXQJ$I3Ur=qElb zQcu|XMfAKP%?Z0Z9i(T-q7Z7Y&B-^_P43nZ1zpO@8!my z&EE*T-vj=Uhe(Forv(S7ZiB$`#-g|EgPLYoA)&l$?n*k7D_cboE|x5Hg#(JDe2Lsg zooCBHA3zn9u>1OIuk=*gLk~4IXf~%mu0EmG__qODVPSrJdat@D9%&Z(kh2SbX?Sq& zc9GHR4+YkEYm@FcZ|-V_{C2FOqoEe&xF)H-Y{MGFe&ngmU6^v>1*59kFu*3&mRXwv zLo0>6sk2j8+efA>Dp7|(!|TM8L4rP|ZeAlzY~Fh`#mSew=d6Xonn+9Q0&o9ceN}n^ z7Fc55_{-$QII;Uv^!KXT5Z{9RVjVwHfoS%gYJ~sI2MI3T=|}xz_C}o9uCyE{6>t>v zdAk{@fVzrsqtR(Wn)wOClcAZXvMZIXW;dnj>#_nI_S-*pjG@Ld9&K&%P`o<=X-u2k5*u&n z7yj&zWCsHTMFT^UEmZ}3C`SFBv_4HVow^9i8k>p)!5_;FRw(2XlI^JHazw)%ww|;g zBC?UpNUzo(rBk`b7c9Q?WH2+eCzVx(Y27xT*mCJ5n|&D z+*>8{iq7v><>alHD=oeFBAS#stJkJ92V#}T>`XW{2#6+j<_O}3T>G-y;*VKN9~qL| zBu#}HMQSD)DFFgg5XOmvWt7Tv6H~rV?d%kMK^gjLr5*GNb?)2QY;yZ5dd1q(le>e) z>b3h}Wf3xH=pJHn-0J0K8yeA*#r|Knx|#4<^?ih0^#tpZ5@yKM8OTQNl_=k7 z^olh2Zp8$f95}k}mM$EMlXRNp8}}aC?w@~{3UY1{MG_jYoj#&zx@DUm`Crq!{+qm6 z>`+0+&*0ET@6F*?nR$|Bo$Mm?mfG=|HBD?UiW4sq?TOo{6;!xx!B4e1isou@*e4(o z(JJ*v*_gG9la5;TVu9cIv11`H_Y1|pTPk^z^G-cifyAEO$~x`Ow`CniCOdrZBx!(pVn$|B zo2{bZ`Y9=?S}M3?0NW>hneqRcwK&^!-@G`du=RkM`2JBlV1rQj$w!Vt+d>S&lM@qY zadh?t2Azfe$CUoJfwv}jXKNm9q7jHm6sP7T=C<@ zb@?=`3Y$K;Vn^w*GiW4ju5m03L)Cff8KuGx= zfA(oIV*LA3^W5l;(YS8TlkuVe2aD?*4}|$VGVY_+%uqGn@2V=i6#tQrC&I( zs$PF5ujby8I4sczxl9)kPqp~i6#{SSJ;~0vo+u~feRI$@ashzaccaM%WUBW{d*anr z;+pT_?_$Y|b)T=VXrSrtChj){evH0We!U~)y3Svvy*nmcj#wRRU2tg`obVvM+_QgN(yw5}H8`PSFF*>-|mxoI)~x<>{WTfR0*|b1OEp%V~4T%d$sGPUH`_3 zJfmKBVEQ~ZC(lUzO9j`;)!)v<3-rp>n7gKyPkxTQHNao;oI5lhp<;3={I{CZ}CsBF01(e zNl=KH8TasBBAHpLF_ydc^&iYQ5)GNvC+pM@?4!uKzaPZ2*PBykLfT|LWxL*9<7Aj4 zVm5Tug!X>$B;)Qqn(*>3m!yifS$|rY&vn8B5YD3X{lY=>ce(fQ+O2zxBTp{U@9=0z~{;wh+8 zAb)R0u0J@1<*-Y{Ap3Qw;o2dNmf>@@>Kn>W!`dM*<5vEQ^kHVKp1N{0UOm92Jb)xt z{$~f65p`Jk=-}Hquh3o^C<$MdmT8|09B^dINY#`6vmus7nAZNYK8G#a=$vaj(PoM; zi(7}od)tlT_0>ev+V{_>X5%Aily|zqq>{RH>8VXGuAL<){JZVB{@326@zTg%+nis)AF2H(t6TEXR+QB170SoOF0)h?)VJQ?bN^@VK{AZ4<;-8mJnmh8|`l zgs{oi99&ph&0d(t3qx`{dH|ACu)mxqCtBH@G~i~f?6!&)v9Qh(gSUw?2lROQ9eGU! z=jDz};SU|8waGDVY|6XvO638+p|@kJ1~uND?IT}G@tgr=H}71+$(Z`z+?8Z^+pV&; zd<{gL+}7MjuTB0*ejQQGc7tFQx9)XRpfR*>BEpx{-LhcJ;?dFUxVW`3__S=IUd6XA zH`lm?-GW`|OMdK3!@HD3vk;;cf7V|Ha%7XKdOc*;I-iSXsq6Q|Nf&;|ElFq_EN#P^`uGk+MJg%iqkM={M&=( ze&P5Q8uU4@hd{*^t{xNiq_*U0=IMW_ll;BrD+%7&Sg13m8=%x!X z7A8EHcYP`bD=(9TvNTt6AD3{V>n(^`S-zD$1}2T$l7c(`+H&ysXv||zszTf1o`T(g zb4hFDX4)&5Oldo{-=-}*^0m!IA!SZ2+Mr)nsjcU(<8Pao z>?zpO^mp?}9$`*Qga^*qfF1E!zj~*ckj424sjRqObOIFH`O5J#jMmrra85e?Mt}q{NqgY`m=RTTD#0r?c{ z)LEOd>)fzhs$_-mfoF}C77(0@m5?;OMt^e*XPbt}&WKtC8fYq8Q~g3TikDkr*HWYl zm2EXkR|qgyJeBfFS*N7%A<5eRI!l8!sKUuoV6E`$>%@3%D-5_I!UBOB)U_b@xO3Lt zoxL{i%rj@MDmB`{#`~L4c&%rickwC>r-i`NCsUjd?emzJ|7$R`%XOK&-Na9As<>*N zd!3M+`8M`5W2}9{A?Dyq`kuNd+~EJ z-F)13k~CU_Hq7=T#>IR$()e=mz%$*BN!VhbS{#vQDI6KB5mf|O2i3QGk0`5Aqv{NB z4uZ#b1dN;vGCVS8g}VNxS?#!USm?6GRg*_#5wTww#K`vFH=4D5QKZJFR&>VA{W>x9 zb%FMc?5Pd~0>Xmc%=7$8*yhW0rC;eRzdh;Lz18NNzFzj_HQ6(;5!h&#`b$C}cj|uI zJ~mSLT#!#jw$I{Q!KT!b`H#3XhI|>keD<7)&;aZ^P;U00k~lU|6*UwR9Q*f_mC4A8p9Y}=*EbHRt?3BhgaoP7 z-2LSw@4BT($e3+9hsSrpEwyfecxZH%l8ENe)2$?J%4EfhITWvS5(~lSIP7$fAA~i2cXQLxj=J zoLIW)OgdGjG)z6Gny2m{%gQba% zk{SoyhL&bq>5*h9i4*oLR$Z7QU|EB)*4Yusx4Q8;kImUD%V^h7F#Lj{(Bc?k`6ggU#89JGPUCMToPn22e9?h#O2G2U%<(|EylcU#G zzu_1d2gk?OhgD=6)$MNErjv65C(PZ>L)m;UsNI)O{nX}%r!ivm&740#4@Rp3Z;RcU z3*c1l{*=7_*?sTv)X~;6&wUGkeXb@uQhbm40bP<7IjzDT-{1{{Q;15}m)MZX3fU+M z{}}%_54ULFxI)h>BELg}USB*`5O&~Rr@j!mQ9_9C@$FF5PUqW?IN!<+Ye9z(r(oDY z9_yKz3IOPLiUv$)RnkyAs)IXS5{0_~(9_G@$x{q5+W$PROv_O2N_J9yQ3-rTfR(>y z30xID{3i}};j3Q)s-2%cjsLAjRgps<*hw(?&>}c$Kev69SXaAA_VDUF*z2VgRgdWu z5;W>~4^56ZuvvSXC44370 zl+ZMK>2oe4a@tH(I*EK<|Gwm;1o@kyQOrx5JlMY7waFK|?0NH)XMRCSCG1~n!5RGa z+v9cP>o?kWv=1tZ!%%yChePSj-|B)V-ZoC&xmj4>b$77!jTtkH zf9!m@Dk&zblF$rUeKV0O5n0o`^r=1b>0|Yk%;ANC|DotQ1Cn6d?)zxwW2RPSTCSFn znVFUok*8_y%ABaU_uiW<%gmMG7PNAu;@;abx8epB#XSL05fKsi`hMSk?tj<0&biJx z9WkkiTVcIw8bzG~P{lc~dCP-?$x$A=vY725A^~wzre+Zv>&8G=d!9`D>*idrEMw|{ zoy~(66bMVe=ikr6rla6JQz$z-ynsbcO}O>&q*&C-`%|J{)t`d!07V z3>~tvM|$|)g4#!sIi0W_Rd-ehNgNwwi^eL}N^jHz$u0d3gl4ly%!Z|rvfEX?uN}X% zpW^1?6xNRK|Jmg7-{Rwdz3e?628UB=N=XpvacO7wb-Vh>H=(D|78YN@b|cm@H9Jjb zRO(;T@KHle2M&H9Lq zS+5D4PniDuGYJ*Ne8){F!I6t!BBiT$E$gz1!@8`(BK{v>ASNa=G!Yt3DM(F3BpTXO#bK%gR%!tI(RSpM@EacpaK&M2f=65cZ}6Y0voXb zt|x7K3!Y#vJ#`J6SLu#k5HqftoJ-SVs$ZH`%`03yByd$Vu%DLXxWTB>Y;3Go`4>)X z92V&Nuor^iSlK!Lx**4=E~8dVl3E$;{L2ibx>+(zt)YM5cz4oMYFiJW%O_iF9dOv0gVA#$DZ81DsFkyrkJtYL{LsclG_^?j zJ*F6kG#^H9M(bu<+6#;sqMmVcik{WD(>Lnw z!Fp6um{`L&uW&GXrkfNMzXm<%-#k~Dm?Wuimu}`p2xg%lySXSwj{7& z!fQzfb+BujZ6b7*4^ghT)N{G`k{LE^p!SM?FHc8Uy+ykYnM+?>ncG64&}>US&QrjE zJ*QnED|e;RX`{4bWjw~8KzScPeqcAB@6dueP)MLiWKt%zljKf9X-K=Wcy8DxJe`Z< z@?hYaX}!#EZDreDVWul8)vAgW>Z$Kj_|?{4ihEOY`&?tiHshpN_-?y1ae>$cJDP zaE0O4AdA<&khofI4BhK*EL?M5G}U`sv;an~!SYz9_z`|t6`(X{j@6GeWE-eId&~C7 zn!>y%Y9~IwQ&W=4Vs3=7LyFLPfa{}6yNxGHm@YZB&4|c>f9(@3QuD@!?V{rShBe=a zzr4147iCQ*14ur;CqlD+gcD}~YGz2jj(r{#1GPH3Q#|2Mu<6wkcGQTgl|Y|qsS2L% z`ENZ2P36|XSB-;NzLR?bN|3Ewk(kp77Ip;&3^i-rGP?p;cEnNTJuI?P^!Hmzu0ZAZ zgK{sVsolra9HXC<%D8J2H>=suOb3}g&I++TKMTAlPuw^kYFx&OjmM5kP z>U3z4zTZaOwNhOgU@jmL*W@bnRDJ%0Hab%l?W-zM+b4J0fgvAK+lh}?%VXb(zdmTo z8$rD!gy28)$53*r(2mdtJ3a|H#>ZLDII}ydrbGR04}e|uZEm2SSu#lqErthS4N4}0}Ndx3b64aIn&~?8 z_x&tBBP}hi@_@K!s-EpUO9YqUs=U$(tVdfq`sMi*%L?txq0naLiYEq)s2Lj6;_PHa7Za3n;Bf0Ivn?!|6baj!A<&znkV;9&zrd zNtXh}pUYuSAaB5Sy?iCD$McWMPz_z=H?SI-53sO*ueE5swqXUWRN>UV7bi~Jdgmjh z(J!Rpx;C3k-j4w8_{Y>OxVJ1lmGfw7f)0Mf)JRBTfR#7YDnr^TH;6w60>j_Ljhk_G z@V+stbNe+VPURF(OKb3w$kxW>?(Gb|d{?B}_BB)%-rvd!q?`4=!c2T5*l~6jO@Ft^ zRd3yPKcjwI+^EVpaX5{Y>tS$jRno!J2~%!eJc{r`@~BbKvs zS*{!GqYZulaE;&HY!+;!W;?qVD{EZO4mMwWG`K8jd3MC5_aKwO#9t{u*^Ib$s8%Dy zL~uJX{X8RQ%p@1?>M+;cooptN_o>hjh$lW>m5jNEtI5%)?fKaQz z$$w9mgeiSbMn?F>eYq5|P?@uip1q{7kL>JuUb8T-SdS1_Ufe1Yl@Y46TL_IHZT>Fh z$Zkew(Z;@Q40-GmlkDrXth;7V@f}LS)hi$8S^oA-qn&xV7veM;M+f_- z8oIyPhQX>=9N$10vn$n6z%~IeI~>2#$r<5vG!q-c0cin`%qYU#zp>wGi}xyu2&+`y ziP-eDL5b6RPKt46YpC}fe7fcqw?Vr{7MnTQswIa1>S~Cv4dJ^VKSfQoBIa(6H zYT&#FDCD7nUFYCapl66Br2sUcmf>ggrlF*Qn`6IvIgSz`YO21LGzu*};1!Ecer}>J z;?k?>O8C%Qhhl275E6dG6smpTd#CL7u-^U&fX@=#&;-T40a%w@xxe)tFw*AFB#9Dap)ZC404mf zGcDgCI3Ik%q%kVaoZMGeXvHT9lBN$TKIS%Ix)PZy?#Y5H=-F#EKz=^AFSfl>xsoS4 zaXFa5$FhIR88UBpQ{!GevD-o4pw2!pNx2@rx|xw=(nhNjjP!X%ATEKq4xJb)b99+C z&agm_oWtk%w!pd50aG&_?Pt}m7}}gH5*#>#itvc@Y=#z)(M=7z6nM6##!Eb;)`PSa z=sYInC)r?e?!^d%@N)0u1_#5b`{d5Q-6#g(wA?17F@aiZ@_kR}&{Ur&%zM3Fslo=S zQhBhI3l8wRo`+}H%q-sr;L>`*1wkOj8P@sKs=V!-&JWFgcg}({=3~tcpC8e+^7=w8 z0z+rNq;%FOaK?dVg6Bey@wscBcjv4@B&v zti9gzlCX6YVxIjZokO6VUH{B>m^&&=vdxg(bs+*(%A9S-fKhSOTJ)R*s8FbeU)QT` zE4p?rRDovxcnbKjuEKna&j=5Z46cl{C{AM@ni?3(m5n()U+Kh~QbwmguuU1-C$Ifk z2DuOeWA|2Hd+dL15YCPkd8mCUvfH_f=`UW0E`J zKb^hZdP63Mfm42}mwaThc-w0>>wTUdvGIFb;hVNdz%F5U0{KyEg<~Td#ac=NhjQx| zK9lpH_FQ!%?=vsUv0YqxCqs+S^*H|I6_SuFR9hdt4VANCz-fo`U{J} z5UoK(pQa|H4%#@B$hf;*Mi(VaHyQuIgst%Xnq~CIs~OnB5eUD}L4%u?BWMlAD1@vn zj#?O8Pq&i`Uy7bY(ra~Ar3eiy8=$uM%;n%c>*vQO0nPrASFVLP)QrT&CdSqVlt`javrSnrnZ=06J_YE~&5{Iw87_5W>B zESh$f{#(r~-@7IjFsQo3-KSHbQ_M0di-XqBz8>JYOpf=m+=c*reUX)wSuvzfJ4v{- zNE_UFq*aouVQp)Gx|9Lsm9X2X9OAsZLK)TwxJq_;rl5nat|ytx`!M<(jn@fsGIxX| zxD?Vy%v3gao&i@6NnZv5GG=B0Y&FYo3|SdJW_q14X7k)p0pzMvF1;{y=KFDR&Qv*f zPKxP*lqUX+aiQ@8ml=U)Wo}-I5(_~K@rC6ClFLK{(#554n1;iM1 z>XiYY9C1{(>3W5cvsarlEYNl^vbzkm(96H2;1``rrz$-`H1!TN`p>o-mxnGM7d(*$ zfi15zvfSWY!UJ#Ws(Va`_`RBjvY{CnKQFNdGG8U2>#+T_xPh@q++tE@E6wZ2TZe`EM+wLiW9+4-|Br}`|o@Y}^t967e(@?+9FOh;dN4m*!B3Lcmd>OjH+8 zSPD7+1+eSXZIl#w>|u7b=2>Dr?*{d3e(t>!j|4Xe=a#^_{@-0$jnMtU0jhUsC+3@*negPd{4Nr6W z@Gzn2QgHLm@qtoYLFYtUFI%IkAstTho_8@6W_tjbBGjWI9ZDIHIv7Z7k{$1_v96;m z#S9YjBQ<+^yD>idsfdWhnsi0v+bFI)V zO4BP53}U#u0JqA%XodLoUZ5Wo73e@E9RYl1&By6NI>h{-^i&u9!Sr1^^t08?4P8W3OGHkb6^$XK62k}W5{J@>m6Pd;)3(T#bTkT z+nEh-kt+H0lv+H1W9#{_{ET^S_0ANTautQP0!r^p7UblX`|X#7KMb)uaMEv=RESYx zQ^f&G_R|^COvlsXkhPGawLcP32chpN4#twKD(n@*n796l z{P;|ZaO-?{1v6oi28pmkw%4xq*YTYa=nfm69zBcz79S_s9YzP7QQxFeN0ifzYT%!5 zfnD1e!Hd7#xS?Px_ZVtJQ_bLT0F6VyTr2zZI>>gzaL}ZMhW*yb%*%v;LxSLNa@fl5 z&dLfgyno?C1d?2Uyr`<@60~Nq_$yy#Mp*C6bdyG3i%0i%}BRUr1(nMO5RO z89fbL;LH%2-Gil1EZm4_;Nz-(zXPZivj%tXnBW^PzHsNR#|V~*3Y$oC0-uM?G}HR& zHnJ;9kfdCp&4;r`W3iHDt#$uNpX<}r@!eMhxd7fUKnl_U)yS7CheOYED}pW@Vv*Q@ z0~Hy!x|3Jw<9evD(ASvZqdPNyD(4e)te{K~ zT*wG=S#VU{8zc*5Ta=;%J4&%99~Bqs^yOsQjz&}R^iP)GA)!0r_JF&`fFkHBSb<(jA8g5f&KL4VL4l`r4_yV3+_pjXjF1AkrpQ*=N7 zVbmS{r^i#6I3?K1>JIj3cRRt4BeQ{IyL1c*L|ts&7#D~(n1qCRqvwLkL_waXFdz}} zJ`}xX%6$X`{BAfhIKGWV&f##wcfcJ8c!gE*9vRx^i7c&Rglr#^0$DM`{A8t1;n0F- z1~y&xfY4*cwPdJ318>Ul2~ksL!<3QzJl1#58F)b7{Zm?ZX4Yn^205i!$X#uY!nOY# z>Fe|gY8SI)BDFDlwovh-!!n~l=kmPV$%W>rJ505wJy1Qus6E+nV!L`{!C z$-IAZT*vSTx*9RnR;QD5Vq@8Vtg^7(A5aEXzT5t~$*|A|ehbN~byX~cY{cYTd3mD_t_Z51*-v6+xii}XesBGGQXTO9 zX%&lUKd&F(H^#)jtq){nn0vK?cB9~B`{m@E3)iF}^wBGfnu=P%6#za`*|#O2G_0*& zFRj&P{%)h^>R00+6xa)jf$86MaH_D3uzk#S4g4cHd>Bz1y&LuGhA2vX-*WL$c!ggL z=^Ez!t!IhK3vY&x?d-Ekz?QDPplx`ObH#k%@^nHm*Viye1K^Ug=~15&&>*#llBVcv z4+@*4jeNJ{#f2QQ=9=6}%>ymqvo^Sn;O(bMY8PR^lMLi^D2y&)X@~f-@$AZE);N0h zLfuWN%BIxRs)OlDqYBMw8yG%2Xl2a<*xtWuQ&DGW4?rNwgv)0p(X7af6INwZ;1EjR zC?~MVnmOgcnm*5kHtQ8odLc1&=p9oZ*`6o#bg{;~iWaH_0TFM?_9#Cv?=PhbKV#y%_(kyr zl@rwsNb;aS2RQYT(}vs35fk=YtTs=dkCDGK+`yU1sgZO1!yz9D^Y6$Bgxv+4=reJbsC&nCMyFz@k37u-H1)=_9O3wB}5)d*NOD+L~aW4`~U@I~KMe#ci;x(wFG$g5xo^*== zp>YegUaIne?O%`UXkM086we2MzhLNpvdisQ!hdgJ_YS`*+ucts~BknMg{Q@8F zBZ#;}&xR-H98%$Fd;^6YL0ax2%@Z@0a*a-&s<`%ERIhTc1dXbtK5wk%nuhD zmuA*G^We{D9wjd)c+TTHOznHF5PtohmX)MEAf6;4m7(OhZZwn>cP~j;wpv1*X>liG zqY3#fO)VlHEeK>cOr>4&P+gdPRhf8)A1D$0@MyWUvaO7W=-t==+27?;N9^`iTX^Z$ zj1My9DcQN)ke0#?bk0`_fbfeFSfMC8r|f-@!t%99ii$$b2yqK%TXfCN)ROTx0!&n5 zD?Co_c%h-F7~aPvAKr9G1Z!vz9Jjyq1hRy4HyEM1Ho%4MJTq`^j>pDd2d~n4uFr4A zybCmlWkk^GYG3nW)^($eLKbRnF)7^Dk<5s#A-%tZ69B~Xh{a~&4ynKs!Z|1Z|LS`E#gVdOP4i10{wSh*vN@m7kWX{ZC)$!jC(nvgY=I>6tctQ zdn*Km>Kx3CEVb*g1t z;FJ9O4oS$Z%Hwgb50Y4fZNUgxy+S|uJg}s0(!!-`Fw}Rk^N;4Ug<48bB*V-;(L4(B zpgnWSWz-g+p6|_gYYg(nc`F&Y&-RB-2BSOep2(bDpXX&B{2B}*MvE7_rAcbog480& z+2oCp20`=(qkFh}xo{7o-^{Lc5YRPviX1aVwgH`#bpk?a>iiy~zY5hlSHhoYL5zHK z^Xp{&Y24y*7PZxw`p^)@O0j-l|CrOEGI}H0QMTl~Qe%eT;Z}_*r>CT!@2NDjjyq`d z@B2#kqo@v~m?hBdS%Mep3Dw}^rF%yG1D}H19xbJXAnu1?`m$B)%`Ycuv`klaasg10 zv6a!c*M)(~CFIkHBtnhW=_ZLdxo zP`eZ1j`4GIN);M|<$mI`e-8G9l?*6H=}z~O8h?~3^~iOsM_cukgxw$e*FqKymNga; z3vLb8OK&HAG!L`w7IYcJ*GU8&@0k!9NjDSeO+9IncL9HgEp4#t<1|l%z(cWP7kPkB z9G>yMNdU3=-504o&_m3C>G@zc>XA%A&)OBwc?Og9v$$a2JhP?Q zSf_cffNIkd_F%T0-p`+X#jKP>lR-J!4^ik|79$NX&Q^9Yl~^;qnedj@;GaLp1)a_&BADhe}pz;FYyfTN|me zBYv5dtSJ)<801fuBhg!C$~(8dFqnc@;%lGO>mF^N6gWhovK&r-;5WAs{uxyJ`jJac zsyq?9%cG-nIGYY@%R5-f<4gM|;uXh+yR5xWwXSp+UnFF3&wz&aN75C-m^%)vuBl37-djH@SW=4q|um0); zg%qQiXWIG=;OH`nOrKpr-7!&vo}oPh7}?<2ooYry8Roeu(QBs)I-|2?+yp9SCZyFD8(e94zA7tx zn-Ift(?EZzp~k3{?r@p#eoHIvscK9>QrI68^P+8q~%xsK=82n9U9-r=#j_Bn2XIFj~MqQ@?5SQxb=qlM!mX(Nz}U zQz%c*6)Em{s?XJ7Mh7#pi2O~6YtlZjX!til+{Jh8*fdzd*JeD zD_BG9Bs=w0@X_l2Nc!{=pu@lbnY?;&E4-4`m@5o{2@f86L~KMG+CNJu7;H(cv9iF= zd9SAM>AzSqcT#wM;CrTMm~n(v%SK$61TgR0^cxZmYG%CrD^%tNoEe0Q0W2TlI#@l0 zJ&Wu@=TdjY=?UkfxU+f1a#>5BMSUgaVK{$}wTossrZkEq)bF!cW5>>mevE#yB&-9j zxZ4h|pl~joaQ@OV%VOCugL>7jiKn46E=9}^oGk;nz8sl~8gu}p!1S&d?539Lnx@Jj zkJQhbvpy&yu0IGW8@fyxRkI8;Y6kh&8rsiOH8hj%`hP}^_%ZBDZ|IM&FOwbhBi1j< zsC&N)?(zRTr{Tr2Y?F`TR0|kF;s#C*2(O`aZ2Z9AjQW)lFPhpjR(bF-DJ=+beIZB8 zclWnhk+a1VQc_8mZ%)SmU0nm_SKo5-9}mstRrlQcIXPM^LyG>g5QdSR$4Sq+yt68H zG9;LPm<1!;gK^CZ#TW<}-L`a5mZjxd(1S$&Sc7*ILk=1Z7x%#Lx9<-5|DH_imE%|6 zxp?p3@3z6|gdz)IAq|u_fk4~0N|DoaTwtFlkIor!)s2gph8w!LCPCt+OdhP_K2`ei zLY0y0uk-q2Uy=srA}DF%VQj~w_oc8i>N#qQh!j3}T-c ze8EUkS^=GtMQ zxT=$1_y-&W#Tb0vZx)AD>8}Yyv)7oef$d3tgw&m3p4ND+83?ImN`HXDbajNk*Y|DU zocRj7Mb164GyUDlY5x#T_z9&KK*?sB7M?bxth${<&Vxm^A4ffjgl zST%dy01u+O0Z>#jA!GIFCnn_bOg+wq9%IPc0O zapZ2j_?}b8wV*;;&{tt8CFJaUEZ~_HopOG!5-Taf{s*qQr%VD<`xWzebf(0Qsl7Zi&_~z z<>4^`$drGIHRE<cR_B69*JLh=3=KR0jFlfX4Ja`Cj7_STU>y`8SI>BnL z4z+wh==pP!n~)=^zdf(8|KE|fn-}{5lEJQ=Ni^+vNNZE5Ax3MPOyq!SIW8MK@0egl zDg&#bM;G_O0C$J?CoPa#o;i`%Dt;8qhW-F|5cirZ1G1Ptm^4{J)~H;$4n|3;5?bNv z&e*PEaB&p(%smA0HL4kI;qeFC(*i=So_QtE5_px86}TBf=aD7$(O0;8sKLWpnXjE9 zHh*O6541i80Eb$2wwa2{*@PNoO2r*;anO%4!=pd#6L|_!0EZn^rwK^_qKqAy_goX* z%8!>nxFalJsw|fjEW7KT0N>XtOOa-T&VE47O$>D_lRkwiFr7AHw@%Q+06s*DGB`5O zpF|G)?HE9xe&xiOK2Cnxw&$1Hndhg|o*8+)19|@-X4xz4v+q}EQd7k8Q<4Qo50`Q*CeCwNiW?rVb@!1?0J0+*t^?4FOX>; zVRoK#l2Eqx(^u)}bAGzXxtPox=u6mWF`%5m!?DEXP*f{f&{{E^G!&m9XfmDUs{ZQB= zpbfofZRPYqgJ5Bjj_e7;3P4XdT~?I4;HG&myB%g|HXj<&yfpu@A~)bQUT%K_oFt@i9R>cy<)AIc)b4WbK4AZbwC!%rVhjeMk`c@ zwY%rH(LSCWJ0{zYdoJ;I=$ftOAE4^;=BSgM@ZvfhX`2up<}(bvd`}AyvdTrM!E~32 zR_PGtrol>I>NbM$X}lFG(!-JNEuaM*LeE|mSACi4=4#G^Lwb);9AXs;g}(Gto@szY zrZviG*WTcFO!_P?vYuZX*$}l|>-X9!dvr7Rq=?kmN;q@z>4rAb>Q9%t2?=ri9WeLvYpM8aBU+hz3!_BTOOkkops zL(3-LL3|Qt(jcsj1?Ol$yz}9TjB8ZGYi(uQIp#zWPz0@n7I$^KX;sN2@m4>2B+${G z&I5TuIHgc?N6X~cHL!Mf9|wn#{tg|(!%`tHB(J@mn0M{Ao&s7{2G27LQ`=3tK%4hc z6hslk$2kbS;SLV9N8a|y#|gxYbvJaY8>EfpS22u*4Uag8jE3c6y1=7&u=z5uXAmcYn zf1t4?d1_foW+j_9Q}KC1<%*g+9_eJ7^IvF1fA$^KxPFZ{*0K2F8q3AO@4OxCmJIPR zV4ant?b#*Dpc4yE;RV`~u!MA&ghVic{BkpHQ$l~Awfd7c!I1PR!g%6=q;mDjS_h)c zK|)RXNn{S=xSASWan>A3VH%$*za7LlvIjT`KeJ$?%?Q5~3Jk`vO^DU6ZifZ_t;;tT znkDFj!w%x7r7DL=rc5q!H!L}ilB<~)a{o{S_#iPLSDXS#IG%!H(`S>bOpg{!Sv!o9ywhXQTUUGCXR^;b`rb3I;> zW%iox9a${~MGXD@!~#Q%@;zrwlzBt*P|aUv6N#U3At@3@TrAt$_PSL7v($?QKD;w|HDQ#_EqS_j zw^)DRXIQWDq4Sf~6Cas{!B%U{v=>m;)}mtI3(D_0UiDC|*_P^i?UQzZ{LK$#F{By& z{qfWfI4W~OEAuF%zQB(MW;HG4rw9Gup{j#FRWlF`g4{qy zjQ$0W)eOvm!-C4Q{Pj3U^s3L}0l5zBb`Kx$q{}O*wwTf)*sO3;a;{?}R*kYjoB83T z78^?yC_g#Q5$|AW585s%BL6fhPO$n!|Fth*POUWj16GM`TTZ$sqS&J~)X_|<%y0%C z=e=%_6K2e5U8WzR`8Q%_{t7Oiw_n!P!KXt*X+j$g04XJ(q1o^LCt!;s$D_Nz>Jxjq zvF3GfxLIKTZ%`@mO9)dq@@U5?@WtYEqa1`BbRo;WB?oL?@7FfNx>39(7uaVayS?i6 zs#pJ0iSPF#GB0rs3TpPwpLS?Y+;3OT7{dN)fk%dnVfpDlC&~k8uzHGT<1|&W8hq1T z%$VS<8v&$EgM?;HV*~Mq3&RFsdF`omd>Sosb1w=R_WEhUaaif8L$m`I6yyK zj_QL4t}ox;vn}G8uJ+QvztwYP5jkFnQ&Qf$D@i3YO`*@G zU@G}V@KILgoBXkRfVJ#fO@VyuA+#Usy?K^KpB7!cs{i1%5OFRyLVR9vS$Ay!%bcDZkR&l9uG|+qTOXXBEwnLNCpMw>A%x3&JXg*V+(8) znmK8^>*&=IfILEvCtr6TnkA+G`KaX8-uX9i1MasOX~Q8p2U}UQ0>VrQP*TFxwa&JAzdimMY zv}23<23hB$2yOyQ^L*6uvmk0YB&lCLzPtd9x2-SvF#Cajk z>;Im4_s_Yf=eW3joxStz+c%L{kKVoOQ+@aEg^Mr$^F7zE>J*q$#-_AO9ugvF$|y1= z&5q@n>-jtFEq8hG*`_?7NB^NZ!!Ue>^51s)YtrmWxT|lFh9L%wCXdIh3dmmLDAc^Yx4{JWn~hl8l*KeDR(0 zf=@hh?Z<=r?#QPGDDCUGjvoq%%J%=WOoOH1xVimXto<0H$`ZOK-jC8JYA#!U3N8w0 zGT@ro){6gEW$Uy1LMmB94)c5qR%3o}Q&Mfp2V`ijC%VGI$glUE2nw1{?(^rXX_|GI zK$|xOllhAtVsU%Moj2tjMe`fqoz-v=fRaIpd#6leERQ@+4JE4*=MzONU z2p%}t{kk^_UN3#f`{asR7_%v_=w`uAYhKIOxkFwNFN~3Ls1L_RFnhz2$o&?s(SE}? z=)8y<#!T|3 zcZa1LbNRl1$Pu5G2JO1{M&HAykrkac@pCI-MJM8ysSHRPH2KqB#?T92)Zq=1lVS&Q z)A!F)vTH0Wz@_)8&gu}LG+jRN7s;BuIFa!bu~Z!w`oCLQ=oCA2U;Oz@{Wm(8COabj zUYtQ2;x>8Vn`6h5gxKU}j&%L!)?)75ZPGxoK5dqUXOQiJS9n8WQzsch7uZhpaZW(^R z{_<+^=UDmIGu`5$asA}&_jZPtZ-HQCP zIoi>yk$JVZ_%#7X&&>bR|Xe`;?d)Xf$^zxH__ooJ>3l&IKZ0zMi?iq=}-=C(0 zzdmrceVd&{z<#&ekdFmM2f ziQV!UzpXk${O1dtU&mQYEh;uPAsg6nxA&Bqxz2BW`+q{P1#uey`lA<%iGwAd_`lF( zvWQ{7-uir*bn@33N}h?}je$;0*xBx9jmX?xuqGG)r zUDY8di&jBK$4Uny+2cD&?s3rhb4sI+OcV{er)qOK?raCUaXx_h7u=ItDN5(KwF1k3 zl<&k1=Rw#zp9$xAto*T}dToM}d2mj%^pAn=tK3Qqiq}Bi$^TS|92X#eKJZ}u48+ycdHnJA?+*i{b>p!Qvf2HwzO8+Qb$LYc z!Tp=7eShxXb_oBJ-7-k`Q@~D0{m5zSxe?Y$h(4P87YWZR402-^3B3xL2*N~) z_ci$amKX`M`QZ_Vap}M5Pvq3Gf4p!*clq%$^$_CUDK+-vVDuB8pg5`n=0*^{RlT6S zbG}!~|Cv9Mt@qPkS#WX;h1K+lz6tTeBOZW6}00C_czF$oIYvz$U;K z@v;Zlfr<={jcWGjll0dQZ?qlKXHwcSCiOhyCdBzOHp6@4DpXD6$51j4Vb_E1k5SbG zD6%r~9*gB>`d``Yg)|xu`zemQZb0`(|q)f>Ex5ESPW^*+Vhlg z;yNhV`Jc5A+kYjORNuy8awjzC?tFc225NF~+po%LRZQg+kNaa3q|@C#K@r^QGwH`R zr*4@K`tK~v3oW)jnJ#9@{;&2l*({(DY;B8;16F5SNH4ITnVc+NyXEQ|$#z+~Jm6-m zYKH`hWGCqkbTj4T$D&2N-aAMS?6rT>3QG9T?au-K9?9PJ!V+1L3kkDUIRC4Ku>5dr zDo;T?NK%=@3{)F?>z%kA`yj1h1M8rBeZ5rMM71C>hkpMFqPo}q-n|kB86LYy4NIUo znlWaSSFAmZKQj=5Z#_;P;j@bqv$>rJOXx=oy?&RJ!VFcwRh30BGhX-;@b87rJ8L$> zD?b*B6|nZ=ieU!Db~`QE4TUi_7sit5Pfi|n_TLa1DR4Y?v4L}FF#g~7M7 zw4BVk{({lSRx>j4dsLMslsHeFXMvfSEILS!E*4?(6R?4YsdL%O;W1ONPNUA$xQ>t=>vD-EhC&m)!A zBKLYc)LiyORX*AOU@IW8OXQVF)3X<@|GUwaptr%uiHa=>alR#l;V8%S$WNxc2J|6z~H*~D8mViJOYmY-|)$4c}HX9Nl}oF!9*M8_m#jWV?;A;P?f!+~xb zeq2P;%!DKL!v8i{4cUmqe=**d0Mz22Y2iupf>)nbV80a2A zw!Zia@CdHC-i~1XkhJqwca>iM`PL^ob9XJq+2&0dBj`1y!uCSg6OJ_i-`m>Y2kGE> z-WR!yiycl78{4R?`m$w%`&Cy(bEb1H>fY^)%UeF@5CPOMyViEdDcIuAeuDf*=oaPe zA$o8lCE$7!VsD7_h$2F3^VJcF5(uk`!I-8_wKEGK7()0+cEY_?z}CP z?0Z+Yj1EbXW-g=!SgiI;hekJ#VC7z6Y?{)KG3oMSW;302nj81Wk~6TNe|}_a`f9nv z^umUi+vYdpG@bL7(lPQVA-1ocn$(EXmw0?qbx_0<+X?`KYQl zK6NQ~cG=?c#N-pa zCQ;A7^F_#mgTjv~0nOI5GDGFt3t|2S>s{~F;{8n?RxaU(*tNDFd>@GxD9bAl{<4IW z+dbo9)Jrjg`b2|>L1QOeQ=7|AC{LUD8acB z)2!_LH|W`$8O@DFFE0Rl zE=n%1i-r6HS0bMei&#?Jg{F#rS{LUm=Mc&tmcG!q!%YqIk2qZ(CgDpG9(k1r$V1-I z=Dbl<(ce0M=g+A)uZfL{*LLvvkuwy!ts39$)r57UToFNWDfhP(yFJD!@v8; z0iJRC$p1Z#W`3p)L)Z}s$xu(y#Bkp4z90#OGM5tD)g$hn!yfSceW?d#^MYlY2S__a z4aQZ3j>pyv2>U*&&eGdu@eI($6VIIzRKLG(***)Gvy6fOqgQxcgiJr0QF^%#?)w&6 z`Y;`g?Z%_&g0cT}iYfzKK8`<}EX4IW(b)!nfd+XXA`SKDlhvpz{LYHX2f z^No1FCImu0miP3X6s{z=wczAmTWZpH*7!=G0cHNK=NhhBH4eUxF?|juC zx!Ckx-Z__LYq)oL#TSyzxB0iuyav~3S8=UJYX99&hn|FcZBjd}rm@{3A3ddWKGpV$ zXk;^$lo(uJC(&3Ap0=y~JIAiq{Dqs}vzByP);?*}I0iL)Q7_r4UqpmHKLAvp3zP2J zL5no|89g>_qCDt%0m-)+pG_urTtzEt{05J{ zrD7sO=mKjR;`g@89lAFY@B1z|zv}w5J$>#;mgMtVEYhGu{Ot+PU&tM@a`#sK`nbj~ zAD%zC5dYEnK}}`831935Us@K7 zH-g3De3ZlsSJs-jCZAG=DmU-M`>F14GQVmsDrBBXI{H8M&Y~#}CQjqHhTt%`22Ws+ z;1b*+xa;5+Jh=M+LvVKw5G=R^9o(J4-Q8v1Z?RQdd+F*+SN*FGedtGi*|BU3$PSsT z(sU!eRg&ZVj-XHN9BYLtXwRbWH92(w8jrT+X+h4wH(^Bs}*4P77ajzgkw=cpmU#i>cSTYQ(n zoWj`d*D}Rpw{hnUk>0ZP{s45@N`rZW`UEjjRiqUL&Phr^^6k!{2#uI7o-*fS)RKuk z>`}^PxL!&@J?APu9M5$X>pyM?9-m1+Cv*imD<_}%@Rg#UEZ4Gcy?dKY+wybC7;_;6 zG!E!A4)Ecq=Z=yKEg_qg&$RK{{0u5T?-u3G!t6s(=N3f9fKbrpj@=XkK`SMn@Bj7s zKX)GQ54PuLLZ(=gc;cc!Ogiv^j<{uE?iNLMO5^0q^d? zs*xG=C3p{(OMG4^vP$B0L)A)Y8!*Y^zAyY61^k-bp7(sxIy$V;ZzjQ0&cBj0lDpuW z&cz19;t2t%THzsz(O5Wlln{&wif%OQ59@{>ExSWstNG>yR(nHJQ2h0|$V?GX<8HFk zR0l)7NoBKE$I?_+er^bF-&;>Cw{e-_0#xfIso%-s|#HvK#V(k2zWV4dLN>J2eo8<$`^3;_ysgb4u$2g((2a((SOy@z zYf)}d8M<;}c>Q)Vsn#OGMa*%+J-Op~UTnbIh55D64X{sB-CGtxf~Hsex%lPRx72%9 zz57Q{UN6}xslah5@U5y|Q)WyUn6H%r`Ey>7rdpbey6-T2+p4xe`ZWjT!yWr@wd-6S zELr6tIo*b%wKwctny!>g!ClMPpt^o@PSbeqUmfeFH%BC&{W6JxTXr=P z&WltH=4DLSe%ukb@99`OzkI&!#l;frqCTwOA+M;AnumHiTZOx(afsIse^S6#QWc8% zJ^MEHFOpdPIMEG><84CjUi7RrQf@MmUZ^Q&@ZmUh^LF%^zI+hUb(A8)Mk>` z()n;UzXPONWOMz@2fy$#yZUIq6UEGjK{+y>)P!EbOU8AA5{whbFmb1Cf`q+h!_i;d%b$;@1JUh=$Rd7I$gJQCJ6Vd(bj$vd zSbA;9Puf@9U;CcC4?yS7@NrupxR9ZZbw0kLPT@OFVzG_Im2}Up+ZF@{%P-;1XK*w~ z^BJ<%(BPr9;*MbYg#ltUk@%laKdEvPw4GQF#^Dm=1`PdfK&c(9n0?!igR4!HG8kII zbg+t6_VbYRdgbe89nj1IXT~d3`d_^?c+byRS=NVB`dE7t_s0`Xq=oV@ z*qI@3dk{h!UhCoJyDsgo5kRx+kM%D9{TNXz_P>Xp)e|T#>;K8U7KL-V)X0gSDLO7C zZeLgkFMeF1p^w9JEV#f#a`SfzLrK!HZAV8F%4w0~>}3Lzj7;vEIT=V-N~XSFBq>=G zwKh{gh3LcuxD|O4t*wqEoyCxB&tYF!sQVm0VTbhEgiO@4$xDT*$yD?H;&82aZJMF9!pwCp{tc2XN6SYm2gzoiTA_^dI|H$AeJ^zTjVvoA$qOfB@)X> zNAR|?!DN)ilkSzH5`ip)O|@@$<=**W6i(dpf{5P|PL_c@M$9E2d-s*JRQQX2PO4sbB4@D=1xtMf^jKB;o?L;m79Ju`z*I-;U zcImuoZO-RVlPpw9jn^12lq_jZU38&O)V<^Y50_bjLW=!HD<-l5zWY1Z*cZJ#-}iXq zq33{uQ4hVEt42Yxi^q3+ zT&-e!fDkv**ffyR3>=?vU5zbC5@-!)w>#-IlJ6k?<9c!B3FfM>`l|oJH9(T(Ii+KQ zcmSjPDtUKIOSTzh)r%{$(ONi|PH7m9`+599H|{3CN?25c%itWvY$KP-W1v56c@t0x zH=@!_&a~XcK!#F)9p=(-pVJq4Y#ZQNEd(cbOKvFsBnqxj<@MgVqW#fPc!=Kkcd*P> z`M}c052;qh1<^w$H)Jcbr|<~(4G`;%l{i8t4PIQ=4_-Av6(~w~C~5n6$W6)j)FPbd z;+M?XK9dSmdx8r4g?{W)KFW&`NP2uT+jBjW<5J|%l+HM;AQ?!%6c3Y5$J}l&!lV0Q zf@JU=-35Jl0Z7`1=~#TBe&jGgWubU`Fu#|Ny6%DpY>WgZGEMOUY$hsmhe+(LOW7+U zwYYepX%$q}D+Q=}uIF7>SlssO`u1c#5@od%|7}!Z!)s*CL%wT5gw^f?HGYbMlGXNNpx<(v z%iw2rKtW=r3I1yz=CvU5;xx|D(mTZS&-I1l^NOq%C5qb~vP8tK@Y-5`XwHRW}^e+RifK%-}i4{r0z#A8e2fj<x{QO1h`bAfh7U2^>(7XIPq)uEs^1-)0k?~>-^> z<^Eu2P4>SED{!I2oR##RrKaTZ=O>(Z*u8)!5g#T13Eb2zSrz585%_!y_q4lD$fRE7 z@aJc5O1_4xWw(A~Q(t#*J93YF(yLOOdG_Cr1N}WmByOe@y6p|(s0&&gAAJH^X)~2T zmrHt^&M-WjFGL@Ozd#~uMQw9$AWm?bR{z!CTRKsfX((kX(~31TQf*mk52lf#%=BpX zax7KBBqg~`>>#n(nlBXo!o{CjtW^e*O|sKH**xay*xuh%c>yD zNJv5{802}h{>K28OR#zh9W_$)_{YDm9`_}3h}j{=`OKuQmwg=f#?w^sXt<8}*S)!7 zA8a##uG3UeOHfqw`#kz#T%=)fs;3xRpz*ch5E8`lOGHG`gZas6gM`KBWWXSV{gXen z*#bQ|)t?~XNK+Qh1)9_fUAaG?;l`BKP|7Eh`-Wpc!GqCZ@Cdm1P$t>&%pi7W^!8(W z5*JJm=k~i*->he@<}0B`m8sltMv_jLvWYR z_(HtDXsh)$;fwnhqQs8P!=A1-L;$JlqGek#$8^yTl?DkChUX(k7sw$BgAenTbKnOH ze^9fxSjiJWx~`r2fm)OxQFwI3|H2o8W7_xa<)YWgggZt!uY5Z6i^qb4+!fALbxlak z6>NoL$xo?PPI)BF{-ITUJ%)b6#aKKQ-$G1GO3<_QSixBr)ZqS#T;Fk-Jym^;taK&^8Ic zGeMo(Dwki+{9SH;h(;7g>*UQ{ffqOQlcTO%F*UWp{Vyf?8HCWvQ)79BGgNd~uQ@^@ z5U}4b#_#l`qnrQk?B7Yp(dAEr~yQ#zNs*5fyxLs!$#T=~I$4_YJ!o)Ffq%w~( z74#!xSRU#sFPz28=SkB=oPAcezn9}CA^k}YiM!+%P=HZ!_4b)eab7v?%h>~=qSk#W zEEg6{{Zd5xG+`v4jpw~TBezZtW-G^k_^E&p+9mKhu6S?M z&pJiWAA%u47(9$lFPlI1t?aYhb4mstY&$!%5DSB|&u+{IpGmum*8W{EJR=U+C|ony zilOz<2VdQE^_y`gOBJ56U;+jmXqMuFgm7N#YUI{sFhYGW1B8;f`;o?lMEM6wv2ToO z>MNam5EF&`au7`aW6yS&c_Yqt2lxfmq);KRU6Yr(7t#eP;3J`|7a>%37Ws${VQ#+> zhEBVzgj^sS!Qgqn7quDoJH1K?Uc~jkL=}KNw=8xW65AQN_h}waLNy8M7SjoMyd8KU z9EcS?%y!^o2~Tr#O?m-QrG41eLQq3=gB_v1+I+nF@JfN#XDi&w`mr%gjo~NOtf%Kp(fO!Z!Oj%XH7Z3KQvlI7mknl)+zp+A zk#8Jo&t~>dgUX6f`H*X52jpo731QQK&*bI=`KDj_@dMGLkN_t!>^*6n{wq~e>cP%> zfO4g0MLLC8l!SBJt%d_@1Q&Uhpl8uRZP5o8hjJ3h4&sj*mM&)P z+l+5U!#wb<8i}E$*zZAn`52b-0j!ol zcaL$H&{mbY)a!aAlg7-I>y;;*ANu}yOZg`XUR}D^B6T8YA8OIRaU*~0hEz4Al-`m8 z{szAXTQOX`H5rS~k7@vmv^9&BAB9DKNZY0g5n-a_e`$hOCw2$%DbmXp zK9QmMkQ261f5R-2kQq}vzX|!{KM^RR%?7D{y4PK3glC?cnYw4@8W0%$T6HN#s++S) zYULphuNL|yBaw&*V0=YaSXSqqEG~VbC=N*n%uHrgus*X>?=H$I4h<~r%#TD1892?;>Av4c*&E@1(>PCB;8OdP% zPVKAW_>{!LQ}^^Fp%kA8ls3{Yb{hWEkH+}de-&tmxY;sBkG0oDzKxfB2gD zn8|+4JlJ%7y|+gSiF=h6*o`;Yu|Q5GjsSv7xEbq@GYA+S4|C@(0-c}aRrHYQMOY39J;kmX*rNac-E zrFso$&-w@bxsEhBPX0N(uhPD_n?@PAjT;@xlDp)OY*AMEj}Tqfuss z>4jZM&8zK$so8@>aJLT+Lu zV&fQS5Su7I_1N77{?HxokNRnf!?Bq~$+!>_2+AA!8%x5O`upNK>;|Qd<-C;tl@BKr zTlx?3PJ@ZtF`1xIyaUJbdJs$s{9Q~y;r;dblL$ER?@5zQ8f$C)F{}7}2@)M+i##*t zNWrVz-_mtgm7vMethwXWM8HU3blkTXLddGb0dd5;soEa-@cPI{IGX_fCWMSgwj2!K zTe@f3x<8%lVM89I+aRBbtmdB|_(~$#jb%?!fkVp>)C=6m^$*!x_A^r!f-7u99>Y({ zd?u&cE{JqwIwy6gqG0vTu^0PHTAZ3eWXBG|K7c#XLUne9Ok3ir$8*k8%-c$e^=!XMN}vALE~W3#$>R}pPQ$LCUhfc#sZS{8oCzL0T(LcAPIjGtd=rDC zMm|1P{zM-C*gnwg2l_{0zGsNXCK7xayv53YZ5myZfZjh$)(+fQ^h|qe6S40g~zno&JY$J*8*M|S(09% zXKo&qjwYO-NgfkQnStea3#5vEgox^|hNYFP98ag8CT_gp#`d>2PwokV5?ZC+O$YH! zT;?`>%tSu+t@95N7R(()PD`VW1moeX6e{XJmT}Zt$??>G89ESHcI)OjLKtuxuMQ~4 zp-Sidc)}##ld=6kP3>-1VMlP0B{G`HkeJY++ zmV?4kYCFkFuwu5A<03HpqQK=a+2}@xY6c0{xgDk*>fk17nOWB3*{!a|3S_T$C5r&Q(r) zH;YhRkB=#eVa-@l^PtL7Q$zGa4~p@B z5>JNKo;fHIu+<@^d@#-oz!rgXRh6jz&+Jv?=ic z+*GFX!)T$Ti7t5gA>9qhe9C$^)((KM*G)ZD%m~6rZrzyG z^Lf5Rg55hh1mJ@A^S~}hz$V_r+L#G)5b_00D;E(ZVph$f6;X_=fN@Lp5kS%q?YJM? zg1kCb^jOs|A+~C=iJ`L;7PtNeAk<<4fB7rR?*JGw8=bvx;?MdV^r_`YmvaW*&7kda zZYzSvt*4V?%Xp5E8d}rN?0ixR?s8aybs6O&j^K9Dmw*DB=Pa*!m=InEJB&}Sjdb(I z?`7Bc-|;E_qeA*8{eshim)Pc@QU*zYZH1K=D6Q`&{Eqnn!Esd1b|9N^M23slDlVGxwe+Mj3sg{*4KMp}%_XWMdK2GMOJCb6?v z>}%cHC+gNx8${WbBZFMPLg19X>1*ygo3lVizP7CA?{Fa4Gy9(em(jTfOxnIn_i)i z6XTo{<{^e56}1nl$9(%B8-fz>-p|zuC`>oLc;#&u8|H2np(2(m$MDU(bI{zEFR(Z zJ}tKKmq5dA(`nM!4IeCQ#up9bE2*nDg5-RxvEu029@+C&w!4+FHz$RyKEF9 zRfyb%?6aMFQ;3H}H!5B?fInG;_k!ANSm6EjXr3G_1R~)u_`E$7_?tjFz-tmLvauM3WkQ-$UX z;JdDr?h!_B7 zyQ0s$I|RYjE?%TB5J9sBJ3^@^jwbd_hrUf$AJPote;y>^iWT{TEI7}xq$!Lq;`gw3 zy|^%+jz058e)+nd#P?e`{!6%M`zV`JptGzxHuu8WxA}rqKKpfSxa^d|d)#Ufx2Sc& zS5vS5aEboXl4#^&ytZaUYv){W((#L zNEP$Yv5>3WK}A%<9}8EwI^n6=)M<0>bMCr6YoMehLz(1U;VV2b zwo_#?T%+`3KN;@zf#p&BjEwG5f?g@FqA_P(przBx(01b|fVxSRju1ZLk5ngKaLn-c zCokGx^P4f@Y&Z7L<-e%~D}6pDsUQ_p{*%c2W4KAy`CTMJF*eMFgx_ zviPkqE$*S_9ab3@pFiTFGT#0-Ps-yfyC0K;-PG??fv-{{%wqemdUgx>0@Ospe5RJE z@zpN#tRij8;A*`E`vW!w1YnPoM!o^QA>nEZVu#4>DuqB|OX)F@&~2Q^kdq4K6JB2WGdW6sfsxVoIKd3?*$uQOdqql}+7| zaah%SBQ^Diy#=+=9RZw>pBha8kD?N(&5PJV%Y(sRj0zmL%DJRG)^^W%W2lhuMvO}k zcatvvFqD46T|jXoTelb=XF;4>ma3ka(o9>N2yv3O;t4uTur*q}U^K~@UaD4xVS~?q zjTehhh*hYp4%=XU?~l|at)#KHv%a)d?>p4cZs-pEryl44w^&2GXjFnNiTPlFe=Vw( zv#)^SCX+J?gOBjOUWuks=-T8jLvbS4t#_Q6GW z+FcIe5nvTWV+Aj5sK~glfg5AbY$zO1Vn)%DXjIeiuMGUs1f1n z>XZxJC22be8&>(_3o2bZCdYk;I_fI6KcSy43nJ<0gD)n`!sK^r&HI@=Kp=tP7rpHR zVwx_y!J5uU|04F4qDp6m_?jT_I$*tAJUUpQf1iR8y9^)}kxKpxoZ=0kO=SpXm*9+iNY%4`!4s?0 zX;gL>5sEk4G!AIC9ECep^KBM?<$n4Ux&7j_iA>#g3j3XK+hgszYM_~+74%Foboe_7 zc_27UrJ2P}_3rvTr$Q=>I$j^mpp@7=us`J-aNovF5PN3YGLaMCK8cTRDIgU1JI7?a z?vhf9@w&XFU(jNz!oa=4udN1N2Dbi^dM}!Tofx2D(7Xl%G%k#-@z8!vtSo&`$#Eb+32gh}av8&hNBSyG|2LaWr%%8ny;?9k@6;W0$g? zOYAiyrdOs`>24O0bVP^Bp61BN+{497rQUuc(%VS%zReR9NU5vFnIm-=4%*%t!T*lP zU_H0$RbAaAlOQ^Fn&s%3Hnm`Dd3TS>I{Eyr=%Y5R!0~IXT#`4%(ycW7GAF@BD025b z4-;Ki0>x^LjnKZMm(QasQ+5~|?u-eTNT=lh-B)wiTNw6%WUHr!VcNzcZd;y?h4(2H z$^W1WU+4c2>6u|hk-9bcAVvs7h*!6E$bQpn&cmBdg|w+0q6ZVPwe$wx^MJNhAQzJa zyPb}_ldQ&{rLWZ`?@yDQO=$n}&f38!3(s2y$Pvxp3Nu;UtQG%L>pNe`1UE{*=>}X2 zm1M?3=Kq$lqFb`bDh*dfR9h(l;}rTwV_%a2*VN|Hr*%@IqLAGP zNawL$giv_j6vbZpKo^jBDl0r_gyJa~@UP74paaJupv9T>u6M-?wC{A*a#Z|n*lyPj zSfq>-7_9E~`O^Uu?)`@f`0JC?YdTQ(kd5K$G4dSTkHev5Y>UaZgq4*LK}j&uF7qxG z3j5#$L-2N#?rE&!J>WgrnPnfcH}PCK-YE4WHeZveSh{?OjpB*?RKY)=pDUfX{)<#| zY8oxIG}Rr=`d@in>HfA0M+v2_<1&@yLGVj@(7=*}b+Ods!}l5o@Z!GKP#bwj0$<2l z9lpa>L=*n%&He|yzukijf{ro3eS*=qug<;Jp3JZmwRhAL zXQ3g{HPxfe_W7XYHw$iFX>pNW6VliA&wq`zEHKNm_-$kKz+!>NKpKP~!4U{k=Mo(5tI> znW!Izzu&t*f>#r~TB!>ilhH*g%y6XWeENxgYCm9|dkPb`#ikPd?8}V{cvp--)yfq+ z$|X<_%157W$1&T!j6kQWnprt74^SDvCYzKnx65}E z(sEU9u=f%C=m-qybuury8(n8o`lr+ob^7C_Yq1SFx_p z5k|zi_PviQTfD(l?3^o-0jTEjB{klJn~{Gdiv(N=d}#D}va}<@auFV3U7ETnw%{~I z-0)T*$<9ok0jdX`Ab+ks@B;xTvexic-goFA*%=iW&qCjt&6Hf$q@{UWcrR87E zq^NsIPb0Z76h(Wjl-yBM@#oxweMKt^XM={bx>5?#jUtx$Z2$E^>Mp_UW-gN3^)}&= zM9n+g;M+!Kjv!Arz3!JX;7&~Se6g_iNLtX}F#jJ=fg9XN5q%#;jWsyk(x>^)J#yfI zx5wpqV-GOPuGxofWC){};d4!oK z$$Y1;YU|!plsfUAKX>v*uds{i4QNm~*$w^LS7^6I47r^<9m`y>V|pGV%|AmY8h$%7 z?|Z{MfpWx5??`H)a~|KDJY|A<(`DrvM!Oxs;KxCPMITIISyU(m+ElT6>gr{eQ@G%`PE~`;z0`LEBlY1 z_!q%V`nA_I^gd~SdVlaV!N6pc{Hwzniv%XmDzAGLMM5$@ExC|P^7~4`=NWbBXKbkH zPUfr+rL_yfp@N#()DkA(@y>pq3C3@d^np_m^q=3l1Lzk%mtnS372-Tmai%`s;+W~q z?g=KSBC!uwz1tSY(pK+VL$o7z-+7DEb302oTim!E5&N=iE!^-uGJ(YC_)yxf2g*3{ zXpmTN6EGem_yXmI&}#3c82SH83(yd*b~J)_`I~J~|AqyA>C<4@%))^5e97dKRxnen zP2E-|_YON#WO-ndi=xS-s0nt7ftTIY)HXu)?Sn={vgdm9nHz{6eD?6x3DN+LeDiEV zi&R;~o!7wxH>#FxDvdkl#l9Z^+`Ta926q%4NmAfiY~`eWb0Of98kM_=u*N0Wl~ynE zpNy}-sSD$#P&I4svnD4)M?oN4d4Ce8{O~=w5b$J3jkCO$%*IDPrBZ*y(R0sLfHnW zo)ZbEBDr0PzuIpz0{0p@V@9Kh2fK6;tz%$G4@3)N9>9cU2n;-Yj+I;Ms z347rC-7e+&Zq1dcK|nmeCy8i}4$8;7cx-=!5{bd1zkdaT{EIu8K+2@}pVuSO+%Re0 z1q6uC60p_o?CxP3h3IU0cnJYy(i+1{0$RjHw-{?xJ1ugs_b<#;%oJ@(Rblu~)YX$GA-kJ*>k_!#vnm3s67IgWJ6PVxgis)9A9do!WZ z>StA$qr-r?bMv@!9#)_pf94D%IjA=Jg|P>a5ZlYj{F}(Ui|{g1{ue3~8#|&M->m1S zgp*qb3x=e2l9PrfBcquc>z;4s{n1p*68E}SWhr3$Mk2!vMGJJX+P7JOPjb)|R;0dF zdWU%3i;7vo-bHuXcJ8;~!YrA=j*U{MS4=Yo?7__)DGdtmY~}D8LpMD8j%gj$604lt z?k&lwB=**dRlk2=)W4YWlg`tcesc3A^>5<3Xg=9yM$UB{rq~F1i#;c8NNvBcH@^ro zs-z43%)Mu;*p5=ptmr$=qjMFEPI!}!^XYJI%SaC$tjv2?O*Ml7886of&3G_rU<}}1 zd%c1~^U^|9VQ@87lyLt%P!$7csf(XWUL#D0lAGqgu*uH@p)bL0CdW?-9KL7B#Rab~WHjKd%cI`VQ=L~p3QM&LnrCu|sCm~3e zdj42Nh4tb`XjpXXDR9mmZZG@H{4fw_23E+epkqXWjD_I#?>~HkP&%Vz7bPwo=vDi) z-ZqJWwqG9Or@GRn8e8eC>Z|!^$>Lj+GBuVaag~m-eGcyFoDPheqYn;g_5Bfo54 z&0Y4c&d(h|k4iUSq^mLzhyRR|BF6AV$Yi^>L*(!SlP>{ydvDOz)4~CATtbk-<18lT z5DVDSSi&>im${A*?>cL>BYG62=)ec1U7{r6K|fVgi)P0Dqjs{{V8Uacekxh}euAj! zr-JNwMvgVcoAS^K8NjJhl*Bh#27z@FmJ$LtOx-`<<3|2-v=6)*U{DQ?(q2pGhpm$+ zGF`9VUEQ9GLoZm*OBsG(9Xe@WNVHc6j%FG5a8o zC%^N`IsM%V3Uf;KmKecb4X@_%>xZq?wHB`hq`eyIB$9QT)hN^7hhKM}XHovhHBxvm z8;F_JdMd65T+qdgmuvbNvr~uM^KGGy#w7F~m2)BcO-#NuDhPm?8xIz{j!SY3Z>$y6 zp8W%S{M-sx`>xx%y&9ZMa_7-Sof}*?4fBoL1(tp8=CO}C;ENo@ZEqE#=<9AwCOXZe zjf;2!zIy!$BY7P(8IFmjC@ILNecM0(^p-1{4=t{CmDM&>l29^~c!pGmzXj;Ium0%8 z;{E^S|Nb2`y|0sD5egC#wsFtBD}1wjUq3>}rWE>lvGFei^?zY={l4^)+m(toeR}St zj)%iBYU8P-$`5C0Yl*4V`}UeNu+VYU_9OK~|8f z+U-KngC%I!V(S$Dz7rvWir4B)k&{zYkg$#Qn&ii72oqO7UAMEVKS!(=#gf_YCVK8w zxPpecra$JV&?23j$_Qk}x@1GI7r_KlUG#0cX338T1>c`wX^*+xmxoLggZE!S-#KCX z4~yv3sZg2zSX+d-#wkZV1WZodZ64wJYPi{i?HW?{Lz>J_sTCbcQV=DZ&1?FU0BA|1 z&u})KA#^Xo(YB<(9TYO-L*M)dfk3n~YWRBra_a_jjGP`2XQds~-ZF@n-g4j)=jlYs zX$Ab+=N5{Up3s2fZ<`gQ$djG+()svNcckrBf~5d{3HY2v&~AsTQ&r4UTU`pB9>G@% zw>pIk#VdN_RpvvLjii-oQQ_0IisOn!m?B20h)KRym5P_x`sG}EcZ!#c!&U>h)iKi_ zi>*ib=wnB@-<14L>%gNn4k|mR72hn+ScBL3qDNhQ&;l?lzRk}#u|qW+Xt~VfCX995 zX}LHyOhJ>eqvXG*E9|`DN6+vt+wY1|)8D9B2h)s#M;H4p^f0t&>ha4KCcqKUldaC$ z+v6AL$z5`1uT4J`Y>M~&QI?0ZNl@P_MUWspLG@j4!{Os^^empCR6NKbR*Rdlu|G(< zI^f+jR!LqIz5Vb%46OITMejy8Dd?#OVquX;$LZ>1#&`|F;UrmnWxUj-=A`h`(SGHl zt~Ht$HWsHJ5K+tLoEo3UJOHDx8cyE^o?ZRG(}b?A zE=+GU@Z{cIxcqt=+Ed{HWU(j=Y zN}TbW3HZ3vQW`Vy98H5?MK~`cPzNPYM{Zi0tAt9vXG0%iVbG#4O-@JOEJ43i_A#!q zvag#jFpBjFTKo~sv5baN$j}H)`iqxlw;%i4r>15ePT+f(U%Ex6i~sSH zh|@dmik*NUyEi)xet;7tg#Qsvm(!*s--ZwIgD#2M^UpK7#Ix*tLV4am20~Lz^we(* zVQq0Nv(E_by&_2(gs5isxz1U&BdFo`nZ*a%D0x?!u)$e(IP01Z21^${=(8>JG;E?) zf!*_~_or$&lV6}8&xD;nTJwGEe9XGS(CsRiVK_4g9_2E^S>AC#Sf*ncsN5W&1_eBU zcK9g%=s$bSIzJ#X@&**7%wmDArbE_$b(sa`uqy6*+?(S(RmP`0z9awIyS@Y0S9bTb zxa%*BMZ3T{ed}w@ky&6#&%Wrqfk)x{SwRi77QT3&6hX;B@(~+?UVq`sXx&oMNhkt$ zHE&+)6IY=cA*rlZ$y}V8X+kL=GQ?lW1;i&GZ{ajcEy{w8gYi)RgBjHCBa{XaD-}72wZu}O15Xv?_(3kc6TGjm8hwJKcEm>XarWxc)oD=0l@Mfvh(ss5TFuEp}X+ z(}A}fquAvQZ2Z)NE!p|1=b6p45K+Fow6|vmnt!|Pz4y7N7Q}}S`NvHa#CBC1U2}`9 z>O*SF_j3*LXe%iive1;keCOBeq|F7%UWQPmFt||nm`nk_o5h^GRKC~`TodT5q)ZOw zt#otB*~w*gBD>YZx3z1k!@p6iM}Ek`@3wiegX(G!FQpTK4$gC>J@AtF9t9pD*^S-2 z`DREa4&~l^2@~iVu;};p9p=FW*u{^Da1_;8KQ+EqD1yd_=HQmXOwdSr@kxm~iE7&{ zYL~$>g6+bsAo=~rpy^bsM$Gj_v8SI0-wl`f^X{6LO`$76ip(%>5zQ~kkG|uMQ`aA; zb6tk60$F%e`vpP18Cz0BJuJp#f^pI{7k~R7@}CJBfLl=fTbZ6?#Re)3|HY)k%|88i`aSXdF7j+dmhUx znm#cbtZ5V6L zsborskh2ohQ1h0InN9Bikhzs9H>~Q4q=nz+O6^U7%S%IJCOUX!#!A%l_-0VE?c{2! zLmF!O!M)$N9vC}RNer6DH$mXhE$r+#PEf(!n&0LT;cz?s1#rQSyo}q16%$_}LGCKk zX!C)XTG*{Z_vV-+&oiLgh@>V&vF^3J_LnEZVbuPpcYY3Z)af6_-D|3tZ+44!Mpf)I zwyWPx;b#a&l5Yw>0A9$8B7S|4K9Eh6ETp4{1VHYbiR!i^OK+YO~66v8o9|8tLJa|nY z1%ak=@i2AYJnwQX-%yIK9zh@+PE&#}6 znQR4w+B9SBr3SL`9vA20`x02-4z)iDrorC><(hg54HLcKDM73$jvYD=bjkl zMw~gza&sI?KT~kae}_%owZ`JBeJxPj<&+v-FAe$3t-qr&H4L=8pF#T>j@Nk-uStvc9Lq?7Lc!s=^z=m5M#+{H4L1QD|Ltnxdi?+4m6Ntg(?K z7HkXTMe6kwLjenS^rGx`M<>G%D6MbYF{qZSc^o-f zlYY8zq36N(iH+8v<8DK2OAy3`Ler{-ss-C^{r%izY>HE$QmE?TSWIP}qygQqpu5p zThntc-3`%zPEr5Jc&B#2UXbONOpQ1-`^7Hi}pLU&aA zO=^5mWMN%MWa&v>RnL#lDzCTHR!%&f+Xo+AqAS9{Qaxfq^^g`pA=DkbqSWfG_^5YW`!z{p?ha;$Eq_Z2=8!zk7xMy zqtJeTs|;D$3P|{tGYrxwv!x~<#Lz9{SK)9QO?%1}tf$q(4pZo1#i|PD!dVRAehWd}5@FhQBr|(&Qy=mb^r9H>TthrT3yq{p( zTvh0MxA9se=aOHIN#M7gB}51a8BHJ5*z!N6rkZ%vQFettSgtu)zkT^IiA(wYww(;6 zXZg96Du$?==R}yv?->)vKrjS@d=WB)mOD5$(jmrZ0H7H1$M~-!|CeEdi_ClsUdXRD zD@d61(CT$|((#lAL4L+iwBpyqw)J6t_#pDa?~`UBiXwzBr%F%Hm@LdtP#og@7{0$4WbvpCHW;e7wnC<-eS}??a%Q}~qZQo#5=P(3t{v!6 zq!r~9LRl7FypP9zjj-u|kRJU{6k0@^h->N$NcY`m@4|H@NWUSj(p2t5!80~nua4TM zb8TNb70Dfg3_aDX?|GF3DSt9$NIUx=d+qVs4*Ln%;s0_w>zE}`@-y{`UOtDqo?+zL zLN@{0U!R+%?aZX z>;}4>&Z6B*f%e$QPEikfj-X>V9#JJP zLR|OO^JA()2Onog&)gBdndDM#U;9v+stkkFo>JT@g~Xor_*)y-IsDd{5V}3kO3c0@ z{pbdp26T9;8{EFmzJVJcODOp>a@6p0x!KPOfVBSE%2t1wZpNjhr|`ojHJ{(?gRLDN z8pDhmT9!6qO6>+_gwk2aM%7!&tnhi@E^d_uTLpHH{c6ltNnwzoyM4y3EsaG$YSr&M zWui!*@$%c*;s4FrOutiN;tX)eves%eh<**RHTG*x=`lB z6qzP79c1>%hLRq<-#B3$X4O=^f^c0jdtRa9u_kU^znLy5bWJEHc(f zoe|%GtTRuyfPMpkl$-H_zebVCH$3)fCw7o^813-NT1A$X9wq%b)c}ON6(M(j1z-$w zU*)I%M;WKP$?)L>h&MdUb?{ULt#-A+s`|6Y6o2QH2ff#TC9W>cb94w;r{&I^JJJYr za*;&2JTYKT{KfU@XC1OLo-p~BB^BAwde!?m-5uzjZai|h{0v#;V)MUy`W?`$A0aQ* z=7TJYLLS{LrM}rQH8-WtMGv5v$uzzpy zX}DnsoRbE_d|si*5}$RP-l`4=zqfH>H8i8FD+b@0?-(Jow=09U0)m0{?tQ-3*_=SF zka4%~NFPvJK8Xmu!vjnvd^NV0;(**sml)USj?BsTZ2i9W7}@%M%0t=rEf^R1nLkza zL8g67lD-M_Bhx1=ool9O!1(fC&0B}PfzD~jdYbJs(E0PZI&LX}rYDzt(9~ZrzGpf8 z;My5v^&PD$-NIf_*G7~UI_xrvot2pBDztHr!lqC95;GR5OcU`z^de)HHLjHzuN^m~~ilOb}*<);8z zG|wM~SNDLChH$!+LJKhV-!P)rl>n+3)&714%DQjETFYV7i|lkd*|O-2gJB)1`Qcr5 zF!+&wRA4aysP53vyqcay)(_qD%{SKoBk{J>Z#O9Y?xgapK^th(CgsS4$%1}B#M!99 zQ$Y1PBc_~k-c*lc+lQkWK+|nXu5cm=sFdcn540Ty>S_j;+-EJInafP-znqFpH?Ph_ z9i{At?md;YeLM+9Z3Zu_ZPUSMN^1VQzbj~kX1(NUe~s*{Yw1qjT?efLI$iS?l>61S z6CyVF5RB}Wcdo3Yfzh>-LD*AoWb2pr4+};)Ffxy=ki2jl+4_-u>GrlivY@U?+Oa!_ zY!%)$w)$29T71jdhq;E3&DPB+J>qZBnm=36Am9X=JYD$gOSF#GaaF98_dy^*hCOh8ubPI&*@SOz5} z4X$y^YoK)Em+b|yEKriUr^EkS7?i?1;^c~yk)`?#1D^x?kquKl-QmqxWL1J&$Z<3S zl%^dXif_CJ-MR9vpXx`DwTnBu@Bc{xx(xmA7q)wVa`BjxZgzm;`+J-oUNWH2>Vq6< zoI}=sGrNlQ??%=!Zjy`}Wj;LEyp$Ch1d6`8=el=&0oH@^BDMpjpeO%*(^w`P^fX*w zy(RwueV1b^Zhx$RJ++=RV#kE6Z9StMf1?fx1Ky%%!d`&lrQ4%sPF|o}&YdFoQxJ4N z+uakCRR=}??#1L=OQ1NMo`2U&8kBUtcJ=hUM5Y^3UlIf5K(Sdabt%sh^rRAV= z^>}9L8C5XJ^>L1hFh(XnQY+NArUJR(@Lc8EKH!wP)umPJ4eZRkV-t$h$SmDpR=uY& zGF3b)Y%zERNHXH(DcVcOf(^Csvy-(TrDx*INaX`m=U2L@FSQ}_&BA;hoN_=->%nx~ zd<9wHGW;tVsDRAAa$-37hcbV?xSMu(qXC#-#Gc-IFbWJjm_(X0Js5^^a3(%!0&R_` zsn~^=pl*53vQpy*7&M&rC;xehOwJb8N5^x5R?lPgqUUbN`m$-z+Vu%!>+PMrjjx4} zSxe9R8}i>M`>rQe`y43sku3d%hMF75?4ow~jp}94Z+oNz>5YnebszL6>~-!RIEyToaJF7K#D^@6@utK-??7gYWG=qNz5`E_>uOs47>L-O z4*OPc5t&BC1pmx$BU`(KFG*aX%+D|R3r22kBD1Hz%(MIQf_&{q`A0)?U@$Vtv(M`< zvMS&8mOn)b)UE_x{$5%K^goaE4_h&TdgQiBg7*%{{jgSDowx&{7KEp#$f{rvdqhl7 zBN|97KH()Bl>Gq9E2_LBa=`L}$bXa87PyRm#{GDJfWp_~TG#n9LH}4qZYI}R;8h$e z7{1a8G6(eJo^lo;v)`$0^EmU6sV=7E+7rA$B`lG{CiNZ^N@wyoh|Qoajyd&v@*+$7 zg>|&fc!PAw1iw?72W3BXasRF87-asqOGiR%AhOgmJx$KQy6C$AFD)u=#&>OUz-m4ddJAmdY zT8Nrm22EYJ?Db`7WZNP*_e4+w=)U%{(@xd^-Jr`dIZBnFJLH-A=SDFQ+m0yNrWJth z>GfZSc(_4ETXdpQyAPz_y5H$ga6%@hZ-&)aTmk9z7TT{BR6rwNZJ^(;1ssn%N;b~O zBily5X6_dGg6`dK_KR;nf{xJEw=L$2KqJHNFXFi$^muQs9gl1PongaT*;RRDQftMf ztD+pJ>q;%$cGrV=NMZEV5)|3#+iRHe%@(Mx=w1AMFddmZdcrAS)E?RHQ1c%&cSoj9 zp^;Z_>HvFK#JBII?~zFhy5eWozJb9b7D?$+7oeh5us<|^4j3-~q>+o;1-g%lR_pZ) zK--^CAFQ7a+7EV1#(tpGgW`gDbU(;I`{BsJ{(adXK66kY&ZP#}E3>9M`C34de(*KkjynJPO2UqbjfbWspR3wLaPDfZm#(>>Lv(=*^63ynbW>?6ry0k6zsZ-ELJj z*AtDPi=94wuqzH(I{Hz{G&>FSWfRcGG$YW)wBz^1mLrp-rTa#o`#_&{St=+`9CS2Z zdWJu*0)6dPJ`4SE(6_EGAGFZ}gWF9vPGv9y!!xlz2Vd2K&dH}b?Wco~NroF&=j^mW zC&=&EK1ylGyr({SoOS_NEc8a@|6W2izR(w-d&-b`UL-2&>RTY4Rc<0*76Yl=eadxC z*FbikkASHCUu2@@IqSS}3CPI#?wfl+3slaxyA#<>fx5``WBKGUko<5@L*$|-Fj&a3 z@Ff)keT;SQAr>8EIg;S^&-*+`R>)z^y+X}rWoj-(Wf=L>IvLMw4u+7NFaH6 zLiEbZBS59Eu=IeQ28gXn;_5GRfbVHHYlDCuGRZ?IGfEZzXug%Fi_JzZekbG`J2 z+!!*+xnp|yYbmh5`0nT{{~pM3&sQRnZvb7KA>)a@9%v243PSz^Nm&cA*c%(o;;ReflS?5kp2`m1+?PFcYnKqfIR0{!_tyqkZdX+ zwv(O(c>%@X$`eG8yCZm+=xqsxYeQM?C&NM7*sD{%;1Y;hnxR^=j3B2O{aWITJ{SjE zMJL&)A`=tme9sGVAX6CzvZZ0HAeX_Z(I`RL_Y>lOU{8cWT4C$k)VUkTB-Mr&KSw>t z-TnN;PUIo7wCqtaklF_lg}dWkl?zawQ#FU;Z?(vV#h;@wLXp5`cs{uN+6BsZBGNv- zDgvBE;V)>~6@lKqteAsq7}yu}WD|>qL5#VR!6HQ&*tQBe!Zy{wT;P^}qts*I(9IwF zw^|PNYZeP}danUVe!F(}MI{i$kW6`>STH?7y8lGN5;zQ8IF~Eaf#bE))PW8gWcsmJ z#771mAiiKN5p&1|wgpuWM%m9`vRq#EJ39qf{5ONwEoYEPpHs6BwYZR}Mk4j~L46=S zW*wEVxC>hPV%mXRvA`69rS`JSBa_`U_X>0hfgz!2%q!CTV@+pHzK#V!mZ}5a?g)bU#rV$KJ1C_tdp(fZqXC5Go0t6CKLN9i4Ks&6 zCup8Qgnpk51&P6)_>(Iepuk=uTH48mOx^OrCepuw_GanjyNN1bI_7NhKKe5<^Xi`8 z;^%c_#*Ce|bK(t{uCeHNic#uFs(AERW+pJD9!p+Ud;&~INiJhAIl=S>Z&Ta5LS+8w zvL?NJ5i(sX;2Zie2Mh){`=#xUgWiic+kO1fzWh^7+ODWo#t-K0Y0W)1LMV0IGT|>x z4w(OXN&mk8#((~O_t6E-^T@27stIx6CfMuN=P~ZP2Nr}cr@3aYf|lRV*~8SQfb|j5 zWZO1GIsVlJHx4nd%%h=64)mndOHy=uzaIsI53g;wm|MY6h(<`(J{ZJm4{h+d@B{PD zANMm{*ujwZ!Jj)(+8}1wSbK8+Am#hPfg{?K^U@648&rdWA-BJX_GcPkd{y3Rm6i-* zKGQ)r9~uDN`<4s;4xzwU`}Ho9TMsZ-6g~8$=K_lb^R=A^1t7C*m;33k1{mZ-h#eWQ zMkXUuFqSehGMULRo8P$vjOOka>HM<5py^O$w((>s&~ z69J%iVcg7C_bkOH&aBi+f%q{(xZ+n^kPGo{XiBBjb(txaS>yCz_BH5b=XHM&nUgF! zh=e0cEj)Q^(z3`h|A{41uPB)7e0Y4gApy9b24A&P_#NR`Bf8oL2^Dy;r9?EhZ04 z5iispH>H8&fX4mKqr+f!WYCxi6+kvSf2bbME&=<*i5aI+9k4$@9UB{OjchFMKTaxk z01-yDUki^)z^unXHAjXA$P>+V!WLz~14+r8$9q7i@y3afiasFs)?C>A)*1whuVx56 zRRy!4BQ93_Zor6zoD(bI0TKI&Cq{>^f?0+BIB#_on0Y$ZgdV;NW>O2ovGvSgAF3#* zM)igAz4wbu=ej|z%S|z$hY#7T{8b#HHjiwcs*&uTFa`OrZzh$ZUx`~I%O77UJ!E|WVnX!2dX&1r z0v++4p4$Ts6srg~1+1>v~|uZAjZL67y`bn%69VC6}?c8v4{L=)L> z2QFCx=g_Vo5qW)N>56yLnTo66ptR@r(?Bm|(yRH{so^$|JLb0ZjPoIwU9{)c$#6rq z6#r;{-@FfII`U)HfsSA%<|<4dC`369&pikO8N{d?oCi5*z~TJUmyMTr!Qr`L;ja)! za8Mk%et_Q-9IjS4|82_#i^IE!2^Qx8bCvD- zGeiN*o@90CCXaxbyy(EaePqz*$Plcjiv_)*#j3jQe_(cM;?ilgDiC9TeDOr|4lrL4 z(EjxDA2^6Mv`z=vgM-VoM|5!k<@w5D*$$-zGqc@Cgh`a@d(PKP_;nCSmKi#CZL<8s+B%1iD)*xw>hT=Z_qd_5B{0UB7ga z#da^UWR+X_S@t}zf8-vU;FSUU7fIpsjlRfI`Lfv_6KYVI(c+#uR19V-!GRVr2f@Ce zk*?5+9Yn;B2=@i>f&E8%>la?aU_UI}9R!6S7x}4h)-@mOn|(EhGlan`TH5E2&KnS+ zN#A3TUj$HDzfgZF6L!f@dorXgf`$J4D2KuznAsW>^Z5IK$g?jES*78iOP=OSu`UF2 z1_Pa^WBy=Pc`iTKTOJ6<4msz9QO=L~+TQGX8E|M(q$jLW1JnH8?l|)npt&4aVkcJ* z5M*9Vr!)fgt6YTiZbM|TlsBv};UkDYh*0%EQwlP=N3F^;s(|e_VKc*b4Wv%@3T1Uj z15M>H%N_#(U_W$&q;>}b;T`6xWCs$m0?*Fp{kn}T_VUQZPsIa^N^SXlo^v3*8u)AA z$YHR19o2VX;tZ&c2-r>DJO-Kt9~=!ntbnk_turA$ULZntsVlVQ23GeU-c%REKy=Te z@dDR&5Z&9dc5=9uay_%^qC3h#yu0;Dgo`MMXiE3_NWK9%RA@9=suF}x@mGqac!I`a zN}R(^8EAdyI{2KKQm1B@u4@~;1qV;%rgEj@!1Cbk`#K2~kn4ZV^ULogPd@#HqJPv~P4T)Payqw(E-NCP4T$hz+X8%mis7+ct;Q z^Bq#bjN3SD%}fEz#YRsb9&4i`!ym}0ipEc7}IZs9U( zdT_9%VQ=wF0OWz1d*duMh>T0s-7qWxBGmwGBJESk??V}=O@o0W*v~<)SO?6iznu0| zJOU0U7#YrtIc7&<}!^i86~v&I4Kb)y&~yTaetcjp&fr5Bh6!nlBTFK#7TQ(Z5R_ z^v6P~a&kw2t!af?mM;Pfo@xigs67Gw;a^P_btAVyS9muLn6@ohb4 zY%FCTWOL8xc4{?9DBj8~+ja!v@ynvjx|DH`>7|TxJpz@EbGMnBp(Q8OFqIz3bbz}?$K15q4ZPvhhQ0k?DR^UkZ(%@ zhpEavfI{^nE}q-j!0AZ$n=geO^mxjypGQ?dLbB{` zO4~T-`IJ2uJVsf+-<}qAdEyG{Sv7Y8e*6H3)1uCX`p%$I_(yJ1>It%}xkqzjkREiK zP7B}iwg65gw7w-v4p}?7ncpq!2)gfn)XvgmfJ&uC=l2`~AfNOc&r}FNrl+!zVA?5S?Ivq z07Tdbg>)T6mIto{`xd?hPPU^BV#jAePq@`4)#?=3U8|P($ioJfCtLT|vYi6^4{LMI z{vu#?O4zlT>M@9S>f7e$mjbzHs`(1#ys6Wz?h=ddfW{kBfeuw0WQk;YC%51mGD+K^ zC`h#phBWEQn)_}8Ipy)5x@k&%8jEoqdfE?2lc4hM zX)w^Um=6x;bY6|W4(SHo)O1hS+J-@+0B5f>Oz!jOfCgu*dlZ53FV@ze?afjc8)LzgM{&ZF7q#a$R@L#lZ&q!=pEU1{As%c>O;#G zIXt$=dfiJjzG)dmN61=+1(P5;5W3)YeI7KgrB3_5Tm<#8AF`t#aIUH0|IwPAtPd}}Z4UM+IX|42H9*S#M@Yh07ie#PR8{yI0JOFfjb)`@ zKxdSv@`e#3IOif04-@Vqi~o*T`rFTdo)T+7eZvkiL%7n?>KO}G6Xzm6F@=J)(7jTJ z3M!C)!_o~BJmBobi7}4+0INBbJH%rJLl-}*Meavf}b z^gMjz+y%1rx_f&(4}!C9ZrifkUZ8ssdvTq34{Qn_<{cb!0$O>YKaKwwz~+m5>OQ&~ zz;QJ&?XM*>DBks2QGcETHemDQ&>9yg86ZB+_i{n|>AU%Ger0fe8JMM4NvW&e2TxhZ z{DeJUT7Lxm_5!W_gVN_N5`gv3_U!Z*YEa)`_UQho2R8oSZRINOg8G`Zx`9L&*f{m+ zxCRM;O=OIER?Iz6Z)lrPM3+Fj{Hgc3=~9rNmTL<>Bmy#p9+67#CqQax{rN`lYmkt- zmr6qE!Q%9hQe(ArP>a?}+qau?-8fm37nxpx+U=ew!(V?uzCo(0HKi8R-qv><|K|vb zf-?)pQa*wW+qWDe^Y388E@UFc`IB<}FZYdqdxb1p_?8v%-GE&~NXW^7gTU16owQCv zxlTHLDS0-n|J&~BmA@A325MI()g~`J1F1)47uxB*f|}>CM&-9bASXEXYxT7tIAtY# z5vDZ(XOi`P=?e;=ZaKy~)kp(2wB4t~Q`5lNfKGXOq#xATwU}7SR=|3CDM&GS5_Sdr zK4VltSvPOa%pEDA)awjw)%OV5p!RtrJ;vt+$c2-x8?vq->*skE(6cDmFv=CO{rL>C z9}T+SpJ)T6ovtgUc3YtK$8!7XEam-?bDqXGj)4u!n?uZ|{lMxYtY1In^f<*)tWEcUSooCK&)5L)15aIeAol?qvg2`xX=nDuLGuloGz68(i-Pu2l z6hUh$@CD6f2e6^`=4G@N2J0s3kn@i+!MY>BHt4B7sEswFX9o9yg~zG;X@e8Mw>cyq zkRSop&73ND{xw`OIzz~y5~Zq~_6WPR8i^SH|gGNq$S62b~#9rWa3Q``&CnR_wukjDp1 z{~}(3>_Y$3S`2bu+1arKb%BA8&dN4m?eKE$c`F^VlxMfJ8AbqS_s7?U_CEv7WqwVT z%1N+36l|!@xdoPz++U?<1%dV|gL`4GD(uR7NM-o953G+-`zl=61)k5k-d^h61J zf_l`jcf-FsAZ0ldvaiz?ob_K)lRIC5wb9v*eS9^b-ctFf;Hec@w}_MuL{WiNjmu%L zqenod8KauFD*Zib#Cy>rU=^s* zo!74q>dl48!lxrahIrk{djA|)wF}Zp=vGsXH!?)-!xl1m3wfTIAPdeGCXI}C=HR&< zu6c@O90*I>^Qf^LsD=74Dm=RZvTHe+&YP*=>}DH2{*VJ4KUJ3e@VX3p$k}=BN&P@7 z>YKT7L zGLwVseiyb29LBTjVrX@N%wliC#-0bxF3jD+H7=k*_otF?!3NYUv{|c2-@$tF@by8% zD`2%EByK5p1+2wX7*9I4g6v~v>4Td0!KzkcJa_XI$a+T_eAX`q4X8aGX!IL2PW6t9 z!x_+cmXVQlVG-z9pE##2#DYe_K~v&M29R@3VEDJ_gv?~{{GpWC;KY;B>a`OOPKo1@ zOV49K&Hc0SgLNHbdj0W5A*HWituiGiC0<0SuXbH}NPPyhX1{9Q_!9#fn-8y(${awW zWxV&3+i4(bHFcgla0FOoyOk|ds(_elaXQO&4P>^SWd#>WfQCUBuPt34Xhc41_;r;L zcGn8m=UNPajkx5ahI4CR8Rf#-^5-Ea#PXGS?dX82ieo{;;Yo0~?{PMIEDtoap4VqZ zEQ7|Gsqdx+YGCzAgu9|Q8?5=iD(imT4_0+*pLPDagUlYkzrVBEz}mJz<|la-tfT7l z4X+-6J%vBDg^yZ;To+4{;S*i3R^O!?a3dG2_m&Hn|2&CI(p+Q-s(l9<{_#&jFQ$R~ z$*IGFwjVu5v+gLy8+-wICNGlX#BCrp^#nb-R141J%W99dkANC$sq-B}X;3|Uch}?1 zFJRO1T1G8+0&LpoeTB6>z~+13C*AP#VAI^Rs-H0iwt_lSwXxc;$Ex@0(v|aIt4F>~ z<01w&^=Yh0RL{V9xw@4@(+lk04%}^_)J?YRSG5`690!@ZH)JgD?0}k+d->6SpD6wR zb04@#dB2=n=fHUpuxaqpsf_vv)@LsIow;BH)<>hYzD(@}>u7Iw29*JjdD8fgeDydIw62D-DbVCqQzLJTC2B2S%R+y4V%p zgKE<|dozU=P!{k}@Tp?~$+988JzCd6vLgD<+Xy>Qr8~FhndA*nWpj*5iCo)?t|Y~T2r zJQv&nmCX`OCHF6vrR#HYp43B~J zl9qthI5Rkp86@h^eg)^1K|;@wP|EoUc7DuX1s9&Wn-`1=z=c(`;&Mwk=&X6%v%4P( zc098y@u5oK!tiC-XW}pA^LMW|($*qV8cWIftB*lU)JS+$h8@JWPgK*N6aj-LF6!jy zliu+s;v)=%z_>;Hn*uZNQdP$qECIq#oN9}VjA=k(E(5ipGt`1Q!$7t}2OP&e(q z4>~XDH8?U0fhOUpsIrmnG|<_JNbe?DH;4y|DZlM&?1SGXKqa zCrX_(epp7ZSOwHdQ`L1ui$EOS=LKptH8vJ$$~x@Mrun<^7Vv7R3P(%K{m*N4i@T9$mnuMA=$^6C z$7)a~h%H3OC4o8>b9U{j53)L>t!Gax19i#nxpk@{aD8O(Ktj0))Mehc9&ETysgGw} z4PU(l9o8?qU3)V?iengaI!am3RPFy9R*wa$!m-cq8;t;BNn82uZ$Vwd_t-z}W$CQD!D{OKNy}JTx z_ue?XQeFVH%jiejk2|0fr)t?OtOTxa(nauIsjvGLeWWCLL)*q;M_C^iVJ^&hpwvxA!vLK#7|GbXw47mTZ%nv$xA84Ik zENYe$K`YeL#ayX1!0#@xw63RE0p_mOSm|v5BO<15qC*1PzbpyXP>DL3fSfEp;`OujbR}Eyuv5Wa0Iw6WY_tLM}cvu zIZeXcGa&BiJJ)RW1T?*@-d(k8*V=3S=}8ih*Mg#Z3)Pm(19D3BnkT}U;ZpnPOXb!MsGIyN zxN9i|!)&@ZyXP{fn@JblN$Q2WU+SgG1YbZ{%yQk7BD9N|Eqdc`!pY6wgMDlj6b~#{ z6-je~lmpKI4katTlF^%R@N8G$(^4N5pJ9|kNyO!<=j%S>lL6CH|t`t_A=DHqweO$M#AXv z?e$wBpF|@V z-u=ar7+nHM)u@BFR9zvtX-dK-As8l+V^R$g%U~q&{nzV%`4FF}z4Q1f2P6_6i+Em9 zflYH?gw~^Hu%RXW2<5#EyzE<9lo$jXzB@mkEy==$&Ch+kmMs}M5*YpFY6*nY4bRqe z{e+%zZ!M6^=`kC<0p2YXiZ)C!> zb35$qzg>N2TndTdt=HDQ;34CW_iowHc9P-eN{!Q{hmc(6(03uX5!NN{yEZGohjp2- zOzYJ-*oQ9)h&R24jlhrI$ zPuYXoD*!nz9S^MUBO@Em8B1-VL5fL7jz5kM`k78Q^=%eGAzw}zHhRKyCbgue{UW5x zS|!R3jzGKbZ}{on1~Pn4>Qcqt2v7@WH!JfUfS0k$Vr*t6$dIs5*rq;6T+p)}_{tB7 z?+X_832X@Z`nD>of)4F~!c&YG0GuysRFw*`<{#hj{eRB^(#{iI@(jJdKOjdh&>Bcj zWRH6E3poF_^VM`3aI%a`JggkZj&>b7$KZ#_vnTGKFaq1AL*6g*#({Tx^QD(h5xO!} z`;Pkz@OJFo+rSJ)|HnP=MoYEZ}E? z#qTHV(42hO%UzcagVk!w9?Kksdc+XxRUi^}$9!YidN}4=x}FVr51nd0ue;6+ycS>8_F~!@;v?6jMedEj zFI(cni$l%OGkfD;`av9eHmb*xg0^kR4%lC>aXmMzkvvj8Lo6UAJFPjRK-Rs zU?t`8xyQ^3R)V3UJkQJ-enoGDr%b?dNBhmd=aR6@J%95e{8Y zQyFMPH!71THdD=m4Q1`5PsxR`24E^5! z!q8n0;Z2g3C%yAQEov{jWmy0x@>+h`9}|cScP0D&aE6#(!zK1pH{cYzCus9VKQbkC zEr_mfLZ$>;AF{}OB6D5g-7E@rP>FZFRqrwj3FWiuLn>KhBCj)|M0u7>ip?b#CkBy) zz_g$H;_DE!rKd{7u>rwReK+^{bdll53%YmPa>*pu)BN>rg0GQ9oeKK4ElNdJAg`FH+V z2!G~z-JaYBkuTq$pI-e6B4W3_9vuG*yJ8!`5Z~bPN^RNXQ~z>7oHdUSu?91a-i`PIE`Ua5^*g=`nH9Zw6D~tbVQXJTn|d zFD`X|GoFN2`P0BtN8HKK8Hf9oP6wf_IPWK=D}wMKpN;q5o`I6KlkuZSXV_0Ze)HF9 z2dF>lvBdoo>$4u@%5K+%VE$-+!7X+ev*zgj3!8&r-K!OSKbzsS!{^71L@{hqP2Ja| zw-+J;*~w$u82g&F7!<930+J1CD{s2@k>U1)E!DHG5XpYsHGcLMXn(%sIek=s<4;<9 zv*a~6N$Pm%e>n~-&wTg%-^nnTc8I>n6A6P>BmQ}l3oyFS_%61y3nFS?M@CRDE9j%9eV!@PM7_w4w`%g)l8s3%FDO{R7Vbf)!c3HGJF13 z*kv1>MdP>aOE~twt%%8y=NoTX!tw3a*we$=Sgv_f>m2_ZcyqoYx%652xd15Wd(!K3eLT5=ON zbXd0qmRPrfa)*nr_}e(N#ok${<}*QCeA1iaUmJ8&4he3i6+$buOj5bw0kn>X*U#(j zgQb{B%H4b_nFv|`?o78esIq&au08K%@L!#7XO68{{WkxC-N+I=EyR{nGc7>1sB$XZ z{TZGO3J-ru_`u^ve>UwKJ17yUa)u&-u+5AQJ@Q)?nvOO7YC8@=er?}9P5F<|YA_N`?D+i!aCSm;Q*b*#WX>%{*YKJ^)^la41b1-}^A2az+8cbc+ zST~4vL!*sl=anA~V4Bb9(cGK>Ta(hit1WIr{cYqawXdPDRZKs>&$5P5=fd5^I|ZOs zSA1n;b0!qG)W7(m*bZN|CAFy!cc9VzvwAUC5L&(wHk+qcLYp(^qZVHwC}}=T!T$t7 zNeVsR8^;X%w#{4qOCDO?`pJ)C&VX_?dsPCX4*Qzl@_I&Pfg?MY*_P}g_{!^0QbezU zQqj^+PFw)x(4cQr<{?m8`U8GwMZ)7&_1U7rgK&H_{J?Leu|UUu2x~-{4yhn6*I8=^kjYdyG>AdVzYObts-p3I)79Y z2SBgMP;ZFYAG96vD|Y|^$cRR8 zcou{IdH%j+~P*$;R0)g=9S}pUEb3F`? zzoiNU_<7sp`}dBMVEQx0t1Iy=K&=`iBtrYD04tgUL~ddc8p z1}$xRH%FV`lUcv*K=~C$E0MRDS`5wJxfEtTLx}5UJS$hzh3Dc4tBE zZ24^h1-4v+M80t3-a$b!EBSXu)F^>WJQVyF*!l^QU&K!CJD~vYXZwBy`U}82-^3wu z%`j}`gJ(NGdP3WRc{2-b0eXI{C47%%VCy{k^w0<&Y+aN%wFKTludCx}we z;`RgXa&}Q1tN5XPx_{A4y8-s|{rh=y8=y8Zamlb_2Pl=N-n_Hwhe7Zu=3N43;V^#u zxb-Pbi1IXtMudt()8I~)Qo~j_=xlIvQQ(Db9FxuQX$Jr9KO-l2wg@t;+>b6-go1Ke z!LsGyS7=Xe3wdz=14v=9LZe^_%sXcCZOnhcNNkOcv;GaJ?7R6@xJ3%4E`xh!BiLXw z;4|#qoJdAI#9xc=7{iJeW@i2MN1+{?V{)n42rC~pl^Zm!1T$^y>^b{M=!Wm`v$VYe zW6}1}Z`|c%@u`9Fpe>yYZJs)?!SW&0GyVOysM29vURC_bfT8!|2JdI`@q>9f;%(xx zRWPWy9ox`cN9ZEQ6PK&_L6~ES|LiAg*lSyDalcRvdt*_PX4h#b=~Qo*jA(`KJtiTF z$VFINi&6^grpU-D-JL=Uu8{vybak>I0oJ^K&r&*tpy>azRjsTC20s0+Dh~_6qEz*8 ze{4J~QVk;abG(L>X{YPel3&o5ci=JK6A4p`*cu_KA?yp;6yM$#$BH)BbIqroKx*^W z44t#uaPF<{k^cG(8i%!cg5Ilwdy`*H(wcU#@Y0@a8T$y4lVU6GXtlyKEV_Jc_&j`> z9u|(@ON7TzitFj0l~|ef%I${{57_!%vT3aL18v}ILBmud3=R4C9krH$!~Ius(s~1^ zoW8qC>!UPGbp6&u{*uCqy{avaOG%L1vp;ArcQ53gj-T5%y#w+mKj_>&`3@^GYah-r zaX|f3o!&7iD;TbtGF9cQh2i0gucg?dVLABLpJT)jM*cf&nnNN%eKrGDqf8iB#Fg`G zo`&U0-W!vr#>hx^*M`$K89eY_eeV&H4jXQt3f7SdGXCE0rLaXGY;NZFywYh0O;cUP z`J)nSj$eLI=oAN9wbkHFrc0oGjM{E-@F`R`TdfG;)&=YOh7iFv572fO?in}O4E;Zs zEvF+1m}ef0Yz^Rm;vv>IWxgAs%U-?m6Z<@z-%^DI6uO~Lvyh#c;)hM$T)?c&8Hh$Z z-F)*(5F4bvpTEC>8|)!+ay|~Pp&DPvT=MG?Bt$EIfBIMeTb^x)*^e@qs;-=I+5eeHf9!FzYWx?DTQW$<>0=g^Bx-VOmg0@dzM{U*5YQGWU%0R&FXbWEw=fOWcFL)`{W0$G$`F8uiFFZZ?L_qlmF>`3&!?b+dmLpwBi_S`oPmHszsP;6x zsJ|908!Pr@?%f9#w&fiSqg~h_{l~3{y&A0Q{uSDy+sKG+4IRNh!M6UnkX@!2*mngL zrD*+u(4cWqhhQgcnHA?s0%c)fE>PLP^bw%1jBY!aCo3VG;mek zdV5O;kY-_S@Qsbk2i2}iI(rCgHxApoKVpK-0~K9TI}K?wZG)et5fEF@T^@R+7POt$ zi48IDT8Q}ZM(7<(V(Q7d3j@s!tSKL|pndFU zqhHh^D2Rn?Nvqh95l2~p&(GD#c&_w_vF0%-)+(i!Slf}&_c9fu$DTkhitoi-mI!FY zFK4GV8AFfnTq1Ye1wyYn#C`GG5Ae>qO}DxBL3c&i3XzaVG7=#Zr#W0i#s)aH^=pkn zqj~ze0?a|VuUl|H!4XP6FSJ7f@4!Xf&HYy6ZwQ?{U!ApjfWaSGx}s&jlVMTc*EM>d zAe$z@Y{k<6Y04RT@#X8Fr0%(v8!86b<%X-~RxCo61pSLPzYUdTBc{pbKA=`EX{5&o zKyka=#<=7mC@Hh5i&VQpB5&b)%JyNL=ot8uJy5? z7?45d*x{q{wNP4IB`X3|(ER0QxzZ*;4c*1ES=0#BVDq#O$Niyo)$6s^={PcWS!(a2 zBVXZCAz8uyagNZRamgPSktTHJJnL=O3dx}0r)N?J50K%TWe232RzR#%>Hpej?$B&q8uPv$QqfI#ts`4>Sey$}rYl z%MXE?$T9yq@dK#oS{g%}1T%Irn#jKNB?+qp)=wunm z`K2Ay$294e*M&fB-XfKo(h4R2W9r$H8yR&(>RAp?BN=q5KS9eG2K8pVk*6bL{`FJq zTu$u-_2~oO>rt{~%=Gsap6?gPn8mSXzmJKaepqYB-TQa7&=`EeXAOGh$GQ=Mu zsLgK=NsZo_xL>AF%#EgGN4tVOGVZbZ!~j%E*aDyP$3v;W&#kiPD3}`)a_@I|h!=cqZ##KO|p=RXcNSg~}E8)QfNQpc4HfuVnB8nHcYXT>UT$s;$d!-rJr^ zCIc7xSLC%pmC3@~<6t>-QZ~i)Hv!_!`*n}*1at;{?Zmx~Lp~z7H}60##IsMd{o2LA zOQ{VV&+M0j`sk#iPxKlv2lelzl>|V1(@1(k9Rn{KjB|A*$bavXh390u9+`qfSk0=A z1A(C3Re3C0Qb5M6vLqV!hk|zd%&8YohRGnbj85prf_g8k$EM&RL^8^JSC%W0nOp3( z{Z(van&VS5P3AIGN{#N<>RCX!f|EC)@fw4V1XN9gm62gpzxoy}ZRm^*a=+eX4arW< z-el~8Y1}cpP?0rcnzK=@u|ERJ_v4>@4Jv{;M>9u^aSn{U2b>hNSCNV4^F7C3o*)ZP zXS`%`BA~U9BugoM3KgZFA$NJHP*&d}phZ_IqC#x|ZC%l@?j(dJ-+bLSe;+!WeP1>meE~Th!wXUM8k#PmS z=bS}N5H|YAHl`L&CN5{E{ut|pirejXX9`5YX_LpjwRb1PC>3h_k>|*i*sT165P)(Z z#<1ef88SYXV9|Z_JaitOcQY3ING2$&F4f)UgyNMIWxrkwkja>wr={<+AQtrLMDd-_@6#XM?~7#4YvZ_$n<{kIaM3+~orS({PL|BRRnV0i zG0V~~hR*XD!-O(+=xPY<4EB^I^f6oJ(3e(F(0}zNeIyICiJ!xQE6PC|-1Cop?^7^W zr*BTZGYREa-BUk;#2LIhMXawR12P|C0x~?Pplx~aYR^G)uqB^H)4UOA+jn`(Z_Ops z%lopKBa0zjW_%`Ggbf-|_0)fNq@l3m3{!W3BXqdHFZpyI85lTfEcSp=*ETdaPTzV4 z3oe<&#Yl5Ve`fuhd_WCKu8wY26I-C8+IY~n{wx`IktbK+R|t`Cg?ksR?n6&3)NZ)% z2Za5Kgo>pj$C8C5dEM4Py5&jMW~zLnQpyrIa;+svgl23_T)yG_hlWHR`n z-u*cZh)2)cY)st>rMpFxt_}CeKkPAVF;673mYQ2_{60aK78ED7)Cq+utNOO%3y{0h z?|oJ<3StuHwSIB^BeUmtEDTpoLW$FjlP1+nrr-P)9Lu~%{+(QLr^-GQf?`4~z6S5e zSbMJXrR%X!d~l&IyQ&$2V}j>01sU(<8+02zzzu3z+?wl_d!cw!qW<`&EQnVn&p&(f z4pdEPJEPV&U=}Uc%eqnkox`p-KJ5~OV&Jb+N5m_jQ~Y&AvAT>b3hZ|;q}_tX)u%bD z0vi6eRrlT->67k4Cf;#-jVY&rVw3v)w>l@3cu%H38G8b|r|)%w-v~j7zR$Yok`i>< zdi$*%g&5~v=^)4MLttvFAGu?t1&ti>zG%5yShgZTZvmxT5ugL*B!eVe^#+-;Nf)hVMC0N__l}E58-A%FYd<5A?y*IHb;UvWUUA1D7m0cR=ak(1+s7 z8^OHXDs^hZMzWL?Zu(HG6FN86-#$`#8OjB&UFXeKL5F6vy6tZ%8Rc(1kWn=Nq0LF# zd894Ct&`ojX^%UhUx_aAT$2tF*NKNXqXWwLYM-U|AuyHyD|=RSj^XEN=iEX(8EE=y zbL`|R2(UyNorrq|wxt9~4$G(Dzuo@aW6lO_5zC8$XqnL2yF07rLk{`Zvh&a93*HdW zf32u@xSUMCpDKUqYYP5@-(0WSUj^UWjzQUOOENIkeO0|A@ju^3?Kb>H0smIf4MK(9 zU`bw0E#qT`0OyUxk1w{8ftPDOXZWvyjpu@V>UN$)z#uQ`yCn)U8n5;coaQIuq9Re$O;} z!x^?=Vnyy}FvZuMSGmLvcK$!=pY*uFMrEp-?$83;o93Mz=QE&K%`E(N@GAJ7tKz0R zkAu0c^hx5m8Zb*|D!9M$f{@&es5skdGSDm;c=E$_GTHQR+mOOG$nO`?9o)_aofYX< z)}A>6=CJzYVXIqU>a5eF=$rW6tqja-1ztE$J;y^MVmJJxS*u#V8TKcWguorCs+GEb*%kIBG< zMNiuSdou8|NvlD>3hK3X?tgr)Fl-ld8~+nfzuLMy;AI^d_-N2HqH+*QGb`n~{PVzK zZL6@wW)>#tmtTxaJt1@{8%x`y%VhYj%bm+YSE1K#8S1QZ1i~w?9-n4f3GTT*mh_8y zWT54B-ojs1u!qO>_BRYe)w(G6<@`VBx7k(=T3A7T>T1=^xsQ-2oLhT_q`>5}@oVE) zZHRJg<)knPfI{tIwzTCUQ=d%!)TkOmu_k1%-GKn8Z(G#6G{2V&{PwafHZUUtLxOF& zo8~}SQ1+iqQ6luCYZR(?c@es_xI;r@JJjM6cYj@Ro{W%neYKA>$Uql!zMSZ9G9l6X zE+dwm3=Ed?D(p;zo#9UAR0j(}U%TVj{r#&5o!awU{dWpf_iX4Cva%v{GXJ# zRdj@ossxaM_S_FsS7O1o@b_d6=OCHRDNSJUH6e5_tt^ZD-OxPJ_k3T>DMHuPeXMr3 zmCzk~|1zKBf}+@$`vKph$iUnDKeN8q2tCFu`JpElq5GEwH$;44oS%iA`-M+1&QH6L zfQuIyH*@q;e#X#8zM`**=nSL&GmTalt|LRCXDVivj*_8mKQ>zLS7zjAir8x3VnTlu zzIHx)h73HnJ9o2jkum=#lZ!6v!N0#T=@tDpq1&(BBpN`2z>J1`hbhzlJFX*7#LCkZ zYC>&~rOb>Vv2Gh{j>B6*-@V3Cf5?^LUn%lvvm&7%Or;1@1E8BY+2wC*%HThVwf8?s zlY#fvOZQhY?t2Nesx~K)f!9gA^1FDTmpy_n*Pa^$Ne1Cp$Ch%9gj! z%Qr$(=GpE5&oDywJHI)j;R+cJ-CXaBhShS5x)U->Zyq*}QEV>L4vwowu~}A_LdF&-8uk zWz=7x>~zn5GL*a7P1EfVq37wag$G@Qrri3B&BAM-^eD}A?b{Hrk<+d7hp$7gF)8Zc z$96K%cjKbpAe9UZr)Z1IkC1^MDy|8&T#$L_BW$veLg)f27bHqEARSk7Iqahjp=YXd zor^z6=v8M!r;FIgz_{odxso4*KBIk!&X>p7FI$lp<+2bEm{N2V(1M7e)UpOSV>0{E zbHIOhE3^}n=)8?q(EFOdzTQWf5ufq0v*vtI(yxkpSz7@?a%IJ4&C3iNnDV^wpq-4o z-FD<+7c2ORbsHa;W-<0z@}|dg2Z%q>yq`SR0>z70$|lAc`_$w7%_PE)VN>&Z!{mWd z@8=!17G&b=-`8F_cgc8;vqxx7CG09&ZT)^607aLwsXurO26dbb>|QI$__}H}W4}L; zjxqVM$FYck2Q3c$@#0VqzTEw)lL94iseom}Yv9b>V+8-7gf6IPm7>MKHP&wm!<+TV zKu?UIaPT!K?oF#$xc8SV>9hr~CN3j%o1=eyt`;zKviMj))H&!Fl>Ep(ag+>v;_3MC zhEYGOVoN9Uc0>5t@0#Z=tc1>^Bl}rn6_jFcFV^jfhUgIW^L1$k?v`)}xUwAtML2Zw zuTMAQJYUV@BaHaT?4SJ<^a_GZFI`3)PXpWLEoH7U{3c2$dc3_$=uvxy?b!Z6WB&RE z-`NnDCw|{hdnO!|jWOmIO)SZRVP#8%2sarXO*Ljsjv#|?SAa1jt|4HWtrb;rGQ_~R`!Z@#8miyUd?vg>a{4$n@sjy=oxIWcU zO9p3)t6aKQ5W2?$%KHjM@Mrz5Z?&n1BGXvj?CfPSn9w<6FT4yK3svV@Qsv3ug-b`Y zGX%-t?utuO#pz_Qi{oX>mSYe|@#2fCxB}U^-A>)eA|oEM84jBSq5fz-`_-W}Wd6p% z*z<}z2%Z1)=gUiZP%V2Qw`z+n84(VsOD-#6)VC*hDt0C?@T&OS*_z+R1^uBRwo%|O=YX=*+PaxK0G_!{1KccN_&k1zrgm?TwKP@VKOcd z_ORG|2GX8CLN)dq5jxy+3M=@)z43{m@7Ybz_*HVk^Y3d=3O=(8$i4@AgxhTy;j3im z!L?T(wp@hvvO@KhWq{tk{o`bBKSV5^9MT;qhF~~qC3URHXo2d+_s!=(Gj}tO^lgLS z_wHlD+6ADH9k;!%u7E)s|C;&ZZn6Ie_;%hf&Spv&7?L{X@NUiF?{ufICbyYg`* z?rK1vuYG8Xh!m)~Y8@_qS&*!;4Pc`!GV=1wB>nd?P?!asBlpXZfj@hcmY6?*wL3e- zsk7<7EgAed-n!%4d+3Whe4l$70^O*thx>0dL0=}~eZVFz$oCJOe{Vo1BPvJY+l*wO z-xcN4V4n|4Qi6QIhn=7uXOdGg9fsa|!M;?-8mN5te^F)-LgrHjPNe?a4}Gtpcq)|< zcOB^n&&HF``ypkz(MuD0W7M8NV@c4Gd?|mPc0j2_Q+d)T1o|fHxbK<2gTC35i~y^Z zpvB!i;2(IFp$GT%28AYpe?;=UPf$L%k5%#CX;Xsk#vh5x{MJM7$zRTlg)`8M z-OvqHT`%*N!4qHAU+gb5fvO_a>J|GYh=nJ8lTN<|z1e#-W2+44Z}({&pgxC`cRVkQ#>`TU?my7<+y#;=)vZ*fi$!#uwm;MbZVxCfzYYTjna4dZ8Sl5= z7vX$a5t?OQew4~%kg}~BHi$k0rLWnV?O6=`Ryx)C`NC5&ds^=6{gaPC<#YGw6?TQp z3(fdzXVb{6F?U4i(p3n4zRev{b_e?Rn$}sXjFLfnd&TPaI&d&hwVW5^+Bti{QzXw**^{5&<+RHj|(Hb z5@b9jAWS;K4NANG!mb|4CiDu;+HN`n7i68P3T(L{8TPWv`cNS#Tc0*nOBN7%&tmtA zx5;qe7%rL)WrIfC&q(!_mC%#(x_-pg7y2<9FBOOhk%=ym`N2+BGW=v^XM7isttQB*+^`Jxlo1#%B-F^EzSDCrn|O8CdWwKp2wz2i7HBUQ0&w-ki9)V;Q(K zC1VVpGS2I3U-`)*ZiqhCIQU}qAEe{!-;DHNCd6!JrkvwnKq9Z_ zz=cN#p_ul8%i!8m=*Pv_CQ6(mQ{A&BYAZs(lJWf4uYIc^TX*OOZNDVMTZfljNGk*N z#%TB+flNX_a*Pr>%SDEsP5sD=mxsc6trzUa?!v63!gEk98RDZ7j17^MMS8CgM2)#S4{kn@QSod3M2~GS3OC2rY zbg&wvcqtL87t9yY@c4Xwau+?Xq8(<0; zru!BdfYW(T<>QD;poE>tuC9pz*RLAhOWQY+kxe(ueu9sh$MqpDk8*@sv!MouaPX8N=X3=MP`5-wqzf{v(Q&JHUS8@P@$G{$TgI((Iiy z2HqS_io|Y4K6E~C**wYuh2q~SVvlZvW4T$5(9LOZFQ_h@kYt942Q|#JeLZ-*R!3{7 z+L3?RvdgnG^2toeXQ6QV8Zwh+9BNZNL>2{&-$_fxl3^EDP4DHhSfLZ6z3akii0v(- zmoGJv#iF2d{|ugxr7(-v>y1~DDGB4*)Zj;Cv0Ap?qreZ^Z{FQ9tnee_!`BYi=m(Nf z6AgY{&nI9H<2~%CuLI#5x1WlaZzZ$-%pt0h56RU2hy;1+J4kfjQ}q+G1DkH=mfMM& z$q0|CL~xoCm_FoO&xx>t_{li#_n(a*Ik>oVi4kAcEAKpBZ<;3yF?R&=oh!j7(J#_2 zoetjbrz-yBRFdJEqsfV(Z^?8IKkvpaGw>Fd^pr*UkkLo1UYotXfT=E~?aQSbV7_Ix zCgpM(p_3GTm!AYuEfK}}p}}P4PUR(Hu1AJLi(km# z=4sVRZev1sIj?H{#)2#+KDMcpsDWyOfOC$zCqobEPVydD52n)P?ovq`$mqnr^Py|5 zf?4{(frcNPWYoS{ZGK6KOtxe!aP00TgD<~5W7ksy7yp(z|Fu(O;J}#=^wJ#LUFT z#KmaYql2PCw=*sT89$!?b_Q?R78CV925W50_W!x$2n&r44G0hUACv3= 1): -Option 1. When `screen_method = "nonSNP_PIP"`, we select regions with the total PIP from -molecular traits (non-SNPs) greater than 0.5 (by default). -It returns `screened_region_data` which contains a subset of region data for selected regions. -It also returns the values of non-SNP PIPs for all the regions. We could check the non-SNP PIPs and further adjust the screened regions by choosing different `min_nonSNP_PIP` cutoffs. - -```{r screen_regions_nonSNP_PIP, eval=FALSE} +```{r screen_regions_cs, eval=FALSE} res <- screen_regions(region_data, use_LD = TRUE, LD_info = LD_info, snp_info = snp_info, weights = weights, - L = 5, - group_prior = group_prior, - group_prior_var = group_prior_var, - screen_method = "nonSNP_PIP", - min_nonSNP_PIP = 0.5, - ncore = ncore, - verbose = FALSE) + ncore = 6) screened_region_data <- res$screened_region_data -nonSNP_PIPs <- res$nonSNP_PIPs +L <- res$L ``` -Option 2. When `screen_method = "cs"`, we select regions with at least one credible set. -It returns `screened_region_data`, as well as estimated L (number of credible sets) for all the regions. Note: this option is not available for the "no-LD" version. +It returns `screened_region_data` which contains a subset of region data for selected regions. It also returns the estimated L (number of credible sets) for all the regions. -```{r screen_regions_cs, eval=FALSE} +*Note*: it is possible to screen regions without LD. In the "no-LD" version, we run fine-mapping with `L=1` using estimated priors. Then, it computes the total PIP from genes (non-SNPs) and selects regions non-SNP PIPs with greater than 0.5 (by default). +```{r screen_regions_nonSNP_PIP, eval=FALSE} res <- screen_regions(region_data, - use_LD = TRUE, - LD_info = LD_info, - snp_info = snp_info, - weights = weights, - L = 5, - screen_method = "cs", + use_LD = FALSE, + group_prior = group_prior, + group_prior_var = group_prior_var, min_nonSNP_PIP = 0.5, - ncore = ncore, - verbose = FALSE) + ncore = 6) screened_region_data <- res$screened_region_data -estimated_L <- res$estimated_L +nonSNP_PIPs <- res$nonSNP_PIPs ``` -*Note*: it is possible to run without LD. To do this, simply set `use_LD = FALSE`, and `L = 1`. +It returns `screened_region_data`, as well as the non-SNP PIPs for all the regions. We could check the non-SNP PIPs and adjust the screened regions by choosing different `min_nonSNP_PIP` cutoffs. - -After selecting regions, we expand the screened regions with -full sets of SNPs for final fine mapping. If the number of SNPs in a region is more than `maxSNP`, It will trim the SNPs in the region (according to the z-scores, by default). - -The following step is not needed if `thin = 1`. +After selecting regions, we expand the screened regions with full sets of SNPs for final fine-mapping. ```{r expand_region_data, eval=FALSE} -if (thin < 1){ - screened_region_data <- expand_region_data(screened_region_data, - snp_info, - z_snp, - z_gene, - trim_by = "z", - maxSNP = 20000, - ncore = ncore) -} - -# saveRDS(screened_region_data, file.path(outputdir, paste0(outname, ".screened_region_data.RDS"))) +screened_region_data <- expand_region_data(screened_region_data, + snp_info, + z_snp, + z_gene, + maxSNP = 20000, + ncore = ncore) ``` +If the number of SNPs in a region is more than `maxSNP`, It will trim the SNPs in the region (according to the Z-scores, by default). + +This step is not needed when `thin = 1`. + ## Fine-mapping screened regions We now perform the fine-mapping step for each of the screened regions with the estimated parameters. The `finemap_regions()` function performs fine-mapping for the regions in the `screened_region_data`. -It first computes correlation matrices among all variables included in fine-mapping (SNPs, and molecular traits), -and then runs fine-mapping with estimated parameters (`group_prior` and `group_prior_var`), -and returns a data frame with fine-mapping results for the regions. - -If `save_cor = TRUE`, it will save the computed correlation matrices to `cor_dir` for future access. - -By default, we use `L=5` in finemapping: - -```{r finemap_regions_L5, eval=FALSE} -# directory for computed correlation matrices -cor_dir <- file.path(outputdir, "cor_matrix") - -# finemap screened regions +It first computes correlation matrices among all variables included in fine-mapping (SNPs, and genes), then runs fine-mapping with estimated parameters (`group_prior` and `group_prior_var`) and estimated number of causal signals (`L`), and returns a data frame with fine-mapping results for the regions. +```{r finemap_regions, eval=FALSE} finemap_res <- finemap_regions(screened_region_data, use_LD = TRUE, LD_info = LD_info, @@ -231,62 +168,46 @@ finemap_res <- finemap_regions(screened_region_data, weights = weights, group_prior = group_prior, group_prior_var = group_prior_var, - L = 5, + L = L, save_cor = TRUE, - cor_dir = cor_dir, - ncore = ncore, - verbose = TRUE) - -# saveRDS(finemap_res, file.path(outputdir, paste0(outname, ".ctwas.res.RDS"))) + cor_dir = "./cor_matrix", + ncore = ncore) ``` -If we run screening regions by the number of credible sets `screen_method = "cs"`, -we could use the estimated L for each region: -```{r finemap_regions_estimated_L, eval=FALSE} -# directory for computed correlation matrices -cor_dir <- file.path(outputdir, "cor_matrix") +We can save the computed correlation matrices to the directory `cor_dir` for future access, by setting `save_cor = TRUE`. -# finemap screened regions +*Note*: it is also possible to run fine-mapping without LD: +```{r finemap_regions_noLD, eval=FALSE} finemap_res <- finemap_regions(screened_region_data, - use_LD = TRUE, - LD_info = LD_info, - snp_info = snp_info, - weights = weights, + use_LD = FALSE, group_prior = group_prior, group_prior_var = group_prior_var, - L = estimated_L, - save_cor = TRUE, - cor_dir = cor_dir, - ncore = ncore, - verbose = TRUE) - -# saveRDS(finemap_res, file.path(outputdir, paste0(outname, ".ctwas.res.RDS"))) + ncore = 6) ``` -*Note*: it is possible to run without LD. To do this, simply set `use_LD = FALSE`, and `L = 1`. ## Fine-mapping a single region -We could also run fine mapping for a single region (or multiple regions) with precomputed parameters. +We could also run fine-mapping for a single region (or several regions of interest) with precomputed parameters. Here we show an example for fine mapping a single region "16_71020125_72901251". -We will first expand the region of interest with a full set of SNPs using the assembled `region_data`. +We first expand the region of interest with a full set of SNPs using the assembled `region_data`, and select the number of causal signals (`L`) for this region. ```{r example_region, eval=FALSE} -region_id <- "16_71020125_72901251" +region_id <-16_71020125_72901251" selected_region_data <- expand_region_data(region_data[region_id], snp_info, z_snp, z_gene, - maxSNP = 20000, - trim_by = "z") + maxSNP = 20000) +selected_region_L <- L[region_id] ``` -We could then perform fine mapping for this region of interest: +We could then perform fine-mapping for this region of interest: ```{r finemap_single_region, eval=FALSE} finemap_region_res <- finemap_regions(selected_region_data, use_LD = TRUE, @@ -295,16 +216,9 @@ finemap_region_res <- finemap_regions(selected_region_data, weights = weights, group_prior = group_prior, group_prior_var = group_prior_var, - L = 5, + L = selected_region_L, save_cor = TRUE, - cor_dir = cor_dir, - verbose = TRUE) + cor_dir = "./cor_matrix") ``` -## Running cTWAS without LD matrices - -In case you do not have LD matrices, it is possible to run cTWAS without LD matrices. - -To do this, simply set `use_LD = FALSE`, and `L = 1` in `screen_regions()` and `finemap_regions()`. Then it will run these steps with `L = 1` without LD matrices. - diff --git a/vignettes/minimal_ctwas_tutorial.Rmd b/vignettes/minimal_ctwas_tutorial.Rmd index ce5d4c07..ab83cc77 100644 --- a/vignettes/minimal_ctwas_tutorial.Rmd +++ b/vignettes/minimal_ctwas_tutorial.Rmd @@ -23,7 +23,7 @@ knitr::opts_chunk$set(message = FALSE, This document demonstrates how to run a simple analysis using cTWAS with GWAS summary statistics. -One complication of using cTWAS is the use of LD matrices. The LD data are often large, creating challenges in data storage and sharing, and in running the analysis (e.g. not all data can be fit into memory). It is possible to run the entire cTWAS analysis without using LD data if one assumes at most one causal signal per LD region. The downside is that cTWAS may miss causal genes/molecular traits in some regions. Nevertheless, this may be a worthy trade-off, especially in the first pass. This tutorial shows how to perform a simpler analysis without LD. +One complication of using cTWAS is the use of LD matrices. The LD data are often large, creating challenges in data storage and sharing, and in running the analysis (e.g. not all data can be fit into memory). It is possible to run the entire cTWAS analysis without using LD data if one assumes at most one causal signal per LD region. The downside is that cTWAS may miss causal genes in some regions. Nevertheless, this may be a worthy trade-off, especially in the first pass. This tutorial shows how to perform a simpler analysis without LD. When running cTWAS in this way, it is similar to colocalization analysis. In the Supplementary Notes of the cTWAS paper, we explained that in the situation where a region contains a single gene with a single eQTL, cTWAS is similar to coloc, and equivalent to Enloc, a version of colocalization analysis. Still, compared to Enloc and coloc, running cTWAS without LD has several benefits: the important parameters of colocalization are estimated from data; cTWAS jointly analyzes multiple genes; and each gene is allowed to have multiple QTL variants. @@ -31,7 +31,7 @@ First, install `ctwas` if you haven't already: ```{r install_ctwas, eval=FALSE} # install.packages("remotes") -remotes::install_github("xinhe-lab/ctwas",ref = "multigroup_test") +remotes::install_github("xinhe-lab/ctwas",ref = "single_group") ``` @@ -50,7 +50,7 @@ library(EnsDb.Hsapiens.v86) ## Preparing the input data -The "no-LD" version of cTWAS requires several input datasets: (i) GWAS summary statistics; (ii) prediction models (called "weights" in our documentation) of molecular traits and the reference data, including information about all the variants; and (iii) the definition of the genomic regions. We describe below the steps for data preparation and preprocessing before running cTWAS. More details about data preparation, including preparing the LD data, which is required for full cTWAS, are given in the tutorial, "Preparing cTWAS input data". +The "no-LD" version of cTWAS requires several input datasets: (i) GWAS summary statistics; (ii) prediction models (called "weights" in our documentation) and the reference data, including information about all the variants; and (iii) the definition of the genomic regions. We describe below the steps for data preparation and preprocessing before running cTWAS. More details about data preparation, including preparing the LD data, which is required for full cTWAS, are given in the tutorial, "Preparing cTWAS input data". For GWAS summary statistics, we use sample GWAS data of LDL on chromosome 16. These data are provided in the R package. @@ -81,11 +81,11 @@ region_info <- readRDS(region_file) ``` -Additionally, we require `snp_info` as a “reference” of the variants. `snp_info` is a list, and each element of the list contains a data frame with the positions and allele information of the variants in a region. We provide a function, `preprocess_region_LD_snp_info()` to preprocess region definitions and map variants from the reference to regions, which takes input of `region_info` and a data frame `ref_snp_info` with variant information from the entire genome, and returns processed `region_info` and `snp_info`. +Additionally, we require `ref_snp_info` as a “reference” of the variants. `ref_snp_info` is a list, and each element of the list contains a data frame with the positions and allele information of the variants in a region. We provide a function, `preprocess_region_LD_snp_info()` to preprocess region definitions and map variants from the reference to regions, which takes input of `region_info` and a data frame `ref_snp_info` with variant information from the entire genome, and returns processed `region_info` and `snp_info`. We have the lists of reference variant information from all the LD regions in the genome in [hg38](https://uchicago.box.com/s/t089or92dkovv0epkrjvxq8r9db9ys99) and [hg19](https://uchicago.box.com/s/ufko2gjagcb693dob4khccqubuztb9pz). They are also available on the University of Chicago RCC cluster at `/project2/xinhe/shared_data/multigroup_ctwas/LD_reference/`. -In this example, we use variant information for chromosome 16 (in hg38). +In this tutorial, we use chr16 as an example, so we specify `chrom = 16`. In real cTWAS analysis, you should run the entire genome. ```{r snp_info} example_chrom <- 16 @@ -113,34 +113,39 @@ z_snp <- preprocess_z_snp(z_snp, snp_info) ``` -### Prediction models +### Prediction models of genes We need a list `weights` which contains preprocessed weights from PredictDB or FUSION prediction models. We provide a function `preprocess_weight()` to harmonize the prediction models and LD reference. -In this example, we use liver gene expression models. ```{r preprocess_weights} weight_file <- system.file("extdata/sample_data", "expression_Liver.db", package = "ctwas") -weights <- preprocess_weights(weight_file, region_info, z_snp$id, snp_info) +weights <- preprocess_weights(weight_file, + region_info, + z_snp$id, + snp_info) ``` ## Running the cTWAS analysis -We use the function `ctwas_sumstats_noLD()` to run cTWAS analysis without LD. It takes as input preprocessed GWAS z-scores (`z_snp`), `weights`, `region_info`, `snp_info`, and will perform the main steps: (i) computing Z-scores of molecular traits; (i) estimating model parameters; (iii) screening regions potentially containing causal molecular traits, and (iv) fine-mapping those regions. - -Note that if the Z-scores of the molecular traits, `z_gene`, have already been computed when preparing cTWAS input data, we could specify the `z_gene` argument in the function, then it will skip computing z-scores of molecular traits. +We use the function `ctwas_sumstats_noLD()` to run cTWAS analysis without LD. It takes as input preprocessed GWAS z-scores (`z_snp`), `weights`, `region_info`, `snp_info`, and will perform the main steps: (i) computing Z-scores of genes; (i) estimating model parameters; (iii) screening regions potentially containing causal genes, and (iv) fine-mapping those regions. ```{r ctwas_sumstats_noLD, eval=FALSE} -ctwas_res <- ctwas_sumstats_noLD(z_snp, weights, region_info, snp_info) +ctwas_res <- ctwas_sumstats_noLD(z_snp, + weights, + region_info, + snp_info) ``` +Note that if the Z-scores of the genes, `z_gene`, have already been computed when preparing cTWAS input data, we could specify the `z_gene` argument in the function, then it will skip computing z-scores of genes. + ```{r, include=FALSE} ctwas_res <- readRDS(system.file("extdata/sample_data", "LDL_example.ctwas_sumstats_noLD_res.RDS", package = "ctwas")) ``` -```{r} +```{r ctwas_sumstats_noLD_part2} param <- ctwas_res$param finemap_res <- ctwas_res$finemap_res boundary_genes <- ctwas_res$boundary_genes @@ -149,76 +154,51 @@ region_data <- ctwas_res$region_data screened_region_data <- ctwas_res$screened_region_data ``` +ctwas_sumstats_noLD() returns several things, including: -The function returns a list including the following: -`param`: the estimated parameters, - `finemap_res`: fine-mapping results, - `boundary_genes`: the molecular traits that cross region boundaries (we will need these in the post-processing step), -`z_gene`: the Z-scores of molecular traits, -`region_data`: assembled region data for all the regions, -`screened_region_data`: region data with full SNPs for the screened regions. ++ `param`: the estimated parameters, + ++ `finemap_res`: fine-mapping results, + ++ `boundary_genes`: the genes that cross region boundaries (we will need these in the post-processing step), -## Summarizing and visualizing cTWAS result ++ `z_gene`: the Z-scores of genes, -### Assessing parameters and computing PVE ++ `region_data`: assembled region data for all the regions, -We could use the function `summarize_param()` to assess the convergence of the estimated parameters and to compute the proportion of variance explained (PVE) by variants and genes. These PVE values allow us to compute how heritability of the phenotype is partitioned among variants and genes. -```{r summarize_param} -ctwas_parameters <- summarize_param(param, gwas_n) -``` ++ `screened_region_data`: region data with full SNPs for the screened regions. +*Note:*, we used sample data (from chr16) in this example, however, in real data analysis, we should use data from the *entire genome* to estimate parameters. -Estimated prior inclusion probability -```{r} -ctwas_parameters$group_prior -``` - -Estimated prior effect size -```{r} -ctwas_parameters$group_prior_var -``` +## Summarizing and visualizing the cTWAS results -Estimated enrichment of genes over variants -```{r} -ctwas_parameters$enrichment -``` +### Assessing parameters and computing PVE -PVE explained by genes and variants -```{r} -ctwas_parameters$group_pve -``` +First, we use `summarize_param()` to assess the convergence of the estimated parameters and to compute the proportion of variance explained (PVE) by variants and genes. These PVE values allow us to compute how heritability of the phenotype is partitioned among variants and groups of genes. -Total heritability (sum of PVE) -```{r} -ctwas_parameters$total_pve +```{r summarize_param} +ctwas_parameters <- summarize_param(param, gwas_n) ``` -Attributable heritability -```{r} -ctwas_parameters$attributable_pve -``` +Make plots of how estimated parameters converge during the execution of the program: -Make a plot of how estimated parameters converge during the execution of the program ```{r convergence_plot, fig.width = 7, fig.height=5} make_convergence_plots(param, gwas_n) ``` +These plots show the estimated prior inclusion probability, and prior effect size, enrichment and PVE over the iterations of parameter estimation. In particular, genes have higher prior inclusion probability, enrichment and PVE than variants. ### Interpreting cTWAS results -To interpret and visualize the results, we first need to add gene annotations -(gene names and gene types) to the fine-mapping result. +To interpret and visualize the results, we first need to add gene annotations (gene names and gene types) to the fine-mapping result. -We need a user defined data frame `gene_annot` that should contain the following columns: -'chrom', 'start', 'end', 'gene_id', 'gene_name' and 'gene_type' for each gene or molecular trait. +We need a user defined data frame `gene_annot` that should contain the following columns: 'chrom', 'start', 'end', 'gene_id', 'gene_name' and 'gene_type' for each gene. -If you only used gene expression data, you can use the following function -to obtain the gene annotations from the Ensembl database. -We use `EnsDb.Hsapiens.v86` for the example data, please choose the Ensembl database for your specific data. +If you only used gene expression data, you can use the following function to obtain the gene annotations from the Ensembl database. We use `EnsDb.Hsapiens.v86` for the example data, please choose the Ensembl database for your specific data. -```{r load_ens_db, message=FALSE} +```{r get_gene_annot_from_ens_db, message=FALSE} ens_db <- EnsDb.Hsapiens.v86 -finemap_gene_res <- finemap_res[finemap_res$type!="SNP",] +finemap_gene_res <- subset(finemap_res, group == "gene") gene_ids <- unique(sapply(strsplit(finemap_gene_res$id, split = "[|]"), "[[", 1)) gene_annot <- get_gene_annot_from_ens_db(ens_db, gene_ids) ``` @@ -227,22 +207,22 @@ gene_annot <- get_gene_annot_from_ens_db(ens_db, gene_ids) We then use the `anno_finemap_res()` function to add gene annotations to the fine-mapping result, and use gene mid-points to represent gene positions. ```{r add_gene_names} -finemap_res <- anno_finemap_res(finemap_res, snp_info, gene_annot) +finemap_res <- anno_finemap_res(finemap_res, snp_info, gene_annot, + use_gene_pos = "mid") ``` -We list all molecular traits with PIP > 0.8, which is the threshold we used in the paper. -```{r interpret_res} -sig_finemap_res <- finemap_res[finemap_res$type!="SNP" & finemap_res$susie_pip > 0.8,] +We list genes with PIP > 0.8, which is the threshold we used in the paper. -head(sig_finemap_res[order(-sig_finemap_res$susie_pip),]) +```{r interpret_res} +subset(finemap_res, group == "gene" & susie_pip > 0.8) ``` +In this example, HPR gene is the only gene that satisfies this threshold. -### Visualizing cTWAS results +### Visualizing the cTWAS results -We could make locus plots for regions of interest, -and zoom in the region by specifying the `locus_range`. +We could make locus plots for regions of interest, and zoom in the region by specifying the `locus_range`. For example: ```{r locus_plot, fig.width=10, fig.height=8} make_locusplot(finemap_res, @@ -250,10 +230,7 @@ make_locusplot(finemap_res, weights = weights, ens_db = ens_db, locus_range = c(71.6e6,72.4e6), - highlight_pip = 0.8, - legend.position = "top") + highlight_pip = 0.8) ``` - - - +The locus plot shows several tracks in the example region in chromosome 6, including GWAS -log10(p-value), PIPs, QTLs of the focus gene (by default, the focus gene is the gene with the highest PIP), and the gene track. In this example, we zoomed in a locus ( 71.6 - 72.4 Mb) around the HPR gene, which has the highest PIP (PIP = 1), and the red dotted line shows the PIP cutoff of 0.8. diff --git a/vignettes/postprocessing_ctwas_results.Rmd b/vignettes/postprocessing_ctwas_results.Rmd index 4b9a09e2..40fba130 100644 --- a/vignettes/postprocessing_ctwas_results.Rmd +++ b/vignettes/postprocessing_ctwas_results.Rmd @@ -22,31 +22,31 @@ knitr::opts_chunk$set(message = FALSE, -In this document, we will perform post-processing of cTWAS results. +In this tutorial, we will show how to perform post-processing of cTWAS results. -Load the `ctwas` package. +Load the packages. ```{r load_packages, message=FALSE} library(ctwas) ``` + Load the cTWAS results ```{r load_ctwas_res} -# GWAS sample size gwas_n <- 343621 -# Load region_info, snp_info and LD_info region_info <- readRDS(system.file("extdata/sample_data", "LDL_example.region_info.RDS", package = "ctwas")) + snp_info <- readRDS(system.file("extdata/sample_data", "LDL_example.snp_info.RDS", package = "ctwas")) + LD_info <- readRDS(system.file("extdata/sample_data", "LDL_example.LD_info.RDS", package = "ctwas")) -# Load z_snp + z_snp <- readRDS(system.file("extdata/sample_data", "LDL_example.preprocessed.z_snp.RDS", package = "ctwas")) -# Load weights weights <- readRDS(system.file("extdata/sample_data", "LDL_example.preprocessed.weights.RDS", package = "ctwas")) -# Load cTWAS results + ctwas_res <- readRDS(system.file("extdata/sample_data", "LDL_example.ctwas_sumstats_res.RDS", package = "ctwas")) param <- ctwas_res$param finemap_res <- ctwas_res$finemap_res @@ -55,110 +55,95 @@ z_gene <- ctwas_res$z_gene region_data <- ctwas_res$region_data screened_region_data <- ctwas_res$screened_region_data -# directory with computed correlation matrices cor_dir <- system.file("extdata/sample_data", "cor_matrix", package = "ctwas") ``` +## Merging regions -Let's set the cTWAS output directory and output name, and we use 6 cores to parallelize computing over regions. -```{r settings, eval=FALSE} -outputdir <- "./ctwas_output" -dir.create(outputdir, showWarnings=F, recursive=T) - -outname <- "LDL_example" - -ncore <- 6 -``` - +If the variants in the prediction model span two (or more) regions, it would be unclear to cTWAS what region the gene should be assigned to. If this happens, cTWAS will attempt to assign the gene to one of the two regions that contain most of the weights. Nevertheless, there will be a risk of cross-region LD which can lead to problems. -## Merging regions +We address this “cross-region” problem by performing “region merging” as a post-processing step. If any gene has variants in the weights that span two or more regions ("cross-boundary"), those regions will be merged, and cTWAS will rerun the analysis using the merged regions. -We first select cross-boundary genes to merge (e.g. we could select boundary genes with high PIPs or select all boundary genes). We then identify overlapping regions from these cross-boundary genes, and create `merged_region_data` and `merged_region_info` for these merged regions. +We first select the genes to merge (e.g. we could select the "cross-boundary" genes with high PIPs or select all "cross-boundary" genes). -We then rerun finemapping for the merged regions. +We use the function `merge_finemap_regions()` to perform region merging and fine-mapping the merged regions. It first identifies overlapping regions from these selected genes, and creates `merged_region_data` and `merged_region_info` for these merged regions. It then reruns fine-mapping for the merged regions. ```{r merge_regions, eval=FALSE} -high_PIP_genes <- unique(finemap_res[finemap_res$type != "SNP" & finemap_res$susie_pip >= 0.5, "id"]) +high_PIP_genes <- unique(finemap_res[finemap_res$group == "gene" & finemap_res$susie_pip >= 0.5, "id"]) + selected_boundary_genes <- boundary_genes[boundary_genes$id %in% high_PIP_genes, , drop=FALSE] if (nrow(selected_boundary_genes) > 0){ res <- merge_finemap_regions(selected_boundary_genes, - region_data = region_data, - region_info = region_info, + region_data, + region_info, + z_snp, + z_gene, + weights, use_LD = TRUE, LD_info = LD_info, snp_info = snp_info, - z_snp = z_snp, - z_gene = z_gene, - weights = weights, - expand = TRUE, group_prior = group_prior, group_prior_var = group_prior_var, L = 5, save_cor = TRUE, - cor_dir = cor_dir, - ncore = 1, - verbose = TRUE) + cor_dir = cor_dir) finemap_merged_regions_res <- res$finemap_merged_regions_res merged_region_data <- res$merged_region_data merged_region_info <- res$merged_region_info merged_LD_info <- res$merged_LD_info merged_snp_info <- res$merged_snp_info - merge_region_id_list <- res$merge_region_id_list } ``` +This function returns a data frame `finemap_merged_regions_res` with finemapping results for merged regions, as well as new region data, region info, LD info and SNP info for merged regions. + ## Rerun finemapping with L = 1 for regions with LD mismatches -LD mismatches between GWAS data and the reference LD could lead to false positives in fine-mapping. Diagnostic tools including [SuSiE-RSS][susierss_diagnostic], and [DENTIST][DENTIST], have been developed to check possible LD mismatch. +LD mismatches between GWAS data and the reference LD could lead to false positives in fine-mapping. Diagnostic tools including [SuSiE-RSS][susierss_diagnostic], and [DENTIST][DENTIST], have been developed to check possible LD mismatch. Because it is very time consuming to run the LD mismatch diagnosis for all the regions across the genome, we will perform LD mismatch diagnosis and adjustment only for selected regions with high PIP signals in the post-processing. Here, we perform the LD mismatch diagnosis using [SuSiE-RSS][susierss_diagnostic] for selected regions. It is an optional step that users could run after finishing the main cTWAS analysis. Users could choose regions of interest, to be confident of the results. -e.g. we could select regions with total non-SNP PIPs > 0.8 to run LD mismatch diagnosis: +For example, we could use the function `compute_region_nonSNP_PIPs()` to compute total non-SNP PIPs for the regions and select regions with total non-SNP PIPs > 0.8 to run LD mismatch diagnosis: + ```{r select_regions_to_rerun, eval=FALSE} nonSNP_PIPs <- compute_region_nonSNP_PIPs(finemap_res) -selected_region_ids <- names(nonSNP_PIPs[nonSNP_PIPs >= 0.8]) +selected_region_ids <- names(nonSNP_PIPs[nonSNP_PIPs > 0.8]) ``` -We perform LD mismatch diagnosis for these regions using SuSiE RSS to detect problematic SNPs. This returns a list of problematic SNPs (diagnostic test `p-value < 5e-8` by default), flipped SNPs, and the test statistics of `kriging_rss` function from SuSiE-RSS. +We use the function `diagnose_ld_mismatch_susie()` to perform LD mismatch diagnosis for these regions using SuSiE RSS to detect problematic SNPs. + ```{r diagnose_LD_mismatch, eval=FALSE} -# perform LD mismatch diagnosis for the selected regions res <- diagnose_ld_mismatch_susie(z_snp, selected_region_ids, - LD_info, snp_info, gwas_n, - ncore = ncore, p_diff_thresh = 5e-8) + LD_info, snp_info, gwas_n) problematic_snps <- res$problematic_snps flipped_snps <- res$flipped_snps condz_stats <- res$condz_stats -# saveRDS(res, file.path(outputdir, paste0(outname, ".ld_mismatch_res.RDS"))) ``` -We choose large effect genes (or molecular traits) (abs(z) > 3) with problematic SNPs in their weights, and then select regions containing the problematic genes (molecular traits). +This returns a list of problematic SNPs (diagnostic test `p-value < 5e-8`, by default), flipped SNPs, and the test statistics of `kriging_rss` function from SuSiE-RSS. + +We choose large effect genes (or genes) (abs(Z-score) > 3, by default) with problematic SNPs in their weights, and then select regions containing problematic genes. ```{r problematic_region_ids, eval=FALSE} problematic_genes <- get_problematic_genes(problematic_snps, weights, - z_gene, - z_thresh = 3) + z_gene) problematic_region_ids <- unique(finemap_res[finemap_res$id %in% problematic_genes, "region_id"]) ``` We then rerun the fine-mapping without LD using `L = 1` for the problematic regions. ```{r rerun_finemap_problematic_regions, eval=FALSE} -if (length(problematic_region_ids) > 0){ - # rerun the fine-mapping with L = 1 for the problematic regions +if (any(problematic_region_ids)){ rerun_region_data <- screened_region_data[problematic_region_ids] - finemap_rerun_region_res <- finemap_regions(rerun_region_data, - weights = weights, - use_LD = FALSE, - L = 1, + finemap_rerun_region_res <- finemap_regions(rerun_region_data, + use_LD = FALSE, group_prior = group_prior, - group_prior_var = group_prior_var, - ncore = ncore, - verbose = TRUE) + group_prior_var = group_prior_var) } ``` diff --git a/vignettes/preparing_ctwas_input_data.Rmd b/vignettes/preparing_ctwas_input_data.Rmd index 806d1002..7b26209e 100644 --- a/vignettes/preparing_ctwas_input_data.Rmd +++ b/vignettes/preparing_ctwas_input_data.Rmd @@ -22,17 +22,7 @@ knitr::opts_chunk$set(message = FALSE, -This tutorial shows how to prepare for the input data for running the summary statistics version of cTWAS. - - -## Installation - -Install `cTWAS` package - -```{r install_package, eval=FALSE} -remotes::install_github("xinhe-lab/ctwas", ref = "multigroup_test") -``` - +This tutorial shows how to prepare for the input data for running cTWAS with summary statistics. Load the package ```{r load_package} @@ -43,74 +33,52 @@ library(ctwas) We use `VariantAnnotation` and `gwasvcf` packages to read the VCF files, so please install those packages before running this tutorial. The inputs for the summary statistics version of cTWAS include GWAS summary statistics, -prediction models of molecular traits (referred to as “weights” in the documents), and LD reference. Additionally, cTWAS partitions the genome into disjoint regions, assuming no LD across regions. So a user also needs to prepare the data defining these regions. +prediction models (referred to as “weights” in the documents), and LD reference. Additionally, cTWAS partitions the genome into disjoint regions, assuming no LD across regions. So a user also needs to prepare the data defining these regions. Note: It is possible to run cTWAS without using the LD reference. In that case, there is no need to prepare LD reference data. In the tutorial, we use sample data provided in the package. -Let's set the cTWAS output directory and output name: -```{r, eval=FALSE, include=FALSE} -# setwd("/project2/xinhe/shared_data/multigroup_ctwas/tutorial/LDL_multitissue_tutorial/") -``` - - -```{r settings, eval=FALSE} -outputdir <- "./ctwas_output" -dir.create(outputdir, showWarnings=F, recursive=T) - -outname <- "LDL_example" -``` - - -## GWAS z-scores +## GWAS Z-scores We will use summary statistics from a GWAS of LDL cholesterol in the UK Biobank. We can download the VCF from the IEU Open GWAS Project. ```{bash, eval=FALSE} -# download the summary statistics wget https://gwas.mrcieu.ac.uk/files/ukb-d-30780_irnt/ukb-d-30780_irnt.vcf.gz ``` -cTWAS needs a data frame `z_snp` as input, with columns "id", "A1", "A2", "z". Each row has data of a variant. We use rsIDs as the variant IDs. The naming of variant IDs should be consistent between summary statistics, LD reference and weights. `A1` is the alternate allele, and `A2` is the reference allele. `z` is the z-scores of variants in GWAS. In the GWAS summary statistics VCF file we downloaded here, z-scores are not available, so we need to compute them from the effect sizes and standard errors. +cTWAS needs a data frame `z_snp` as input, with columns "id", "A1", "A2", "z". Each row has data of a variant. We use rsIDs as the variant IDs. The naming of variant IDs should be consistent between summary statistics, LD reference and weights. `A1` is the alternate allele, and `A2` is the reference allele. `z` is the Z-scores of variants in GWAS. In the GWAS summary statistics VCF file we downloaded here, Z-scores are not available, so we need to compute them from the effect sizes and standard errors. Let us read the GWAS summary statistics and convert it to the `z_snp` data frame. We will also collect the sample size, which will be used later. ```{r read_sumstats, eval=FALSE} -# read the VCF data using the VariantAnnotation and gwasvcf packages z_snp <- VariantAnnotation::readVcf("ukb-d-30780_irnt.vcf.gz") gwas_name <- "ukb-d-30780_irnt" z_snp <- as.data.frame(gwasvcf::vcf_to_tibble(z_snp)) - -# compute the z-scores z_snp$Z <- z_snp$ES/z_snp$SE -# collect sample size (most frequent sample size for all variants) gwas_n <- as.numeric(names(sort(table(z_snp$SS),decreasing=TRUE)[1])) cat("gwas_n=", gwas_n, "\n") -# subset the columns and format the column names + z_snp <- z_snp[,c("rsid", "ALT", "REF", "Z")] colnames(z_snp) <- c("id", "A1", "A2", "z") - -# saveRDS(z_snp, file.path(outputdir, paste0(gwas_name, ".z_snp.RDS"))) ``` -If your GWAS summary statistics is of other formats, please extract relevant columns (rsid, alternative allele, reference allele, z score) and convert it to the format shown above. +If your GWAS summary statistics is of other formats, please extract relevant columns (rsid, alternative allele, reference allele, Z-score) and convert it to the format shown above. ## Prediction models -Prediction models can be specified in either PredictDB or FUSION format. -PredictDB is the recommended format, as it has information about correlations among SNPs included in the model. In this user guide, we focus on the PredictDB format, but will provide information on using FUSION format. We also note that it is possible to run cTWAS with only a list of QTLs, which is often the case in practice. See the section "Running cTWAS without prediction models" below. - -In terms of the choice of prediction models, cTWAS performs best when prediction models are sparse, i.e. they have relatively few variants per molecular trait. Dense weights may lead to a “cross-region” problem. Basically, if the variants in the prediction model span two regions, it would be unclear to cTWAS what region the gene should be assigned to. The problem becomes worse with dense weights. If this happens, cTWAS will attempt to assign the molecular trait to one of the two regions that contain most of the weights. Nevertheless, there will be a risk of cross-region LD which can lead to problems.. Given this consideration, we recommend choosing sparse prediction models such as Lasso. If using dense prediction models, we recommend removing variants with weights below a threshold from the prediction models. +Prediction models can be specified in either PredictDB or FUSION format. +PredictDB is the recommended format, as it has information about correlations among variants included in the model, and also has extra information about the genes in the model. In this user guide, we focus on the PredictDB format, but will provide information on using FUSION format. We also note that it is possible to run cTWAS with only a list of QTLs, which is often the case in practice. See the section "Running cTWAS without prediction models" below. -Note: cTWAS has implemented an option to address the “cross-region” problem. It performs “region merging” as a post-processing step. If any molecular trait has variants in the weights that span two regions, those regions will be merged, and cTWAS will rerun the analysis using the merged regions. See the tutorial "Post-processing cTWAS results" for details. +In terms of the choice of prediction models, cTWAS performs best when prediction models are sparse, i.e. they have relatively few variants per gene. Dense weights may lead to a “cross-region” problem. Basically, if the variants in the prediction model span two regions, it would be unclear to cTWAS what region the gene should be assigned to. The problem becomes worse with dense weights. If this happens, cTWAS will attempt to assign the gene to one of the two regions that contain most of the weights. Nevertheless, there will be a risk of cross-region LD which can lead to problems.. Given this consideration, we recommend choosing sparse prediction models such as Lasso. If using dense prediction models, we recommend removing variants with weights below a threshold from the prediction models. +Note: we implemented an option to address the “cross-region” problem. It performs “region merging” as a post-processing step. If any gene has variants in the weights that span two regions, those regions will be merged, and cTWAS will rerun the analysis using the merged regions. See the tutorial "Post-processing cTWAS results" for details. Please check [PredictDB](http://predictdb.org/) for the format of PredictDB weights. To specify weights in PredictDB format, provide the path to the `.db` file when running cTWAS. @@ -118,19 +86,19 @@ Please check [FUSION/TWAS](http://gusevlab.org/projects/fusion/#compute-your-own ### Running cTWAS without prediction models -Running cTWAS without prediction models: often a researcher may have a list of significant eQTLs, but have not explicitly built prediction models. In such cases, it is possible to run cTWAS. One just needs to use the top eQTL per molecular trait as the prediction models. Running cTWAS in this way would be similar to colocalization analysis. Nevertheless, it still has some advantages over colocalization: it estimates the important parameters from data, and it analyzes multiple genes in a region simultaneously. See the Supplementary Note of the cTWAS paper for details. +Running cTWAS without prediction models: often a researcher may have a list of significant eQTLs, but have not explicitly built prediction models. In such cases, it is possible to run cTWAS. One just needs to use the top eQTL per gene as the prediction models. Running cTWAS in this way would be similar to colocalization analysis. Nevertheless, it still has some advantages over colocalization: it estimates the important parameters from data, and it analyzes multiple genes in a region simultaneously. See the Supplementary Note of the cTWAS paper for details. -We could convert top QTLs into PredictDB format weights using the function below: +We could convert top QTLs into PredictDB format weights using the function `make_predictdb_from_QTLs()`. It will create PredictDB format weights, and save to the `outputdir` directory with `outname` as the file name. -`weight_table` is a data frame with genes, top QTLs and weights for each gene. When `use_top_QTL=TRUE`, we assume only one SNP for each gene. If `weight_table` has more than one SNP per gene, it will select the top SNP with the largest abs(weight). -`extra_table` is an optional data frame with information of the genes ("gene", "genename", "gene_type", etc.) -It will create PredictDB format weights and save to the `outputdir` directory. ```{r make_predictdb_from_QTLs, eval=FALSE} -make_predictdb_from_QTLs(weight_table, extra_table, use_top_QTL = TRUE, +make_predictdb_from_QTLs(weight_table, extra_table, + use_top_QTL = TRUE, outputdir, outname) ``` +`weight_table` is a data frame with genes, top QTLs and weights for each gene. When `use_top_QTL=TRUE`, we assume only one top SNP for each gene. If `weight_table` has more than one SNP per gene, it will select the top SNP with the largest `abs(weight)`. +`extra_table` is an optional data frame with information of the genes ("gene", "genename", "gene_type", etc.). ## Reference data @@ -156,7 +124,7 @@ To run standard cTWAS analysis with LD, we need reference LD matrices for the pr We have precomputed reference LD matrices (`.RDS` files) and the variant information tables (`.Rvar` files) accompanying the LD matrices for both UKB and 1000G European Phase3 references. The complete LD matrices of European individuals from UK Biobank can be downloaded [here](https://uchicago.box.com/s/jqocacd2fulskmhoqnasrknbt59x3xkn). On the University of Chicago RCC cluster, the b38 reference is available at `/project2/mstephens/wcrouse/UKB_LDR_0.1/` and the b37 reference is available at `/project2/mstephens/wcrouse/UKB_LDR_0.1_b37/`. -The columns of the `.Rvar` file include information on chromosome, variant name, position in base pairs, and the alternative and reference alleles. The `variance` column is the variance of each variant prior to standardization; this is required for PredictDB weights but not FUSION weights. PredictDB weights should be scaled by the variance before imputing gene expression. This is because PredictDB weights assume that variant genotypes are not standardized before imputation, but our implementation assumes standardized variant genotypes. If variance information is missing, or if weights are in PredictDB format but are already on the standardized scale (e.g. if they were converted from FUSION to PredictDB format), this scaling can be turned off using the option `scale_predictdb_weights=FALSE` in the `preprocess_weights()` function. We've also included information on allele frequency in the variant info, but this is optional. +The columns of the `.Rvar` file include information on chromosome, variant name, position in base pairs, and the alternative and reference alleles. The `variance` column is the variance of each variant prior to standardization; this is required for PredictDB weights but not FUSION weights. PredictDB weights should be scaled by the variance before computing gene expression. This is because PredictDB weights assume that variant genotypes are not standardized before imputation, but our implementation assumes standardized variant genotypes. If variance information is missing, or if weights are in PredictDB format but are already on the standardized scale (e.g. if they were converted from FUSION to PredictDB format), this scaling can be turned off by setting `scale_predictdb_weights=FALSE` in the `preprocess_weights()` function. We've also included information on allele frequency in the variant info, but this is optional. The naming convention for the LD matrices is `[filestem]_chr[#].R_snp.[start]_[end].RDS`. Each variant should be uniquely assigned to a region, and the regions should be left closed and right open, i.e. [start, stop). The positions of the LD matrices must match exactly the positions specified by the region file. Do not include invariant or multiallelic variants in the LD reference. @@ -164,16 +132,13 @@ Each variant should be uniquely assigned to a region, and the regions should be ### Region info -We require a data frame `region_info` containing the region definitions, with the following columns: "chrom", "start", "stop", for the genomic coordinates of the regions, and "region_id" for the IDs of the regions (by default, we use [chrom_start_stop] as region IDs). +We require a data frame `region_info` containing the region definitions, with the following columns: "chrom", "start", "stop", for the genomic coordinates of the regions, and "region_id" for the IDs of the regions (by default, we use the convention for the region IDs). -Here we use the b38 European LDetect blocks, included in the package. +Here we use the b38 European LDetect blocks, which is again included in the R package. -```{r region_info, eval=FALSE} -region_file <- system.file("extdata/ldetect", "EUR.b38.bed", package = "ctwas") -region_info <- read.table(region_file, header = TRUE) -colnames(region_info)[1:3] <- c("chrom", "start", "stop") -region_info$chrom <- as.numeric(gsub("chr", "", region_info$chrom)) -region_info$region_id <- paste0(region_info$chr, "_", region_info$start, "_", region_info$stop) +```{r region_info} +region_file <- system.file("extdata/ldetect", "EUR.b38.ldetect.regions.RDS", package = "ctwas") +region_info <- readRDS(region_file) ``` @@ -182,15 +147,11 @@ When running with LD, we will also need to add the column: "LD_matrix": The following paths were from the University of Chicago RCC cluster, please replace them with your own paths. ```{r add_LD_paths, eval=FALSE} -# please change this to your own directory for LD matrices. LD_dir <- "/project2/mstephens/wcrouse/UKB_LDR_0.1/" - genome_version <- "b38" LD_filestem <- sprintf("ukb_%s_0.1_chr%s.R_snp.%s_%s", genome_version, region_info$chrom, region_info$start, region_info$stop) region_info$LD_matrix <- file.path(LD_dir, paste0(LD_filestem, ".RDS")) region_info$SNP_info <- file.path(LD_dir, paste0(LD_filestem, ".Rvar")) -# stopifnot(all(file.exists(region_info$LD_matrix))) -# stopifnot(all(file.exists(region_info$SNP_info))) ``` @@ -198,68 +159,37 @@ region_info$SNP_info <- file.path(LD_dir, paste0(LD_filestem, ".Rvar")) cTWAS provides a function, `convert_geno_to_LD_matrix()` to convert genotype files and definition of LD regions to the LD matrices (`.RDS` files) and corresponding variant information (`.Rvar` files). -Below is an example of generating LD matrices in hg38 using UKB genotype data. - -```{r convert_geno_to_LD_matrix, eval=FALSE} -# specify LD reference -ldref_dir <- "/gpfs/data/xhe-lab/ukb_LDR/genotype_data_0.1" -genotype_files <- file.path(ldref_dir, paste0("ukb_chr", 1:22, ".pgen")) -# the output Rvar files use the positions and allele information in varinfo_files -varinfo_files <- "/gpfs/data/xhe-lab/ukb_LDR/neale_lab/neale_variants_hg38.bim" -# prepare a data frame region_info for LD regions with columns "chr", "start", and "stop" -# the positions should match those in varinfo_files -region_file <- system.file("extdata/ldetect", "EUR.b38.bed", package = "ctwas") -region_info <- read.table(region_file, header = TRUE, stringsAsFactors = FALSE) - -# generate LD matrices (.RDS) and variant info (.Rvar) -# update region info with paths of LD matrices and variant info files -# we could run this for one chromosome or multiple chromosomes -updated_region_info <- convert_geno_to_LD_matrix(region_info, - genotype_files, - varinfo_files, - chrom = 1:22, - outputdir = "/gpfs/data/xhe-lab/ctwas/LDR/UKB_b38/", - outname = "ukb_b38_0.1") -``` - - -We provide example R scripts to generate the LD matrices for UKB or 1KG European genotype files [here](https://github.com/xinhe-lab/ctwas/tree/multigroup_test/inst/extdata/scripts) +We provide example R scripts to generate the LD matrices for UKB or 1000 genomes genotype files [here](https://github.com/xinhe-lab/ctwas/tree/multigroup_test/inst/extdata/scripts) ### LD and SNP info -We need a list of variants as a “reference”, with the positions and allele information of these variants. We provide a function, `preprocess_region_LD_snp_info()` to preprocess region definitions in `region_info`, and map variants from the LD reference to regions. +We need a list of variants as a “reference”, with the positions and allele information of these variants. We provide a function, `preprocess_region_LD_snp_info()` to preprocess region definitions, and map variants from the LD reference to regions. -When running with LD, it takes input of `region_info`, which contains paths of LD matrices and variant information for each region, and returns processed `region_info`, `snp_info` and `LD_info`. +When running with LD, it takes input of `region_info`, which contains region definitions, as well as paths to LD matrices and variant information for each region, and returns processed `region_info`, `snp_info` and `LD_info`. In this tutorial, we use chr16 as an example, so we specify `chrom = 16`. In real cTWAS analysis, you should run the entire genome. ```{r snp_info_with_LD, eval=FALSE} res <- preprocess_region_LD_snp_info(region_info, chrom = 16, - use_LD = TRUE, - ncore = 6) + use_LD = TRUE) region_info <- res$region_info snp_info <- res$snp_info LD_info <- res$LD_info - -# saveRDS(region_info, file.path(outputdir, paste0(outname, ".region_info.RDS"))) -# saveRDS(snp_info, file.path(outputdir, paste0(outname, ".snp_info.RDS"))) -# saveRDS(LD_info, file.path(outputdir, paste0(outname, ".LD_info.RDS"))) - ``` -When running without LD, it takes input of `region_info` and a data frame `ref_snp_info` of variant information from the entire genome, and returns processed `region_info`, `snp_info`. +When running without LD, it takes input of `region_info` and a data frame `ref_snp_info`, which contains variant information from the entire genome, and returns processed `region_info`, `snp_info`. + We have the lists of reference variant information from all the LD matrices in the genome in [hg38](https://uchicago.box.com/s/t089or92dkovv0epkrjvxq8r9db9ys99) and [hg19](https://uchicago.box.com/s/ufko2gjagcb693dob4khccqubuztb9pz). ```{r snp_info_noLD, eval=FALSE} -# use reference SNP information for chr16 in hg38 ref_snp_info_file <- system.file("extdata/sample_data", "ukb_b38_0.1_chr16_var_info.Rvar.gz", package = "ctwas") -ref_snp_info <- data.table::fread(ref_snp_info_file, sep = "\t") +ref_snp_info <- fread(ref_snp_info_file, sep = "\t") +class(ref_snp_info) <- "data.frame" res <- preprocess_region_LD_snp_info(region_info, - ref_snp_info = ref_snp_info, + ref_snp_info, chrom = 16, - use_LD = FALSE, - ncore = 6) + use_LD = FALSE) region_info <- res$region_info snp_info <- res$snp_info ``` @@ -277,94 +207,62 @@ Another potential problem is the LD of the GWAS data (in-sample LD) do not match ### Harmonizing GWAS z-scores and the reference data The `preprocess_z_snp()` function harmonizes GWAS z-scores and the reference data based on the included allele information. -If `drop_multiallelic = TRUE`, it will drop multiallelic variants. -If `drop_strand_ambig = TRUE`, it will drop strand ambiguous variants. ```{r preprocess_z_snp, eval=FALSE} -# load z_snp sample data z_snp <- readRDS(system.file("extdata/sample_data", "LDL_example.z_snp.RDS", package = "ctwas")) -# preprocessing z_snp -z_snp <- preprocess_z_snp(z_snp, - snp_info, +z_snp <- preprocess_z_snp(z_snp, snp_info, drop_multiallelic = TRUE, drop_strand_ambig = TRUE) -# saveRDS(z_snp, file.path(outputdir, paste0(outname, ".preprocessed.z_snp.RDS"))) ``` -### Harmonizing prediction models and the reference data - -The `preprocess_weight()` function harmonizes the PredictDB/FUSION prediction models and LD reference. +By default, it will drop multiallelic variants and strand ambiguous variants. You can change these by setting the options `drop_multiallelic` and `drop_strand_ambig`. -In this example, we use liver and subcutaneous adipose gene expression models. We preprocess each weight file separately and combine them in the end. -We specify PredictDB or FUSION format in `weight_format` and the specific method in FUSION by `method_FUSION`. -For the expression models, we limit to protein coding genes, by setting `filter_protein_coding_genes = TRUE`. -We set `scale_predictdb_weights=TRUE` for PredictDB weights because PredictDB assumes the unstandardized genotypes (see the discussion in "Region info" section) -For PredictDB models, we set `load_predictdb_LD=TRUE` to load pre-computed correlations (`.txt.gz` file) between weight variants. -For FUSION models, we set `load_predictdb_LD=FALSE` to calculate correlations between weight variants with LD reference. It takes a while to calculate correlations using LD reference, and we could set `ncore` to parallelize the computation. +### Harmonizing prediction models and the reference data -In this version of cTWAS, we allow the joint analysis of multiple groups of molecular traits. This could be: eQTL of multiple tissues; or eQTLs, splicing QTLs and other types of QTLs in a single tissue. In a more complex setting, multiple types of QTL data from multiple tissues/cell types. Each group is defined by its “type” (kind of molecular traits), and “context” (tissue, cell type, condition, etc.). So, we specify the `type` and `context` arguments for each weight file. +We use the `preprocess_weight()` function to harmonize the PredictDB/FUSION prediction models and LD reference. -In this example, we will use liver and subcutaneous adipose gene expression models trained on GTEx v8 in the PredictDB format. -We can download both the prediction models (.db) and the covariances between variants in the prediction models (.txt.gz) from the link below. The covariances can optionally be used for computing LD. +In this example, we will use liver gene expression models trained on GTEx v8 in the PredictDB format. We can download both the prediction models (.db) and the covariances between variants in the prediction models (.txt.gz) from the link below. The covariances can optionally be used for computing LD. ```{bash, eval=FALSE} -# download the weights files wget https://zenodo.org/record/3518299/files/mashr_eqtl.tar tar -xvf mashr_eqtl.tar ``` -We get a list of processed weights for each weight file, and then we simply concatenate them to get a list of all processed weights from different types or contexts. ```{r preprocess_weights, eval=FALSE} -weight_liver_file <- system.file("extdata/sample_data", -"expression_Liver.db", package = "ctwas") -weights_liver <- preprocess_weights(weight_liver_file, +weight_file <- system.file("extdata/sample_data", "expression_Liver.db", package = "ctwas") +weights <- preprocess_weights(weight_file, region_info, z_snp$id, - snp_info, - type = "expression", - context = "liver", + snp_info, weight_format = "PredictDB", ncore = 6, drop_strand_ambig = TRUE, scale_predictdb_weights = TRUE, load_predictdb_LD = TRUE, filter_protein_coding_genes = TRUE) - -weight_adipose_file <- system.file("extdata/sample_data", "expression_Adipose_Subcutaneous.db", package = "ctwas") -weights_adipose <- preprocess_weights(weight_adipose_file, - region_info, - z_snp$id, - snp_info, - type = "expression", - context = "adipose", - weight_format = "PredictDB", - ncore = 6, - drop_strand_ambig = TRUE, - scale_predictdb_weights = TRUE, - load_predictdb_LD = TRUE, - filter_protein_coding_genes = TRUE) - -# concatenate weights from different types or contexts -weights <- c(weights_liver, weights_adipose) -# saveRDS(weights, file.path(outputdir, paste0(outname, ".preprocessed.weights.RDS"))) ``` +We get a list of processed weights for each weight file, and then we simply concatenate them to get a list of all processed weights from different types or contexts. -## Computing z-scores of molecular traits +We specify PredictDB or FUSION format in `weight_format` and the specific method in FUSION by `method_FUSION`. +By default, it will drop strand ambiguous variants and limit to protein coding genes. +We set `scale_predictdb_weights=TRUE` for PredictDB weights because PredictDB assumes the unstandardized genotypes (see the discussion in "Region info" section) +For PredictDB models, we set `load_predictdb_LD=TRUE` to load pre-computed correlations (`.txt.gz` file) between weight variants. +For FUSION models, we set `load_predictdb_LD=FALSE` to calculate correlations between weight variants with LD reference. +We could set `ncore` to parallelize the computation. -After we have done data preprocessing, we compute z-scores of molecular traits. This step basically performs Transcriptome-wide association studies (TWAS) on each molecular trait with a prediction model. The underlying calculation is based on [S-PrediXcan][S-PrediXcan]. +## Computing Z-scores of genes -The `compute_gene_z` function computes z-scores for the molecular traits using preprocessed SNP z-scores (`z_snp`) and the preprocessed weights. We could set `ncore` to parallelize the computation. +After we have done data preprocessing, we compute Z-scores of genes. This step basically performs Transcriptome-wide association studies (TWAS) on each gene with a prediction model. The underlying calculation is based on [S-PrediXcan][S-PrediXcan]. -*Note: the software uses the term “gene”, but they could be any molecular trait.* +The `compute_gene_z` function computes Z-scores for the genes using preprocessed SNP Z-scores (`z_snp`) and the preprocessed weights (`weights`). ```{r compute_gene_z, eval=FALSE} -z_gene <- compute_gene_z(z_snp, weights, ncore = 6) -# saveRDS(z_gene, file = file.path(outputdir, paste0(outname, ".z_gene.RDS"))) +z_gene <- compute_gene_z(z_snp, weights) ``` diff --git a/vignettes/summarizing_ctwas_results.Rmd b/vignettes/summarizing_ctwas_results.Rmd index 852e8618..d03eecd1 100644 --- a/vignettes/summarizing_ctwas_results.Rmd +++ b/vignettes/summarizing_ctwas_results.Rmd @@ -21,7 +21,7 @@ knitr::opts_chunk$set(message = FALSE, ``` -In this document, we will show how to summarize and visualize cTWAS results. +In this tutorial, we will show how to summarize and visualize cTWAS results. Load the packages. ```{r load_packages, message=FALSE} @@ -29,7 +29,7 @@ library(ctwas) ``` -We also need gene annotations from the Ensembl database. We use `EnsDb.Hsapiens.v86` for the example data, please choose the Ensembl database for your specific data. +We need the Ensembl database `EnsDb.Hsapiens.v86` here in this tutorial, please choose the version of Ensembl database for your specific data. ```{r load_ens_db, message=FALSE} library(EnsDb.Hsapiens.v86) ``` @@ -37,18 +37,17 @@ library(EnsDb.Hsapiens.v86) Load the cTWAS results ```{r load_ctwas_res} -# GWAS sample size gwas_n <- 343621 -# Load weights + weights <- readRDS(system.file("extdata/sample_data", "LDL_example.preprocessed.weights.RDS", package = "ctwas")) -# Load region_info, snp_info and LD_info region_info <- readRDS(system.file("extdata/sample_data", "LDL_example.region_info.RDS", package = "ctwas")) + snp_info <- readRDS(system.file("extdata/sample_data", "LDL_example.snp_info.RDS", package = "ctwas")) + LD_info <- readRDS(system.file("extdata/sample_data", "LDL_example.LD_info.RDS", package = "ctwas")) -# Load cTWAS results ctwas_res <- readRDS(system.file("extdata/sample_data", "LDL_example.ctwas_sumstats_res.RDS", package = "ctwas")) param <- ctwas_res$param finemap_res <- ctwas_res$finemap_res @@ -57,37 +56,56 @@ z_gene <- ctwas_res$z_gene region_data <- ctwas_res$region_data screened_region_data <- ctwas_res$screened_region_data -# directory with computed correlation matrices cor_dir <- system.file("extdata/sample_data", "cor_matrix", package = "ctwas") ``` ### Assessing parameters and computing PVE -We could use the `summarize_param()` function to assess the convergence of the estimated parameters and -to compute the proportion of variance explained (PVE) by variants and genes. +We could use the `summarize_param()` function to assess the convergence of the estimated parameters and to compute the proportion of variance explained (PVE) by genes and variants. ```{r summarize_param} ctwas_parameters <- summarize_param(param, gwas_n) ``` +Estimated prior inclusion probability: + ```{r} -# Estimated prior inclusion probability ctwas_parameters$group_prior +``` -# Estimated prior effect size: + +Estimated prior effect size for genes and variants: + +```{r} ctwas_parameters$group_prior_var +``` -# Estimated enrichment of molecular traits over variants: + +Estimated enrichment of genes over variants: + +```{r} ctwas_parameters$enrichment +``` -# PVE explained by molecular traits and variants: + +PVE explained by genes and variants: + +```{r} ctwas_parameters$group_pve +``` -# Total heritability (sum of PVE) + +Total heritability (sum of PVE): + +```{r} ctwas_parameters$total_pve +``` -# Attributable heritability + +Attributable heritability to variants in genes: + +```{r} ctwas_parameters$attributable_pve ``` @@ -98,80 +116,50 @@ make_convergence_plots(param, gwas_n) ``` +These plots show the estimated prior inclusion probability, and prior effect size, enrichment and PVE over the iterations of parameter estimation. In particular, genes have higher prior inclusion probability, enrichment and PVE than variants. + ## Summarizing cTWAS results +To interpret and visualize the results, we first need to add gene annotations (gene names and gene types) to the fine-mapping result. -To interpret and visualize the results, we first need to add gene annotations -(gene names and gene types) to the finemapping result. +We need a user defined data frame `gene_annot` that should contain the following columns: 'chrom', 'start', 'end', 'gene_id', 'gene_name' and 'gene_type' for each gene. -We need a user defined data frame `gene_annot` that should contain the following columns: -'chrom', 'start', 'end', 'gene_id', 'gene_name' and 'gene_type' for each gene or molecular trait. +If you only used gene expression data, you can use the following function to obtain the gene annotations from the Ensembl database. We use `EnsDb.Hsapiens.v86` for the example data, please choose the Ensembl database for your specific data. -If you only used gene expression data, you can use the following function -to obtain the gene annotations from the Ensembl database. -```{r get_gene_annot, message=FALSE} +```{r get_gene_annot_from_ens_db, message=FALSE} ens_db <- EnsDb.Hsapiens.v86 -finemap_gene_res <- finemap_res[finemap_res$type!="SNP",] +finemap_gene_res <- subset(finemap_res, group == "gene") gene_ids <- unique(sapply(strsplit(finemap_gene_res$id, split = "[|]"), "[[", 1)) -gene_annot <- get_gene_annot_from_ens_db(ens_db = ens_db, gene_ids) +gene_annot <- get_gene_annot_from_ens_db(ens_db, gene_ids) ``` -We then use `anno_finemap_res()` function to add gene annotations to the finemapping result, -and use gene mid-points to represent gene positions. +We then use the `anno_finemap_res()` function to add gene annotations to the finemapping result, and use gene mid-points to represent gene positions. ```{r anno_finemap_res} -# add gene annotations and use gene mid-points to represent gene positions -finemap_res <- anno_finemap_res(finemap_res, - snp_info, - gene_annot, +finemap_res <- anno_finemap_res(finemap_res, snp_info, gene_annot, use_gene_pos = "mid") ``` ## Interpreting cTWAS results -We list all molecular traits with PIP > 0.8, which is the threshold we used in the paper. -```{r select_sig_res} -# select gene PIP > 0.8 -sig_finemap_res <- finemap_res[finemap_res$type!="SNP" & finemap_res$susie_pip > 0.8,] -head(sig_finemap_res[order(-sig_finemap_res$susie_pip),]) -``` - - -When we have multiple contexts (e.g. tissues) or types (e.g. expression, splicing), it is useful to evaluate the combined PIP for each gene by contexts, types, or groups. -We could further limit results in credible sets and/or protein coding genes. -```{r combine_pips} -# combine gene PIPs by context -combined_pip_by_context <- combine_gene_pips(finemap_res, - by = "context", - filter_protein_coding_genes = TRUE) -head(combined_pip_by_context) - -# combine gene PIPs by type -combined_pip_by_type <- combine_gene_pips(finemap_res, - by = "type", - filter_protein_coding_genes = TRUE) +We list all genes with PIP > 0.8, which is the threshold we used in the paper. -# combine gene PIPs by group -combined_pip_by_group <- combine_gene_pips(finemap_res, - by = "group", - filter_protein_coding_genes = TRUE) - -# select combined PIP > 0.8 -sig_gene_pips_df <- combined_pip_by_context[combined_pip_by_context$combined_pip > 0.8, ] -head(sig_gene_pips_df) +```{r interpret_res} +sig_finemap_res <- subset(finemap_res, group == "gene" & susie_pip > 0.8) +sig_finemap_res[order(sig_finemap_res$susie_pip,decreasing = TRUE),] ``` +In this example, HPR is the only gene that satisfies this threshold. ## Visualizing cTWAS results -We could make locus plots for regions of interest. -We could zoom in the region by specifying the `locus_range`. +We could make locus plots for regions of interest, and zoom in the region by specifying the `locus_range`. -If we want to show the correlations with the focus gene, we need correlation matrices of the region. -Because we have computed and saved correlation matrices during fine-mapping step, +If we want to show the correlations with the focus gene, we need correlation matrices of the region. Because we have computed and saved correlation matrices during fine-mapping step, we could simply load those precomputed correlation matrices. + ```{r get_region_cor_res} cor_res <- get_region_cor(region_id = "16_71020125_72901251", cor_dir = cor_dir) ``` @@ -189,3 +177,6 @@ make_locusplot(finemap_res, legend.position = "top") ``` + +The locus plot shows several tracks in the example region in chromosome 6, including GWAS -log10(p-value), PIPs, QTLs of the focus gene (by default, the focus gene is the gene with the highest PIP), and the gene track. In this example, we zoomed in a locus (71.6 - 72.4 Mb) around the HPR gene, which has the highest PIP (PIP = 1), and the red dotted line shows the PIP cutoff of 0.8. +