From e57ab42a344f1ce9d8b5626f145bec90264229c2 Mon Sep 17 00:00:00 2001 From: Armin Toepfer Date: Tue, 5 Feb 2019 13:11:13 +0100 Subject: [PATCH] Version 1.9.0 --- README.md | 70 +++++++++++++++++++++--- img/plots/summary_bad_adapter_ratio.png | Bin 0 -> 44599 bytes scripts/r/report_detail.R | 2 +- scripts/r/report_summary.R | 10 +++- 4 files changed, 72 insertions(+), 10 deletions(-) create mode 100644 img/plots/summary_bad_adapter_ratio.png mode change 100644 => 100755 scripts/r/report_detail.R mode change 100644 => 100755 scripts/r/report_summary.R diff --git a/README.md b/README.md index 185b4d5..324060b 100644 --- a/README.md +++ b/README.md @@ -76,13 +76,7 @@ The sort order is defined by the barcode indices, lowest first. * Enhanced filtering options to remove ambiguous calls ## Latest Version - -Version **1.8.0**: - * Add clip lengths as `bx` tag - * Enable single-barcode samples - * Implicitly call `--dump-clips` in `--isoseq` mode - -[Full changelog here](#full-changelog) +Version **1.9.0**: [Full changelog here](#full-changelog) ## Execution @@ -196,6 +190,7 @@ how ZMWs many are *same/different*, and how many reads have been filtered. Below min score lead : 11656 (32%) Below min ref span : 3124 (8%) Without adapter : 25094 (68%) + With bad adapter : 10349 (28%) <- Only with --bad-adapter-ratio Undesired hybrids : xxx (xx%) <- Only with --peek-guess Undesired same barcode pairs : xxx (xx%) <- Only with --different Undesired diff barcode pairs : xxx (xx%) <- Only with --same @@ -212,6 +207,8 @@ how ZMWs many are *same/different*, and how many reads have been filtered. ZMWs for (A): Allow diff barcode pair : 157264 (74%) Allow same barcode pair : 188026 (88%) + Bad adapter yield loss : 10112 (5%) <- Only with --bad-adapter-ratio + Bad adapter impurity : 10348 (5%) <- Only without --bad-adapter-ratio Reads for (B): Above length : 1278461 (100%) @@ -509,6 +506,8 @@ A `prefix.lima.guess` file shows the decision process. Example: ### `--guess-min-count` The minimum ZMW abundance to whitelist a barcode. This filter is `AND`ed with the minimum barcode score provided by `--guess`. The default is 0. +If there are in total less barcoded ZMWs than the provided threshold, +the guess feature is automatically deactivated. ### `--peek` && `--guess` The optimal way is to use both advanced options in combination, e.g., @@ -522,6 +521,9 @@ Equivalent to the `Infer Barcodes Used parameter` option in SMRT Link. Sets the following options: `--peek 50000 --guess 45 --guess-min-count 10`. +If used in combination with `--isoseq`: +`--peek 50000 --guess 75 --guess-min-count 100`. + ### `--single-side` Identify barcodes in molecules that only have barcodes adjacent to one adapter. This approach makes no assumption about an alternating pattern of barcoded and @@ -678,6 +680,10 @@ HQ length distribution per barcode: +Bad adapter ratio histogram: + + + ## Implementation details ### Smith-Waterman 101 Barcode score and clipping position are computed by a Smith-Waterman algorithm. @@ -828,6 +834,49 @@ For asymmetric designs with *different* barcodes in a pair, at least a single full-pass read is required; this can be two adapters, two half-adapters, or a combination. +### What are bad adapters? +In the subreads.bam file, each subread has a context flag `cx`. +It annotates, among other things, if a subread has flanking adapters, +before and/or after. Adapter finding has been improved and can also find +molecularly missing adapters or those obscured by a local decrease in accuracy. +This may lead to missing or obscured bases in the flanking barcode. +Such adapters are called "bad", since they don't align with the adapter reference +sequence(s). +Regions flanking those bad adapters are problematic, because they can fully or +partially miss the barcode bases, leading to wrong classification of the +molecule. +*Lima* can handle those adapters, by ignoring regions flanking +bad adapters. For this, *lima* computes the ratio of +number of bad adapters divided by number of all adapters. + +By default, `--bad-adapter-ratio` is set to `0` and does not perform any filtering. +In this mode, bad adapters are handled just like good adapters. +But the `*.lima.summary` file contains one row with the number of +ZMWs that have at least 25% bad adapters, but otherwise pass all other filters. +This metric can be used as a diagnostic to assess the library prep. + +If `--bad-adapter-ratio` is set non-zero positive `(0,1]`, +bad adapter flanking barcode regions are treated as missing. +If a ZMW has a higher ratio of bad adapters than provided, the ZMW +is being filtered and consequently removed from the output. +The `*.lima.summary` file contains two additional rows. + + With bad adapter : 10349 (28%) + Bad adapter yield loss : 10112 (5%) + +The first row counts the number of ZMWs that have too high bad adapter ratios +and the percentage is with respected to the number of all ZMW not passing. +The second row counts the number of ZMWs that only get removed because of +too high bad adapter ratios and the percentage is with respect the number of all +input ZMWs and consequently is the effective yield loss caused by bad adapters. + +If a ZMW has ~50% bad adapters, one side of the molecule is molecularly missing +an adapter. For 100% bad adapter, both sides are missing adapters. +A lower than ~40% percentage indicates decreased local accuracy during +sequencing leading to adapter sequences not being found. If a high percentage +of ZMWs is molecularly missing adapters, you should improve library prep. + + ### Why are *different* barcode pair hits reported in --same mode? *Lima* tries all barcode combinations and `--same` only filters BAM output. Sequences flanked by *different* barcodes are still reported, but are not @@ -878,7 +927,7 @@ Even if you only want to remove IsoSeq primers, *lima* is the tool of choice. AAGCAGTGGTATCAACGCAGAGTACACACACAGACTGTGAG ``` -3) Use the `--isoseq` mode; cannot be combined with `--guess`. +3) Use the `--isoseq` mode. Run in combination with `--peek-guess` to remove spurious false positive. 4) Output will be only different pairs with a `5p` and `3p` combination: ``` @@ -939,6 +988,11 @@ false positives. ## Full Changelog + * **1.9.0**: + * Add `--bad-adapter-ratio` to remove ZMWs with molecularly missing adapters + * Fix rare case, where a read only matches one barcode and not a single alternative + * Fix `--no-bam` to automatically omit pbi + * Allow combination of `--isoseq` with `--peek-guess` * 1.8.0: * Add clip lengths as `bx` tag * Enable single-barcode samples diff --git a/img/plots/summary_bad_adapter_ratio.png b/img/plots/summary_bad_adapter_ratio.png new file mode 100644 index 0000000000000000000000000000000000000000..22904ba82ab85858399e230a810d28089358ca74 GIT binary patch literal 44599 zcmeEvby!r}+qPn1AaE!J3B^K6lx{^rDFNvck!Bcyp-V(W2~kmc2#KM)Q^gq+BxLA; zK?jMUM!Ml!oAW#8J%{oA@BOap{jSgd(Aj(Kwby#qllT4XXSdZ9sg5uo*|TR4)lH=v zntS%_NAB6P&xP_J_@Ayl!RNt$oV1aXyKN_@C1t;y&am6(R1rFY$t6{OgM-2jkvE0=&y3$+nrV&x+$wst zdNWyF{&Fc`tpYy7Z<5<(k1+Fj-2c7f$XJus``+xV2z2~=*N^2w26J-~``J7WF!S;bXv4zhs4LOJjts z!hu^f`%c<2kKFX*dzEtBkbl@1mOkF5e2ht?lU^oZYmUVF!m9e}0WYKSQ(5JI?;lrB?8b{2ViuzhEq>59d3skP~%a z&30Aiq*~aGoP~gepc0IKDPBr6BZuq0^F3g9HpML2&NCYf zc`mXAtg<7Ibbk83Rr>x=^Q-1d>o7`sZt7CD%!hC9mR#g#f2c&im(zcbzvL&0%=+1> zmX_Mnm2pV3WtsID=cO5ZU*Gn|R^s+R-yjKzWGz{0)4D>-A^UL8UW)ycG=KkbQRD7| zXE2-KUH=V!mp#paC_4RLSN?U&Yd#LdCFGree~$9!Z`Uv2elZM@;ySA1rjXL^Cv$+| zss67HpEWC7;5yx3V%wEv?^)rwI8kER63a6YW}usA6r){R+4scpxK?M3R&L?os3zRL zSgeb=ta8e^K^SGwr&BBvI?Ztb7jM#YBiD>2>QuTyE%dvF+i@+vQsJ@yqlvJ>@UbGJ zvVaUBduI>T*0QQv<=)XE){jHFow%m62wVC52KMTE=-`Eqbm|EAHe8$CJ*NsAJ6@er7O!Y>N#ivrxtMLx;5th{Pl~s{FTcm~DYLu4xZ;7?p^!Z8 zX(B4IHLJTIY*X*J)|8!d7W7c*kTLv$ql`ZaPNsrM48Lli=l5jy+n7Cdix;maZ13zY zcM9(7)OVbaNK1t~g)c#`DAtN~GW!`(;dF-XOB?!SXu#D<6>4cEM=qo19VFF%T$ax@ ztM*=v5P8ImgW3^~Ybr^w$MNG&JHC0xH30cMIr`a~y~!yP7-7yT#|zby+a2%UrK^i`9J&Pg&LI1uHoWvw~a6m#HvZ?JpYaH#Svg4gYA5ZNX>fs)L#siPlcqu$Z{~=O; zssmX+Fqj{m^FM6Oe-6C?J&kxN<@opAYGW_OCKzl}^~>qqN`S^%13dpsQV`v4y+OIa zL?izDWq$9>Z!z0_XstDp7?ZJBbr~d^=Xf$BSGSb+*2bz$ORUhOrDQ{{MRQ~tD#3H& zN&Q^H!PVK<+Jn|w z!cPkLZEr5~&xE|N>CVw7t`BYddj5QW2>X3S*NJN4j+Dqc`;UK4kKHe?Mx|W9-H|`li-MeXqk`-pSv> zb%o}WFrsQQuew8QEIxwG{hNOYn$b{_bJ1scFbWfTGD3ReJGIR4m2k;9rK%M|2sZTO z)eeEiQ(8Rbjci`2D|!RxqKsZW2a^{x`|$YQ)|Z&8*pzVTSiYK{9|Ro+DiSQHC8i1x zXa-$+|2toJyejX>`;t$B_R~_l8p(IL8cRTjU{}Wy6HR_TVv^)^d#qMHU52}gXLtLO zASJ@5Bnj460@^EFnpk zptmluIkQ{){y6gjHjW@XyTGW(I9=Lf(q_#%W*9twRW(6U%-Fd_i1_7d8_)GD_N|{m zlTzm&Wk%nt#ADSYre8lfcJkHeml&-t&*?&CR)0inMutcX7tT7?i_T?Djd=F_5Z*0u z9EbevWlWYM%WEtm>gu;QNg?p$>bWNFknnG4#UE$x2Jf-0-(2Fpp7B255Yvz>Q;%ap zdYzWHVE!zXOW0&?IlhW?53RUz&yt7knzo7sudbl*JN8+P_m~n82hd01=b4vX2rlto ze9*xRHecV$BlYmbI$8|F3}pcasV=o@EakY|zpC50GWtoB-@yB0DHeqp(m!4*Y@gD7 zsjRrDG0QAYIrr;!Z_tP)oUe^q!%8JvH{a4w>Gg7oSsb?Pg;s9Z%V3}Dl844Fj}&d& z>S1I&oeO<{1TtH=O}UsS(j1W9S{c^wRcL`?iP-~sV)o^;RRSo324Q>%zfr|y2iFDq zvP$5hoBJes4M&(5!b-ayh80sG4|X)Tp9%gp_HG}wisSt|-bV)qDxE{-1qKmVRbp0l z-mM^ZU-wxm1gGNjfp6@7cQJ8KqY0lfqbVv$owDH;%_e!ID&ZxsVO0x~g|J}$@+J~d zU8v!vsoLhss~x8Myi8mFjzBt?)n;2JKHk>C5N+C1CG6bK=#>S*x5$*M2EW7 zQ9_L}tvIZ>MHnrH^yKfUU{y4<2eOLJFEmmS$R~A1mybjmGETp#G#YY&O)V>g6@JOm z&JH#hw$443E~?n`^Ls;x^6^C-ZF7Lt5-M_6%QB;6kpJ;VLn;d5GdhRAcoV7o>!a`eXxHoQGv*Sj5Gvv1Nh1>>up{PuYg`+|{sC(SW1ayJ+C&3F8khnFPG9 zWfigwa~&q}EKwn)(L^z)3?Y*?LFqo)Hj6NIIHCyNJ9@`X4>K8CHF{7Vl2|m8j)ibfNATZ;-+8Dvt^Gg6|B>b$xn+H{%iS=xUpvZ@sH*Za_I$Y=R&@O1f5UBz<`9~7p>{W_OdW!+?dFNsqjbiA^9_ImAd zO@f5PPlXz#a4OL#&+IcfZ?_A1KOPOb)X#@IZ5t*6o82EJSF2%Hqtu+%9bXpc9&l3Y zc$tPBC!#4C-)lK6e}k_SGcQ+T)0T|SxVqKqlqqDg{$tLKaUGW-atKH1!ji1ji6bU! zO*Ifo?T?ogl5ACG+e*}^74G@5K(X<4l)K!(6iYLvLNcxg3NSHCG>e!tEQ(U$3_aM3 zpP`=4z2RtQedsy*r{1*(+l_78B?Z=F*f+-&i85QmGQ#fnHNVOC(qk(ehl3|Y54)3y z>V9ilvv0klj9s}b8zbb;diz9$NSK>Q5>?|x3*9xjoGbA{m8!dq;HPUz-N{%;H7%HK zFPakLhBih$&4FaPrFzVdrui)Op`!iX?L9<)b8E53v>xvSw7;?ArK-ApQ!fH!bLCJy^;?TJK(>{mOgwEIpTUdQ z@ym5l(o2QmX{@fjr&(P(HP9!}VsYu}i6-%S8%520%F9E3epPD|+1MQY;tpEHQ!kkn z#nB%ek7tP9Slh?GB3m`zDuFe28+;?Z@}8RUdnYtvf&jl)Wx|s|GLQH^$L+~>E((<2 zdbyu+#Cg?QrW*t%>{ms0WUKW0R0+X?dN>{B~Rp5T>IIM2{ zTGzr8k($%2ZUc^!RKqK?!`K5%Mz`bdlWIvm@HLd*8j7%WwN*4_Y}ivE1-K2|ANY6s z@&3n|GSJUHa^Ch#=Tc*48*>Sc;A~Blkx`AloIJ+)pq8EP!dx(C@ihcfJZ_IUOkqFy**K2Y!dBTBD(*O|p1o{EolDW3_$n!{Bd1n$y!oCHzjE z6nc^ANHM3%Yq;bTXOOr~mU@_G%vfxk><_csB%7}T1C6-L!En$(j$hT#Y$^3umHL9(F?L2NT9KDS0Dq8Ybk$M2lZOD7v@@Dc}NqWfhJ ze0S6(8we6FM;wL7h8oua^DTc~<4YDNT>>$e_S0+|vZ2bse5+i3EXM72x^5?fSG;bd zN;XuLHi*S?IrE3e*6|2~F>O+m`&x$o{q^zs>tUb0tatvmTcmx4j`2Mx$5j>C$G!{rOGu zxjH58R9PVA#n>GjiGp)mscLwh{k!^-9iO_FNbjwC&8;TmD@MwO4DT@Eh>^y^a{y zd3zXu@}p;~a{~5iegs@{_9LB5GjF*HNHR?QxWb*EyL=+DwCF;j_Z2Eu@sdoWoc^Rk z88XAAugJpY<5OB2kdbj%`rfuS7sAMUX=`5E@=LUp&tOa5<)uUFGMkA0dm+3XWjF`i zR72w8NEDC$D=58J)@K>FbJYV-3FYEzctWL&gV*~MQ~PyCnN}fef}34&u;rZm>OKY# zz~g3KoR%wchbp1H{GqH8Wqoluj-#I;3dWeF^EQD1?vy;M%I+{=^x6U)#C75&MlPsXjkCQ+fCG(q9ZKAKnz=V8ZYtD>#15IVf8eAe)RTaFavjR&k>Nw zYp^Unf)tL+yUO0XF0VT5I$U*?+eBBpv#%m24(C&zbs}FOE^Jaw2iAF_Ck%%)f~S<7 zvW0biDig3BF)rq9D}pOi;i8AuU57@*r5Vl9F-;kAXCwPG`zq~TRN;cN-hx}fyzf zQu(b`erw@>CvSh-li%j|-+95`-pc<T^ zby!htQ55dS5#C)IpOVJ9_+7hje*H*SzCr0Fq-(#8&HG0Oh<2vy2K7s4?#Bw*>8x!n z*C(=h|9mfa@8eShz!ab&)TPC!SgZ2Z_Zjo)JZO_zPF3cw6&!SG=CxU!8?hmLd`1M% zRMkpX(T8tZ@Ryqm`;x(l5CyW`isHFOl~<5%cn=~d$Ga;Kvh(YzD3NEZtH$4kbl~ty zWeZ4>_XP+WFezo3)Qr~9MkV@su3;0m?v%`eYE%_v(V34`q`Gy_eo&|Xbs?r;=62$- z;Avw5RJ&;fkbwnIzD?yV>wUXyELG0#@#7_87R`h@dDSx2{JB_qbo~xUpaS<#kbq$O z$XoxIt5)fX>^_3V?%xA#1}Yq-U~I0vX2f9byo(f>mm;}!ejM8Eb#u`9+lM1Om#c@6 zjQl3?tuosi3y+v;?gk7LTeYbvNAal!opSo*xy%m0A!t&Mi5JXKs(s|%cNT|=)2c}E z(hCoQJyI0HZei0^l|}IG=Csz@h5rpYVEU!ng#`C|o|}`Vrl&^e71?3dw%!A}uRPIpZHY61rkR%p|(3&xgC-$*g6hYo@Vu7gXp~#JXZb zU{loh>R&{__!dX+PZh%12H|}I=9mTVd84b% z@^Xgo+QoyokeIJJ`wax%lp1~Q-Q#CJj$05GggrWPX{nVyu_1|&kPhI#%ubRxojv#g?%52tzIuA@O&%rk?QmQy zfokc~OHS;Qt?H649mB<%WSC{8o|Sie9*KS$O%Io3floCX{0Soy_`Lmv0M`- zRW`URpDUn)++xqpLW2E?2WW|2S93$_{T14A+NBSL{J&my$Lc-y_L9p>s>qJ#lZoTkCi&-_S6Krnh$?d7nHWGJP zC`Qx^SwGT~T%6qNBbr@V(#84#=vr{8vNn-DcF?z2s)I{7PgNCXU)*x_Jo06~0}w2E zx)vTL>fYM~+z!ql-8%ovNHFH$%CK0Qhg#b!T3`7yeO@MQNLKgZXBQHbzB64j6H(U= z)%E5PG?}s=&!|@`er~|I$}Vv4jnUB zdVi(kiq?ReSj=oK{&T2OPpH`M-Ft41MO z1a6{J`ov0dS=jlkTs4Nn9zep&g52|0R7@Xt@ikKAiYy(TT>~;|F2k5;8N5t`l$KiR zHRn2R@Hx2~EK0cYO@Lg!cEZ!tBi)xt51$cID82a z2G|k!wC{!8>Kjo7Rl7lkg<;FliG&A=aREi3sWEo#ON$1Ygm5ul$!oogDePm-ihA?M zXLKPcE=@*xZqp8A6-6hS43}{PRD}-;H1XS9oD4&|^cd@UKX-*-Wuw%HO1+i8OYu+j zXTQm8KMlV)0FW;>P&qz`J+juw+&GV{fGS4W3BV|A%W8?zQpPS_+Qk2^)=q_}Jkj1j zEmN~%v~@`Or$1s0>6TtSQ^oD`Gk~4vdd@N^eoL8~EICP0RbTtJmr6f;%OgL&fadp|J()Wk_n4rRO8JJ(cjbEtKn1{ma$f>UHRs~>bng6` zD^hcxIME^u=$D!-u}vlpa+8gp21YKe zBKlvm-#Ob8yzO_!Aw;Sm9ADGIC^j5=wFAP8w)oZ@BFlPG{x}U??K;)#3ott`)gAZW zQIQKeuG;1^WMuwC^pTPws?t)*D@u7(kEe zT6ZWXQ>6QCFSToqF9fh}B;6I8sTdY+nok$0Y|ZM4rhYGE&IYXdi~ z@7CEIw`XPuB)*ym(9bH4Y$zr`*5FdR;tOggsa#s= z3#Qd#_It!+USh3HnlZX7G=;pAANb#mG+V7D@t<>2?%ypF;f z<3Rv~W%HW)>J^D2AKLglkWs|@15Bwv2f2kF;7a^_!g8v72qVmuv^3~8)Z&}-s{_Ib9=`B0-1}@|Lgs6-HvnAC09=_R zHyJjFCL>qYCK=1`x7xh$Ox^p;W`oRzTyXeTUI3GYJQ*Sc=nI+bKQeKG{6%E^)1^&% zNMHRD3K7H(L=zV(6^fy{vUQhviU6El(5f|_p8fOsl$A_eUhPua=nV${R{Etx3YfCLx{z6HZ zLXykAz0w%?1wiv&l=a(YQ(Z!DQ2VVkpnZSy}FM(|{Vd z?6Zy=;^PYo`lv|K$u-d4s9NkcMCJP~X4`aTXoOXhJ6e>TdP!-V$PX4Ea^+IH=sOx? zl|JZ-+7+$QLoPC&M@tZ?S}DnZ<#5GkA%&fn zk5;$^DIpTZDj^(pO01~xsU5jW{ImE=Y(EU#@I#+b6@$-Y^Ct{qSx2?$MhzFl~FsyifDfrRn9?vBoZi81J#&J+|6P^=rF5_pcUT) zLa#Cg8U{S{Vq-6-6Q_!9enfWA?ID#Qy^u_uh7Msu)o)GK^NH7hbU|8itGL}eq+eII ziO@F-UUGWxi+o^quppkQFi{6`*Y+!gv&SD#5#Y7mF*OYm2xoE)O0)Wc^b^Q+#EyLc z&|MMi&3x-Zi0Q`Lm5jmOLUXK8k5L4E1tY?^w-Y2mj0Va- z+w$y=t7tfGAN7xvL*b9is=7i!6;_T~>UAB6LDGOV7_?dW`JGstxM}Z6R5p!GS3|~b zZ4PZqNdbOEE&C`LqvD}#%l=zbNM>QVC&zBC&s4)bTkK3CmWKSc0VAiKkSEGK0U>*L z%v_9BLQHp_G0y#0^2vO7;~dqcOu1~G+)K|N?67#s>lSlMfO>s90u+T9tIxt$8^D1f zu1)5vo34Cd>SzR@=LLYcJ`Z9-i*vcgS9c57UJB%y+qcFMkX@vVt=n46N3G6}G5N3v z+`AzTJRqzDe)c)yoByHo_W|2Zx4XDaW`hwQ4qt$bXlgTmy(AJeyRZRZ?kRd|7og;J zDVf(Cz%mPZElnZ52sAZPM%X*RT-Ax(ukwgwGMnermlmQ%(jtPhN(m)B0i~`wfmDO#lKL zMe-h#Ac7pJs-vO@Gw2BNXD@)cpO)z>3KGA#Q|bV?ldz$3o=jKlrSLhh<>X4Lb0&vI z15JMbfPjy>v)$UY8Z;aQW_y#dPmI5u+r3aU*4F^0mm(pE+`YBelN~QHV{#_<2@Pnh zFCo@ICNvf>x$qw*B!6MsS$4dHg>;cLklo{}opH%mPTZWb!UAVukV;{VTxqnYzH7IP z1e$dkTjMd`0oX={MubA(Fe~NzOxuMtUpvJ5`rgT3td`qRK?d>7wVrylWxN->ikKeP z5$?N?LRd=HHv$;_ZO_x4fMXO(aO&2#O4CSDNZWT@G-*sJgQ@Okps`X;^KSX}pP{oB zM@%RSdd&Txj{-!mstUWAfl7HwAt03WQ4ocT=(jcI7rYKZ^C2&P5%z(!dz8HtLUh(* zh>6L<;pYj_!imlvvD)oymM`NsBOUVPNba`2d@c6(K04i^^y{R$U%+GYxNra*dNF z%T;HrgW-805n*xUhO|l->;e(!TQA6F0kqU@&(F#x9duVz$uliqDQv$hx*bc;r?0hwu^}e=XWx*@H8!HR_J#}TG z6GX_Uwt2KVGsCKXVURWb-?uERYbCDTsRFuShiHQMD$kTwzg}IhtQ^U6Gs?*6rualg zN(W%FbeMQezmZ=18lD~puM>l;UCbyiSz`l^+jRcv&k?o@Vs*$FP=l~?)CaNC+GllN z&3Ge8waegal%qg$i5IF(K+zvpT`iv)A6wMTi1R|R56DOX9JkG(@%mph0B&RlFq;#I zr9P`20wv^>_>6Fy-c)4sxX&!*iKrTlwd<(HFXllKi@idlkbb|`2v05M1J{J z+2-Q_#ho2sQc=EJ^NA?#yoxiB~`GNME+LRnDnpi=ZE4EP*|!LHCvHd9*BTUUdmPPm7j!R1nm{H1*vvfO)f~d z3!!dE&H$zS0dBv7SRIK%K?I^}0SD8-7{!#)XYL2jivYHs&Z^z`Uh)IL1$kd^FC4A| zxX2;lU#5Z8Nuy7qB~$MHYlneELWXUTuAP(o=MLy3<<+fCK$f?K?~L>-5s<2YauW^) z-!f4N6Z+PR6I~F^H}g*>7^dkBcZhQ!+U3ubcEb;-hro0)7M9W-O{h8&QnWsl0TmM9 zJVBl*twu{`VVm86dskEtCJs`{TMo{BjD$)NKCz83d87c^L)Lp$Q*e4v}S(lIOk!$}DyHbY%lh%@Cnhu&Vg;Bsjy6ZscnvmVDMpH1h|nBf0T~ zW5O0KGn;SI9jhIUyT!9?Q5sgmRauJ&&uZ<3`mdJ5Rah$K#gW_a&Jz*@fTgRzxIpC* zZXDiFf3^HSbVkGJGG^R$6r4J2jk#wec zs@K6eO9Vn$K7l!!U-boBBh(+xd48$iE?*G9CNVjCjnxw*J!G5pt>+H_?*}{oZKOF0 zwL^iiAigc2L;}{`R_hZgmdHTZ#@H25ze?#|8SMrD&=r!w269BxZrHlz)P!-Fy}n~A zD6e4b0^P60$khUzBkXh#m~JoQyC4(KW9eC~ukQ+$gys*fL{l%G5JuJv7zUXlI`GeqxM3Cxp zTx($8deJIDs)1An(4sj2-JsnW6 zxB%35k*Fg9(A?r6pWPJ@ITr>=zU9ku1LS5UBNJpkwBU)LY%GKXrFHpF3wds}*(;dm zPz*%0l3CX$p?`g!){uz^O#n!xUGy1d8M??=7HStd+Jx?>Q5@sMUoZ7X+9Yx zy>91vv{OUaT=v^`n_v}Gjk$7;)$t=O8ssHhKgZu%CHG7AG4xKg6T|&B%+=c$psklu z6s#SW=btznOOR9}cU&QV=n^7U!p#9(4+=&l7p#L*qVa3p<&0)n5Qp~p1beEGyR8XO zAPLz{kgm`byzEH$n2$3LH_L8Gkc-+0lK>*2C?+Tk;%gG2Y~)&Bn{RjLFc0~X zyg)%q0;TavfZKuUsTma_g6s(EgK&@qVo#@2mlQAtZ2x_}p$Etv_c;cz3ATv&(Kxso zxJb8Za|y>%Y((Y*Q__sn_U5$P0FZ?Wwba7xolXH8bIo4NzsPSrl%=%S*+~&|Y)%ht zft~^gn;pOBkj%<^(07aZE%`d`KLV{dl)zfZh5K+Lotn7hIYH^o#B*@cB#Mzo@8KJ= zU3&O#l$A*D07bKjt6!`-(xOaa6+%zE+=-8MYX}TPkb1?C$3E>N6mc%8TJYi(_dzI~ z3ks6hcroV=;CE}CI?&~H@)ud}D&?e+cCHHOygJk_cJBnDe`Q{L2y}@`c0GC1e5$qN z8+&p8t;5K-X@I+(4w^;eiM9c+*DBso!E-)tC?3EuLdiIivPsSRZyiN`$aPFWHS<;H zh_>w{3!tC*Nj4Yg~dWZjP1i%U(6ID@?K!Ws-KQVzm1 zf}}c=(SvI{sU(eRRhgzgcHIT0Pzi-(Fg2@TNp)1bQ`4wqr#ZPRd4=KP;tfI!^F=ZO zs%+AI5XgQ`isd8=SCwv#-{&R^_h~K~117C_=0S3if2>wUW<~b+cL2U>T7>NjvgMK; zzt8e}_}_~IWahWz{61EH24TOI!f%cFAEWR;rW8u{c7TL83MB9Bza=bM@fw!#fNoe? zlWha<-~hbpuN^JF+Z`Ui=%380t$F}G097XBKnt#ehR;{mN%ek?rjJ#MEn9AZ>c}!S zfLb~l@{5j3uY5zyeZGK7mg}*q0>z#kkIA8OB(emw`4j}cGy-%bO$Q0}obrUCD8qcm zda^{ixZd4O!L$I-a;RelbkbT_z^bE&SQH8*Pjh>|KwLrJdeFTjq2Z6UJx4$OEt?ILfBa-!6ZJJQyjcnX!KZXr$UfdJ!9~N8^C8N2FUC!=wq_pPOMv> zQ3j6qO{j%SqN6C=Zp^Hn11O#<%@$TD^HrXO@_05M9v?<;kO5|!{{fhd)?ac$2XBx+ z@yzXYXL0#)s7DG6A!^x{uoDg(7yI#IF;W7!VIf2LpjMo)Dl+guE>_&dR{a~jTPU;8 zB`6&^=sp}Kf16>x$%L<8kgUGsL<=Bp^ki+bRK^nheAm8H`&HmqZd_*2;_30)sl)(I zV+Zg~gZmsG&4PSE(U6m)m89!P(;K4yAvPro{+{ad7N0JIjZ` zpNG!BMJpOf0JLmjv{Oc@8B3{7A@Yrsj@WbqGw#|{LvLjDIYE2fe3O zRi*P7#7z2REoxO@ztV$Bnz;D zPQFACr17;@(G%p(WV0=O9 zt2qbcmNSH&cQ=_C)VMn`OMAOJeon?`*-lA4vB>2tYZHy-FKbgvjO5YR^u8lh7DyQY zwru`A`Z-8-angsw=j`{7Ag5Rq@1wu+0?_Ir=4Bp>xNR9Wewh8NV_&zfB6yZzU=I$Vv0rr9G357X+Dd! znwl0$wZ=HN){cEq8H!2po%Co$#vJd`hA$#xw$ii2HD9}RcnyC<-&i=FF)Z;Zqb(&X z=OuhYCS>`^O?3e&MY<;DA*-yHth4$$Vfg6#Lj^ER{e~-bJqttHK7WJ_OP#J8;?uvf zdOF~uU5&^Mk74F}9a^DN&tB@OUbOR0K)VIv+wkaw`Ut;+T@Nt!pq&{}&8^08L@uIi^}UlqJ*&6Od8Zd`vva$RjeRso#@)tZ z{^|r}&-iz$TyW#ZE}b&1fh&aDiue@Qr$rcfL+^J&$u%0?E3dj=7YP%Mo*w<8q*@SG zq;m~*{u{U-P^`6E3v~vr01(R)15S_QbO!C9Njr z1P%Q8nsBt0KKMJ)>Z`xmiJs5fHg}(HhYf*fu>+h+39Gi$cOiGRQo}^>6_S{<-)e~< z=f!Q>L~<>!euJJhqou9en!abj6GB}Jib6Tw8KF+q_$wY{t(>eMy}NU02XeGf81BOZ zmLPyxhq59^t-Pwhk)!MpaWYAb&nMA~9W}B@S>V&8fCs zA65k57u^(WVSR zV6c&f4lxk&?_@-69f2WV39_L(jr5>uBmiRgT$Gx9rmg}&fa{@_1NS3TJ67Vbz z+2z-EfGnDRH#_M1IQZSO@NOrr<>F(e-J!YvVW#M{XCKWwen$6jomZ@Z7%!|d4}xZB zD)42`p4Gf<)ZCv45k{_jUGa8n=(UW4b)cjdGoeIiV$!6StiDSI8J5rat z&sG@4BQ+sH6IrsRU1#q#R9RV#C)rrG^fZVpTHS84B`*W>EiE1CB1?1JP64kNR#8Nj z1uoGB^UdYW$CHi43|`Up;~Syf3z?J)<_m*)ERv1Yd68rM4rg8t*)qZyXb`q`&DrFw z&lJEbGCNw3E#vPk`Mo9oiV$RfpOW9F8$+EdR%&% z=J+DEfOjn-LjzJ10CcUHu*^4^7ys@8Qtg6#Oq?dZigGf;uNPy=>Ji?bEDn#;G$~7w z-nSKF8h=z(j5^j8FZ@-wH8i~@NNBnS77Q!#(}L;hnFadQRjPg~p4iq|6!n*%V$Afc z4?CcG5guRs(zDr^QRSPAwtkAWSv4WJWjv9lOjVRA>pWc3BTln>VgPo_W7x?0mVQe7 zOV90av{%DD=^mX=RES9lv(HU3s!XEW8}o4{fd)<68)Gu{k=~P^79hL;YhTkY=>D;D z1pt$kF4V1O{@!DpSaLVZ(!S#J`8dQ;&&|c29k>+?baWrIh}ih?N}2B_08-c*fX4S1 zDx*uIa?7a{T?#<7!o#aybYBC0JVkuz-u$y76SqMZ8;~DS1`)&3JvY-CoTwEWM)~|W%f<@!6TOnwoFOu&0(wgKnUn~nZE#L>#i#v^OFEn1~ zi-c$86!O#!W4s#9qE>viKOE(^zoXx9Rl1pQ71p8Ei`fZV$K<>Zg^=z=BoLtJUN|jh z7EjZ2+>uIv$R*8Ah%G;CQTzP~OOsn52>Gn*Oj*06NQKKG%)9Tbs|+SapSVOCOqsrD%(?w#Nm?? zrrfw|;q9URRRP_tzPtVXy%grOGzjLS)d6G`;B-UTcAF}k?)1O(EfhulzUB9d{T7(r zDgNK*<+r-{e@t8IU5QK)zxbE;^zQ9`Kz$6Cht=nVNg&m#R!8HSEtMlfx46H9j~bNc zoYDJF4o#)<-FUs3idgi-U$rRA=?62FW6y^Rp8rM0JE8aPXO%y3ikyaVIq_i<*855% zS$>Q17!Sv|;AJE$Ste=!C>;$VRvqK?&ronDQQ+?n@>K`Zi^`SSB_j$;gA<&%A0~5F_E6Z<`S?oI>`LER4&eg}bfK;7L8YCBZOCVJ( z!iA6jlG(qr=XaHW?3~^^{fl!2IgDk`&XJ3#9GD{2VQV*X5$ywKq*#@Qtd2o*LK6(8 zBwJ0cCB87~UGaK(ZupOinX1Vhl2SLo?>aKPcl_%uIBCIxNh!Qie{&Osy(u!}?QBqV zf5Fr=9DEyV;xgNxm-*8k3Ek52dK_q{tYd|A(E+^vh#gwjABcX?-=_n7r6m!f@s-hJ zT=?^BU$&rt2(7QXlwAIFiT%+>PMDPI71F}}T7-Haa>&JgfIy%T zzT4|V;VaffWN1Sk_PtWb$ng0ZMKOr_3!y7^_=Hg9%f1|HA_4VtAp#1E`(dS)>S6{h zuFDhODgEgiXmvmkSwaYoV7Gc+grNz(IGNFGxAO7MR#eQ|t*+Swn(f9NAO?tzRBKFy zQZF4bXe`KjcIdY3e>Ooj0BBCkfIHyq=wdf#-0B_PaRnD+nC(&leM=ewcq`?5c4c>f z?>~u;YKm{n(O*7g=3y(*ow0v{GW5pnvlPc&_AA8FKU_$-oK<{B`tMzDhQz=lWcd{7 zZ3y2jfbZ6?YBrwfC|>AB74lGMIu;g!&oQ+c$P7I?k*Av$=TMco@oKO-$3Wy+_}7M! zCoB%pAK!9eUsKksp!wN6RL@-g!=|&gNWr8%-i+a6W=an(5yZ9P!A6D&Dbh@ra2Fr2 zUMuS7xUce1miEM=TiVKTowQF9c()5;+KCn^t{Ew<+WMR%5RAGKWc-F3BhudpHL-|@ zLZ2m5?z^y>_!HD?xp5UUF$?tDJ%sbRpe(z%&C_Hkm3V)ssqC;tSX=Q$E7Dj4AS1*R z;K?`JZIU&G?A-uKC(3B7%5c!YE?2u-8^^Qo(-|hyQxrH_UQ-?tAzmk|w10xe+Td>L z@qHfnOW%whoTuSv=Pm&3=Eb`Bl`*Ix1+RsDhZ|;31^5V=onr63BN9F-LChy zl@nAi+3i^{sp+A(@Mw!NTPAc3xiV?o?H7O3lT~9u-1F`YR%!r8c&Q^~#7SVOH9WuB)12Y{ zT0|4)HBA4QcDup#W0F)Kdb@C9TfdF)%$s>6r#0~W0wN14<0>?WWY&N=Y1ya5&l)Em z2H*6WEDGChNKtud+S@&gbQGw6!k?<94c7AJHk-~Hkyc##>r$>VP)2@V1zkNXP!}%p zUBaVRI5Wpy?(RaSC-vl}x8BrhHVQkAvv@C-2T0!*xXQPKT}9~Du@VqLyN0u0%YCY} z^k>2aV!OrE+b(wLvU?B}U=KidyQt)cw3kDbNU5k3n3L?ERJrVKYz<3BXES|E)i-rS zGn~8J0V>%1I92)A zs4r!o5+)og*rY%5%#6*EFLce>s@q?ofTIcYQ;3Q@`dwAH@BV(V?l}1x_K$9X6BzDA zHdw{LuwuhX9u1fM1*u+|y{L-9mh?CWaU*n81K#oUDXlviGSR$Dly3 z@2z`c)T`ykP?K#(TlDEBU7x=2knXKJyxC9s7W0?MN{W|111)QAj~s{5d>yWwQK4@p zM8$#1xqaw4g!OKn{O&HeIG{4jJ|A>Vs)XES1+Ip%x+v2rJMrj?AsO z4ICluKg>P7K^GYl%qj`gkcpvJTa~O5JS_X+^_xwpFaY0Lj=nWg@VBn`%Z7aip30MX z7Ee;N5kQ=W^U(P=PJ2Sc{#RqCgHFAIkPyl9QDoG&yQe5-Ohk?^)|oqi&%Y6OsFFWF z9BoF1$38%%*Ijb|%NzYWR*OCj3|@eF4S(!|qu*wah>7oV+uG^(HcZBLXf~1EY5fA! zQUjQkHU6?Pf4_t($e`HXQs(-PnfNb6dD1nolvK7n@?;=>L3vPr<>F0?|9gb`ja=x=zfJAbfCCiAhP1+lb z+jYL5pB(va*xu3~26^gfDAm1{DD+;qiP_JdX@%X2Hcq>&P?=7;Eo%;FaMeOs$8u~I zq{&M%dxSyS)5om?t12>hr|yD^r+#0z*Lk!D&1Vk@&#|WN(!Oz6_Hv+ET@R&rSe@|? zvY61i_T3LRE*l@Zy?-Cnzy1&eAz7veEd3+dNqNQuV#1a4ntwSoe|xI3^N`s0Sl%Mz zHbuLGGxQ98<^Wmf_kZUx@i+VdMTgHA5E3KwN)%obA&&^K=9=GG)Iw}%)2+EE|K-F{ z)}X>#=+m`O_a1R%s~{bfQTa7k!&Rf-Oo6Dio zV7*^MOhWqb4tswyd#f?aoB8i9mqDzwHQdlvWFCNv;d!LCWx$Wnw^`{E13SJ^!e3W$Ssph9zHv*0$lV5@E`W-0h;d=*UOH;Ja@`XUDkJ!+r;31E%bi|5v^iq~ zKFF$UXkUoF5-v5(a?B1lIS(41oPjAw+!h{qUHHeP7o z@BcKReq)4ReOCF#{#ds8=*Cp&5j(Cv)K~Ns^oI7;6Y%Jbc4fy&dZy>(hA2c4IGiCf4(K)~?y#fPU zuv4X(C`9>;TrR5>p2eT?AFsF=%)WixqQ~?XnPlF6q&9~7V?xdDdT;e)EL<)t+=jmTT({Z+k8~EUYk(EJKM3StO;a3n z(ra3NvZ9feYjXle`T8UEInYEazXNa`H|93ofokHF$-=k^Bw8Jn$Tni<`C^k}D6f9u z*&ct3qcxz(#gycf8TF*mgZBpp!{+7H#x%k`?~FkmiK?0Goy`I_e|%j|GPX0XKBP^;HvpP8hujwgVoByMy0 z(a2Nov4u6N{7*`yzCGM0E+#to+4NMky(3$r7y#ooxPQkp8`Adyxt6iBrvwh&vZV{Ll)bf$*hQR)q+x_Z!VOBHH|Q*AF$ z!{p<|GCIi0O-OOWh3>KOkO6~n^UsNS5PWS{ym_k*c zT4q3?i%yMcm#FSn;N0`d&qrP0sN9!Et4Fvy#}01t(J+k)^!Xf{ z555`g3`f28l791@0L8){`LKDKG_y9`*7nM5Edi17+bvCyjHc3$8!;;%{X1r+V!!CYtW)s=&0VbyMhW6qQn z($V`YD42M?m+*&qYR=H=nNV~gCkP>Ce!*h<&uHTo@|V|i;#RIe2C$VtSs&PoHo6VW zR{=|RnM_Y~sL+fW-6$+U$@_KL4q_GROGR&b*?cRRQmIso%8q64OUh^r`u3LqT4!6O z@*{Uhr~qFRdC*YDyfs;3kf=7YDdw;TfLOiHK%{>nXvX1Zml;%;`^eRd#hI`H;wT=b zI04583!-6oHGyy)Lr5%QM!v6)NVYaGju{MQS)>#|JueIg~fh zMo8Fp0?h?%5~+Na?f*1)o?%U9UB6cb1`rhyr3gd>rB@XY2q2?~0)m3{4n{#uC=nqE z2&mN19J&UaLQ|U3LI>&6LJ>jEg~!rsSo3aqR_ znrNmsqksmJn6jbnc}(TSo(Pkt(v!i>H6@yh%ZSSoK?Cu(5L`-eoCgrGAwzK%2D z;<`SSbl8dsAXH|OUbDGuG|R^LA=6Y<@@<5ukX-_3zZ;Rn;+N~6-`~ZIeC<>e1eKa6 zmrmGnDc#_!9~PUar`PoBDE-3v7;`3 zVcF-hPLTZ5#7v(kdy^8bNF?1QNlEV@ z>`i@??YUCs!P`BL^&EJz5t=&{*goP&4^fkGuHL*lrBrYC+<8}Iu%dK{^a)#B+~80~ zSiR3A1g2TJp(W-|{IN&PR1x=N3Cerbga%5{GH>=ekOcMcrV*z3u|;}FI}p6 z?+~vyO~Tjc(VC8XT%+l5(%vB5>1La;M%oSp%Q0stYQ{OCzt1ld*a~d4K-{0^i zdwEJ4CI6msBq|I2@(W$;*Awc9K}}XmsnucjX4-`gVls5DK{o@UzXgIjb&5)2rzp7*umbbK}aeLf2PtQh3w!kVC~1#)cm9> ze#8ZdQ4T`p-=~cO7?l0YUAGl4dU%qn{?qr&Pq}KV`KE{iQHk2i?veTicOYP@W{vQ) z*R!YRdS!=fMhXBR90J)ulqe0!u9(ajP$!`U2C`Y0+V|j^fcnnJQPt0cd2i}pYj{eM zF?JKfTVi&{m8(5^Dl11RP&L1~A5_iXTK!o@^DKuCvc-ethBBmrB8b1W$!zrMGD)Y4 zg1zv&Zz$;Whv-$Uoj>eSgPRIvN57yrzzrj}9{Llu%;dhCfSHap?{YJzGko;3)ana3 zE1gCfKIL=J47|^>ZfB7K@47*DD5Y=0Or^>UfE9d-TG0w8YzJ}5iJ;VEYB5RrnSI`c zy~@FosjW{Xo`4pq$sEt|Vl*x#8XMmSNh8X*fztV2*%8ADi7(PBKY`8AT?BM}Mp-kG zHhx_3tAQ!9M?@zw!2(1Vvq)1Ju7VnE#f1oMl|AZn;26C?hg`57-Cr;50*z9ZFV|-C zmLB?E-~9u|dbfeLr&1W@0<-+;x>;?B!S!LPq`bir(|!wK7L))?OB-U7nJQhdtKO6< zk$$heJMaCpw ziw%1bM#^TGspshd{$OePoa}7R9SrEri?TS`!ByC5-%g>CekOpGzxCgBap7?LFI*-eqx}5EJh7gW!NV+=0!@%#w9?un=B-TxPF5J>EMk zO}R02R?OGCexAKzj1C&A{&d@O!Mj-{Xu*s=cuS4X==9OI^hz_TXk=#=f@!Z}&Qr(V z%I>mdkq>7wfshSbl{Dai-%dA4ru&YWt+@Oj*#R1i{@S-{E;18Aj?NL^~(U^XPS2}_T&vh9?>I%!+;Y&5Y z=^ZzI0)9`{SFrxxG{8?5x%nbm)~JQim10V>cXXa8BEgG%#gocnu4JG4YAWuq16`|+ z;9Bi2J*Rf4N5jIF`JYCn$Qe`N!RIpv;(ZHFlx~UzJ;iV0i$KmgI?3K$J&HjKUvq4W z@6#3W$otGR&6>OoChDbzuI8o{Zq(cU%kxfti8G9R!P;I6Gh#bH?8&}l4EY{C(f5-9 zT=>UxCv6q42Mv<;SR1~I-iuOVHEYidtXt!6R<0ew3A52&w&aK&M(QP+wJ1w?8G($hxJfA}5^Tw`l8ZX$vmWbO`|W5h6WQx)W;WCCl&uY`g^!VeK-aLUHJ4n_m_V)*iTFr;A>RT_4dNq*#y70!ooTLbQ4 z6>Q}s&)r7P=i~9^z(Y}hJP~3@fxA7_JS{#r%{$sD%eA9GZs?D@eiV9xNf&6a3Nixl z4@uBMlnSPSu!3M15H%;u929!~MReckZ|sE7##FGUT~njHg&2OHOU!3pff*5!N|P6u z??93NK%bZH}9@h7pb%$FqZWl7ZrfPI%?^Ap7U&;U~%!~BITB^w0Bjvhq zp+0q$?T#X_MWs99KwupQv15Vz)A}yf${wjqOS)V89uqy5b8Kq#L=z6SuoF zGazN(up&_Dh`Rc*IO5PNZx9Ff)*_puJ9IbL6<o+ zz%E1{eqk(%0)4ud-=Ol6vCy#}UA%Lo6k!a4n zAYfM^HDrFYNAmZSMlPxvZK1nW~zD9)K2GK#vV zv6{eD@R)iA&+FQ9K@chts-c2vg1eqv*%zRXb=8cQhUu-%%JM!mog}RhJ}qrO5iF+- z3CKWeXVAs3-L%RdIdZXQou1qAI_Kh<`~#4CN@rMsj~!+HU=VGMD-e@l!0DNBD{l}; zcUg7^K%=7W@#ylMVM-z-aDsPJq)A5~ob(_L4(F9JOks}i;lLOC9x0_MwSkAjMnUS-sm_o#wey z>Zu8pLC?+>6?xEuTW>I51=x58NH$>QL{)3Aou`0rsnrO^ z>aSiSfE`dfB1u}iYP;ty2^qW2LtDO4z-J9WcyyP_`SbuII9{T9cg6Xndh9p_7=mal z@M_FugID*!d>1GAQrmv}c)gBQ$EvY7wCxg_9~JL7x-(MBF}fV$DkgWb41AIA_jjrf zh^0T^P-&94J*(ZfA-%xAQehbyN$+k$@qqD@O6O2-?J6nJa8Zb9npl9W;x~!Sktbd0 z{?3k0a4#78extjoaMvT96+=$C|3lDn($$$1{-;irtr;13iq3+}JoD3XeK*Iqv~3dH z)d__XxIfP|ecFpIjWP?;W~01^6{Ea9XtpS3yaUZU?_ zd#WJl0s1ksC`zEAg3W$4Us{`LK#PEB}nQ9unkbCz9{RY zJUON-IGUbE%z1YkE>IVG1`)D1(eLzfBZDAtI&?jBVS+nBRDWl3tRtw$I`TK2siA-a zrY)3P0>D5kU&{YR0hJ{4@Rhfrzb1({UpQM8!JZQ=C9RK;n- zYLcWnW1V4pZ*@kiJ9aIyl2=ptkm3iT#73u3 zv$I3wV&U8!&vxMrvAM~SD~KJB^&gq6joG@i{#p2qGDpkXS+BHNeg%%Htwjrx7k9zJ z!U=AC5>^d@9No`_ia+15YRK~mMf!EO8XZ6g!W=yQ2ZQzlVVxR_D0L};jbB-8t}Caq zwDqGL7nl{)!qT1hz^d*VmK@+O3H8$P16AJ#znzAUyfm`7cd&6p<;IO9K3-yKyv5Kp(FeS6?te(a*=~6Rr$XJP&qI zR##+$ysUFnKqXN5&Ki9}9_T>OUM%ibdf{5Tn|6$MZKM(vlcX1MX^}gIe0&}@%|mfh zmw%Bn2xg5xia!rv8qVOQ*Fmt9cuh_rL6sEDAm!tl$euec7D|XVF@{kLE04Bua$yc9 zF9d^V5hP_xbu)z$_H2Tly1x3Rsb;WCyfGj+&IxQYx zkrb7g2*0lpd#UPo9j>P|{Q{QgOXJb=zdrVEsI32voz8Zk_RgI+`2ce>G+`#-wD2E z1$KJy{h{f|P!y|M#W3^0zDhQ9Re0|F(Qoj&oPmm+7(uk;Jb#&fI_OrXWx&1;>)Cc| zLCZwK#4(F$i{txAdO4fDrB(VR_z}siJfC6YJcGg-Jwi16gQq8nc&#scL=F>@iz18% zpeGDMRNI0@!s7hG+oXa)+FG{uMczg%Cl~IoD|NFKs)*)4T72tc5t_a5?x&hW)%aE3L3de_NL2!vuf!$yNUwtoTmt@eB-X^YC1Pn zEE@(gZqb;8?NWS&&CAr7F44A+F_JkN9ld*G+hPfNs%}dc>#CNxc)gi1#>gU8QDpG} z$3k_{iiCJRBDcAA`_Z%0tGfCN=h4U4GOqQd{(|aGnv@Eg;*EQ-#Z$(dI4DKxJqzd);?CQ1A8s86Dpodh&djNoRb3p8dIGtV z=DFABdH2xZelv`jF`&$Gl5FGPa`y;!D&Q1D-h zUwl!`WY}yuD-~cnne2S&6bq~3LG{7LfiBFXU+MYK3DH;2SYiUbjA#v$2H*+rpN1Fa9MhOiv=ZKrin1AG(AF2;c$uBp<`zV(t|#CC#-+_j zc@)nu=JAusQs_%9L|=foKdVTE?v-mgJ9=lY;u)uRy3AA!$hWFHbMx7;s0!fROR{c~ zEs(QdjERJ&w7h5HYl$hsTodjOEh)*rbqNWy=jg=?@j9UZ)owzBdV~fNQ<3;zP5qK$xJI1YKSg^L`^n43tu;rEE5om;4^z+tm4%7&}Mb3CO?2w^R7q2Rfuym?Kn4)g-fuWc z*ltL&zq2XRryb?T3)BDU-^xI!x1SImuR; zrVQ-Aj6PsurTbaa4y5{J-h?9Zu1k*dnVZjh?jWJF-K*vSLj|DXYJ}p3CN#`#_Tv~* zzzJGmns|#>WLkXg#I^ds7=A(B9wxy=$&AOkJDduymsic;r6&0i?u}_;4Rrl%^lSNd`|3HXjtmUn;v8Z_TIZrg+gm6sAO7jO zT8w$dD(tJMUpY|DIVFfGT(}H(06=DOmj!dVzE!)|AM}n#n;1xD@9#myil?6ijbOIg zaC_O!@@zA#-S{Z zfiT0jDafp%{)rQvFMA5cm%V;FDHQHw*lbf~)-9@Am2F+@HX=I;xX^l9DLme_-$VzJ zG?w848~Y&f6Af(T>ed)7;y}iMG|lfUrvhQF_`5w)#2H3bq=I1}TDlXO8q1lj17EH) zMyIs)ThU6!bJZqrEg#7H_NaMX8NK6p3t!j%O{Bl5NdDWSUz2o)_n{`QXF=HE3OnLG zB^HLG)w?Y2D zQ#)}!E1Z(m`yCCBTmxbFT+p8LsSGiaYv4Qp80XJ&1IHkN+irGJ{-pxu@?xRH25jC8 zw5ncuf9+-Bb@y8^bmceK7SWxCb@@}h!w-`i@jPG=_vaTjDS0XoGH;UsXn7F~HNKDl zZZufBmb}rKIG3)@0eE~Z-`rcjIu~d@z;lR%rVSe~j+z=X72L3!pLS_irH`C)YS};2 z%=aGNo$IMj1|B`uINEQe)ufMT1=1)IgKElz$6AyuZ5WMFui73sZyI&m-03T{VH*A4 zmG+v-;b&YLlp}l0$ae8^PVsQ|{veNWwZw;2In6iDCUEgzllw}A{GE&r3#f2K5h@@M z+-G(GEGYyg-cz2JJ5B@IPC@~_Sn0gHi6G^Da;R^W)`n@WC|MioW55l zsNC_rlA&p(%qygNh(9e{a88`v5%cHh>9lW^`sKyVtm`jUxV45oR_?_XqVgb-2Yj8L z@8Ne{U-LKxF=*WM;(aB0FgTe8gg z++ry=4bK+CCBy+J?=fxFvfl!@{BVU`ZlQ)l)ZksTQ@;e)aaG^~A zcgcIM9{INv?f!s`{H(okX_`0|caxa|qKIVZ;qk=neKzI8)?LfL=k}cwI)zS)2VoFm zF*U{=POr`O03Q>iusB!Tp6n$3YoV`t-v!F^t30(xr{i&?<*HxGg~Su0%=+ZTR{6K< z?tTGx@V0r$Wvjz0x0#eiYz8F-7(=(hSJbeaSPoXd#`bdVCF4iMz*my~B{h93@Ta!C zS8?|RUo?-9uXg>LmF9Of?s*>UpT*|(LHShPufi||Ndd}KF>_KlBuU(&nJ6!{<9xahjodyn;ixEir1%Gqq6KG?81 zqOP5Q{;bt$x|Z3i0l6t!&@2A)(vSf4W#X{tTl)RMp{w_Q{;1|6=qb5tlP&bGc|A!# zp{iSl>-GO?uhL2Ze)YwYz<% filter(Barcoded, PassedFilters) %>% filter(BarcodePa titration_pass$Filter = "PASS" readLengthsUnnested = report %>% select(ReadLengths, Barcoded) %>% unnest(ReadLengths) -scoresCombinedUnnested = report %>% select(ScoresCombined, PassedFilters, BarcodePair) %>% filter(BarcodePair %in% unique_bps$BarcodePair) %>% mutate(Filter = PassedFilters) %>% mutate(Filter=ifelse(Filter==0,"NONE","PASS")) %>% unnest(ScoresCombined) +scoresCombinedUnnested = report %>% select(ScoresCombined, PassedFilters, BarcodePair) %>% filter(BarcodePair %in% unique_bps$BarcodePair) %>% mutate(Filter = PassedFilters) %>% mutate(Filter=ifelse(Filter==0,"NONE","PASS")) %>% unnest(ScoresCombined) %>% filter(ScoresCombined >= 0) refSpans = report %>% select(NumBarcodedRegionsPassed, PassedFilters, BarcodePair) %>% filter(BarcodePair %in% unique_bps$BarcodePair) %>% mutate(Filter = PassedFilters) %>% mutate(Filter=ifelse(Filter==0,"NONE","PASS")) dpi = 150 diff --git a/scripts/r/report_summary.R b/scripts/r/report_summary.R old mode 100644 new mode 100755 index 0e81bea..3699b7f --- a/scripts/r/report_summary.R +++ b/scripts/r/report_summary.R @@ -32,6 +32,14 @@ zmwYieldVsMeanScore1 = report %>% filter(Barcoded) %>% filter(BarcodePair %in% u zmwYieldVsMeanScore2 = report %>% filter(Barcoded) %>% filter(BarcodePair %in% unique_bps$BarcodePair) %>% select(BarcodePair, ZMW, ScoreCombined) %>% group_by(BarcodePair) %>% count(BarcodePair) zmwYieldVsMeanScore = full_join(zmwYieldVsMeanScore1,zmwYieldVsMeanScore2,by="BarcodePair") %>% rename(NumZMWs=n) +g = ggplot(report) + + geom_histogram(aes(NumBadAdapterRatio), binwidth = 0.05) + + scale_x_continuous(breaks=seq(0,1,0.1),labels=seq(0,1,0.1),limits=c(-.05,1.05)) + + theme_minimal() + + xlab("Ratio #BadAdapters/#Adapters")+ + ylab("ZMW yield") +ggsave("summary_bad_adapter_ratio.png",g,width=20,height=15,dpi=150,units="cm") + g = ggplot(zmwYieldVsMeanScore) + geom_jitter(aes(MeanScore, NumZMWs)) + coord_cartesian(xlim = c(0, 100), ylim = c(0, max(zmwYieldVsMeanScore$NumZMWs)*1.1)) + @@ -98,4 +106,4 @@ g = ggplot(report %>% filter(Barcoded) %>% filter(IdxFirst == IdxCombined)) + theme(axis.text.x=element_blank(),axis.ticks.x=element_blank(),plot.title = element_text(hjust = 0.5))+ xlab("Barcoded Samples") + ylab("HQ Length") -ggsave("summary_hq_length_hist_2d.png",width=50,height=15,dpi=150,units="cm") +ggsave("summary_hq_length_hist_2d.png",width=50,height=15,dpi=150,units="cm") \ No newline at end of file