From 3e5273114446c2b57f4f52f9a9cd919335113c86 Mon Sep 17 00:00:00 2001 From: Weissbrod Date: Tue, 21 Jul 2020 15:51:22 -0400 Subject: [PATCH] added support for genome-wide fine-mapping --- aggregate_finemapper_results.py | 85 +++++++++++++++++++++++++++ create_finemapper_jobs.py | 101 ++++++++++++++++++++++++++++++++ polyfun_utils.py | 30 +++++++--- ukb_regions.tsv.gz | Bin 0 -> 19390 bytes 4 files changed, 209 insertions(+), 7 deletions(-) create mode 100644 aggregate_finemapper_results.py create mode 100644 create_finemapper_jobs.py create mode 100644 ukb_regions.tsv.gz diff --git a/aggregate_finemapper_results.py b/aggregate_finemapper_results.py new file mode 100644 index 0000000..397fe1c --- /dev/null +++ b/aggregate_finemapper_results.py @@ -0,0 +1,85 @@ +import numpy as np; np.set_printoptions(precision=4, linewidth=200) +import pandas as pd; pd.set_option('display.width', 200) +import os +import logging +from tqdm import tqdm +from polyfun import configure_logger, check_package_versions +from polyfun_utils import set_snpid_index +from pyarrow import ArrowIOError +from pyarrow.lib import ArrowInvalid +from polyfun_utils import DEFAULT_REGIONS_FILE + + +def main(args): + + #read sumstats file + try: + df_sumstats = pd.read_parquet(args.sumstats) + except (ArrowIOError, ArrowInvalid): + df_sumstats = pd.read_table(args.sumstats, delim_whitespace=True) + + #read regions file + df_regions = pd.read_table(args.regions_file) + df_regions = df_regions.loc[df_regions.apply(lambda r: np.any((df_sumstats['CHR']==r['CHR']) & (df_sumstats['BP'].between(r['START'], r['END']))), axis=1)] + + #aggregate outputs + df_sumstats_list = [] + logging.info('Aggregating results...') + for _, r in tqdm(df_regions.iterrows()): + chr_num, start, end, url_prefix = r['CHR'], r['START'], r['END'], r['URL_PREFIX'] + output_file_r = '%s.chr%s.%s_%s.gz'%(args.out_prefix, chr_num, start, end) + if not os.path.exists(output_file_r): + err_msg = 'output file for chromosome %d bp %d-%d doesn\'t exist'%(chr_num, start, end) + if args.allow_missing_jobs: + logging.warning(err_msg) + continue + else: + raise IOError(err_msg) + df_sumstats_r = pd.read_table(output_file_r) + + #mark distance from center + middle = (start+end)//2 + df_sumstats_r['DISTANCE_FROM_CENTER'] = np.abs(df_sumstats_r['BP'] - middle) + df_sumstats_list.append(df_sumstats_r) + + #keep only the most central result for each SNP + df_sumstats = pd.concat(df_sumstats_list, axis=0) + df_sumstats.sort_values('DISTANCE_FROM_CENTER', inplace=True, ascending=True) + df_sumstats = set_snpid_index(df_sumstats, allow_duplicates=True) + df_sumstats = df_sumstats.loc[~df_sumstats.index.duplicated(keep='first')] + del df_sumstats['DISTANCE_FROM_CENTER'] + df_sumstats.sort_values(['CHR', 'BP'], inplace=True, ascending=True) + + #write output file + df_sumstats.to_csv(args.out, sep='\t', index=False) + logging.info('Wrote aggregated results to %s'%(args.out)) + + + +if __name__ == '__main__': + import argparse + parser = argparse.ArgumentParser() + + #general parameters + parser.add_argument('--sumstats', required=True, help='Name of sumstats file') + parser.add_argument('--out-prefix', required=True, help='prefix of output files') + parser.add_argument('--out', required=True, help='name of the aggregated output files') + parser.add_argument('--allow-missing-jobs', default=False, action='store_true', help='whether to allow missing jobs') + parser.add_argument('--regions-file', default=DEFAULT_REGIONS_FILE, help='name of file of regions and their URLs') + + #check package versions + check_package_versions() + + #extract args + args = parser.parse_args() + + #check that the output directory exists + if len(os.path.dirname(args.out))>0 and not os.path.exists(os.path.dirname(args.out)): + raise ValueError('output directory %s doesn\'t exist'%(os.path.dirname(args.out))) + + #configure logger + configure_logger(args.out_prefix) + + #invoke main function + main(args) + \ No newline at end of file diff --git a/create_finemapper_jobs.py b/create_finemapper_jobs.py new file mode 100644 index 0000000..58bc41a --- /dev/null +++ b/create_finemapper_jobs.py @@ -0,0 +1,101 @@ +import numpy as np; np.set_printoptions(precision=4, linewidth=200) +import pandas as pd; pd.set_option('display.width', 200) +import os +import logging +from polyfun import configure_logger, check_package_versions +from pyarrow import ArrowIOError +from pyarrow.lib import ArrowInvalid +from polyfun_utils import DEFAULT_REGIONS_FILE + + +FINEMAPPER_SCRIPT = 'finemapper.py' + + +def create_finemapper_cmd(args, chr_num, start, end, url_prefix): + + output_file = '%s.chr%s.%s_%s.gz'%(args.out_prefix, chr_num, start, end) + cmd = '%s %s --chr %s --start %s --end %s --out %s'%(args.python, FINEMAPPER_SCRIPT, chr_num, start, end, output_file) + if args.max_num_causal>1: + cmd += ' --ld %s'%(url_prefix) + + #add command line arguments + for key, value in vars(args).items(): + if key in ['python', 'regions_file', 'out_prefix', 'jobs_file']: continue + key = key.replace('_', '-') + if type(value)==bool: + if value: + cmd += ' --%s'%(key) + elif value is not None: + cmd += ' --%s %s'%(key, value) + + return cmd + + +def main(args): + + #read sumstats file + try: + df_sumstats = pd.read_parquet(args.sumstats) + except (ArrowIOError, ArrowInvalid): + df_sumstats = pd.read_table(args.sumstats, delim_whitespace=True) + + #read regions file + df_regions = pd.read_table(args.regions_file) + df_regions = df_regions.loc[df_regions.apply(lambda r: np.any((df_sumstats['CHR']==r['CHR']) & (df_sumstats['BP'].between(r['START'], r['END']))), axis=1)] + + #create jobs + with open(args.jobs_file, 'w') as f: + for _, r in df_regions.iterrows(): + chr_num, start, end, url_prefix = r['CHR'], r['START'], r['END'], r['URL_PREFIX'] + cmd = create_finemapper_cmd(args, chr_num, start, end, url_prefix) + f.write(cmd + '\n') + + logging.info('Wrote fine-mapping commands to %s'%(args.jobs_file)) + + + +if __name__ == '__main__': + import argparse + parser = argparse.ArgumentParser() + + #general parameters + parser.add_argument('--method', required=True, help='Fine-mapping method (currently susie and finemap are supported)') + parser.add_argument('--sumstats', required=True, help='Name of sumstats file') + parser.add_argument('--n', required=True, type=int, help='Sample size') + + #LDstore related parameters + parser.add_argument('--finemap-exe', default=None, help='Path to FINEMAP v1.4 executable file') + parser.add_argument('--memory', type=int, default=1, help='Maximum amount of memory in GB to allocate to LDStore') + parser.add_argument('--threads', type=int, default=None, help='The number of CPU cores LDstore will use (if not specified, LDstore will use the max number of CPU cores available') + + parser.add_argument('--max-num-causal', required=True, type=int, help='Number of causal SNPs') + parser.add_argument('--non-funct', action='store_true', default=False, help='Perform non-functionally informed fine-mapping') + parser.add_argument('--hess', action='store_true', default=False, help='If specified, estimate causal effect variance via HESS') + parser.add_argument('--verbose', action='store_true', default=False, help='If specified, show verbose output') + parser.add_argument('--allow-missing', default=False, action='store_true', help='If specified, SNPs with sumstats that are not \ + found in the LD panel will be omitted. This is not recommended, because the omitted SNPs may be causal,\ + which could lead to false positive results') + + parser.add_argument('--regions-file', default=DEFAULT_REGIONS_FILE, help='name of file of regions and their URLs') + parser.add_argument('--python', default='python3', help='python3 executable') + parser.add_argument('--out-prefix', required=True, help='prefix of the output files') + parser.add_argument('--jobs-file', required=True, help='name of file with fine-mapping commands') + + #check package versions + check_package_versions() + + #extract args + args = parser.parse_args() + + #check that the output directory exists + if len(os.path.dirname(args.out_prefix))>0 and not os.path.exists(os.path.dirname(args.out_prefix)): + raise ValueError('output directory %s doesn\'t exist'%(os.path.dirname(args.out_prefix))) + if len(os.path.dirname(args.jobs_file))>0 and not os.path.exists(os.path.dirname(args.jobs_file)): + raise ValueError('output directory %s doesn\'t exist'%(os.path.dirname(args.jobs_file))) + + #configure logger + configure_logger(args.out_prefix) + + #invoke main function + main(args) + \ No newline at end of file diff --git a/polyfun_utils.py b/polyfun_utils.py index d1c39e8..45312dc 100644 --- a/polyfun_utils.py +++ b/polyfun_utils.py @@ -7,9 +7,24 @@ SNP_COLUMNS = ['CHR', 'SNP', 'BP', 'A1', 'A2'] +LONG_RANGE_LD_REGIONS = [] +LONG_RANGE_LD_REGIONS.append({'chr':6, 'start':25500000, 'end':33500000}) +LONG_RANGE_LD_REGIONS.append({'chr':8, 'start':8000000, 'end':12000000}) +LONG_RANGE_LD_REGIONS.append({'chr':11, 'start':46000000, 'end':57000000}) +DEFAULT_REGIONS_FILE = 'ukb_regions.tsv.gz' +class TqdmUpTo(tqdm): + """ + taken from: https://github.com/tqdm/tqdm/blob/master/examples/tqdm_wget.py + """ + + def update_to(self, b=1, bsize=1, tsize=None): + if tsize is not None: self.total = tsize + self.update(b * bsize - self.n) + + ''' Logger class (for compatability with LDSC code)''' class Logger(object): @@ -39,7 +54,7 @@ def check_package_versions(): raise ValueError('\n\nPlease install the python package pandas_plink (using either "pip install pandas-plink" or "conda install -c conda-forge pandas-plink")\n\n') -def set_snpid_index(df, copy=False): +def set_snpid_index(df, copy=False, allow_duplicates=False): if copy: df = df.copy() is_indel = (df['A1'].str.len()>1) | (df['A2'].str.len()>1) @@ -53,12 +68,13 @@ def set_snpid_index(df, copy=False): df.drop(columns=['A1_first', 'A1s', 'A2s'], inplace=True) #check for duplicate SNPs - is_duplicate_snp = df.index.duplicated() - if np.any(is_duplicate_snp): - df_dup_snps = df.loc[is_duplicate_snp] - df_dup_snps = df_dup_snps.loc[~df_dup_snps.index.duplicated(), ['SNP', 'CHR', 'BP', 'A1', 'A2']] - error_msg = 'Duplicate SNPs were found in the input data:\n%s'%(df_dup_snps) - raise ValueError(error_msg) + if not allow_duplicates: + is_duplicate_snp = df.index.duplicated() + if np.any(is_duplicate_snp): + df_dup_snps = df.loc[is_duplicate_snp] + df_dup_snps = df_dup_snps.loc[~df_dup_snps.index.duplicated(), ['SNP', 'CHR', 'BP', 'A1', 'A2']] + error_msg = 'Duplicate SNPs were found in the input data:\n%s'%(df_dup_snps) + raise ValueError(error_msg) return df diff --git a/ukb_regions.tsv.gz b/ukb_regions.tsv.gz new file mode 100644 index 0000000000000000000000000000000000000000..54ba875ee9c7a0eef576adf260f1411cee8631ab GIT binary patch literal 19390 zcmeI2dpMM9+y7OhEh-jSwrxdAQC4;(QA$d(Tb9vOqGT7t5HUqXvMVfV_(V{57Czg z9@nfbwh8}%zW7YxCjnx^DEL0rl2yHqGTy~{E}>DN?#wM2!Od8bd2exuG9m~-n0s)u zXry-!j$cXE7tpzaaU%f>ne|3u_;sP{&@80gTOf$T%|F5k0?2$$_q-sEj1(i25x9AN z_xvk_$1)PUQbzjYaJ=)#Oaz(7;UZjQmW~rt8qG7w^Y6WpIr98$rZR3)mdx+QA>szZ zGdiDjMrglwDb9VxA@%Oh5P66DxwyA_L-PE6VI;jdT3W7CZnyZmJ!Z%~%!K;{Po-{P z5RHlTi3wGH-KD^HZ`LBd>k>k_7Tn%OZ{Q25iPZSWOWD(&0zS0Ftmi>F=1#A zLr-Cdrlrn|+yEgX;S{s`d~c=w)ZQCJRG63`_qe| zTUx# zyemuyvPh`bAsFM{8czjY@bf;n+h1_Q3p6w!mfr- zj=T)NzP}0cg{#lER@ISP^01I$1sVTL=!rl72wuoZVabpGD&+Cr^DHj%#D+rL*;sd4 z>W2|L*tH=TKKpAw4!pkHPhK4T1d_0~*vb;NMW-;k6Bg!bdwS_&^@&bF<&|X!-Jwg%AkF&6BCWJ%_ zA<+=>yOXabZ+w$Eoey-qzpujo!~ktw_~R`RJjutlF*l&N92i0NK=Z7jpT z$ic|IXl)%0yOOz={Znx$;V0SPO8pX7p~N0{u`*_?WMrz!#9uMallGWe^r~7XpP>FR zGg%=PhUNePvr=xk1Nz|dFcg(Hj{W;|BnSz!M65Kh8&BB3hX!mg#_P^rPzcYaaCEZUd%syAmBDK7I9Yyhife+6hu-IY`vx`48BF zDpy5pO9_$rgmW(Vt^RO^Cz|3+|1x#YzkP3aH2|&~(P94@To(c^>roeEr}6Za=*KH) zMb2;M%FpuGm!xBb>=96)!PqHad7eQVYsU@=$yJctCnSdm$zbOQdLhu^&Uv|pBvebGE>NKHlD--HGggOe;%?fq#P`7Iw>NcZ6 zL99Ml(8=n*Y!(~q{iq-6o(pwNP?tIj{q8{><}1|c{pR=nEA$H$`VB(eDn6|99n|T~ zL0uly6;D9jF`?fV=$8Zhuq_{818hQIYYSlm8s=aF9Lb=y?z8 zV;hCK4yel$>Pn!_jR$q!P`7Rx>IQ_38-WFN2@7Jtf;xDxplV^`reQ(XU_s|cVL>a% zu%OQnyeagpaJ;D;#&IUpU7m-!dST;cp)N^yA>&Y|Ez}J_ogdH*XSBeAmh)jjxv(JX zIattf9NdvfSWuAwc8Dzm^ZgCUr*cc}=9i!4pL(Lg=F|hmd`1VF!lY*d#?0Olz$mrm z>99wGF9L@5qbYal=jjiC(S4i)7~Fios7~($#dTZr*nw9JNrt@h`My?AY7*e=a1LEWK^|1H=1nB8MLGKCjDr zdb4@Qo(HIcL~lak(O;f=A1(LW)*Z5qja^cg=wC_Jv(|KLtk&sodwAvRLsV~K2&VDD z@!5gn<$c?DQ;A12$KPi9F1MLlt)6_K|GtTfVjS1e=9fb*@ZI-mBqX8nWd6*a|K1MV zueh;Eptba;PeDJYi%_7650k)_RJ!*AGznkiK)E-Gsjc`#uM;j~gv*2?uqRa-a3kzJ zfg{7V5;#^~gd0fc1dg9<9f4!ON8rdIbOOh!LM(6$_y`@HN`PZ!AqhBg3!H$XQwbab zdk%0MDZu|<$G=&QKs)JlPa$)k5OzF5%M_p^=E;HyM86N2e1yFEHR7Ye?k1~ z5w;~lv=*Rc1?WffWX*Z<@mZYH3{LSYQZbI`4

