From 3b0644fc4319aca017461255b1f690314a30c40c Mon Sep 17 00:00:00 2001 From: AnttonLA Date: Mon, 20 Oct 2025 14:14:19 +0200 Subject: [PATCH] Add GwasTrack, which draws a scatterplot from a .gwas file --- examples/gwas.gwas | 13 +++ pygenometracks/readGwas.py | 80 ++++++++++++++++++ pygenometracks/tests/test_data/gwas.ini | 16 ++++ pygenometracks/tests/test_data/gwas_1.gwas | 12 +++ pygenometracks/tests/test_data/gwas_2.gwas | 13 +++ .../tests/test_data/master_gwas.png | Bin 0 -> 34493 bytes pygenometracks/tests/test_gwasTrack.py | 50 +++++++++++ pygenometracks/tracks/GwasTrack.py | 76 +++++++++++++++++ 8 files changed, 260 insertions(+) create mode 100644 examples/gwas.gwas create mode 100644 pygenometracks/readGwas.py create mode 100644 pygenometracks/tests/test_data/gwas.ini create mode 100644 pygenometracks/tests/test_data/gwas_1.gwas create mode 100644 pygenometracks/tests/test_data/gwas_2.gwas create mode 100644 pygenometracks/tests/test_data/master_gwas.png create mode 100644 pygenometracks/tests/test_gwasTrack.py create mode 100644 pygenometracks/tracks/GwasTrack.py diff --git a/examples/gwas.gwas b/examples/gwas.gwas new file mode 100644 index 00000000..ce7eb4c6 --- /dev/null +++ b/examples/gwas.gwas @@ -0,0 +1,13 @@ +CHR BP SNP P BETA SE MAF +X 3001200 rs100001 2.5e-09 0.62 0.09 0.21 +X 3003450 rs100002 7.8e-08 0.55 0.10 0.23 +X 3007890 rs100003 1.1e-06 0.41 0.08 0.27 +X 3009990 rs100004 4.3e-05 0.33 0.08 0.30 +X 3032100 rs100005 9.2e-10 0.71 0.11 0.18 +X 3045600 rs100006 3.4e-08 0.59 0.10 0.20 +X 3060200 rs100007 6.7e-07 0.48 0.09 0.25 +X 3078500 rs100008 1.9e-05 0.36 0.08 0.29 +X 3100400 rs100009 8.8e-06 -0.40 0.09 0.26 +X 3150000 rs100010 2.2e-07 -0.52 0.10 0.22 +X 3195000 rs100011 1.5e-09 -0.68 0.11 0.19 +X 3199800 rs100012 5.0e-08 -0.60 0.10 0.21 \ No newline at end of file diff --git a/pygenometracks/readGwas.py b/pygenometracks/readGwas.py new file mode 100644 index 00000000..4a2089fd --- /dev/null +++ b/pygenometracks/readGwas.py @@ -0,0 +1,80 @@ +# -*- coding: utf-8 -*- +import sys +import collections +from .utilities import to_string, InputError + + +class ReadGwas(object): + """ + Reads a GWAS file. The expected fields are: + chromosome, position, name, and pvalue. + + Example: + gwas = ReadGwas(open("file.gwas", 'r')) + for record in gwas: + print(record.chromosome, record.position, record.pvalue) + """ + + def __init__(self, file_handle, has_header=False): + """ + :param file_handle: file handle + """ + self.file_handle = file_handle + self.line_number = 0 + + # Define the fields for GWAS + self.fields = ['chromosome', 'position', 'name', 'pvalue'] + self.GwasRecord = collections.namedtuple('GwasRecord', self.fields) + + # Skip the header line if present + if has_header: + header = next(self.file_handle) + self.line_number += 1 + + def __iter__(self): + return self + + def get_no_comment_line(self): + """ + Skips comment lines starting with '#' or empty lines. + :return: a valid line + """ + line = next(self.file_handle) + line = to_string(line) + if line.startswith("#") or line.strip() == '': + line = self.get_no_comment_line() + + self.line_number += 1 + return line + + def __next__(self): + """ + :return: GwasRecord object + """ + line = self.get_no_comment_line() + return self.get_gwas_record(line) + + def get_gwas_record(self, gwas_line): + """ + Processes each line from a GWAS file and returns a namedtuple object. + + :param gwas_line: a single line from the GWAS file + :return: GwasRecord object + """ + line_data = gwas_line.strip() + line_data = to_string(line_data) + line_data = line_data.split("\t") + + if len(line_data) < 4: + raise InputError(f"Line {self.line_number} does not have 4 fields: {gwas_line}." + f"We expect at least 4 field, corresponding to: chromosome, position, name, pvalue.") + + try: + chromosome = line_data[0] + position = int(line_data[1]) + name = line_data[2] + pvalue = float(line_data[3]) + except ValueError as e: + raise InputError(f"Error parsing line {self.line_number}: {gwas_line}\n{e}") + + return self.GwasRecord(chromosome, position, name, pvalue) diff --git a/pygenometracks/tests/test_data/gwas.ini b/pygenometracks/tests/test_data/gwas.ini new file mode 100644 index 00000000..f0aa169e --- /dev/null +++ b/pygenometracks/tests/test_data/gwas.ini @@ -0,0 +1,16 @@ + +[gwas] +file = gwas_1.gwas +height = 4 +title = test_1 + +[spacer] + +[gwas_2] +file = gwas_2.gwas +file_has_header = True +height = 2 +title = test_2 +color = #50E3C2 + +[x-axis] diff --git a/pygenometracks/tests/test_data/gwas_1.gwas b/pygenometracks/tests/test_data/gwas_1.gwas new file mode 100644 index 00000000..dae8e6b2 --- /dev/null +++ b/pygenometracks/tests/test_data/gwas_1.gwas @@ -0,0 +1,12 @@ +X 3002145 rs7823451 4.7e-08 +X 3005892 rs6610293 1.3e-06 +X 3008734 rs9034128 9.2e-05 +X 3031022 rs4420981 6.1e-10 +X 3042567 rs2093847 2.8e-07 +X 3056789 rs7734012 5.5e-09 +X 3072341 rs1182736 3.9e-06 +X 3079880 rs5502918 7.4e-05 +X 3112450 rs3847291 1.8e-08 +X 3145670 rs6029184 9.9e-07 +X 3180230 rs2938471 3.2e-10 +X 3197650 rs8840123 6.3e-08 \ No newline at end of file diff --git a/pygenometracks/tests/test_data/gwas_2.gwas b/pygenometracks/tests/test_data/gwas_2.gwas new file mode 100644 index 00000000..ce7eb4c6 --- /dev/null +++ b/pygenometracks/tests/test_data/gwas_2.gwas @@ -0,0 +1,13 @@ +CHR BP SNP P BETA SE MAF +X 3001200 rs100001 2.5e-09 0.62 0.09 0.21 +X 3003450 rs100002 7.8e-08 0.55 0.10 0.23 +X 3007890 rs100003 1.1e-06 0.41 0.08 0.27 +X 3009990 rs100004 4.3e-05 0.33 0.08 0.30 +X 3032100 rs100005 9.2e-10 0.71 0.11 0.18 +X 3045600 rs100006 3.4e-08 0.59 0.10 0.20 +X 3060200 rs100007 6.7e-07 0.48 0.09 0.25 +X 3078500 rs100008 1.9e-05 0.36 0.08 0.29 +X 3100400 rs100009 8.8e-06 -0.40 0.09 0.26 +X 3150000 rs100010 2.2e-07 -0.52 0.10 0.22 +X 3195000 rs100011 1.5e-09 -0.68 0.11 0.19 +X 3199800 rs100012 5.0e-08 -0.60 0.10 0.21 \ No newline at end of file diff --git a/pygenometracks/tests/test_data/master_gwas.png b/pygenometracks/tests/test_data/master_gwas.png new file mode 100644 index 0000000000000000000000000000000000000000..99eafe05596dbae7ea66fb28bca08cc1b51548a5 GIT binary patch literal 34493 zcmd3P2|U$l|MqE`YNnE@XwjmD3ZaD*S(;KI*(q64QOJ@#dkakzlR_v|_O;TnbBaoZ z$eJxt2-)}TocFp9n)P|+c|Py^|NP(cnMpY3_q%`h@?Ea$`hL5rc=!OzT#mUE3WY`H z;67yvWe$Tvnepm}+4xD67|$1c*<~et#7f1?z{*z7LZ2e9XJvlY%<8P+i8VI*7M6x) zrrQN}3J7gpbH>Wb+)`9f(B#`21k5Z>3vzO=x5h<&G(ULMl0sQbC;!iwb8FXS3dNmQ zX5SuFyE|X&&TIMV(-}QOd;D(B+hn7-tTsHXOzUAmtwpp_bikWyC9|p0eDX<#cbj8{ zx69}6G&DLn@y3gL){c%Sm7Rg9E$(`H9%C2JZM2Z;|K-cX*rf|sXMBEmH}^tiU(d&P zJ=uw4>5^CSUesL}a?(zDAtbMY@l3y-8jskePk%$cDvDpm419ap?aFca_6>!Sxw1j> zci-xGb#wmiTfqavtiOHh>c4(`&Trp7&0k2L`P;W2n|@%-`0d+ni~oz;YJ6&7T$nkh z|LIKDDD4cJxj+1%7^@yHWZ!?rq%NJRJaLUvQfG;0cuwL{9lLwA64}Fz^P~?Ph*gVK zCvQD-Cdz%mD&CAb9fPqD=TVvM1_g?zPK60%eBE|5UfsI0_UO%x3Kh6MuiEdQz4OFY z`O>)Pw>!?h552Mez?%m1D3!>=>CJ)0R^1=>nAB%$XJ==RcIUDRV;ugY_G29>RMWcj zZ9{h@oHRLI#vgGR76uA)>@+`qbz-T2R`KN@7x&Z0y7G;vF}gJ8xj+8+uIE!8)w1QV zw6wHg(;HuIud1Ze-m@31j3R^7!;*V=_;`DJC!5yAJUHMTb2Q07GtDAAA#1RDoAF#0 zmiTV*x^sD}0wrCt`OJR5O8(_XhFu1~9+lp=@5}o~C+!=rZtO3S6MUT)X|8vHTh8yx z2Zt`TPcN@iZM$^b+}zaKwajP}$*b)?&anRaUa7UU_3@F-+=G&X?@tD9KkentvLXT> zio#^R80+op8<-fXk3C%!%>N)-0g9GiKyd^6ZE9}0fp>sZ5v!$(`Myp8FC9A&2 zW&9JT+*dhyc{T6V;%rNnM3B#!RF{QI1PL2QltwBT7X^!Z{PN4DM}q6u@5SumA>R}g z`6nitg&PZ`?{b~q=x1) zc&0PxW7-So+&pX6gy!rh3=&CH?w&eMp?cfr@4YlSRKC4*qoHcCKu5qjeMxTod42vN z=Fb!v+1uvlmWF>vp*)o%>pR+SD5Gp{ZXPgpEcq-?^p;FPqhsuzO9lP7Se9?wXHt{Q zPnG21kx_1Fqqysx>8wq&|8nn&(8-tc`L)wkTZ%#ilj?ODqm@2kEg>#CBO@a(YKlT6 zbu-hSKEJv+R!d~^;JUiMZnfgIe{MeOR5*``MYaKp-#Sg%Jry~Ow;Y?p<=^_ds0ZFZ zFyFtlwcqMBnUY zL&R(kW2^au;iM?pemEa__^znh{X%vPt$TvBab*pybJABk60=8RlvPyNmoAMgkJX4( z2$8_5Ri|6WK7Ra|Jj=vDqD#D0;pAW`9LdAOcB>8x&Zh_YC&jo}#`ZVec9a%4r1JC6 zKL^AHW9lgW$|tuo%rblS%alCU&$nuPes$5<+_`g~oOyS@(WRxiS-o{73`5fRZimSS4|^ukQ)>7v!wIfP@%>XPo=yQi6C zaNT*VroL^s=hNtTU(rC->5$eh?;fO&w#5pcdi4|LSuNUhYPc+w7t@l2qYsQcbGd;V zKQrznTW*LAp2^id{ACW~-rg%QX7$;kmd)2p(_1&hYiFoODuh%dXlL}-)YH|+hPoJC zGEOlz)%Up>dDoUdczJz=f3-_b0`hua){x)_TgIpSe9h zx47b#Hb+e%b?{TV+%cQ5H(Rb%^Uue#3>?1h##rJh9pcz~_5S_)w*muWU`vU)YO9D3 zQcuu|o;7a?TMEm91qoIR@&Fq{9M{PNZqfTWz$8%ESbf~S|I0fp0`DW*xlb=SL`?5+ zwFa4G^L%iokDn(NPSkPup~OVLO5OREyXoV-`PL&%{>BBjx4vj`>V3V@*Vk9AwR*~+ z2pw8*Wv>rg$#d#pGUG}9>(|lY9vGsc?nuR{(^#PD%=1Zp{Kq0be?Hap`E!)h*dR`| z(w;pR!g46tC-b~~I~`g}!r$cQmy(U>JTdB96@dqtyI_I3y1HL=ikYrK(ESdYqtLSm z`JgATpVn<7yDHatNE&9$CEM(54U069bj8*6kJK3~i7PRBaUh+|<;1 zVb1De_}cPlm5uPdYQs*HkK!=;Uu`l*`@cr~`s+{doSV3~_EiqH31rj>L?&IVe>1tk zcz;^F_(AgB_15hry=IhGHnJ(32Q+Hl-oP$8XLbRP9JO?Fgkhj?c@-?cKu;bgW=?1B zym>cs8;VB<+I?@{e6YeSQ%TpcNBU~hAouGNj}G$QiI%kM{;2E}7!mEYO04bUvr8nX z&>sDm;XGkqXH_as2ET6odh2pE4MR9FS$TP18jaT2*oXk=%c8uvBaePgsj~@sdsi${ zJzg^!R?;U79$h^kAi${lNwj+0k$d!sF^yISl_%*^`}gnn^_rvDH$PXg|Ca{QJHZh$ z;=RJjL2M;q-OC-%MVAEexY5ivCLBBwEzVp*MSj}2*z-%3pN+L_Usvqg);rgK)WvwC z%N6`qYM|`QFoPI5#srNKdFvM21gxG>N!Zndv}TSMFJHFV(#GyOI?4(PDJl&!3PtL_ zyj)?X+Sgp@d*jC1*zz`Ie)xdV(hz!>FQ3MnE%Rs1nPXfWy2pOBHL^VEbkR0%89BMq zuT@D@`vId>Vm60-$|B@fu&mf|NMU!+Wmu{&jkfoNa$G`9&xzjUD5E8Q_qe40#cfAn z!A~X1R(Vd&=3y0aHJ^J2Cd1Uw&4(jZVB$2l*d`}`J;?vT*(9-Z%C##=Nx9D9`yRGe zB=%i5XA?r(wj~*L38y0zl&e2KzQkBpb(qKE9!)23%lXV(vn6ndAiHSUX*z3EQ zGXq+2+Fd#J40mm~%#f-pi99WGeQq`L^odzgpR7v>%C4hp$5+l`evtg3D&1;5JR>w1 zzWa7q{OGdbE^!|GRa#*4eEMPyu7%^@ zl8ey8C!1+OTlzSof^cjXEn4)1#A;PHG&3ivb?OGXA8cn^ym;?)Y=uP;_3o~h|ABEK zy=cpP*!@jgwkS@=n^UtR&^KT;gQBSM{_KvQU8bjuV#m#?u=Lw^_)E-{=`fM`)Bmk} z=zlB?yOS6>4^KwT-!gyK$_2A#&3bK+bUuCM^b{!0yFJ26krMxGC90oNQ;r?(cKPzA zb?bg!yLN4QvW9I{hYi+W!T*Qp4;cN5bhNZoTIhQ?U!SmR*XFZod`F4u`m6QX_om=x zqjI>2%6NIQW!NGOu1AZeFG3mkzJV6>6Ozl?XRTI72!-dqk}X|&vbVfeGhIGeC0#$7 z)_Z9^)u3*}(Rh6!66Y~rI&S*NK%)*In77`UHU3f8C-7LfVBY6;bNfJSn%ZM)rrt=2 zVkgOI+4_0SZ!5x9H_c&4wbj-gO&aiA%o!4Ip>dbGV$I~Scm)$P4oh8eFqn>-nJF)vmoO49eX8wVo9 z6EJta9hp`cXCwP74!ekN>o~^@CZ`7F`?o7mgab5tpexojPRMEN2HnE=@S0O%;m7<& z`6PELzit=dTeIfrgxmC#x^le0aiz2^j4h9O&MqUmbY=7_H+o>2X0^lPE#37F%J{2? z){$l7@E_}jr?G+1;Nz3WmGbfMzc#>Ye|!Jxwd1yn@m#Kv%iQnXxwCTC?AdSrv=1G{ zweeh0+B<@-7jsTctE>InywcxH3(h^Uic`PDm#6uUzG7i06e08JwDXq(FD8F)$78?* zFQjJvXM5XK_TlXaL;Lh!+}8*RsruZyaigI@mU(R}cfC!>m{uKDC=FXZ!cikX*KnPF zvkcR-(v50%DPQIDIr8J=U7yPRM9=t5ZW22LGeNB)$~nWC)09qJ$qWQ?sg412ZbfG=?)6A`0I^%^DZ&4@xN!z zK82!lF@Ry3LyiQ6Vtuu=RrTkeFD@y#8xND$auDzYX@|P&@4|cX8No*u{GE!&WvW)dI z!?e+>yZ%4kip-DMvHMLV{E}ab9#Fn+$CydmEL5 z2Y*?&ZXMgAMfZ&zto`dunW)26f*75sauJE{&?ireHw_-Tu})LH*11cj)ni22z#!u6 zj$`p<(}57nb-iUCtG8-hDJ>O>K9K#==H}0V2LsM4YsK^%c+jeCyi^g0(09NK8nJC> zwEuV|cWrd{-Y~wi3`75$H&wn_XsO&rD^Ygw6R8L}Z|&Hnn>2gr;flGBm(Xo1juWV4 zJgJ0D{KQT^Opq&w^suf1iOqO39m?2Wx>UoURVYg{=MIPX;p8C>%wDEK|8z9yFFq=& z@&^ngs4yaGwKltdo8}mIy$mMgcH&|Hk`G0?=CVV+pM&nfBJa#i- zByE&y8;pp1*h3u3dBI&;PuaYx44xvUT4Q1xMi&zj}e1qp7m( za+dYV$Y;mqTD^-qn$(ZPC>GeYs9kqxWIm4`KrqcryGJC=H_IN0csZ7GWOGwz{;OAt z?G^FH<985_z$8^#?g)~WR^?SY&83h=Y9^lz=*;XxWcNMAHGk6gPCom9TcJT|pmor#Uw;jq zWqG>YZ6x}|SVBj+wqW7id$(^dyT}mrc*sflSV#iL$(Iz%W6gMv>(|RuYG^7ro`4H{ z0|KH8YixLdmcE_MsCVqUzIlmfunYaXObI}+Q;&|wI$r0L%znN#!jO0E+C3>B>|=9s zaz1x;B?8Pb`ciz4Or?NM=8?{k!LItmXqCvenAim6n;YcWe)=i2yW4!Av+h`FgnVgR zd8}tmb%)K-%}e<4cAqB(g@MXcpk=WyX&WgNyTy}UfDBYo?(B{r{1=%Y#lBG?xU3_k zKH>A{&x>l39zT9!RL&c(NtHc!?%e3-n+j!53`-Z8&?mZF+5zCa1GFG_RuV88()oyv zniOMPUsgsYx+7%q3M_Z&TXv>Ge{eYBF+_BdZaMOoPg!2N|qEopW=_ct869g3_u z#Mz)#VJr4dYA|LA5v1$V_X0O`kHymX0r?$dtlb1&FmjA=eOb1yQ) zSyM?#iC{B0Q3LN!dbd^^C;D;$`IS51wekj+6TqUj4_OXZPv++43R`_CUewYzFrbD- zX=-mjc3&Ofv2MX{(VFizvk$vvN=g+Zwbu?E36v%e*DIn0Cn5Me%m}0|6MD+WH~K?yDwr8oQ#UY1guHMC!0Zc0JN%TTkn*D z#2HA{&bvF$-QV9fixPTo4JM#$Q|^{q!NEF>rc*$187B~2e&;L~hI2r~(T2)rm*&U~ zShn1?HNu?{?i5QhMilTCc>pA3BYiCdEni`px+g6y&D8NUHnj*sX(=iaGo)t&tyS_! z)kyzAk}(o*h8TN&zZ_4WOu;Q;_fv30sIMb>$%Y|p&CvrAS%p|YMqA8O3@ z-+WZL0|CvMO(c*jKe5c6BoM$!N}@uwtG7$PahAoO9LsOVke`LposNEsc85nsWPE@* zDmK%BV#i=re8R%z{42hs@h|aQ?eKLgGENyEo30Pq(bW{nnOEcnYR!N0vvLdh!`B!$ zxSVG3E7{PJFUcqWZoPZbl&N=JdZy}eI9DrJ4Ryp<(G$`$vbeJP-im2vJC1OV5sVQC zSP#t`M>8`s)1v0Ifx!~NAzjscM#W%0$XltkWVU=?Fym^tvbv<2YD`QWLNx;9 zjkZWKqQ^Qj6ilk0yaPJuIi?^ht88bN3ZycuNMh{J*kFf>jZHF8ARmU1G#N9-A5Ua1NI}td&b{@B_Xw=puEJR)^yZn%$BM%1$e}?jJ79M0_K)A@xFfD z7x5OC{^`@kwW;Ue4UW1^IMGL1(g!NQfMEv#Y)w7opsTB!)$?+N{l^P)WY0=E48|p- zxACspZTGNj>-N)c%islhITdUGb?F*3?*2HMIwcH5a#-_D^($)9MJq=d@;YnlNE zygb=)#G>%_*1dp)ZQVBel~s7`1&d_^;M!)#rpo*8zdzqLB(okeuq@9y1N&d$zExDo>uAMef`?ab(soyG96s?zs2 zgzX}aQugTNfmz>wUsgSYA(-VBxQOUAOS+D3KqQNB9J_rk@9EbJ-8VL}{5k@=6R@~c z-8GZu2-|5UNe`ra|2CMJd{ePG@z62P{{fkEo&BnkWVe}>Bp)I?p?$rhFztB>V=AIE z?EKWpw_x&`EFOk|Z}KhA1-#LHs2Sfoz%4L?Vkh}GfE~ritF*NYhm`$S-w45%x1455 zeFIc{Go6#Z`|i71-A$7hZOghehc#yJl?A>MoprYL<$$=2v)524wbuyFqGs@4Bu=;> z_3~8iG!XF5i4)#$v~UsxIpp)RG?6D zT$q3mxnhFFAF=Kfzv&=w{Nus+qbV}##>=mhH&R$;=hGkT2Zqw$Z?)XLp%v;Eyt?gs zdU{0f=s0%^$`r}nA>UA_C!f_;1(rP+|42Is zM)jlc2tm)MlJByko5#8~447281dA~N0I3c<6Z{JtfXA7m=f<)&hVwPcA%YxVz+5L5 zB&JG+$m#Ihe8|fUZ}tfrn>TQ-zZmDnxGQ`_Ht3EZBTh)#d)eE$%-d5`F9JqWwY!=1 zrwlA}MuQ~p^_B%pyRYbQ@!~};_Jo9%xt}CEfW;e6Vr{M>YsL?R_tfrTUdzaPniI2* zpYZ7lLr(qZvn2VIkojds6x+Ydfk=*d9;}u9_7}QU?7L+IYk0D#rMke#YsKA|nSI;l z92eEOMXQ)uX6WFH=g*@|S`OBwKN(*^d3?8oUv7e+nm_`kK*!LxYca6y>4Am(Tl~9U z1K2sjcIXD4mU1&ukKNb0Iel|1=P++>&WeywJ|dp_%wSH=`e@gqxD<(Qfx z=6blMyTJq~hZxY8g1(VGc(537D^>r1QO_VU{K{PKw~B`j&40$c8}kE{)ETe=WqI8U z|LN_?{JEVMU~8#~<4ILNOiThT}V|&hGGkF<`>)a_(6dTUc z1P(xf1W1lU0DgZ72*#vWRowMBiUOZ2RP;!uQN1nAmC^4u2sK=)9?x(Q?bEnw- z`Y(ZOSCMrjWbY3BS70B7^5GYPEIz%2OzOV9m6}$JafiR-MiPeZ2#TKu7I71Hf@4!y z9j1@jTe*~c8(FiCe&X(I74)W7j*_VN=af7bIe=^W%YOq|{n?RGLf?{TgW`Gx=>T`a z!IHF(2c-r2D-i{KxaoF*@!@7kLAmnhoScu6=cE^X+kl}+$TNQ24qqoF{h#non-)*e z7Or46r$B>S-#iH($W`<$)1^~Fg@OwAf0<5@%>%rdEX?#HQ_A=L@5t`{$7rViFGYs` z)zeO4ARG&(jSE*$VRNg+)RKjMhm0L@(r@sD+ruSSft=Mre=RFF`yuUZxV{P}a) z0b#rD<4CLL(jdpn+U6}XDUWiWzx;l_AWdt=j2UxTSv9PzoXuu zouU(g5l2<-x?Q;fCJ$+rJQ2KU##aOIN;aw4AY=daQC|NLjZqh#tz+NH-SHtSqrhq7 z;XB0|x=!4lFYZGJ5({w*_h1X;5-GHacp~EgR?_Oe!cAg$q>^@!Bq#CJ3Pe2)xLBv5M{D-j19nRmAslc4Wx>~g|LAy>< ztih$}As5*64@LjRM*NvQRE^s2KS677Um@$(;!{R*ae+@i=z~{`wHI_KZBRk_bIIS{ za<7c8pYPWSph6AFU_aK<_Xj#R!A-(1y4aR`$(@ur4s5`GQko_?>>mcLSN_%!qP)0# z`TLBLi!5Z#C-EI!lUL|m zN8<;jM8-)(T13o(Oas#>k8@K>L^JUpC4rR@mD!Ypu&bV?iy_?Lia(mPUnSxY0p2`^ zXwcfFcc!=ZkFAYPIAQr#23d5i`n5kUzNE?1|XdqDKU$Tt;+>Rnu@qewxc0gF69*bnTO zRz{DDi)*YTLShJbh^m+93;^B%@-1&J3?drgluEDOM`R-oDBuWj@oX=m4*r=8P^6X1 z4s&vH3PL#viD$?AM|)4cy!I|k=BChz=f5-tjHo@b>p3CAYg;jd9NJUY_Dk^|sEORt zZMzb1J`JC`&cO30mZhFM8GkJ0F!qpzamS`0y=B#bkL6Bgq*XX;*&611sV;X{%0TGq zYV_zY1hc`!B(X-3V#io?aqA6j?;c$F;54Kye`lwXjEs!^K)FV2StT$cWZ9x4!^2-7 z^U;bYbQL89O-c$^NY&xphX{PnrP5~w!$ss$>&GW4HDG>99kJ4fTxP4 z8Bd=|5ka_*6;avMuo0DAbbd?X;f_=xzwc*LLg}!g1S7Y6cEdyKOJYE5c-F>wGa1Np z@Kee)!$ME#45}ez@=kOJPHe%7R8q)F@?28w_ z=d>hhVQgt1tB+4-P$7O* zr9m;O4!7y+?;i<*h6rNWD|+7@;A|zn2jE&+tcD=@h?s4sDvS{nr4dc}{;@6$dU;IN z<7sV@TIT}^$gHN)?bVqZbfIER-r#dXBejfiOl*?YEto$YM$9Y7nr4r^nK z{Wc1d()8)m18i)_yA@b>ce~D>zf7@lnMMxd{_$+6Xl3Q(EYpI4Nb5GF)E_HN9NPD%;w70DCj;aS+Mn6mxuN+L&VJMNz(4GqMkL{z)708%>0F|GzBW z_^;xsNo@cbFkb&xH6#CD-nKlTn*r506IM*n7K51}brl@q_D@jdp&9M$JleW33y?7R z33ks$7Z(>{kzL?(*ZtR4sx{&PDR-SXH)QKrGn%{ z7)}TwSHp*9JWWDoHOjIw^pnJ_wu!g8pE=9oRD&E3BqYl_s_R9B?BGiDMG=Nn`1J58 zQmcMo08EdCBPI$_A`B%V?Z_~jz`=*7UN0iii&1UrqruL)*oO}v&K(QdVG@XB--3;nO2rX z1{Xh7^UR-5_xvno;k%nly-#-}h)6<@teqT`(apc$hU~p|!;EKhIM0AwA>50mW*~>K zk?IKs@&h6hZoj%QxR@_7Q#sz!RZb+S64#-Yci&-I^Q;Qy1&u?Ol0ir0mq$*E*IWQXU}V#`*DZ@b1=aGi@pgUoc1xzm<_n z^@5cTKkwW%$Ep<5s7Sjob9e=xl*f{4pLXv2ehzfjk+34I@jx> zqQu(13Bnr4`G7vqjklcHy_zY9Eq?4K*D-M^NRC&FXUceg$ZNG~YgSrsT^lJiQ-fdK zfBE9_Cwo_n%-W$=y1=AN|K&awE$Wdx{ub-Ls2jk4f(4diBG|mIFr{f5rlJ;g9CRHY`uC8-yy5Otei*BMS3+niKwf0)->n-!cBO@*Icog>%3{s}$ zaIi(Rle2m<7tfY050hmlQyeLNvI-ye&M6p#+et_WWGcOJM7Tvm9E#M9K(IB}Zj^?Dq`_WCu(U zbRYEJj-Q{)xkB5tTlcwSR6tJG)e7vKT2)MtYdhA6b(N?p(n7Kgpi37a42h`I=o18F z;q~;v6r`6f=ZklT1$adDd6u} z@GKWDTo66ico~uXW7)eyE~Kh~*D?38>5YHytA7lP(eBvv1U2uCah&ja#|5SE2+`*xKYTHXEk?~h!Z z44ho;kB{Dp7Ul4day-z$@ z01@GTGy%s~@lY?9C=KeIhU!RtC~@0r3h#UvDCtS;cZin%lMNf!`7TNV5kI0trtlgn zIr1tCg_4FWnaR+^z1kI&-c@PUJ9>^vKjCvc#gBYR`xA&Su2wUbu^!*X?kh-FeJ>WR8di?4Fmi zPuK%aF5dgv4%mjfV)r1?7LCTl(Iz^~^q09SD=R}TvjOLHA1a!NLW-Eu*LisZv*tps z!X*O@0MG`Rv<_nF4Tc%W&QO0x3Smae6P?vYL{vNhCMUGk$sLfCw`i?WIUp&c>>c%? zKzT}m?K&4Kev^YbN^zZC@DIA_|982;B;fsV?%%?uNX(fs+_ue=qSw_*S?u4K z)u_`UoBmQ?pt$kx#J?X`K}JycV82nPJ_s@qo%)KM?k?33BX4YtY^v5&4HK`H)O3l+Y;>7ujzjM-GZ`o$2N#;8< z5WF1bpebrxB~;PvzATkcvRC|-Zt~x-5 z$o%gqp<0(2ECmdS?J@b+t{PR}yFUJ1P>kV~9f+LoV|L;``SqorFFObo+Q)^Zk%!8$ zMQ{`as|4zA#KM>v!7G@?kdj4M4ja`$rC+Fq6>&#Vj7(g2q#TbvKrehp)GXJB>VQa0 zK%iG~XP4!Z^L+;T`ufpO1wo%(0+r7I3gO!^=n6<{jK*tW-M`25M0^QQ%BmlcOOrb= zGDldYOIna@UT~u3mTr_az@9h3o|)9-j8>i(G;$aqtKO4~oicV?c^`AQ|tUFLBD z$7pgS^V*T1y^(M)%9Ck_(dAh0pZ9gxwTAmcxRYlwRFA@32oE#33}1M#ck)8m4MnWT z8eb=aTQ267GH!CKX1H!V60Z&mZLVw@c$0Viqgq8)tQuRa~?b-_H`cFJMgSY#HA@#o%V}gw)Ocx)~B_HBy zh3JOrEI^^~=7TzX>f}69unFgyJ-MjpjUN}$l^YzTd|VV%_;NEjyNH?+$w@{!L^9Qe zp?!Rq-(FTzo#CNsa|MkfC9iL8v_@pYw`*70^FizsB5MUGrL3iOM`GmjVqp`_Ub7PD z{%F>uJ{hH$+{v6svjHBKLSaUVZNtOD$_8tqZjr7Pi(EvYV3hw0SFU;BAB%0r7_WQ(dPqnJ zgi2PZ7biMyB2Xoo)YBnO>JRo`S4PG*&{rhMY#uRBFq24Rc{gm>4>_&gIZkNxQ~{`1 zBa$V+Z;>P;Ga;i>6gVZHpr8VYmzcWVjIL}I2+g?2K}L!NoL?O?G&2DZzvnen$1@#Yw0)U#m2B?XBZT0C zt>gp#K{V;7>AO2>le6p`FmZu{uh?+XUl~#DAN~x8+heucvK1@%pcB1$^OIjdGf;MA z5G9^7t*{x5UOL}MvutTY5UOTw{&0tBZ4_TZ>U}$1KpFk8J$$^pvP8BDEdfznYCCjL zeZs&}v;zC_Vd@O)GHeCUorMk&T~jn09ns~4lo2~RImwl^9fmG34%%cvnZVuWA0I1X zp}6un?$)N%=KG%B50gnXX-x7sQ>lY6%QOKQ4_726MmR9D7WvQ-v3|PuhbiV@SAZ0+cywL z(L3*KOvI|{D7X(~U5~tuy=}*tV#>U>QJxkOa0!z-S+Wg?j+BYCu-4H zqWZDECsaNC5PPfOEh;c@);Ps4Uc6v(-lnDrE9$FUP~uYO*jHGcVVg);Pg39yqG5G3 zmP2DYNR}#bIc+if_fgfIYTh3yigVCi=RZ+a^R*WQY=`?+mzDK!baW&=gNVO2O$2xL zhN-u~io>(U)?gDvd3~k)bmdw>qm`@jfC)5rm9tfK_JW|p$dJTm^b`f z2WLsOfO|3GXh7sfUR>-F9NmA8#`-qxp+udA5ST}XkRoc(1Q@v}NeRkyW%hXiB5;09 zN8X95NTu;EwAJ7d6iloggbeUS{>0MqUK(}cF)h%dr%`&M)0S51!Oq7eLt4F%N(K%| z$Y7|;mRP8F=2$qhc*2pfutoCwX$$ALFu*0!Sr&NbPx6cg{Ql>0rA8 zP6-85FY7~SEwu7z#ik<8d7yYND?zTvZUUD11RY!^`a&3G$k8Krm{gyE4B^EDh;Sar zgHb_J0yh!?(2GY*>==PPuvg;CmBTz<-(j2&B~T zSj0l+In6DbLJy;NjCg681P#aK9FjEA$~cVpYjqSGMKK66345r?qss;XY9L#ZQAPj^ zN-q*%liK@wKoSY0w}#Io8+TbQc?gd{$aSOU39#p+#Y8Lk>u|4tU6v8xlhyS2n;sq6 z>}#ZuiLz~DW)TnXi~k^$e=1A%j?SA&q=4P#od)t*s>9&AyAtO6FnIDMmP^d*^7|ROwF%X@t)duOV<= zHT2o%pk*aIIiWj=K%!E|As)ppO$c`(rvHLsAc7Xf5nEfpJtPPwpTHZ0zysFyigATP z8J;yQ?i$#mfBFL8gM6kf(RuS%e2Q16Vxa-fQ3VMd_xg>$Hd5Vl__l=V@n>|^w^#wMsdo4h znv=JmfU0xPRQ(LJy5yb0O)H$Mg&o!H0Ch$#eMc$Vib74m$GJYA%AV{Enp+cs)#97u zrDkBTq0)+sl%m*vs8m%)Mc@hPR_)vD+^+V#9Vl6)6O88YRmc}jnrEa5jY7v%z3^c znq?A_e8+rHSbO9c_hgO$^T{}lno+$fm*@kkW_MS-*T%Tw9(EwBAz21|TKc_h6xY~P zzbR&&qA0SpkJD6dnY)c?O#<+Qj53KoBr2|NZcFILc)jH2t-s;$+6Fs()x;Vlq>|Pu zQ&~Jz#_;gw^MV(%PzLtmfUm2jRwj?p6+2ldQ>Pg}(ZDZ)NE;I7D^zdY_?2Vbcg|3v zIAIS+s-?HoiA3@^J0-0%58Og3h%ZlYuz(_|-ENHps;hw>8Z~bpIjB79ZPZILQwE);+^@1rr#d5tA_xy^*XbGDW@=?Z0m6n4G_FuTi{)c?p+> z&F?@6eE}k2PwVsFW=7d$f7l9RzHVH()p!c)B+sq&@l)$}#pdL>zrCVX6`(1A0iCh5 z=<4*Nk)&6e*=|Mu)2?^l-bzt)0?bZf3HlSvWkk5H-`+(Dwf_^_jMIl-KQ&(FzV=G) zJ!o22b56!tl=1I~>=TEJXx9OO{2jl-q>>g+UcYhOZt`W?Rb56Be7ZZ%UiOo_UQ zd;&LoW^u_yrkB+y;)AD78K@|%{Kw|GlM~{aeu>!tKTW>E|F?n`R2rEHL1wk`Hv+ZBCOm=Z}A9p!s_x(z?y_+V^ouAD@yr;1?CL6{uNVuxrRAF z3P+e`(1;6!E#iFz+c|bwgjFY-91E8~ASr|@K(uu<*37Vp&rj;ggTqEPK_n-#va%&t zIGv*sX>?80QNAPMGKV44i5_p4&DL7beuF9q$6U364J((t*BL<5|Z75<9|LPCm36amGv;d^?ba*1wz&c;>WfB$`I(@yN5QE0^fV`4at;>sKJ4<=ozP6drgks$S> z2)K=u%dPj`tz25CA%|VvQ@{WWm*%L2Fdp4zK0`OOaQ*%=_4AldowdK^ZR_B?^pku3 zd@hxAArk&u&~L1y7JFIi%K)W|5tpeUHXRESg? zLUTc4S0EuJsK9$8-pR=rzOn-ZB+7(a(X|q5b(17B=*$&`>e)^Q;`6MFLHx2w>`1Tu zNc-it5Q79S^5yp=9!+9{b-Rtie6Bp_9S=U^RxhDNe^%%l;v4@UgsT;>27v)|`H#Kj zC7oLEI#X;UG@eq^0c4aO4hET&Xve%xi2aPv4)xpTu9~3i*1!4f-B=)_eh#%bj?X4$;KsH}rw)JPIP|tl5dm8?SQmZptpx@x)0rB8p zJvK+gp9{!rnYzpw3OqR=QR>NjbqMB2KSm@*)%1=fKp+=QH@S*=*6jibB4VL;giUmD zVFb1%?-v49M`|3g&F(8mjO+uetqKAS$&SEzlB<+KutDTQfv6~kRJk88b9%JVtrx9U zZ%k8S5$gQbuIUe0K2SntvTbo>pvj`IkJ z34PycXLux%_mi^19cPd9$p7f`1f2)@b+fh5ws@0u?H@Iy`;GQBP3IrGG8B0;XRN{T z2hUnWlpDp}fYOQlAV*-7aHlp!I$@=Ayf$H0@Z}7ahcL!lN)=!k5&(ZI_)oGbNbrM^ zA15U%lw>m&lZ$Z-bopqK{CSR{(Y9mVsyjQcD$VdD;Fx;)`J2 zc;AI2CMFV18esW@J_|S{;(!nnz!pUfR8k_0X2nbc&yg}MH1c{glJmLo`@_Yv%gRdcgIsx-PKm$c5Xse0P3mq0nf=iVD zIL#-bCHJvH#!@6cbO#bJUq8PH66NIO9RSf!DEa}k#y$h_H+5*gk^p+v$RNAR2_J+- z=4St}`DfG*Np-+zgU=uuh;J3NBMiaG9ZIotJ!U|^6u>C|U4IG+CH0A@WMBO;Hn*{+ z5g@T;F=~7wYAQ7F>)NFm{$9?3Va#Jk|FTi0ZDD!PlwgBgQtDCzf_M%OB7XG}>+H5# zu4DY~v=vjsaBV@4uFQ+52!=S!3JNZeVxFK;$4BN* z-^}h6Kt$_$ry=70nSEf&lrA%LP5?~}MlJZ{>*HGVeg%Dr@2vk+lJ#G!h4X)f8ACD7 zF;W5nYZ3`j7^#9GEFxqZ;Y45wj9%1*?no`t4~&4SL!3ykfN|h0ZN>L4Z8&@<3J1V{ zYy{jIsgi=?oCsv`0O{yseCy7guNItuu-|Fh)`bJ8+q7nl6jYHwN-LJ?Ix4x@&Q|F4 z+v?a^y9=0+(P)2bI(@XA^a{>LIzbu*KqWc=Fq}_ZJS&F-4QAwYUC^sJX2Q8a;=?vC z8oZT%-beKG=8JPPJ~z?8e=-%YpMbX>5*ZjBjo1(#l?c_bi5&6Tt?q(ooO@Tyw(YL% zhj$QYv9q&t!$&vvSO8F#9ne@@S&T;^3OJauUm-j~vkkB*F*vCTz@2+n<=o8FpRb#0 zz6V?z>2eak=hEzGs8P{uMdSGK(1h$!Bk}&?y#f=audGw)J$NoC!)(^oTG2%b=zgy`vG2eEjgbj@udQ@Oo%8EtQWt_)15ZlWFk}Ktx>iU! zk8?|2d0>l4 zlNg(j{-zd9kbR48`um?=uwqACEq!8qyZ{ZJwH|wepSz6*G)r?jvUy2g!0F)VhkmDn z&%bd@0_5X#-F{re>dLE@`$rR~#KF9`5&iaRS8V&1Q&6^E>oPIWD4`3oblVuuTn2ityS&bR8!Ub- zd}x%`m6F-bzyfV`a&$Z-V{qii?Gd2Ow~Y1sJ;>DXrhm!xLSbJ$wiu(V?5&v^P4AJuezr#Xw>t^!L!lXMM!# z^H3Q$w+4uyJtcPr4^dj!XhWZ7G*eVAR@I6F=NXHP5D{N7GG;2%voo}NEPzJLE- zHe}|9*z}XbSlVWghS+KajV5RT6Q}xvB1~gnOVLtB+>^6aCA*!+?V@*}HPxYW+GzVF z-}GTqTU`1UcU1Yk-j*A@4xP_M_U_#)A8oD{htqVc%B$TxS9R&jZoyoJT%l<;TSBh2 zu616ZVMT$Rr8LcX;!w2ttT^4^c7Kl?225)Ct6llZO2d-2owbd>Qn3M=-g0alPaKQs z>+2yK2uW1f2($tIH3sD|YA;TPXO?ef&cDa-4{E(1FR?D6*)%PG%V-~{)b-~gN=f|V z1QyW~1YUuZV2uo60?}YwDW*`8-ouKz$b&vQ=bPB^PW+1%Fc__Xv~F;rgPSZ6CZ6&< zkei8+y~vahAp)q0TSh8nw&*`k2cVLua3qGT>BLGl$zSm>%Mio%Y|`_*HZwtlZ-7)N)>$AR;}vjd88dod@u%m(=d>LW>lLFP4J!X zO7-t~V6qoN#z#;!fM}aj3a^WNAQB}s!$PjwSYPCPl{lKv9~xZPp`%zzTY&o4vPZzH zbe#r#xuU_aq7+&%sm>&&PQ7-B$U2M|q^v}SsINu`JB}G}l0b&QX?mRC!qU;S{>|6* z;O|gKu{X`>sTdybs|PIe5}E|G)xM)r6HWY|QrBsDt`IO=mJ6d`+S3m8uB0)>B1q3X$ zBnjwmP>7`7MQlO7md5v-@kzZ#`|nhLM3tBF z8_$4r)NU<#g;YeE?t*p^^qTI{9B#|BEF6VT_nY}5sYyQXtpGwRgC~*LvTi9--z$ta zdVFE*Q!?tIMSVA$<5Z!CXX!5}S^OxeGAj-3MM*tvV7<$VdZM@uIkCjaH!oV)?vG<* zEmKSiw#u8)=^F*n7P|nZh=e9w1)1ScQ@Qo0p9WFzysso)e1=}(rkirh*+EYn=L9||m333H*+jO1#2>x@#q=a?(`-f~q0pC1BA|#CF!q6) z#OYUVWoDc=MJ^0V@D$Af&dp^U0k%&>McNHNErJCb~tsY^aq{)U%Xu)rC{~ z-4U%|-hYzlkV8i>0M`eQ;0J2ff<5>IjVVn2&`qolE@H}t_u?g!$7ofgv=)v?5?&$DP;^xrn7AO;UhPROddg6;%r~gN-%<7Ye5_4W5bw-b5fQlpQ)z zI62S@!qx%|i-abwehZA@d22Ni1HG0LSTXvX7-RO2mqhR*y)QD#P%l@W;*02V~aObImxIfQk{V&ojPu@Vvnzx$ zexSMY)j`Mc)w|yjTpAcvW=SL13%dRoq0AsGLZR)IjMC}5h=|~~2f;vy_={2O2QY(H zM<#)~PrTtn`KNB88K%V-((flqX`?0Mf+H?kxpboKg$qMKDI>pP6G;#1$#5~|osT4z zIRmskVt1&BqQ10Y`Y*7=((IEc{-c2`4+R*~fqS}fk*P~d{vc20m5IkY z74iEdAxx;b`{04$66ZZ>KH5^RMU)w*aoC)Z<)W@1g^-lKH0KW3HVx|(@s~(RB403X zW2>ScLpA$xPA6n|VoE5Z3X`SixD{Gb#RATq#lHRd{Ysqu#AN`pCsdDr9d3bZOME`J z_F+7vf*w9giL~(3%(bhNq919qA|K>q;$tOx>fAOi%tZ{wcF48ALK?YH!dp~B8!=e` z0+xcjO<*;mFj*&G1ZNi% zn@Cv5MkAYkCVAuHh(0DM6kBcNWfEqUe93$$@itH)cnNF0+H~Q%2ovTw+>%1DTq3+t z+{eS9sQVESajF3ery@u;52KjTH2r#VctYV7bTg)qQ4HcX>Oq8q4I&3f^XA}}5Ee2a zMI+bzBN`EuD>rJUKXpy3u_-umKSmLx*#jy_(1xL6NJL&i;UKN&fE0)>1er{<3I|T3 zgXLV2n2+Jb@b2i>_SR-W@*&;4`>6MIEpg?`rA^@+W+QS}3i;bE(S0|vE_N8g{E-n{ z^g3EwGK?Lu&+rl$5_Sj_)?@^S5q}ZikcDW&W3fl+i!a?b-MpIE4R{?*kI~#T(5VqB zX5t?Iq$bsmb4JuzfNv$<6!SY3F@=jWnAc{J<{f~GeeYa>R(IyrWoNXvXKF2*Icd%- zj}_<6O~Xp1wW!IC@>rM0jv2Q0uYDbBe(d{WUw^gjJaIy1JVJ}|>)sGuYpvv`5eaOE z6x|`8cTwD%Holh&lH)ZUt%s8n68m*VB+mVro)ED{PyKdq&D+hCtS?_@XV>T2Z~dTL z^VPSl!Jb>YZ!xfI;K*3L=PUmP-}=Q@GRq4-=y>gYyFo4lussa5ERgjVn-Cj0_VEtaOnPBK%iv4?vlssh*R=#@N>kPt%~^zGsBDh z+8ZYU(iKtOq$yMGQ zAfYV+rgA&ydxGWiwZxuhq_f?ik1ET;1(mBmGFhw9pZ^k>whDO*oNW|t&zJTsG>%C6 zEyrYkX7@*zkIyCR6$c*sARXUDu?D$v9Pcyh-FR(jdt{z)hAZ;cy8e7@v(Io!QhvLu zSiRhNW?mNvX9C&NA$Oc35Y!yBhDnQi0r8$j$!@EdxmG~9=ZZjC1^~ulq=8Y?9iXcCEs$H_Y&h3$7VNps} zrZcmx%h6vf1PPs74$A{CcMG$)G=3{AkZFG??8;w9z3DyNB8ky{;W0sfykr7M-}N|~ zJ}@M^C_M=GGspFZLNAE9KAOnB;d}tHBb*BftW32@x5Is&Ov2o!qTa<|D1htd!S^i& zMYJp6C+?z9tdY2YoNf8Vs}_vS3eeL3D*UH#>(?k;7M+#4(wUQUQ=tzJA? zuPh+o4h+2vfk^x3$HUD4b}$l9M&@FERNqbi)ZUNh*p_8sGzZx|;?4yN61z_6rV%=0c{+tv(&2m*7WXC2vby|i@O)G3|1aB_rsTCySRTN%a?M) zxfw|Zhnba6me0G-nMs@RHsCX8TP%1+l9V1O!1}7I?>nEf8@f-M!Fs#+{uN&#%a0ue zoI_qkGdiKrLLB~_0o9NSFb_pujOLr_#Hy;0v#&Ipjfgk3d8zA%{{DX5kC`xCjer;! z3?B&2I*GD;C%GK?o-+#Fzg^6Sl#j~ghTqL`ClV)V2!jb1%G9io6D}r-1kE zH=cX*8KopAGeFZQLx*k>gHeF|@ln5OGyCDqj%WTU-9`8A-#-E4t=?jGmD}x7RcQih zeuKHk^0l@6WMK}Qep%?Oq!B8nfZFoVnO)W+{qr}mp6Ez$e7T5XiR<0w8Bhzw&$ zpby%*cJ@b_NcOa<6$w{7X){-Ovsc~GvMjFGwqxP z8vAMCT#ZbZlMkX^JwwB-0(DooEeOUv#Wz0;?s!1Xq$AGc!q^j^o?F5;E}Q}?>J#&f z4>J>NTTACQH+6Ls`9ikK%RWK!hqXnPZv(SZw;;mp&X*QUHa~^8M;-A(10d1r|zMf^;6+7SU8;T9$}=nLr~}+wPWWKMA(%K|9Q=tB81FCkJ%5Iv- z(Jw_Q?({B=kg9Sr&Qv(264MFw2WC)Rdj?0}_XMr8oeoo;@0D_Xsl-Z>L9KMvO4yi4 zn4Aa-WqPo27a||mn^bO3zT^g)+1*sdT)5F8gbJ5_5yO>cdM<+X!g@G{Ntfg6$`k2% z>M->1ScOpp*}LrTS^IZ9^U&?%{F(xHklCN_!il5%RwZp^{E#AhE`?BWHw z#=7@bBvw~fYp5jo7F?{>;%l3I;_@y+!Vw8z#+I6GYg36JI7CINhV^Zs)V2o*dVm9Lf}pwHvGtop2MGr zR*>@}{|t#Z(K-8MRSTlgSqK+lUofau0kr)?lJ?#$n=`|knwq>g?v?#c^MEg0y*3hW z7{tBG$A9G8>qD4w2f{?cp(`_A$kUz!s)2UpUJUXKfnZn_hXQy~by~L8`O4(Pr?EhJ zc{>p19Ratb6ci|+!jg!|uT<9OMU+6yoHGz!o6A`!WDxBU zx$-4@DCA)MDeR;!3KW{sMj1Nbf#`zxb+sAlB~3;p3vBa=@4kOwv31x7N6HayMF&)O zU2vh!Ae+^O%+YD`O&$$NH5-(;j?1pV?Hi7aM+bL2+1K)I24t5qkSUx_P$l2r2ZUZ7%VonD1`1HmM2D)*ftd@j4sa+Ne)k!Zmz`?USBO?)*2H z3wu_(q3QKp*jWt3-HF&HfpO}lsJF2nl)9T??3$M5C!1l?9tBfTCEv5G;^=!mrJLMi z`4t#9U0BC7Pu*g84N>bTnBMdjH_Xy)LT#}bgp&yJW|8pbl9>r)ExkG}vJwGm4fLSS zmEGD)+r)apP0qV#@|CwNg)}aOkmf0b4Jhni+-Byft2CS(--N8887~=)Q}_TS{ASR4 zo8S)~5-)Fd+H&z@o0Ci2VfEwDJaXhusx9cba_*cIdx$tYuP!RJ1e;)fmgV}X2}2yC zh|wcQiF~7SKjroMyY+q6hipHcDIFayoWm77ll0g_%M}f7T^3^lHXX5R49gMjp3k|#4c`BbEthbEIsZV-%t@ffws zt(33ucXQovZP~&y){TR^H_(dV0@n2#FY7fw>POf=GX2#E8gGG6vum5{|8BzQk3+&& z?VkB^MvHX9kEqyrx_3nUF5GD`y~z{k!}L0F=jlzw!+s%RDoV)DU%hd6G?y_bUfgQ!w6?e%ajn98h5v^d0vc7IzG#kzVU}KFf#8?%^;Edh|7`@brUOW zbwg%&%lGsNOTg1cJYcGEo%k57&8+$9fcOF1Q%p-Q2gZ#F^c>QFoDz9 zzrs<2t8Huf6?PDwW;=HSmn{wg#Sl_GjJIiVF={Ahk0ks(oQz^_^(QXV`wc%qU-n8# z1p&NK*(z%$MB?3aQUYvtYr6}%ZI5N#LoTeUtD7Q%XTyM5v!-ih4PD8QLYGN30?1PQA z!DVd3#S_BdZ8IdCRG1yV@cH}>*I@AP%I-;??Ew!o1O{OcxWr?(w?ANrg;5z!`AF!o zl`4I{vGf;Ti=mJ*suv@t!vNrX$g(;P=vR>+u$RywhhH#qTs2IQ8)Y87zD{lf zorC(74@5Ps!6+C96R(oQvk3{iF-n$%RE@}a$d-{u!IPgIqdtQ7-q+$r$Moh~2<(nf zDEuh+YJB>63I)6%xSle+xDAS=*?!Z*b6va@LBN(lfJ4+UAzTQcgX2aKg{(JJ#{is> z)!+*9K?lg+kP%rHBZSKK=&v5c`x5c2MKd0!pa4U_AB|eP{o;2%tvkjF0V@0oaH3#v z@hV6sm@!~GFg43i9*cwDYRt4OQCOa8?nR~>ew`moBx7!OYokcs*2W1>E&yAK^~hol zmcHc8;_~P9Y7DYJu*KGp;2@qqF)@+E3v<;*uQqqhRmfXoWG~5=@=m=Zc^5o~5G?a} z)voZet1FYS3cVo*K^S(pd8HDkrk&&%z-*QAtG&4UTj0w>NmL0OK=9`nuO``s*_u7l zp6*+Si`kKQ&qVoFz8!*%#+O%S_ex6#lQpdf+1C6~@>qPb)yF^4by`!ba!&Kzr8`BX zpDlg9dwLyh5<@0TOst%Le!ms(d=H=Yp8^HzHA;1rdARP>+!ZyO`LXQXAn zjONQ@b;XC*Tzl|rzGTDF=%uT{3~#_@|P^~TXG%| z6qxK?RlulHi{@d)Ey%4$G_qDeU>v5SoLKaMxks)=k|3T#Gl*c`USpjBkprv3tn;a> zXA1&J*QDUxLorc?gtEmjXzlyu18QQzw$nK((T}dR9);sm6e1V*>2r$2QlHi_L00C0>Jj5ZRCI)9a%jz zJ2EjJdlxiegari1D@t@iC`kdoi_gh9fGlSVWY@f$lU9YFe2=X@OdB@ey+!Zj3B*z&NUW&{ zp%nlT*&yNuLK+?=f5`&bpei$F4>`Mi-$#%F~*`lEM?m zuRaSQ9akWBqfmN~UW5UERh)r0Ln{IUgQP1Fk(ekRCvUh)eD$bwk|H><*x1^51i$Wx zt;->SDTt@GS65$xJaG_V3^N}~kVpMGm){T(yOOm#bS8E+&fNbl>3)G^mSIig@ZEgV{;qLm#9D z@mcR@qsbD?c})HQajx6)sPBtJntPu=XK1Jao)by<UB;pI+~2HN6zdU;s<61e)9N z_B>gUqbrT1qhz8&?&^F^6e=F->N?SfB;jccdJAHg`hKRmeSXdO;1_%i8bMj*9;6@X zntvY|I(gJ0>Y)AU#~#&^C=LpMLl=U6y%t5xYM-Zqf$#yu|HsHQk;VN4O0mq{eR$}a zzHO@}kC_3iiHL2HC3pQ)R1SA$Ii`A2#9bX)dpju<3lv~{Kxe;9k0uYh&3gP_Q^U3m qKr7{suPOy{tyaJOkNjA`TQ6VhL^W8KyhJ{h=84nl8GrrTum1q_BJwH# literal 0 HcmV?d00001 diff --git a/pygenometracks/tests/test_gwasTrack.py b/pygenometracks/tests/test_gwasTrack.py new file mode 100644 index 00000000..336201ef --- /dev/null +++ b/pygenometracks/tests/test_gwasTrack.py @@ -0,0 +1,50 @@ +import matplotlib as mpl +from matplotlib.testing.compare import compare_images +from tempfile import NamedTemporaryFile +import os.path +import pygenometracks.plotTracks + +mpl.use('agg') + +ROOT = os.path.join(os.path.dirname(os.path.abspath(__file__)), + "test_data") + +tracks = """ +[gwas] +file = gwas_1.gwas +height = 4 +title = test_1 + +[spacer] + +[gwas_2] +file = gwas_2.gwas +file_has_header = True +height = 2 +title = test_2 +color = #50E3C2 + +[x-axis] +""" + +with open(os.path.join(ROOT, "gwas.ini"), 'w') as fh: + fh.write(tracks) + +tolerance = 13 # default matplotlib pixed difference tolerance + + +def test_gwas_track(): + outfile = NamedTemporaryFile(suffix='.png', prefix='gwas_test_', + delete=False) + ini_file = os.path.join(ROOT, "gwas.ini") + region = "X:3000000-3200000" + expected_file = os.path.join(ROOT, 'master_gwas.png') + args = f"--tracks {ini_file} --region {region} " \ + "--trackLabelFraction 0.2 --dpi 130 " \ + f"--outFileName {outfile.name}".split() + pygenometracks.plotTracks.main(args) + res = compare_images(expected_file, + outfile.name, tolerance) + assert res is None, res + + os.remove(outfile.name) diff --git a/pygenometracks/tracks/GwasTrack.py b/pygenometracks/tracks/GwasTrack.py new file mode 100644 index 00000000..06834b1c --- /dev/null +++ b/pygenometracks/tracks/GwasTrack.py @@ -0,0 +1,76 @@ +from .GenomeTrack import GenomeTrack +import numpy as np +from .. readGwas import ReadGwas + +DEFAULT_GWAS_COLOR = '#ff7f00' + + +class GwasTrack(GenomeTrack): + SUPPORTED_ENDINGS = ['.gwas', '.linear', '.logistic', '.assoc', '.qassoc'] # this is used by make_tracks_file to guess the type of track based on file name + TRACK_TYPE = 'gwas' + OPTIONS_TXT = GenomeTrack.OPTIONS_TXT + f""" +# Title of the track. Usually displayed to the right as a label +title = +# Height of the track +height = +# File containing the data. We expect an IGV .gwas format file with the columns: CHR, BP, SNP and P. Optionally, extra +# annotation columns can be added. +file = +# Y label text +ylabel = +# Fontsize of the labels +fontsize = +# Color +color = +# Optional. If not given is guessed from the file ending. +file_type = {TRACK_TYPE} + """ + + DEFAULTS_PROPERTIES = {'fontsize': 12, + 'orientation': None, + 'color': DEFAULT_GWAS_COLOR, + 'border_color': 'black', + 'labels': True, + 'line_width': 0.5, + 'max_labels': 60, + 'max_value': 1, + 'min_value': 0, + 'fontstyle': 'normal', + 'y_axis_max_val': None, + 'id_fontsize': 12, + 'marker_size': 45, + 'file_has_header': False} + + NECESSARY_PROPERTIES = ['file'] + SYNONYMOUS_PROPERTIES = {} + POSSIBLE_PROPERTIES = {} + BOOLEAN_PROPERTIES = ['file_has_header'] + STRING_PROPERTIES = ['title', 'file_type', 'file', 'color'] + FLOAT_PROPERTIES = {'height': [0, np.inf], 'fontsize': [0, np.inf], 'id_fontsize': [0, np.inf], + 'marker_size': [0, np.inf], 'y_axis_max_val': [0, np.inf]} + INTEGER_PROPERTIES = {} + + def plot(self, ax, chrom, region_start, region_end): + """ + Plot a scatter plot for the GWAS data. + The p-values are transformed as -log10(pvalue), so the y-axis will show the exponents of the p-values. + + :param ax: matplotlib axis + :param chrom: chromosome name + :param region_start: start position of the region + :param region_end: end position of the region + :return: None + """ + gwas_reader = ReadGwas(open(self.properties['file'], 'r'), has_header=self.properties['file_has_header']) + + # Fill in the position and pvalues lists with data from the GWAS file + position = [] + pvalues = [] + for record in gwas_reader: + if record.chromosome == chrom and region_start <= record.position <= region_end: + position.append(record.position) + pvalues.append(-np.log10(record.pvalue) if record.pvalue > 0 else 0) # Notice the -log10 transformation + + # Plot the scatterplot + ax.scatter(position, pvalues, s=self.properties['marker_size'], color=self.properties['color'], marker='o', + edgecolors='black', linewidths=.66)