ih<_8pZbaBA2rUVrX(2=ngy<|l z+Y8W?c`{+1918@q3xYpEmzu#&&*9R+XnYo zK#TLq1TbWCI42%XaRRBBLi7g_f6l*2z{i45^D6jhUP>f@yRhIRxb;?DgI_fQeg$u& zkspe{ExR%Td>g@fpLGH^XQv`~NB*7-Zo?zazrB&th07YkWfwg7(7NvcZ}5BE!Cjrv z48K3{9N_)?CKY@de0nEzdjcfSP2#nq-reQlphp zO0HpIuHjCz%W7tqPsfg*h#h~S(4Vc)e+>xq8kb-aBYYPGe}XQBgPrzE6lbYPf?;a5 zQfjkO%E~o-mTS1*?D7S(%U5H^TVuyX6#KKm#OXJh7&I>NON_uKiVvzua?~VSt&}BK@oRTEKG1gD>`F3sA(=OMBsF>@rF5t#b*O6< z+o~7ao^>s>aV>n^WS74n_%n2=4~)iJy9O{C_waYYXn^^-Kr;6vnYVf*b$TRab*Ptg zs2dmCUM#i^b1nSfS}6L~F8{6F@_YD`_wj3YJ3jbvL9oyz&)~?R5mohi6B|oOD`Tac zprFVV<5C~ro7luhJkC(c2{SaEoRE^~Q!RKD@mNjLYPqD9l~T@HrJNi?)6IsaH!eq- zUyfYyRVry*N@hS+ydosw4DB(2UU<(yQ?d1+|+U_tO_=u$5j zjlHI|U^Idvqrhl@`59DI=cuZ;n%H!j*knaKeird~zoeC$q}4U0oHnJL#d}R3>@^h& ziZl<7TrnY)G_@f3w!3tZE`JA0i`;#9V_Sy&ouzK3oz@%MRALSXUo@4^TjORKb2xnW z22TDOxA;SOdv!L%v_w zt%>;TGo+gjJyx$S*4gdiejx>)z0J|&mB;F19rE%ElK+N|op-;`Ms3w~G)W-ctONV< zOlNm1we=Z3`x@96&(*8Jxu^SuEPS?`qsbD_)yL@a@)z7MbW&UQJDR*F-CWu!e~+{v z`Cr^)L#<}Ot^?~;Fft2_$Li|{)D$_}osc2GHd7#eVW^@41 ztXj_n8aZu-%n1bY1K+Lw$}`gEnijPyWX#9LGthKShGv(gEv_i&2UYprRiW_zGx3- zy~%PQuQ(e^7g{)1}oo*_fsoE>W?n0=wR?)upxMO;$MAN~fiYe(U`>3zC0f z$DX`h9%r!CP~~`_q2A|3J};RK%OxY%DOny5+@CfY==*X>?{!MTIwh&Ch8bz2S-PK9 zbw6KwM>2UwT9WG-k>e?Th%RxEE)h*qj3Oz@c^XK28tm-Es&!&dKjWNu#(8mpo_&G7 zAo$~VspkS+)s0>QMx%`s4n`x-Q~WSp;uu{bo}_q>q^RO)u-nsMe<$`rC-&+yPU|y{ zh#NiIjlQUj1in8k$@7eO_-7MfzTb5}pvmJvO2?J6uIpYl4XoRBU&{5z%#@DnJBw#u zH*L9|`gW`9kE%@``?QLWoN={2?V2#B*SIY-*7u~5nO<#+)$vvF_gBQhY01QCX(>6YD>C~`In=YkIPwqm0Q=R^>9FIOTDwjTW5=SO3FP-ib|Z8PMnsp zlJmtsoM08N(`%4UtIEMi%hZ1S`K{pQ&OZWZnWIGhifGrDOM1rEE8@s(dG z+k!(XK@j&j?kPA}I3ko$Gbe;{=`j&;o-6(paW?$Wj%W!SIuI?TpqTI^J)c;6?XnIf)Xj2ZV_vl6O@Y#aWhwDhzT%@ zKJ*yL8fwh57C(N4(+vD&rbX<;Z3H!1GzLIQ2mtk=ao!=GduD0fsO|!o_Xg(8g1q}M z?<>gL2lEDjVT=G6wgE$ZJ}}IOhI7F1X9O5d0>fHpI1CI6`2a+RARf(VEty2;vomjX z55d)Lu=?FRz`ug{8Gw(4_;G;$0Ri{{fDaM?{5ycx;RAdD#Loe|3j**{0M8Zx5Ep_3 zG}#UsmWc~V2q~NgtY*lX11vmbO+gj{SYwcd0M-y>2>`1fvILNY02TsRi$z$A2e&fmnx{J{a6qq*$31&J4scDlU`*rYu=Sv^pZso0Ybbv9)-X32@WLyr%D4@UK7!4?vc+>&qJ6!ebU) z_efuv(L^m52_rR&SzctLM9w&}!nT(BYhHbU$J8Wr>C5gPq?Y$Ic%`pY31p*2^Gd0& z1hgNkOHNIZN9KWWaTySX0ih}i2wwrA3IPavfl!VHgg!tR4}_b_o5iFW{WSxhqJaMZ z@OL8se+=+9qybC?z~s2QTtL$g(^kBz9q(Q?fu-FFJ*{1-T4cMoNfIS7+ zWsn^V*gdfBYQP=`Y;81a@rpV9J@07_*SofR_tJ=a;9>l`+%6*a9@!;`D! zrSybE_l8)E?mINLZ%d$ZZlLnI=W%DA$8BR!b}}d+SB~opa$~2p8hCOPUP@m`G&96v zd|w?%zpTs*R8Dvv=kq-72!o=_plsqgD{!6fepBgb;YXTlben57u}&+nPT$4sjK%D1 z(37L+rSvI8GZid+CF*=6w&;m7b;T3X!hQZA=_;(#_b@x-F+1Pt$<^zn3@AkRDGV=R zod#V{))n8VC(e7JC)c3YF5X;oxVdKEC5c165>wYOJMUtCl}ro2lNKJ*^G!)+_J#Fy zONl;`5|v#RX1y$IrD)AI(HgO^9_g^2#m1Lb7+<=kqI*N7Y`>I(i4;n8S*+?Zf|Mvl zO4NQ?7-m`6QPCPb(VER+JqlqxYmF~S7+<=jq8p)7W-65i(l;z?R$ZpKMwBHbn*2?r z^BUG98|n3SS;m&lH>_cW`o8YBjMbe^*rPGaInfm-u z(#5S5Xlk+HYO(Ip6k}(u}w3vOA9a5jVRI;Ul4qwUHZvK!pTSC*ItE+UWL=OdVaNfV4JTz))jZ~iItRL zsc1BhW%nw80oZ4g7G9)FiHO30379$IBLPN3S-siVE{$pOJ4NLi>>D@O`yM1D9whWd zp|;&dHOQ2o{I2}2hVO&zzMJj)%pLm73H3=&>yuD-G+gdzEC{~QF4b(XAKG9ac98Jl zAVKstD*ra>$am%ACU-O%HrR)%`EGXTYm_`l*mN5;n_ZuD?v93$bot%wzWos!?C%~V z%%A8p2P2>yg{qx2GoLbxV#g&8$4S$aw{<9QBN&~0YIG8X%R7h5D-#@U5FGa5-?_`b z(>u3uGjHR9;G65xV&29Fyp3X0X6DmoD~96|N8+SAl=CI|cQ)`gPG%dO1pB;i!mQJq z@nnD;=`Ph`SZ*#Uy6>cPZld51CB$5Paq< zdbI_bO&JaS;U`=ar<*!?zMWrv>TE9O`(U@2YW`SBr#HQH>ZJT9E@$)Z&g$~y=E2Z_ z9jWi8I8`=>rg_19JNDG->d^ej&^{lc`@26KqLvjnCfJt#mnAZpTQv5?6wd5+7;{6_K5B@vWUB8w zW_Ly2e$4G@@B6YS%LJLsI@>Gvie`2bdf%U`edO~el}zxZs&Fgc8OY=WtjP&bH|1_H z<=R9t4n{KarK)aARavVB$*KlDIxwMmU?TR$pwo>(#g*?WSG+U$E++tJG)=kMrd(7c z!zGeYE>-nLstThTl&l)W0s@nXW+1pZm<|N1-WkZ|1g!t>7U-cZlB-jiXU|ql)cQ5c zzJIXeA*MP-Gej%>Ypvgx+h5XlJUoiFU|rF2nmen{y8We1_IOB&i+$|$AL!6(*F zORb;I(jHOO9x-<5GC9{3-dGkK zdVpZpNDUCUc7?wwi)bnnPxg^a^^w$MDrqy7P<4CH)$J{d4z7<5zO;TienId*(xnXE zqLKOUdXB5VNZ;@x-R7#n!K()OA^}k%0VfxYNG}@6Tims6an}*CGHcPYbt^ADT5)MB z$d>6jzWqfy(4ejwxLh?T7YS$-3GiJslDKH34+xfYH2}fVvIHPlb!n??&++v=$2DK1 zYya5=1tK&u37yMl3*F&T78u%%Mc8z*z?I9e54^&PK1q?33LO&cawg zED{9H4AH-#ssxNqT+AsPCM5X(oA}${+h|$}SmW$F2j964C}yqpxp)M+wLMGy{80z1nZ* z=gV5ZI%RRJDn)+k?8Twl4`s5}2|FIvRi)g!qNOuZ`=Rf4zR!+_ThJCU!CE@AXD@c$ z&To*lW&%ldiv0B1i{pO~iT?ER{_%B<<(6v6nJ1b<4wTPd_Z?Z+*uOOqlX*hZC_9aJ z-FJ3*f7jN;H04q5F literal 0 HcmV?d00001