From 63a7932d69d168459a562a25a8302485b11784c9 Mon Sep 17 00:00:00 2001 From: David Shiga Date: Fri, 24 Aug 2018 12:40:28 -0400 Subject: [PATCH] Add ability to sort bams by arbitrary list of tags plus queryname, and to verify sorting. (170) (#46) --- setup.py | 2 + src/sctools/bam.py | 110 ++++++++-- src/sctools/platform.py | 96 +++++++- .../data/cell-gene-umi-queryname-sorted.bam | Bin 0 -> 31074 bytes src/sctools/test/data/unsorted.bam | Bin 0 -> 32984 bytes src/sctools/test/test_bam.py | 206 ++++++++++++++++++ src/sctools/test/test_count.py | 34 ++- src/sctools/test/test_entrypoints.py | 89 +++++++- 8 files changed, 508 insertions(+), 29 deletions(-) create mode 100644 src/sctools/test/data/cell-gene-umi-queryname-sorted.bam create mode 100644 src/sctools/test/data/unsorted.bam diff --git a/setup.py b/setup.py index 77ce3be..306e26b 100644 --- a/setup.py +++ b/setup.py @@ -44,6 +44,8 @@ 'MergeCellMetrics = sctools.platform:GenericPlatform.merge_cell_metrics', 'CreateCountMatrix = sctools.platform:GenericPlatform.bam_to_count_matrix', 'MergeCountMatrices = sctools.platform:GenericPlatform.merge_count_matrices', + 'TagSortBam = sctools.platform:GenericPlatform.tag_sort_bam', + 'VerifyBamSort = sctools.platform:GenericPlatform.verify_bam_sort' ] }, classifiers=CLASSIFIERS, diff --git a/src/sctools/bam.py b/src/sctools/bam.py index 9c6194a..fedbd5b 100644 --- a/src/sctools/bam.py +++ b/src/sctools/bam.py @@ -17,6 +17,8 @@ iter_cell_barcodes wrapper for iter_tag_groups that iterates over cell barcode tags iter_genes wrapper for iter_tag_groups that iterates over gene tags iter_molecules wrapper for iter_tag_groups that iterates over molecule tags +sort_by_tags_and_queryname sort bam by given list of zero or more tags, followed by query name +verify_sort verifies whether bam is correctly sorted by given list of tags, then query name Classes ------- @@ -24,7 +26,8 @@ Tagger class to add tags to sam/bam records from paired fastq records AlignmentSortOrder abstract class to represent alignment sort orders QueryNameSortOrder alignment sort order by query name -CellMoleculeGeneQueryNameSortOrder alignment sort order hierarchically cell > molecule > gene > query name +TagSortableRecord class to facilitate sorting of pysam.AlignedSegments +SortError error raised when sorting is incorrect References ---------- @@ -32,12 +35,13 @@ """ +import functools import math import os import warnings from abc import abstractmethod from itertools import cycle -from typing import Iterator, Generator, List, Dict, Union, Tuple, Callable, Any, Optional +from typing import Iterator, Iterable, Generator, List, Dict, Union, Tuple, Callable, Any, Optional import pysam @@ -464,29 +468,87 @@ def __repr__(self) -> str: return 'query_name' -class CellMoleculeGeneQueryNameSortOrder(AlignmentSortOrder): - """Hierarchical alignment record sort order (cell barcode >= molecule barcode >= gene name >= query name).""" +@functools.total_ordering +class TagSortableRecord(object): + """Wrapper for pysam.AlignedSegment that facilitates sorting by tags and query name.""" + def __init__( self, - cell_barcode_tag_key: str = consts.CELL_BARCODE_TAG_KEY, - molecule_barcode_tag_key: str = consts.MOLECULE_BARCODE_TAG_KEY, - gene_name_tag_key: str = consts.GENE_NAME_TAG_KEY) -> None: - assert cell_barcode_tag_key, "Cell barcode tag key can not be None" - assert molecule_barcode_tag_key, "Molecule barcode tag key can not be None" - assert gene_name_tag_key, "Gene name tag key can not be None" - self.cell_barcode_tag_key = cell_barcode_tag_key - self.molecule_barcode_tag_key = molecule_barcode_tag_key - self.gene_name_tag_key = gene_name_tag_key - - def _get_sort_key(self, alignment: pysam.AlignedSegment) -> Tuple[str, str, str, str]: - return (get_tag_or_default(alignment, self.cell_barcode_tag_key, default='N'), - get_tag_or_default(alignment, self.molecule_barcode_tag_key, default='N'), - get_tag_or_default(alignment, self.gene_name_tag_key, default='N'), - alignment.query_name) - - @property - def key_generator(self) -> Callable[[pysam.AlignedSegment], Tuple[str, str, str, str]]: - return self._get_sort_key + tag_keys: Iterable[str], + tag_values: Iterable[str], + query_name: str, + record: pysam.AlignedSegment = None) -> None: + self.tag_keys = tag_keys + self.tag_values = tag_values + self.query_name = query_name + self.record = record + + @classmethod + def from_aligned_segment(cls, record: pysam.AlignedSegment, tag_keys: Iterable[str]) -> 'TagSortableRecord': + """Create a TagSortableRecord from a pysam.AlignedSegment and list of tag keys""" + assert record is not None + tag_values = [get_tag_or_default(record, key, '') for key in tag_keys] + query_name = record.query_name + return cls(tag_keys, tag_values, query_name, record) + + def __lt__(self, other: object) -> bool: + if not isinstance(other, TagSortableRecord): + return NotImplemented + self.__verify_tag_keys_match(other) + for (self_tag_value, other_tag_value) in zip(self.tag_values, other.tag_values): + if self_tag_value < other_tag_value: + return True + elif self_tag_value > other_tag_value: + return False + return self.query_name < other.query_name + + def __eq__(self, other: object) -> bool: + # TODO: Add more error checking + if not isinstance(other, TagSortableRecord): + return NotImplemented + self.__verify_tag_keys_match(other) + for (self_tag_value, other_tag_value) in zip(self.tag_values, other.tag_values): + if self_tag_value != other_tag_value: + return False + return self.query_name == other.query_name + + def __verify_tag_keys_match(self, other) -> None: + if self.tag_keys != other.tag_keys: + format_str = 'Cannot compare records using different tag lists: {0}, {1}' + raise ValueError(format_str.format(self.tag_keys, other.tag_keys)) + + def __str__(self) -> str: + return self.__repr__() def __repr__(self) -> str: - return 'hierarchical__cell_molecule_gene_query_name' + format_str = 'TagSortableRecord(tags: {0}, tag_values: {1}, query_name: {2}' + return format_str.format(self.tag_keys, self.tag_values, self.query_name) + + +def sort_by_tags_and_queryname( + records: Iterable[pysam.AlignedSegment], + tag_keys: Iterable[str]) -> Iterable[pysam.AlignedSegment]: + """Sorts the given bam records by the given tags, followed by query name. + If no tags are given, just sorts by query name. + """ + tag_sortable_records = (TagSortableRecord.from_aligned_segment(r, tag_keys) for r in records) + sorted_records = sorted(tag_sortable_records) + aligned_segments = (r.record for r in sorted_records) + return aligned_segments + + +def verify_sort(records: Iterable[TagSortableRecord], tag_keys: Iterable[str]) -> None: + """Raise AssertionError if the given records are not correctly sorted by the given tags and query name""" + # Setting tag values and query name to empty string ensures first record will never be less than old_record + old_record = TagSortableRecord(tag_keys=tag_keys, tag_values=['' for _ in tag_keys], query_name='', record=None) + i = 0 + for record in records: + i += 1 + if not record >= old_record: + msg = 'Records {0} and {1} are not in correct order:\n{1}:{2} \nis less than \n{0}:{3}' + raise SortError(msg.format(i - 1, i, record, old_record)) + old_record = record + + +class SortError(Exception): + pass diff --git a/src/sctools/platform.py b/src/sctools/platform.py index 6b01e81..b0601f3 100644 --- a/src/sctools/platform.py +++ b/src/sctools/platform.py @@ -18,8 +18,10 @@ """ import argparse -from typing import Iterable, List, Dict +from typing import Iterable, List, Dict, Optional, Sequence +from itertools import chain +import pysam from sctools import fastq, bam, metrics, count, consts, gtf @@ -28,6 +30,10 @@ class GenericPlatform: Platform-Agnostic Methods ------------------------- + tag_sort_bam(): + sort a bam file by zero or more tags and then by queryname + verify_bam_sort(): + verifies whether bam file is correctly sorted by given list of zero or more tags, then queryname split_bam() split a bam file into subfiles of equal size calculate_gene_metrics() @@ -44,6 +50,94 @@ class GenericPlatform: merge multiple csr-format count matrices into a single csr matrix """ + @classmethod + def tag_sort_bam(cls, args: Iterable=None) -> int: + """Command line entrypoint for sorting a bam file by zero or more tags, followed by queryname. + + Parameters + ---------- + args : Iterable[str], optional + arguments list, for testing (see test/test_entrypoints.py for example). The default + value of None, when passed to `parser.parse_args` causes the parser to + read `sys.argv` + + Returns + ------- + return_call : 0 + return call if the program completes successfully + + """ + description = 'Sorts bam by list of zero or more tags, followed by query name' + parser = argparse.ArgumentParser(description=description) + parser.add_argument( + '-i', '--input_bam', required=True, + help='input bamfile') + parser.add_argument( + '-o', '--output_bam', required=True, + help='output bamfile') + parser.add_argument('-t', '--tags', nargs='+', action='append', + help='tag(s) to sort by, separated by space, e.g. -t CB GE UB') + if args is not None: + args = parser.parse_args(args) + else: + args = parser.parse_args() + + tags = cls.get_tags(args.tags) + with pysam.AlignmentFile(args.input_bam, 'rb') as f: + header = f.header + records = f.fetch(until_eof=True) + sorted_records = bam.sort_by_tags_and_queryname(records, tags) + with pysam.AlignmentFile(args.output_bam, 'wb', header=header) as f: + for record in sorted_records: + f.write(record) + + return 0 + + @classmethod + def verify_bam_sort(cls, args: Iterable=None) -> int: + """Command line entrypoint for verifying bam is properly sorted by zero or more tags, followed by queryname. + + Parameters + ---------- + args : Iterable[str], optional + arguments list, for testing (see test/test_entrypoints.py for example). The default + value of None, when passed to `parser.parse_args` causes the parser to + read `sys.argv` + + Returns + ------- + return_call : 0 + return call if the program completes successfully + + """ + description = 'Verifies whether bam is sorted by the list of zero or more tags, followed by query name' + parser = argparse.ArgumentParser(description=description) + parser.add_argument( + '-i', '--input_bam', required=True, + help='input bamfile') + parser.add_argument('-t', '--tags', nargs='+', action='append', + help='tag(s) to use to verify sorting, separated by space, e.g. -t CB GE UB') + if args is not None: + args = parser.parse_args(args) + else: + args = parser.parse_args() + + tags = cls.get_tags(args.tags) + with pysam.AlignmentFile(args.input_bam, 'rb') as f: + aligned_segments = f.fetch(until_eof=True) + sortable_records = (bam.TagSortableRecord.from_aligned_segment(r, tags) for r in aligned_segments) + bam.verify_sort(sortable_records, tags) + + print('{0} is correctly sorted by {1} and query name'.format(args.input_bam, tags)) + return 0 + + @classmethod + def get_tags(cls, raw_tags: Optional[Sequence[str]]) -> Iterable[str]: + if raw_tags is None: + raw_tags = [] + # Flattens into single list when tags specified like -t A -t B -t C + return [t for t in chain.from_iterable(raw_tags)] + @classmethod def split_bam(cls, args: Iterable=None) -> int: """Command line entrypoint for splitting a bamfile into subfiles of equal size. diff --git a/src/sctools/test/data/cell-gene-umi-queryname-sorted.bam b/src/sctools/test/data/cell-gene-umi-queryname-sorted.bam new file mode 100644 index 0000000000000000000000000000000000000000..f14155a776980e919c3cedf2e3788b69c0978c8a GIT binary patch literal 31074 zcmV)-K!?8{iwFb&00000{{{d;LjnLX4Bc3JjGaXlpY3+r?p7W}ASIw&|1bpO?tU}h zJXaI?0BW(jg)Wv7G3(vE-?mroW4(Kqwi-|)!H9`OsrbNv*a}8M3?ynm6v7{9D#qXg zf*6I0B*tLW1jJy_Idi}HI_H$)|8BC`o!@WfeCIXi%*>X}(?hp!WNg#>E*{-IGsP#= z=(e3xoo=TXo!vRL(z&+N-QOA8G<(J9?99}{z8)XlF*7BkbxKPi5<(0J!CkI|ONbl@ z?yPW1sg#fdA`@JCCLvTHELXV|XoU`h$!w;aNeL4O<#OqXJRxi#q~f{NxlRZdh}`B* z8g5bpvjLMEXRw8w4;YnMtuZ6YHNvPoOBLe->$%Le%!N&vf%%!qWNtWO2KHy%InFJT z798NR-0LkJh>Ff+=JO>Xt_}-~5JvJ8pAHMO(v{KMc4Y8lCkE-2Ig?S&MUGRB zJ2^3knd4GZPL-S$)@j6vObujm>ob-U&JNZ&p)D3EcwUP1TO+4jR3$e_v`bExpEFu# zl(SK$Pa;D&S8~Sr{e?M^RqJdf2**oKsoW8nLRK=h%L!MyOmNN(<#d#p9XH|q+l9&U zOq9!fxCs*-=A{Ed-(O@4IN2pRc=ZMx^Ra$GpL{t41A&eO? zxzR4eRE!7gD#Aqfc#9|%1MGJYX7Vb$ig+O}eJV6HHLp@o6I*jteTOp9)=UuH^3qn$ z1@U$sQ&tKhO=c%!bm@fL*H6O49PverFfm80A~~WiP4pf^X zmqc@nrAuOhsN4`H8Y<{su+dOKx1-CVLoJc!q6?%Y+tJ0nZAjc)^gvE)RqPnP&k!PR z5#P!Q6Fp!pUL#}K5n5`(#AM^f5hkXTw3I3RDKesgj6DJlQiQ1zM-XYE@+aB8GS(ws z-iS0&y3JX-q^j1WPDxR9SrQys)wWy{X|WnPk(@J1DrOEBifmvhV+nC0qatOw&Nzq5 zrPNVDLzdj#d+nv>?BV%7HNAsV4rw*Af^i5knc zcLy;jR%d?*#Mx-Fs@W*I!Aw@=qdzwhs_HjJTgsKXYi%`ll8Kv4`kRF&6pzhSaJnj% zPors}Ht~WWGec^tiFMzXd2R?5OO*B}MBZ7lT)p3jF{um?)Rb48y zrAo%8;kz+8I+*HCk--6(aM2}Bp13Xs+j367Cgax>4R+O5jIcBS@h0N0)j9Eqh#rwv zkTaHvwS-$gP0@w165GsDq192`_asumOrA(Hn)ifeSfMZCjx)tk1jEiZq(Sq zed9PGb~4^;ct%gBrmDg|WvRbp^oef;#3xe6y)Q|pGBUcvS?^WCMVI9C>||uzJid<7 z98%Tg$66XVmE;L&e@4Z=ZB8f}s`QB#j$B23wpGQ8nypHY!)wRv=Qda7OYd(DOAxkkNTpX(RP{jp7#2X^+#&AUb~pBen!J@`8(Ce%b; zpKEs;&BbQX7~Q%f;I&$9e7v{PneH}3HW4y84D8}2(PJg^nTsIzPj28!rg_VA@+nGEOTAtjw)Ni&| zmM8rNx5Yx9kIygY#<(<%#c|un7RPhFpwxn?%WPgvF4XVn0mU7bkaD`ra{paL|_EpbvDGyS=`5GMZ5D`*_c<_(<0_w_2E6 z>9p%hOMb;w&Hlb-2a`*^X1m@ycxkJ?7glcZsWQ9J?G2cgKNZzdGfY zW92}zXD^$G1H)AN6Qkw8vL_DTy{;TM_MwL_TDPYBXvQXP``6ky;_U5PPN_J-?z`jf z6_>MbpZ4UMIFjt957w)ZVrOkVvl?l3VSj!kjtr~sdu4STS@zDwqt(c<)4Mkemm`s} zwO2h|jX3+m!#`aWM}i$`uN^K;*ux%qeMdCzYW8K@jLbB9a1&$e{owfPIvdZ2Q^I+8+g$6)ubMi8jl%@71%d~p~-AnOGr9_|9{2o~>~Hk-ON zg@SW7Z9ou+e_<9u;MLQSUU1H3{WS;zvp+Y6ATWCmQUXq;Y$KKhs;g;@`o4G+TXUZe!P1hWt0 zg<$q&qX>d8E+HG?(U&*}-UrxO2!h)`yBQe-*M7GDbOeE3Pa##|qR&2tA<*mhD4F09 zm;D>hf^ZCB2vmI(>jkQwfjj{>em0BDhO0h%0?z{Ne}=&jRg_;#r{mexw9E zz_LB5bHd9j+dYOL_~PxW5d`8N97PbA{Th}9X4lst2+STqvck67fMkVj^{o*EfsMZz zMi3}*4H6&9&^1Vjuph5Lvck49*AAr+sCw`RDFjz4cKvH3DFn0Z`Kc^;1!aF-gCOwg zlgLJR_hgr!h9HpjbP9!+PWItPt9FosA%nb0AVXLHo`@5aj9s zw7V{7nQwVm@*m zS_s?31Z49}(>J_vQwErU#AIi|D zzlI>te(W$>2--jU9W)EHf7N%L1On~%-+EpOf%a>+uTG$_74AukRd_CA$9AM&gjX`Q z4S6-3I^O$ZGz+wU=tcyA_V2j|EevP5$I&d%{*EspB|!Vh8xaKB54|IW!a?M~7gF)? zvc+CqOvQ)lyZj7-K>K5QIE6s_jbB8wK>Ndok$BMl`@={H(0=DSBn!0vYzxf-?JsDh z;^8fd^)5jWX#W)QJ-j@zdMi}|UYyu-FCz%F{|`!Ecwu78w;%|#e^IA;!HW``J$s1# z7uQ(}7=$JO03VA81ONa4009360763o0OK>=eM^ih*?AuCof(NX${2NP?1cNn%Fn-426s)P9#MLf_RY!kX2T0iZHmyLII?`Aq3uzMxH40D(dBly z80}msr6Vlq!iR&O=qo90$9C0+R2PqZII{S;@NuT<1&p-+H-TF6X^Z z!L@ZhU&5vNviI|gx3aXi{OfP!KNtH+Fv)&4>AiXV*8FcSKKX3fdwuy<-g~zC{FBuu zKmX*HJ~{hGPt+rN_Vb_b{nF3RK2P8OWB6wN`u?5Xc>S&Qe7T;(&3d1Dy?WRDjn98} z`NmtX_U504AG^N3ey#lIy-{B4ed^P1RbF{#AGBV3(0p$aO}#KmCh$_6Wrm>P^A#m# zjpttjO!8ZV+qmf9&pM)bzf>>4zXmBjte}0(0JL=Kd44eUvM}(#(t^Yz(6+P7v%zeT z&vx5puK_@_d)xr0_~ds+F$Don(uk>)DOZSn?_z55N7Sc8m(r7*?2>+ zNULUO6;$}iA8?7K1xkD0H}cQnkA9F1?I32_(Egv1q4_>~$83_s)6nxG_Kt^e$@9x0 zqPQ4shiCB8crmyb4tC?&Y@82ldf*PjCL0f*;Bur}b$NKjvH)numdaN;^+T^Tx5oh* z0tOfITS3daddF`X9ciBDC2%oKgQ=gzerRTC4xj~Rk`#hmnXsh=kXR1EArTUI@U9RL zt2j2=6ucQga%=}2skANIhW{)9p8+3R;_(P?mJF;YK*I+Ik1?OmbNE+Jpw00C<0HOb z&>DAYT3UJ6eM`G;fc7eamLKts|H;5RV(>U2xF*>&OHwwp*5SwH;2aLV#mMeUG#Lj^ zGYbO;o6jM-g1<|SJ{+OI0=NJ@zzL#D2t6dCXVY=!09s^S2(AEd;d^luGk9=d2=HtD zyfK4s^7R_c`}OKwk9x=JTMg}-3bTKw4DD&hfZvAzG7Exe8ia8&=?&oRA0yxv5J1i- zfLz^<;Oue>PP5;~Lk$hb5`8c-upclWG#*R^95f_^4t=ErBNSv?7GPzA^}z>Rfgb=U zrH$`~m%(u164)O8CM>H2ZVeA}j-SZb;Az4`gjW!_AWw47pU0DVoQFmjzEs8Vr4qyP zrVR`$e*wVoQw_VjuG?Ksw#XRjmkMY3J7&(3F^ofq|D!bEP;%Hc0RYnhT;LbDH=FYS z*Xc_~#VC@|gKIX$CrKLX zTj%)eji{t+Ma%Ob&yF$=i(45x$3YIKRu|;J5klZ$_sOl*4XkNyHl_+3os`E`@s;DK zvKZJZdsp~ZAgq>#maMah5o}MYus!K)YJblFwgaZtA(#Z%qG{$s(iS9MoS2y!xI#20 zrw`8|+HHVn&HEKc;i4l2w!!)t5EjDD1-6G>mXfJ|g#+pE9+=LtbAWbdU6E z@9jz*(lx-$52tC8MbkJ{CtBw!Bt#vS)b^tahgNYS$+EIV>6pO74!x)LnN6H9)TBSq zp@K5tZ9NvqPMm1R%q(v}rbEYX7?zGF{?tn+DR{@s57^GSAqMtVoJA5t5x0uP!DH&+ z$Km%#QI?J$60~QP(6Iw(uhf#Xy?EcW-ao;5f8s}&uSJIS{z~D#IHepRRX;}V!KW`c z>b>Ylo$943ijM>|GfO*B?|n%kar%hMMy zlo3T-B1_aP`{9qYt^xmnGUWAZy1rF__D6%MR-yKh^4 z%d%`~+g$K)u;@=Q1r_xncN5am=u+wp`J2sQM~S1vjx`o1%L|Vo4yiV2I>jLzOExU+ zHyd%t)2^U}09ulU(=ZH@q^KRAoeik`xNU26SitWCIDN~7hYTNpWWO&Y*N7>rrbZ6e zAF)|uK}ODDi^Is;R^c%n>Zhra3gzSla7#Q};6M0hPu;aSTv)I18B*3Pf4GeYla z6}_hlz1r;>=v8hF==~i70)jr<{?Mo{;`R@w2&E!VJ9Nqsl$01?36n`9x@Zi~a1O32&+(DX*VL_9(*-b%d$qXcof?sU z<*uD0e}3iNM*Xpi#i^n(B-xN=2NS<+)-2j1c6ko(-1MuD-jVu{E6@=|NK)xiQ~B!% zdXC}(Nym1D1ICHCEvcA`y)XEO&h2b`fM5=21BhF4Pq;Pj)TkR*?h3el#-N@3dkVL& z|AQxxWIsE9v13d!iKhTt2EnB7^)BJf|MZ1k??3;QUhj)f0lwXEv&he|VRJDYUTj8# z?KmF{;6yu{ZFY;@YPH#IFK$*dI0G+-v)vFoW&OVFTLS%~g{?9N3%}e!gDZmh3KIP! zT*7o5n`$cR#QF&~_{8ubH-K7~5GO3jU!pa@C3vg=Ul`N9QaGn=<;gejHV1CGB~!4VHQ9i`;j{_uZ-hGS(`O#wUMZ$|gysNfaj2AP6F6B6n+!USd<5c|L`8G;OOOs~UV1x%Y#bS|dKW zqp4~C_*Xw>%tZFREK8Az5OO*B1hWqOH1%BIU~B z8~~f$s=DQH8qTq$Q?Mc{4IEoIB@;n#HCa47Xw^fl6KDquCht}ZtvlzRnr+}C({vh6 z(j+-6-W+n!56HoG=k^UG17Qk5qz6$+X4bw@VIV$WOK`lQ zjiqWhhkmPmDXXDn&qkIPjq?$(Zu%MZpYzP`8Xid=Gu7K9Hx1BOEGK-TGK82YX zWc-{NEmctzf?B5P2y+TEk~IBv60{0@kHgN4gAu%@GVJ|B?R!QQE=tzYAj|-~B#V2W zeX&Qimli&8xV;&VXRGbS`EWej*%aFoDmFE004SnJs!}pu!8(bf>Q~62!`h>x<~UqT z@A%x~*(eS+k&}QWPE}IIAqp(R=UC#+6>eC$z_|({|CmvK4&c_f(_S!A+!b(JvmOOX22y1<$oN&R^&{9hG1g0Qvtxl7zkg z^-{0*g;}pR>0w#&1yW_49ANL@T!Z^z=qI~ zYy_teb>*1gL~gl7VLaJj#cQ&$*}4X%8>-cCL|L6Jxk$3cNeXaYONb4Wv&4EK`m>zC z3+8^jjO(@K=u}=)Dxhu;8EW1N!Pqnlf_lu)sW-*HHF!KF(XT2|1FLKz-*$qC0(jpTTJ03PRmF< zsIVYuaaozZv@QIOK9fa3qK5BLha#CJjZnEX{f320n~WA~lQat_%KqC)%+hq8-Ba?` z3AIY{c7^dHL?c(PftSoHsFA9igit zqD;!9Ec>XgZf_HV$U-86Dm3eN)2S#wIH+~*b$`dO*Nr6CEQFL`633aDe{6U?*&NzX z)Ysbg#d}@NyrVvwwCX9jl(Mz0Jx~8wqXc#Fj|VY3U}@dTBoVmAFrJ3K7kg%w<`m-? zZkc*(P!5w3s?Ze#2D$gDtYmKwY-!KTEbSF^ko&0Vl>3Ke%H63>G6|-k7p2oUi6gH6 zcb$;xTf*tveW_^v4$H|DZcsWEh?Fy|)|yLqDXM_$481jI&(qM3LuNv9)&0?LDag9f zxB#+j>Ss7TItfY=*P<)Y)xvojNSgbYXIYp~QSz{KPs%->WFKr9 z+U*hv8nnMOCc8OLHuYH#E&wF^;@G1811AWc(6$8!kzk^h%hhBgn3rN|)v!}rW+eq6l0K956}KE8EY^4L+wkgl zYOw7pcQs7g-cy+U%WwY_c*~#Yy#QZ#Nwpy&@luik^0F-GeG&e+;9tMe`{Ike-WGl1 zV6|8+E{3D+U>WdU{+Vvyty zfI!n77;p-l^kB9RM4eb^M&MC|##S{j+hMk>Bw~X_3Sth2q`-;hO;A$5o9> zL71~~=^9tPZHnJ%v-^Q5^9~Y%D2y^9bSYj*$0)_o>LL^~+8ki)&m2~Bvs+uIQs?g+ znQ}IOj}QOh=X26AqgKG29_0F-`=)z$Y6K~jyBb!LHD*8aF*D4(1dh5anFfBSMuKl| zY@Y#uv*Fobb9Odc0C0=51z4BpLrBrLH9xtcC0U^yOP)LwoX!=%36f2dLK80WE2kJY z)oL>M21SxS8 zRZ{QU#)4OeWI+3j7v*3iETbB4P6>rh5U{pPfu~?~Ey@9qTnWDk>=F=Vne@19q)R%O z)LX59i_&>uhTUGB>kxL$*|!qALmjs|cKg&4T{L;ubb2bDMw8Hkvrg+IcSUW~uT_w> ztx5i(9Mdk8nW!r#$`VScf>D-{i%y-5r)4|Xl6iPf3MJiOy9;OXeN)vf_!Rn)hqD0Q zbdv6;cFt zt6d(@N3sYvu$iJj6)=iANd*-D0GrZ497lye3`;Vi?HPAZRmyG#cWagM2ypt0d(Tio z0r@3ff+}~DMC-|4t`_GvU~!l2O3Qc(+5NtS8lIF$FW9#`ZRRM}{d` zbDX#c4$sC}h%SNY3ahNGI6}6?ELu<2@iLhwW<=!|tJuAG1iPmCMTy;;ts4WS4xWxF z|Hu$Kh)<$4@WJL{CBoHtCf`wBudGbETxArBz@=blJ%dy2=-?^$@G)gmNa0jWd0Vma zx;Wr|=3{0RW|Womyx5uC5&=4fD7kid`i2?<#HgcZE*L(pU)sn+P;#JD zhZdSgBwe8b5<$5(cgShTwA}$WSw?1lvef~1JQpwtOztKR|M)JuT{6fWu>I$6ffx?Z(GY*>vTujD7rb9`A zpmAcXjiJHOJP6Kuqw=ytry~4ZE8K#)mwVnim@G5XL~c(!xf_Oc@nr54&0KIhW1#>?%&F5aIqW41Q@bWE#S*myM z0lg>)ylI#w6T>{*F3VgTInwzj%RM$)`l;xMGniIMF4A>LRL8EyCsbzl|C>dI=Vl-nd3ml;>sLX-``P@0Rxfe&h zyn2i&D|x|kiHg+t`&y}(tg&!^&+0RKciN9p^j-I*LRpCGx&W=FKJ#iLw)|tm8Mp}! zoO%f+1VNUnr5u>ME!ubGNZ$#ZV4swns1VsF_Nzt3Mu!@4s$xQ%6ZI3eTt#y`o-LOn^&Otkk6C4~DTle-Fu70K_FhfV0Apml=9Us)aa;=OY9v zwzV8ILXLu{>n*g*Tb0_~^GEY0`e zQ;?=HWZz-RYl}vgi<%=`ok9*rGWx1dBlvt`W3eTT!y$j%VFf(F8w5B#XOp6+($c4x z5GpU%F^VJF9&6#+Pv^0jn8KbK-XY`~X+Qyv$7{ z={B|Eq@wdVl822^w6fwM@d*;OiYNlKx`@aG%bcqwYcEd@S&c9Y*T6nF9iErJncAtTpyi^^r693MJ2hgKDlgcELgAl?F>UbHm~0k&91nQ2@_ z)kIAisI*Px1(Nq+YKx=cLHtb1pm=T(|IMQv+^?H#+-YxIvACxW4_gu#`ydkIHb$ysH z(}rfM5`+axh#ozn2Hm0(bQ_9gjZrO-a@4U@cETlZ#G&hzMYOhBWkggjg&6J%BSJb^ z8L~1TfuYs*72O#$eX9Myu-)y47>`6WNfo43D%Ya5#o07{qYE6Cq`(mUpe@x*zR)#f zHM`AcTth8XI##uMLN!X^mR03bpm|}Q2KSY*b;0e{jP2_7jLpQ@X!Z?11re&xCrN5% zY&LJQRJ|B>gw;YnbAgUb!Ra?*y;AiFM3LuFgrhU!zS#0eONl zlfIYmwt$XZx6z96k$R$=?ZvLl+Gws;0aEpI9XO^cTQr%^tBr~#_k)5~TQqxD(7sx+ zwC1z8+&L|hbqqCr!K zV@@G?GS0zG!o>B$H5~BAP?tc6auVPaW`Va1qor@wU^`G3>Ea@9t&6;JSg%mK$nTnN zE}BfEB*mqdky20Yx`}U5^NMa&N4cGMERVa!PcrYaEpi_mYtx+aKQ$aQMfK#=OS5SZ zdWvv2Mp@)8lQ-0V=?mV)LAKS`%4sDmj^c@{8WyJ2s@SEFW&KyML9AdZ;X9i7a9U<- z(K_=VU9xtfn)O|?L3|(Aad^1(GRmUR(A4OeTbeY7R1n1&Q5sx=QKRX5hE6xtxa?Rt z=(vm}xnM~xyP|gdkfye76G7%h8~Ye$xai^fx@Kx7nEBx}4p8#~J}Jx)ri%`iMl2tk zaKceSYQP9AY2XSDs!U<>j2te3X}HZwC4zY)#VRpcOy$}*Z)&6ZGPQC@iBf9X$jQ3M zWSQO@OR`H5@-DOy!`*JkNcP}F!=`B9$7Z|Tc=V!(4Y>qqHBYn}Be|KN6-7+Tt!?k% z{@TL_k@s7K)S!LgpT7*Cb!)AuDQ-WQ_yI&EsoL#!c5!}ocD5T}n7G)kX6K{vYOxs& z#^)=*a-1)=gF(KT?QY=Z;o>sKPPXkRqJYf=A4e8uguOKk5nWE3?HXcI^MQ1o-z~){ z_0Jcv2>w!AX68%QPr!{d%Y0r3jr(f;9YCw!X-}!7xa&S>rHfQtpS&4AR%U<^G^L-= zMNtwXe?N_)a0(gr#0)b@d-hN-7aMCtC*~8(adhj9#hnVEN>qxPRXvi_y&Or+@Q25& zE>v<0W~U78x*CX>Vdndktr|uqB+5maPD)SGY6Y6_nGUoa2p_xPb|QSdx=m8@M`q%!9>pL*oKC}7k;PHSx|6K| zWRBfOfcium>62z1QTDkhZMimyHPNym96Y3ajgkT+^3YNgtau^q7|%a$f?AZa;b2@l z*rNIU<=6YD)fSZ=LrujZ3|%2nUhS)0*YoS2Hj=iW{2DAS2&PdO#)?ySa&beQQ-^7F zUs`Iyx=%0^p#}}^mUPt$-?-q>RfCaAOv1lE4i8ydlz75JURN#cK>4v;l=$E7(hJRy zoD~6bY-%0%Bwn&uyd~8xU(AUuiYjbo6icaiJW3tZW}1T9xJ`R;FdK(-nySjIFUiLp z5n!E%2kqHWL%T~`?avIW{irCJMYPmeF?_me6AtMN!nT;8Ma?uMG0 ztXY+2#u>IK8#u)(XUBdFn3;tlM=;aE2!C!mH|j%jjT$l)jkEz+*m(zyG#*t$W_MrO zg*cW?GlVulN2DyJ3`%?i#rfeGQJZ&J(q0^Pl)_;tNehL}?UvrrUng_JQr3Y6O_x*6 zgmA1qEO#?YZ7b(8lAHJ7v zjRg7O3H;xwgd#@chThaqVwC;yG9|j(ITq>LLaVLf*qIbmA}Lh;CCx=K2#%v`A50VH zd<=)Ni&B=XEna*ivB>>`rYC%V%V4EcQ8e~rZyIM{Xg&0h@DmF-ypC&JSKQ}_-RIWA zSfAMVq)6EBlL`sXV&kZWByTBMpoA9!;~QNqJIHfy^fJB`e!z0|yHN$mYKL1qnfu<- zGowPb@{`9S$^8l6y-4zCuKktI8C8!_Nzw}@39Xb>)#IBXfHWALUtW$bXPfbAcD5Ve z4D-=&wpol1+H%Ha|H^sX|^d#QnF1Zh6acTI;uVGP|3@t z@kCAR;>9?Q;3ZR1t$@QE&r&6SdXyy0$^ZCgZ!K{%5I#;CsfmJh7#a;`9Lw6un>9rw z%NrlfwT}pqv_lcYHf>b6U+svSNNN8N=YSmX#U1&hg$~r_WWRb))uy$R4iV-zRnM5HQ`dC9A43f z+!xIBu5mUyjDrYc%2=`M5SYbmGdde=$BTS>KA5fY*>E(QjTYn4c00V-)Cb(a=Q%h~ z-=`sZJeIIZ=mP3~sp5#`~yj(b7xC^Jwm+ z>&Rahhg>)G4k&F^?m1WJRk!ff(0j95;=0O}*}a-nx2Et`sa&_>3ZwgdBlKw4H;R13 zF;gmOoqZ)25EZ!!3&yk`iirqnGn8AD`Xy`q>qJ4-QPjdPOBRwB$}Fj%Y&`~F+2V|Q zHNUz-+>&b&h`By+#$6@w>UEh8{1YQz39v;|KTL4wH}iEtd9xYkqs=Pcfc=d&i;=d! z-4Xk9Xuv|!%r9CGTACyD>ODsf|tjzCqe43)OR1XXvH^hiD$4L6&c@fr3@^A`P5t*c9QDi`%W zsTk?-!aX$%Wzo`7S#-FSvTMB;S7Rz@GINsDOY89Ew)cLKS!i29BJvIW%#P zBX78B17lF)AXwPk!7RcdEd`;4noi*s=h<_FZEglzo|98NHlS^yxSR5E8@v{&bzPS zGwxR4O&m++Ro!pNqXO~Btu^knU*OVrJw!#;l@iUg zyTa|V3|a!$A;v5^I6q;a8AFofA$mz%DVyqe#|hCptK{ zPPCZ}*VAwmdm?I^Vzws*VvtIyghlDJ7!$8ep*0LHGq!ZWa*vE~n8Q_`Y!7K| z&)Qg9=NXBKuH{X?pjKv0tTN@sV`b=RfSw_FtU#6!z=dCQ5eX0ZUIBt;4$9tbs( zh;fLc=#uCNDKD_Nkb)#HNAG#i)U7twc3%>92ViS-ODnLoIoXb_m3_NhE8D5@hE>Vb z52snnb~WDk|fP6!5bt zQYeBTbI2(k-qMb{$I}~$_91(BKN#uOc7y0K^N@|~dBLGz!%8 zZu~4xw4l}*^AJ`o1+^?>jax{Rw4J_Euh#nzsWpl1PmU2bQ?1*=#0EgI9t5(#&%9bRc^bNz%wthu2^B@zaY ziXoQG_|jA&Hn^!gW)sd-Ug78;7fxbbV#(Jk*adh4nLH|b;8ldEG|hCoSF&^Pdl6I!M3Ya zmkQ}>Qrof$i9yLzzbt%)%oK}!C+UOQnwUrQWeQhFl3wUuj#?gyk8ZwHRpm8KQ=snj z{hoG>J2gP7-1YQk=YwAJw{$B&RVQEQsN_%|382yhCg+) zIGf#!^3iy<8J&-}=bORJ;(Ui2R4z7KIP(VJS8;l3w%Kfk=eupbx8<@Cnvlgc)?F6! zqp}gjXx>8wK9$m~aZM=BeWGoU4k(>IOZsl#^rBNMOT8~Y?O9)W+N+wMU7BD0Zq0e_ z^3H(Y=k5o;ZsUp(hAYUX3^V_3_~B5FYtfn;_}sT>369M>IA}#MPjN#~;;MzgYG)WV zxk+c9_k5$4GUgMEo>qra>)jq4zuI7TC;YCiulDf!Rt3M0cf^lYO@ZZwo}Zy!iI=nl zaBu*?EdVHJcnT45P+Xo?YEVFlAHphh2pt|~Ia)D$>{?zov@x ztHecl6)V#JR>Ale283CHF`EK@+0;w?ln=u{g+K9F+BYoHzR_@dGak=Y+rjy8Jll~x z2L4jJpOL1AGV!Av7NLuC7u`5}J@_ZJ~%4(WJVlsaJ9lJa&oMO**I={%a4 zOYL;pwQkV_s}M|j=8Vzex(hx<+*Xe(>*uUNPeEe2r6b63>!O#fv zS*H*Qe!?cPC)^fBZCE@c?J|Dl!ZvLua%8D}`bSEso%G?1jnQEKy|wb&D$LvBl-1A{xaTnainW z^E5U7u~=eFGpdUv9wZ=X0KabrAd&qb1jVcn^7>xyB?i#Je6TBmRs|bp?~CG_dZDF) zXi=+i4nX{P0IU|&?(Rz8QvhCm*`%Pxiy7lAj!?%w?R~Sy_q6g3sw8o6kV1ivaDA4MU!Qi7(@7^3gMH^ z5Ps$`F==0#=ztCWz`zEH-{OG{`V-oW8pXX|RZjGm;g9X;LWk#46|$?SRZUA5J8Q2N zZ#v{!>%3{_aoFIqHsQ=ZfCsoP*LRFZiF%7Iq%26LUK(b>-A+Mi^UiP5^|H7crrJlu zH?wUuo9s%pfm8DVd49nzD%`rKFEz6EV*>DK7;~!=JwjIc?aC?GwW>|x4wwmsVKhzB zY;w0lP)$9$4KM1dghHDi@6&%h8}W3gYorIeL=xH19Qu zZ{DeiS$S8(B>wFRX5E=)9**LrAsl*1syLmUA#)rS@MUcVE|*E0sqENO=j{UyF> zP)B8Td7;al7u!Z5n-H;y@t{OVX;^`HEe=n^BI!kPu~SJoSlmmHQYq#rHWLpLspbQf ztu4i$J_%VEEQ6(&EGK#38G*My`g#DoTB=?H?@a@H(`*+LLN2KT$2db7Wz!(Zd{zo; z&2rI|P&PE)?oFE+E0g99Ykz5HSXxr!Y&$-crx!TtM~Rv8oTWMF+`)F-p~LO&7N_sD z-{{kKZRI32)$Do?wC@;FlqQgvXA{qx27#}{9lFYJ&#>-Ha#4Nk>PVuI5m{PwWF1y& zpybbWw0J-zy{ph+>S$sFYBx@{2LY{ivfTpO-r#GUSH@H31@--H(`W=fWnq@m0OcuK zS$+^dRdmL1uaE5?<@9Jii{|oZJv2$HNUb7(tWIjIqZ36=An`sO%h3zOrQ$p79c@Q3 zi^6qyUqf97U~5oYDzKUGalfHp`>)E@y0-Li*%iDNOufkS64EZ}(bPg+7wqy>lw$WS zT25doIU?LNz&$&fW}wpwEzcIb7@T(DYO7;YL&+k{04+y}2dJ(C&Kk9Kl7Y82+$!A) zyW@V+gp#Yd)bATAD*FC3OHm0Yn%nr93oZ=rlDXroXE5|UqXUY*Z`s(Su_en?Q>z# z96HNE2#ZIKTn<%*)e<=Dn@4jWB1JY?qgJ;s8P)_(S>1KH_C-f0KQx{3 z0)z7c?86qC>RspZF6dFadf5lh_;8=D5NU94`6#Mf!nGT3u|gB6kncfDjVzh?f_18G$8jMXMRbyQOrLtc8dNv8y%CD_^P(yR@q(9L5Ofu9Uh1 zv3}H<>l*0lx>cThAku^L@SxQxVjVzhj+fL7&19z|EiZtJDRx#pPc2XFD$IS3L?L8{ zSt>7+Bpik(8%>GuM2SlzF^WALSp<+o^>axUn8YBk2u7!&sv;w&Q_XihenW!~1|`9| z)&WDiVuk#|(7tE5$9<9@p8+m=k*~#qJ1v-lcUvrdLIDX0>FJg2SYQ@v6_0eap{c;p zt*Y-8ts6y3v6YShvQoq&@QyVKUsvzA$6@%c0`0%6dAx3k8_B7Au^)#JNyd{&Zwqh# z)t6A&ZFIKXY#@rvSKFJ>cD%dUjfSi3=3lnWgb%dPxxBlEqW_OOj5ZA0nYlBI42y2nshb5LLhtuMYX+OD1EkaYPOt zC(m*vYcI}|JYGlZI1iVQ5WpAf$Y?6>C#z`wq(ZZH`$j5P`KuhF6!Cw`1g{`ZxOz2j zpsB~1-!o#?dPpfsCs8yFlQ7Ln72%Gh%1f83AUH{ffxm8oN>Q>+sBu6-9g?+jXmmjo zv1rvAU3RgqY7t01(Y&T$Z>dHcxetfRYdD_Daiv)@aKH*1ciO8j7k4!*pAn(f9w>wgo^$E03VA81ONa4009360763o0O2*= zeM^rd*>zrabA}8O0dT663KVb`zLm@XGDs8^?}$v-!jYNEEVjuu#Yx)afTTv|$^vY7 z;hh&&rxQm00D9q_U>IInq)19bk+cZQf((HH>Npk0mStiWj|55ja71@Y`=_4!YR^ZDmKH_K;X5C*d#NYiWka979Yq41uMX^}p zr^RBk8GLsBRz{wC_+%Ppkw5)-Hu&7?t;OG2e&YFh@R{{n z#o+nPm!I5x;>%Bd`4jVh@I>9J=b!rW;Fmu&|1v%PNAO_rnZql+3D3df0v-=O`x*19 z#hYLLyX!aKx*9A#|5j0~R;y>#$KcKCxxr^Y|5o|JYVcJB-p{`AbnsCC@5$heL3{83 zyKoj}VVKO~I7|of{e!`Hcexzx$HVja5b#;fC-XBv@O)Pk+q0YPZamtZ0ci8fv*Gsa zY`$EM5GO>ywFwsfaPdj6eCQEiUiHxU>>hof{t)=kn86J!;D&9A^&0LMU2nF8MB|eDFTfS*eXe8MHi5P&@E>?V zpV`pxA3xmSMNvE~XpJj149&ReA!tVot!Nn9v%#+^(B1{mypDzjpoM-mn}&WcOM@V0 zJ8R<qtvMQL`t2Pn(V;RQOdXTy1GR- z{YeE`7du1WnB-t=Q5FXT+1+NCOcCJ zC!D1Bs*?`v?Uja%eB<4&vCDA=wlxdUN>4A53h(4M~06SV0pnEKOn7KYO}B+yQK!QJHZN9SNqQ;T?}))+tj1znJ@J0ujIO$9+#U zhz1?Y#xPJgd~ji|27w7#l<*h{2qr16bag-?)mi58fd^06=uU-lm}IGN!r9`{F?-k` zLk?kM z)Gr<%v-{cEDa^jw%Qr&cn9ag8&;3~je^7q;i4|Sdx&AfP@Y~!K4jNI(W_&kl=xpc|zLwGX5`Ya45d1hR-ju`-S5QJRAY4+^u1ZIEH^E?aV*)&fB2pr?w=fJU5R%X~7 zyE1faDzS!B2oMkinIs2EafpuS`LknK79DFb86{kDSb>phaR3}K7mNh%w|k)~{IRGcke+?1)>m{YYIBRJXC zkhCNLhLt7)+m?fa=gkK7` z%^GgjX1&G=E)|O5?!g^gE8P4=xDKX^AS%Le4bHv~gs+^YLN-# za2mHitZ?g=YWuS=&cb{)or1r48*uxN@F#XeZN9ggF(C&h$3ALLONu?w=Zdmk>DIpy z4^San;^x11xIH&cx;k#lJ4z(SBFZE7)+wU&qgfVB@wiLV@E-QYhzKk}!r|qV2)Ul_ zoGDO*&eo0tdH6JUY98O-o`aiQQ~q&_z147A-GbXU4ctD~(NAK85`?oLi*h*X;P0B> z$5DbabI>?!3iKtqOFZ8+ggP6r%Wx6{QFN3J1u$?hvghPporjMp_@r9rE8{3Tr``X( zN2vsRY}gk~A`Dd0>6AmumTj6cApmQzxVt8}u_h_GhKQb0J4w(`mrZ(n=5c9+a&M}u z5Iz)}z(I)kafyV=&k(2x#td#1-M9q=1cVu;#u)p7^Lp5HpnjzWX2w+y(ryNAbq8j; zj{L{H;z|f8r?VhOR|%m>$Q_#&<^Fns%TnhNfu#%fTp)3|Eh@Z&4{_2vs@e*ldx*RC zkdFvjEg854w04H3bGRS(Yi45(7Y0F=&VnR~*wPTOW&4uk_#7jC&jUL_zerFB>|9GB zRY(D;n}hipWo=9(f#^`;)94ain@S}3Js5I2^e)#{HD#=pv;r;PYk?HuYzLpTOUGjPJ)LK-X}D zN*)XcaR_w%SQ^ObvyI7G=6f;pv^m5_%qVi9qY=C&*2S>)(hj)vF@U?nL*8m_U+r1t z3Xm2V-{*aNJR)DKgnlV5XCC& zzZIcFF01W%Rl!3O>7-7@YN3Z)dtf+EXcf@rI@Z6OLwL&fF{6x5jUFnMzwa&xk>I4W)lPP@sWa!PKNnrrMW zdXw8a#}3j>9=vz?+VeKi2=l}1_C^gGqW)C2c&wh)jjol=k{}IdQ5vUl&w3U>M(s~H zkAx)^YZXD{(xayQBL_hnpg!VE4~hwd4s56@iTnklRr)FhBNflzz)S}tKj{^WKB3)!*&=Z<#P!oogkZHmD`V^{U8-mrxv`mf)AGn!o%>uliPq z^{p&SA}p9F3Rh>-bI9~Cn4rwdARaTh+)>iT?GnKj**6^`K$(HNaCZmx6a+*R;|}~tKg0mYW{G|q{Rnz zmyIhO%paClJ?bu3xA$yT^~QFeo{f(+$rOoPQ5f(5Vyoupg%PD_gAraW8*9uG)`IN^*|*mQXWi6~f?1GE z%fZpE#-?KqyBs6bg=CjzeO0;GJ>m*aHv)8qURJjr)7VafcIY5icg7}1N|N1JTOlAf zjnY{V$2xiIT9UM=?yu^R#PjFI4c>0DEa@4kHfD3S*N8Uj&XK#KQgVm<)d7vD$*X5q zG;X6xlGU{DHC897)^(pcGF3{Z;31O~2iuZ3?FX9V)LDi3t}YcH+KuV_P7{7`3{fz= zY-umSvG%+TXxB7!XjI1R+N~1zb9r8KIUNzhg+V< zN&85|te4=>jKxJLh!9=o{Yacd!?K8)P;-r7Q*LfK<3`?5z7*kOd(3Lt_D3Ld zy8;JEcD6b=(Z1ia9O=VZmis<+N7JYuX!@RW35&~(4mDeZOMOmHPq8r;x>j<^)*6QQ zQaHypptY|3|5l^P)|JePC6X}9;Y3TPVb#)J8lATBvW26xE#iJNvC6T9L&H&6xwo<6 zfl-BHb3zCHlJyT9<|<$*#pd<_b?30~rG~vl*=cPFTo^5aMNca7kq*ES^iF1MZyD%O zz@aC9%MQRHdcUjC`)|MWVLwla29M-V;3y+IqH2EeJS(1A<@Q@7jDZKL1CXtwMYKtaYy)4$>u4Q@ zg`VMvlhn^Y9=8PYgC>w4Bp@4CZ%!#zKM4T&+%!XQlLZ`i&$6r4RWlnXGHmq}7;=73 z0lE1&pLefyqyvQkoOfA1i^C|G4u1ODVDM)@7!2NkKNT(V>|(Ot-<$)8XP3MEd^Fq+ zx7+O=;M-qb4!8U5#qve5y`1k}+>9>gFD5(UU*Y}7FdjQQ=s|5)6OL8O4M)-B4UJ7nsy8!BSuymGy&#aU zz|&WfXTaeFo0*dk@CYQ~k--|kVF3}vmAFocp&-1lUL$bvI$lgS$tDOFK~lh1i?D8r zVee{l^V&`3a`ki6tIgCw1;=__R4~F5*F2U3M!0%*jp@X(F@bKK9%P$`NeVE+=ki1x zca3Fi3&|t4YR0)-u#`k2w`pLJED2}xsi5f*`5Y1l>8Rib8|nPAN8jq|ltK+d5<71} z5?G6Z4k`AlYo1qIJTS?td8H>h)MC@ubc|G&xt#F4Bwu*&03r?vi;-2>BxVM zA}@1LfqHXf2RzAGENBjt{ALh<jWRYb+^j%XKSX_$5vROYs|+!_ zUS$L!K6A}X=Ig5!M?6&u@n8B?GKoEuPUCnMrIFHkYAqF{BT{7{Xi?HOPA}J(7NC-* z)Z{!}g{c%!+)_x$srp3ZLvu4yofveeJi z1o4ZrNfpErJU#6!lcpF(g+y zZmG?6H)arZMw_JH+`wV>a0ED$_f-eZ;15yd@>z!ILuTfFt5bc0ul1{Dh9PPfPIDX_ z(5%=S8ved)uhT`WWjLdFp2{P=P)V7CqvA!Zwjj2u22F}U9CXU6%4sS*RTPzSit$mc zZG9577MzCFslLHGJzXRv>=C4VRcEzR`59vrzRNhksPDL>F#<9tNza$FA4;}Rm8nr- zfksy86&%r#4vU;nvq2JJw&@}vG{8yo;x=7tw{dP~4P=G7) z^C^UvX&h@J3K`kFwHTg5@jpr=xXkpz-)&-#ljMR!15->u;atz*TA5%KJhVM>swj?u zFj6ZWk^=k_T&L+HPk|Q1#e=D19Ko$_Y$k5ULQ1Zw9O9;%3OX(SZqF$lKY)ZFjuM1W zNrf_OTYznR0NWg48xw37^WB&O6KSak5L}{#O;sB+JSGca&W0p%cmqLZ6aUrr9Eu*1 z_l2qowyw=kG%1x{2B1wMLS-Is2X0}Jt)q047eTm=*FirF-*000{u3D1V@ealO|2{q z;K8pe2pXKOm+Qlq4v(p9L%^Yo#|LCSF+9%@w4r%uD^cN+Pju8@#*2Z@!Y|39|2 zX%m^;$3gzzo}4dqkiH*IXK|LRh_byno)V-qQ!c59LQ*PS<*MdSO0bU38S4BYeqCn> zM!UdkR7&8V{Dy@5Nyr2Xi6oo#9)MO|so@}vtM&x#Hyhz($N6FE5I`gnVxJ(L4&Ffe z(11e9ahGW1=kTne+9Y)dh?ip1uoRXZX||iFUYgUaCQen(CsgXqd7Toa7s}+Mpl-p_ zR={Nsr?GGZI@8|j(0#Fl?&`YaR~l#Afmj#3@*M@;ulKb#%H5I(pbM2=uO@!o=raz! zES=`f5vOTctcvGQNFE zrz1mM^CW?nrh(s2Bm+Va*Rsf2+DchoID$B{Jo=5L5SKeUij?P4NyI@u80qF=F_XsV z2tEvFN++u6mShWL3Ao#Yg`dzA78<_?^OwylHHDw5tG4o&hZ_Vm3+#`24iNm zXOIS@qyrW78%uTSWLSe)S;d;dJ@z`k#61L5T0qNkY!9bVZB)5egOHDJZKpB2k0N~^ z>8G5JwJ?v8V3to6Pqy=^Mnk0*pJCYlMma)hsZlDgtYx`|7xgN<$_SGB0dC1Ia!6dS z>X$w$X3zV1NK!X%WoO-ZviU4QQhJaF)1Kqs9>gIwZPh3x@pR);U?wNZu{S1~Fm^)f zH6EN4AYhnBnB!Tl(Jy^9F}rom+qZhwBLl2zO?^Dok|~Feo!X>$$Wc4fQl)${42Sfs zJsuqvDvreS4Z8wDP&v{jV{#Z&o7MBMHGnaR%lf05pYr&inKP*OcaECd{QG_986jH> zaW77k=RUJzAx~Q_$1li5E?cpH`X!;-HeyXUD0J;&1ph`T>WN=|?Lwu^WBIKn-xCj?PP-Ol3Hb*nI+9jJuTXp%v zjBQmz@A~R!#{IQJ^tuhmMG%ujh@l^)p+>JG2Tx!UZIM6Zpr5rANIQ*blY)q{sG@-Q z5|mh4acE(OM`<=SZm@GlB(HL`AESS}6?zYAa&^#ZCRYRPsL7cP)8g5I6r~2O?>Yg} z?DBKM&(EuM0$p6Ct&~Ek#Z7wyxJPCy#+Ij4)>yn;O0A*nsC88&7fZnh=Wq^+EI2rB zAX(M$PLq2N{A3HzdiY762>;goh?BSjB9Gu~i;~!%YU1Qr6DKS4TB3Q)j|L?c6GBkf z6}_azTFN$c(vTy9E8HYEZ+O^)OJC>*Z8aJqJwXdr1XY{4H$YNt^beqh8Q#YQV@A|| zTtKwuie$8jq9R|!5A7!PB@bDz`tm&sTe|7ibN5>qq_a4PvRM?x)8o6}mUj*BNDI#* zq-^naMM0bp&!Hm@{hgmL-08?N$5moch!Y>d6l`zsF1-EIuMY;J%lUr3zu0asZ*C@| zJzz(3AN&1qe+gd}+wpLFIoU1Gz?s6c^F6MAvxFuZ+ZNM*NU2;hO9W#>o}=)T*7f=yIrPGWc|gCMjV+S#>` z8y#-X&gM%1aEWn@g-}^&N_d#@IHW}}KdPv<1 zjdQlcR@gn*(9RLy>ifq|0C%|R?*wN^x5R(ujec2PNSlU}5DYF*6-)|J&W}nl%#T6` zhn9<=v0KCnKc0lxyCDnXT8W@lMwFUkN0pY534#o3l~iVNsi7)q^0l0lY;SOZ4`Jt~ zXXJU5?Ky1q z#3_7DptaFtt~p7PUsI~))L6svm(Z?P1@8@dpq7loH?o;yK{vqjQRo}pD#Vv1{Yh)Pb#GyKNk;nsCi9Ds+* zH8Y5u_2D+IGVuA`CI|U|gK2>ZBan%pdpJ6ili1pkBp1uNdz5o~w0cG8J2a#IPtRH6FrFnihZ;sv$R)9EdZU0N$`6s+2@}pe z9#1rfP+heTLJ~v0oraCJUPrZikc3*0>;cYav#tiESP#B|xOt^U8g5+mp4+Q_*pDrQ zxm%nBC==_4YF4L3r(kTBGHqgYVG7{lYBvu1ZC;5(OgxGxg+oQknzE=~@PU(8#KV*jYNR@19FZfn=w>69Pp6-H9aAV{(p*(ZTo zU)7d%J_PqUIy<|$)Q|{Sk-7Bt9JdbN9WK$W5OXVE`x1a=hLQJokPoqQY?|cLSrleL z+34xEbIfLm8I(5^=m&Xt)Dx^QLmWNylB8V1*a$(>AdITik%v&D1U65EbGHD)74hhx zwXn1bG&AwQU7p$I?Nhdz@~Eth z5;#YV;_z}Yhgl4G9ugZ9^%}IU(4&J^D|#qM!eYQE%d<_CtUTm52w6_#V5e0dcZ42|g&^7dCGd{V`(gGGRh@iL%!qH2 zWwaf7ng{(bgFFrEypB)|WLu}rG&k(%rI^w<-gx`EW~U(3NU_wbQ>jKmXkxouEo6N8 zo#SySPH>oye*J`PU{*yYKcRoTSR85g?~kVmpWX#$YLJs#3v4Np3xgS8ZTG zG@e$KXC7PB&C*u{v(>Gb(W7r`%)Wvf;@&)eU-wmvQ!ugF6ywHUeQhxKo&P)-{0jV` z!{KPVJe$856{E>~J3622&bPxC%k%yD?nQC2-2r~X^GPwEEHCHV?e_A#GAObdSC;za z6-2k}7^n$nmsKns{u@C*j}y{Ldo{X9D{dBxql~;el_`UN`E+3a>l4b9&8szzHsi_w z-!Jv3X8i{PzK=0{U5+;72nA`5l=CDF{O$mb3ot(gpao}~%^Un6I6)2q&gM(4d(PEI z*Va^XoGD zq4>oPj#JA3|9*oq)Wry?)rRZ-5^R-1tP%`jc>ABe&J)D@2|#mU5>9TAa5CKQFSp|t z+e;+Q?i6DP2`78nUxm}AEMjOA2NsDX6q-Pw4_G;<3gp;K=@Yh+UZ6-uW1`T-yQ*vZ%oZX*6+O9ka~k zmSG7pVMzf>8J=$3)Q7N`yeFjMsjw-M)B)bS0>9*AkPoWO1NY$?S2_qc%d6fqba8tu zx1|X=Kr4-K2qcfVoztyx3VDR37zNynp=3^#8g=)RBN5s#rx=v7%I7p#oU)T*B2A2t z7^y6WKwbsR42DT&tZZqr6f|b<^qi2xl6iJ zDAjzVc(R60#-n2Pyq}$&h*!RO=o`DKKV`El4igAhl7x%P6&Zcp(l^qE9cgRrWeXQb z(14C^f0LX7R-YYe)z~T|p`Xsc`3n^EApqhG-M zAcEFWR<~j03R;GvLX?G)nI~p*7=EcQfKLuBuW_{lg_`opJpufV0pLd&z)o~LR9uKL zT}YxJxF3FQy>x;%-KjD?so^r+sI0sr#V@Ll$Iq%$5Qj{t6e?bE{Ep74Y6JXQ4+B(t zldu#P&w?lm^MS7b#Dx-7*P!LeCy;KU0u_ysZmr=kk4wX+Sbc^Ms#BY1qNL$SObb|N zR4PgzuRoUM0;M|v45nY|E5gtUz*?$)0>G=Hk{MFFZa<<>C4OGS4Oo=5+bo@D^HK?$*(o0H8SrQm zxJ!gFEKL)OnTTPiq-bN~4;7-G->smLZ{mJY%K?$$;9?E7SAE>L%+EEmt3}kJ`6vGh z^7D^&9s^-&Wt98VS(4;xP4Q$|?9Yeun__-Bn$Jhe$!NE`yx78X=cD0dJm2jnqoLl| zvwJa;wR$u4OU+%ywKMMEq-#8JuzyZ(;pDx+G>(Y(lVBqoz)>llHV!vYzeM~1wu6f` zmGr6$H~A+8X8^Yxd9*}A!}jBT_ zrVk{|d~mJi;M?W;gaeh&t5q^)j>m7 zo`-5}yzYbF>ryd|hV5KZGPNCTay*Z?U}1UQxP!$hy169Jya_kK#$P9!e3KUGX1a;` zOdX7IY3waAb{nwJE9ia{BfeL@)v!Vg@?@Ci>6`BGprpeVHcO$|GBLDWbpBn z-Eqm9Pcc0zQ0-0S=D4HY=kaT0(Z%qb8#;tY?v6d@f*c!WT@RU~A6)6MC8V!<&sx}b z`>93I>bX41lG!v?gK}L(wMXdRh$X2pA@|2JLXMvFcRZmP=~5qaqExqGWm25vQ3B7N zTYAlx0kcNCygg<|Bb2WhcGh!!aS%tfS?`)OMBJd1Ewrcpptp@`!*mgrW43Ku@hRGk{c?+^c|DTe+ws&avD2VrzCF;Az3eXDUDdgW-m?iy`Yxp!Widc@GyHOZ9T3?%OipnRl;_e}@ z&ue~=EgIy8SM)%;wch+^6|tdPOFRbPQlAdAC>Z>qA~t;b2@c2=^U1gvPR_^Yqupe5 zIocPK%bQ^_c`>;ep6w^#J1+}J4@SHBa59;U6i1n~uqlc2D1g`TpcIa)PcGFsEy22F zc1P7Lb;n?XtD|sbFnmF~61ii5VO?>w%X2Sg%E8qR8wUEyb!JNI@$Y}HDd}H9qwC5{ z;e~(E6t^S!NVrP16R@%o&M2Wi9^C0wCrHP?E$4Av&P~S z(~7t?dUT`2z4lHYT&*dTSY7#m-KUIf^{mw)XaLw%b1z6yP&Q8eIJ+N!O0=Huf7SWn z1&uW(gQx*mu{xz&sQ`26;;PnyCzbZz;{#Zyf}I9XV}vSE2P++wV#BIzy8g3|mrd7h z8Qyf3MaWP|gD4sNt*4CNpS};PnUm$k_WbN-vcKHj3@=7Qh+Veu`#JK7CIvoO?B{#* zpT*6-SeDHJr;-qIsVH6R8c0XMDGfVyd16%=BW(kr6kbrmh+D;98vJVFz_;LMKzPvj zDPKBgRGpBaR%jyCS_!S$-7NB21a^~U!&%pqsbf+%d^Kwq42pkr{L zM4oh9Wk!^#;{*?hDca;v^|&hv%hpayl9)Bg?3C6!H2Ec1PRg;{8DJc*vG>{=#RGcu zF%&l^t?-M-;+X0fBsZ$($XC~Fg}S);!@kmONT?0s#8_k2c^HWLZc2HGv>Nb#;Yu~Q zRj2ToL5|yaY^L$tkC0Li-Qhhe=mf$(wi4Umm1gVWLG^+UuGCbvsjk{Hf!E45f8KZa zKb-~B6wkgm%y@$97POW}i_)^VPjtZv0k8{5zJQ?}u}omQ)2$&5$ITbH91CAzX1EgCCn z9H_LLL!Ab_ed4f8)~vv|YCq@W8r$=Fo3*X6-63l8oxaOLAkHNEx}W8$lhBO_J)B#P z#Y@O|*EigDTUaC^g$9H=+jQ7zAyRD1(IR*ff%ZrvpV{v4EQtmq6G_lqGOH%O1elq6 zZMR}pMUV%Wp|os`zW(QZ>FWWUXVG*zOTs9*%Re4(+BrB{WmkQo4U!<^)_B*6gQ3d4 zON}AyZddWX8CXCX=auMGzK#>WnPCOFRB za(p@5P4*Ytv*CESyjd=nyW#n0emPg3&DFXZif>}jwkQ_i>SZZgQ`MKX|CFqNqrQb^ zkOl35SK#69(gmxQFVq@b%!IK$YJj&WJ{{eaVgFmWYj3^tGqpO&pu(16GH^~UKhee1`pwCg8I0aE>$vq0SQ%msRF)eo}*0tp>NLJnO&BVdud=Wg*?&5=s9qcD4 zFWk`BKE&<1xV|px3VYx+e@)@`2OmcDh8|0K;#rhHu9k(NVpn$CqBNWnbCM0jB8HI5 ztNlc|Q&oj2wYbEjgf&I?n9;y(b~Vd$H);>(CJ!O2<<~XH?o$`}NiU`zYis9$kHa9r zluhk48~yBf%54+#&J~Ag^wAIH7ljf6mKEJndr8%jBe#l?e>>2UXdOS0nB_6t>T6w2 zm|E){s(P5>$PpO@`W$&H2o%WEwc;F$gPrz37 zc6229y`CGeXxCgG`m<>kDK*P(b&<(dua}lQUdl<9OZr~6jxV*T2`y`17K*EmDQ{Tf z_X{1(Q2UKALvfOVJ|X*xMnE?E_FN9GOOW5~=_2t|^C>3rA!qAH zvJ^}vcMK#H9)ea_l@CitsE{)RW)B(yoT9tDckax372O!R9OQG0}|Ca7pPC!BRO-B0B#TC*uUic@I5pDOuMvVYEenu|o>GF1ygW>}H68DKbUxXhZ_kL1H&%4KQC+VMPQl^Gaga%A z(Ha8UGAH!PG2Z8Ble8Dl&tZ*N>iuMBRF14!mF|-Q??J9r&9d*+kRcB)!+o=JxB1W&P64=63aC0G5C6wq>}ub`H0KG*)zK zj@SS0#j*^-0L4(^2-hSHzkU;QNkYj34bzlf@+bt~S<9EmC+&ev|#6#TFiTgUkIP?8}I#S1MNR^I- zYgir&$_fYFBo06F25OS+rtcsrHJ3g3lbG&O`-*{W#iGw>W1T%-OR zb-Qt0IO?HelK^!!TwV^30Fj1+AlVwVCy?84*h+O?70|k`ae1e&_!kD6SscycX>wGs z?Zus#(GalI%`}uxkgnc_YZENeQS<${tp!X-@>TrM7HxKbtv+JS2zQM69sE7Rglb=3Rq8`;!)dpi>~0LT z6JskbiIh&xnj)o^tpyh#-jtG+6@(-m2I3(|z4I8+v^mNo@ddkFnyXghx(-Be6NQ4#9xW!+R6-ftpGF=zAA z&B2iB#{fis;VxAMoHjV28vJg*3IC9WSff0dk4UXz=f<#@bA)7#zvT!*HIF-vh% zE}e5IS+LNMEOd$B2YHhjDQ>tnA>oiVA2op-mv*Xxv3hE^nX8mH82VD5gOUqTfw-YX z-#Bc;g-p^X!TU*uSVTo!1jTfn6w$N@{37V52xgfmS{g<3+Rch|^>YAw1BD3-oIW6N_E?e8vuTJh_YZ&Se2wpEw`2@q2(c?#3ey64murb zzH!S=>*Y&L8Dc4^MuwrOwP=EP1WMjcE@@pgBZ>Njm^P!MEA(h9ZT?oHtq)a_ZeHnd zv{hgAuDLqG?OM;(9R!@e-gj9~It!vCOJ;EpB+S;^bp=gTDa#bDWfG^~u4~VXrof>| zTIY#%G^ffG?gs6lVwWSJH8w06pxq~S`Eg%{Wjc#}9|D+k8pW!QVAjQ!cej9W{|6~T zm0Dn~qP^-3tHy4u!)fvAL-H(vqpXp8Lz+c&CuU*3>ACIkFm^ct+w+sQR%?)hlcX9- zQ6}f8#8*{UxcXUErT@o02z!2-VR9aD`=^`y7rfKLWsmC4JFKUt*T>=Tu;XdJX_Q!+ zW8(@>sdIVE5U3G>>paDyo0>^xK51&S`XvLH5WKRWf)5$9D=L^3&j$Yw(mTEy>iz%# zABzYC000000RIL6LPG)oO#_{j--_Et6voG!O|JxLj1Y#o84C-ff`cXhku7|Y8BHZ^ zlHC$pTXsb$NqyB^^c`$O(c|=@^bLBG>pVa&3c1Z2q~95fu9Br`HsjG)=R0GKIX`{p zw;nU&?Sm8c&HLy59{c+AC&cgDLAh-?ema?i(IT0Kv!uwbpG8VtT_wvT3d1m&$I~zg z<0y)m_=+)A>rpLb;i|f})yT|RyQ6xj-DYI>l`XZam9HKm)yK{t@B-n90K*epc!9Zd zy!OWPTFC@~`1c+^G6|%)z#@d>SVqP}gkxFA3&A-TsjA96rJ2fH=4mP=yHr)0m!*V| zmh4m*s6i+9sM=`s~s`NX3c9FA9UP*RQ?0Uu7uJ`Ip{oHFSbaB38ug=xZ z#J?a&H@8pPZYr6VDYTMZZn~$+?e0ast;Q_9uB6Pe?8(hw+Z&BtURRyW?=7>}H_X^w zm_^ehoX+ATTFhrhX1>;dF)XAW*|u-{S{2&14(8tvNB z62uKO0`5`8feh9Q$S@w~0pJ;MTyT5?HC$y;p6n(78#}qx8kBZZ%7h+w;5e)G9F-KW zekw~8wIqJ1VM*mZ@(brb=20las0^2Rv}!@l=8MVl(d;upo_0arGax@Y_9Fv!{sKUL z-vPM>KN)8LxgK{6D2)N3QywAvqXqepeK%!^DfRiv@B+(pQ=kY0Hf8&YppZ>J} z++n@Ct0>HP95OaFP t%IyDI=RW~ep;PY$001A02m}BC000301^_}s0stET0{{R300000005S;&Vc{` literal 0 HcmV?d00001 diff --git a/src/sctools/test/data/unsorted.bam b/src/sctools/test/data/unsorted.bam new file mode 100644 index 0000000000000000000000000000000000000000..746c2e8f557fac658c11c75fd63c4feaed273145 GIT binary patch literal 32984 zcmV)eK&HPRiwFb&00000{{{d;LjnLX4Bc3JjGaXlpY3+r?p7W}ASIw&|1bpO?tU}h zJXaI?0BW(jg)Wv7G3(vE-?mroW4(Kqwi-|)!H9`OsrbNv*a}8M3?ynm6v7{9D#qXg zf*6I0B*tLW1jJy_Idi}HI_H$)|8BC`o!@WfeCIXi%*>X}(?hp!WNg#>E*{-IGsP#= z=(e3xoo=TXo!vRL(z&+N-QOA8G<(J9?99}{z8)XlF*7BkbxKPi5<(0J!CkI|ONbl@ z?yPW1sg#fdA`@JCCLvTHELXV|XoU`h$!w;aNeL4O<#OqXJRxi#q~f{NxlRZdh}`B* z8g5bpvjLMEXRw8w4;YnMtuZ6YHNvPoOBLe->$%Le%!N&vf%%!qWNtWO2KHy%InFJT z798NR-0LkJh>Ff+=JO>Xt_}-~5JvJ8pAHMO(v{KMc4Y8lCkE-2Ig?S&MUGRB zJ2^3knd4GZPL-S$)@j6vObujm>ob-U&JNZ&p)D3EcwUP1TO+4jR3$e_v`bExpEFu# zl(SK$Pa;D&S8~Sr{e?M^RqJdf2**oKsoW8nLRK=h%L!MyOmNN(<#d#p9XH|q+l9&U zOq9!fxCs*-=A{Ed-(O@4IN2pRc=ZMx^Ra$GpL{t41A&eO? zxzR4eRE!7gD#Aqfc#9|%1MGJYX7Vb$ig+O}eJV6HHLp@o6I*jteTOp9)=UuH^3qn$ z1@U$sQ&tKhO=c%!bm@fL*H6O49PverFfm80A~~WiP4pf^X zmqc@nrAuOhsN4`H8Y<{su+dOKx1-CVLoJc!q6?%Y+tJ0nZAjc)^gvE)RqPnP&k!PR z5#P!Q6Fp!pUL#}K5n5`(#AM^f5hkXTw3I3RDKesgj6DJlQiQ1zM-XYE@+aB8GS(ws z-iS0&y3JX-q^j1WPDxR9SrQys)wWy{X|WnPk(@J1DrOEBifmvhV+nC0qatOw&Nzq5 zrPNVDLzdj#d+nv>?BV%7HNAsV4rw*Af^i5knc zcLy;jR%d?*#Mx-Fs@W*I!Aw@=qdzwhs_HjJTgsKXYi%`ll8Kv4`kRF&6pzhSaJnj% zPors}Ht~WWGec^tiFMzXd2R?5OO*B}MBZ7lT)p3jF{um?)Rb48y zrAo%8;kz+8I+*HCk--6(aM2}Bp13Xs+j367Cgax>4R+O5jIcBS@h0N0)j9Eqh#rwv zkTaHvwS-$gP0@w165GsDq192`_asumOrA(Hn)ifeSfMZCjx)tk1jEiZq(Sq zed9PGb~4^;ct%gBrmDg|WvRbp^oef;#3xe6y)Q|pGBUcvS?^WCMVI9C>||uzJid<7 z98%Tg$66XVmE;L&e@4Z=ZB8f}s`QB#j$B23wpGQ8nypHY!)wRv=Qda7OYd(DOAxkkNTpX(RP{jp7#2X^+#&AUb~pBen!J@`8(Ce%b; zpKEs;&BbQX7~Q%f;I&$9e7v{PneH}3HW4y84D8}2(PJg^nTsIzPj28!rg_VA@+nGEOTAtjw)Ni&| zmM8rNx5Yx9kIygY#<(<%#c|un7RPhFpwxn?%WPgvF4XVn0mU7bkaD`ra{paL|_EpbvDGyS=`5GMZ5D`*_c<_(<0_w_2E6 z>9p%hOMb;w&Hlb-2a`*^X1m@ycxkJ?7glcZsWQ9J?G2cgKNZzdGfY zW92}zXD^$G1H)AN6Qkw8vL_DTy{;TM_MwL_TDPYBXvQXP``6ky;_U5PPN_J-?z`jf z6_>MbpZ4UMIFjt957w)ZVrOkVvl?l3VSj!kjtr~sdu4STS@zDwqt(c<)4Mkemm`s} zwO2h|jX3+m!#`aWM}i$`uN^K;*ux%qeMdCzYW8K@jLbB9a1&$e{owfPIvdZ2Q^I+8+g$6)ubMi8jl%@71%d~p~-AnOGr9_|9{2o~>~Hk-ON zg@SW7Z9ou+e_<9u;MLQSUU1H3{WS;zvp+Y6ATWCmQUXq;Y$KKhs;g;@`o4G+TXUZe!P1hWt0 zg<$q&qX>d8E+HG?(U&*}-UrxO2!h)`yBQe-*M7GDbOeE3Pa##|qR&2tA<*mhD4F09 zm;D>hf^ZCB2vmI(>jkQwfjj{>em0BDhO0h%0?z{Ne}=&jRg_;#r{mexw9E zz_LB5bHd9j+dYOL_~PxW5d`8N97PbA{Th}9X4lst2+STqvck67fMkVj^{o*EfsMZz zMi3}*4H6&9&^1Vjuph5Lvck49*AAr+sCw`RDFjz4cKvH3DFn0Z`Kc^;1!aF-gCOwg zlgLJR_hgr!h9HpjbP9!+PWItPt9FosA%nb0AVXLHo`@5aj9s zw7V{7nQwVm@*m zS_s?31Z49}(>J_vQwErU#AIi|D zzlI>te(W$>2--jU9W)EHf7N%L1On~%-+EpOf%a>+uTG$_74AukRd_CA$9AM&gjX`Q z4S6-3I^O$ZGz+wU=tcyA_V2j|EevP5$I&d%{*EspB|!Vh8xaKB54|IW!a?M~7gF)? zvc+CqOvQ)lyZj7-K>K5QIE6s_jbB8wK>Ndok$BMl`@={H(0=DSBn!0vYzxf-?JsDh z;^8fd^)5jWX#W)QJ-j@zdMi}|UYyu-FCz%F{|`!Ecwu78w;%|#e^IA;!HW``J$s1# z7uQ(}7=$JO03VA81ONa4009360763o03|uyeM^jN+j&-NJTS=Cj0cbUF%``kX=)qQ#$@9+Ho`XAr_{lD_0 zM{nQy*wfx$`N3$;5LV zc;^K{IEjGfXf^}h5Udp!X^bu*d{2BXdVVtO~qM&s#bbUEH$ZU%Sr%iZPnF1y-n zr{m4wa-2=a^Wk)}F>QDg(w6WW{zq@{ffQ2O`ZZ|_Q<}E0i#KT>eUtb@EHRtDK3gtl z>ukB6WwR`sE!VyGO<&K}vl;xiTxY%a%wJEF-tw=#p1nKrV&6&M5%%7{czyOa=I?sG z?7g~tJ?lMRefY`hT_1k(Gw+)IgQxIx+>bx)8UN;K@A-Q^-20jLOh1ed)ci;A$?VnU zX|KKdI(`bT!pX0ePtIQZ@Q0SKy?)c1ec<&hTPzmOinrct#cRFyec<)dS3afi``<4k zem{Bnr1whiRQ%$J>w3W?OoL$J`M#I-#G87()o8li4yM=iYdarZZ-?{I?s{-J9c1&- zd^(<9z?Uw!S+==YZMN4VDb2nuENKvM{0`o%;_1U{_?Zv@gDFqItAJISmWjB*LxU%o zEtl{d*UQ-qal_{ZzcP51P2A%9aC@SjP#LoyD9jos)i5ia3YdNH0hs;Ij+nXbB=SRN z;>VE_6K2+g4luhPM!-apkCXzT;0I0bZ0$4m*%rQJ2{JzukfU!m zCDGfvdpI{{)sp}-bcNt6mxx*E)Z>_C`!M^g z!tDEinGKjdZFjeM=o@`MocKu)^#<_GuOMtKEp5~9n}*pp`i6;yg$4wcCk%n^kU+ur zfQ)HN0K$}jgfuMtI>0F`fsPcv3%&qGY5-#BNAV(ruNsIPz=rN=nc>4+F4tfT2pc}S zHJbp!7M};U#z{5zSUMG8`#S|}4AHaO1E!W^+WWME?K=S5+dG(=i?^BM`V&7&9iMD% zFrH6W^TA+qJzfDs^b1}Zt)>HrL9*<6I0wY$qwVklzBHZ>&;trGF`GFb%LQFd+e#mkACHE-&zy7FrZMa=6mK^enpX*pGrq!E2U^TZD=#!Qcda1NH{N@AAhBN z^q4ScUhDyALEuFydOQ;0Z!f^=ro)TD=Hg;HpO18`A&pWrUrXF*=aboZW_YFrPGrC{SYiNy=0JG-uuE9gl{iT{a%B2HE&-yc%5W#$b=bY_uJWw$s6Q zJRYe_3hyS&tl?gg1}1PS@(>U?Y8#XaD3L>R&~eh2FcYw-HiBVWHp3i>yp+1y@FBws zvvb!*8mBd`xBTR4Z3MwODgM`1_sNden**W0b><)n%pYZIN5u?IOz*Xy`>__{z^K&_1gY?KVjp29jWJNtnj&Hal2sLwn+^q>qr~n^p91NATD@Efy5v3}nV4h%e zPKP9hD9wubdm{#lREl8YFJg+@DU3}s*gd1Do4V7 zR)t8EuCh{GDeOk_EeLJCA|ap=rxi!e1lhjgNVyJKEsi`S`u}>j+Og-wlPHN}h)E!n z6EY?ZEU=NEHuV-|%AdD126!m@PYRwQu3 zP|#`{lN^iW0i~e<)uATk6cqi>FoA;K@OLiwuays5pEb>HoK!C!7pFcDv`=*?Xd!5x zJ8|PQaR53uWJ?=tTSg(K=wks&a`-^ye}k(XT#&SMl3i(+){RaytWTD^G}wT|tAB}Z z37(k>9y5Bx=gqV0C)Iy~;(Z$7pYbP63 zJDFnbWI7t`Dzy`g*86?x)fkv0BOnL|7QIYFgE4}59;?WF1HvBxu<=irT#PmhcB}XY zr6IAML)B==01Ue4W_0t=ql72O_uP7oPja~=6oX|jOV{Dr^=59E!ACPM?o>kQmtpMV z%~wDBHFB)>btSMV4UVhBWLEhOV{sYgip)pJw4z}YpHn2o4(#K!bK&eyD!%K zT&I-IbzOikO#A?X%{bulNvj+jB2OW7l*l<8wQ=e-=km>FRjbM65Avizs!Orx!{IZV z5l;s{RCO;!DL8LRS39XzX)c|5>;`<>{WY)8ckIIGDa)CKaW_~y9o|LhBqH?>+VKsXz9Sqq0 z)KPm|EDkxzPj|}gT!>31o}VUIMM)j5e6{Th%+a{8x=a~4_TI=#f?G712TO*75r)*Q zI--v-T}bD774 z7m|S}p7>GdGHbGpf28_KABnE&QDbaK7Zy``EVL7XH7TeyldCTc ztSt8}s#*!#w21ywqHf8&;133A5v>lHuZFE^Uu6r+FF16W7(nvH_P7z`x?2Q!r115* z4qa19y#p+HIS};+p&xb(`)zfPX_H}ZY2tMxFqN|@9JDYQ1N6XU;)mco1$OCd3=#29 zfl0!5Q<{$HH`yK1m^tiUh$E|)Q+HOs_XA|>%I{k2;vY3 zBVAWfCi(WtPIi(UE@I4w0XrnzW}nKJdGdxm9ilD+8W_kMs&i_=D;T5@U(ev?5At4- zXU@?-PJk#w%c#4&=H!+yG~x-*bob2-)=9SvyAg;q4yy-|#lZmEPu~yQExvqv;7KO($WM~pfB>+`Q>j9U@1QV%DR?Q^xtv2KI37#m%b(E$ddlBhYXKr@aPO63crBe^{ zx0I>fX16!l&7wBFp)q@NH+r(`IukE&++-4_K8O6L&LQS)Z)s>hh!!Sr+KHYuUu~wv zev;iMa!!aM#Pc{BMl%Oaa}x>(8M8GN-ZckxrqeWf;mm1T{iH_GI6w77_gJSVe?^BO zAL>L;Mjnvs#z9(|xY=xmyVYp4+s+5O)#W(b?XGsCtKDwAxxUzpujliN?fi1J+tIJ_ zb~-dvxwnrkHw(Kfk{qKd17ohO8gN*+;tx%3*QkySzQEsa6T^=Fgv5R+vL+4A(=ujF zl?)B(9rUlBdo;^n?ZuhD3^Q*XW$qfAs53XJblidyp?qr@%C}PXRy)5|+A95B0Lphi zz-4yKxIHI@n=S}bYz2hfGHxnWp@|{2l~Qd{;oem}q|8eqhu)-4!kH+hiRP;DHKRcb zb<@$@MYcRM)ueV(J>xD;?Kr3S5X>k#fp>B0V~i4YoOdfbIEX2vE@M?${aMXk^@JJ;+1Y`@ z)>k{VMj@K?C!Q0A@F6NDsLZV;d$LcJU+lhMbJW0MBo;m(%2!n`Db=*0G^m+F6%$*Q z6fgvr`~`aw3{=ISvNDbA2-PQr%)!JysU;N~2)@AVEoytMD-G>+hUJRo^1M+<Ik% zE@6wr0x`|hm;exR{!Qws+Dtj1drV++0Ml(tc}HSJVIW5yixW4zTb47O1DtObxOP&5 zCYhhwsTTRH_r(VcODU3alth4@nF&g6b9#RNNq)m5DuN9tAsld8v-PKfu%1b{9vsMYya@)PO2gMGdFYHweKv` zu4IIl${49^pgw$R?%3@>-vtdXQu`E_i5wQMa<5o&YX9D|nG=$da%ksi7-p1NjE0+lw_ zdqc{HDR9+KtIy*XCwElf(wE)WIz}jN0s$p91U%P`6nm@fOpbcOrkM{MT9y83{>Ec* zmaauI^G%r&4%;&=xVwd&#|6tQsag0KDYw= z>^_muN-$u2Z!&V=52<)8ct(*)3^q|3tUR|vfhUVAocVOlQRzE3SqJGd=xlvY%hvbw zfc{-2Kd9*6%{u`3OA5%(bQobI>x(C;8@rP@RMqCaGdip5D{YHc+nO!nGse);)H9=$ z!=au6#+U-7oRC^44O160P3r6z#8}`#0>QH%H1<(cJ@Kb$0v@)%D4w7i6gy4on4$Ng zF#&(+FWq?QWWL+U3STT+;fs`!*UqeiSDX*vy}tu^T4CX5-@&f4-POmu*9a1GNbZ~< z?l~{?df)z9ulIfUXE_}WHUqG@-DsQbhQq;Tx4D|%Wt-u2d$$@5r*~tFakiK6B0xIc z4lL~%ZPi`XzSFAK|EmM1vbYUZtxI2qI-S?U!>T`r=`z)RS(<{?)rirvH#0wKeV>3w z_d7pvKw+VNT0 zP?TN<>PHGa(v3E1shWz-BxM5DBFQKo>L;xReup=pVD-sI?uA!oa@zF5_m#nZ{$Ujr z*Ta4)$}ISvo4Q9TDi?<;Dz?S)TPMA#P{vX%fLy3c+awFbQt^`&Zj95$zC(W*1m|jo zHz&UeA1XPD`rOec_Nq6Kysx@ z4bWSuf9yE>nW@}-jk(vE2fynlHS=+O>H~|HPdiQPV?r=-QzwO+%}bRi-U$u3tKncb zo=(Tv0I(g=IH?|r!tzL7=BLlGHfeB9P^8$ForfEB<038J65UI_0eU;NvA2D_og1Dq zrNU+9QiK1|plNo){nXps(4M~1jo<|@1SufrIer*-!>sJpgsy#*?zr-thQ%llU?%4a zjXTdYm}7+01|phUkB(Vgp6ld(+`XfUzumDP=Y$abqjtF;`KpQ9y5g*Wq3XX;iZ&|p zx{kurf{m&t`<``_cAH})rVdNLu*SF@;uZq^b;qY zgki#lwuMY=cR84@vgvR%osQ9rpEigmItoDz&lyK^MEy`bR16PhtQZ=3Yd`En3u*?B^~37*uHxW&@R&6M zkp^8MD+P~#)@{NRDJHHT2I0g{RUmRIvm|$t7IRA&n?h@gd)g{W5j35sG@%qV!Kf6Y zb0y{y-%`m2Str#=&gUPS)sf0-he~y%OCX}FBj4(zKS`x6FACBCWrNb_VY@kuQVm2j zKFB|*R{Tr-gS<@Ayyyy>`%&A?p*uVO>>X@$Z2(k}WbC7SLIk$gd*@TGZhh;8UT?e^ z&#yL@7pw7ZxLpmdz)!9QB>sw8PvZ<<%y!cq!k4Xf**v!)Q`5RQSGKB-lBN=6j;4YK zzsb^qqL_X2;DJTVT2wfB!&JAbs=bG+TO|k^Qv>;z_U-9t!|_7;dFTW+?OVUPu9nuOa|m8>usDpkq-qaOD&88mL6>cFQ= zI2oYn?kLoQ!nu;d_O!{We!tJ$I_jvQnOTO>M^m$O*OdLHN%dr_l+mFk32*k0DD-3! zk5Ulr0S8iPJWCOVj1um*IUG-S6iebD2U?Oz zB5MZ8B^b#5!dU{mr*NnnAB9d_@RUKsblCCr z_G&a8O(AI;kA@d3xZbYrhLDC2Mwi#uqwDEryqaF@#&^SP1X1ODY?`Icnn_tV3zC$B zL!oMqEVEE^pyjHHU7}j5((~8Bc`~`iS@luI;@qaA%!n;? z14&&oD0xPOF!+jglqE6Ib_@GlcoA@$1Q45q2gOA1?h_N`X?(uEP%a4DWu#ijyH${S zK=Y&quZTdVs~QRrB(|W}uh@#Ef;S+$O=LIZdTPh`UV*a=&os-pkLFd=6D`ray38_u zp{5tqx^NO8Zm25ihiYz^Bxmkcx#Y|>*S=Pz+SjHMYW9gM$%9qo@&xOeTeOmjZY+|f zw+t(vr-C=sDkO@6_M7Ao(Z%v-BxjC#hFZCoxTyFf@!sKT*P1t?@QA)M^OkNpbN$To zmR>gq8%({%xsBRa{&xTf|K@!V-exygTCPi+f3_4`b|^96PX6Mpt z%Ab;)sEH&a5}Y;O5P5i|BX>_8p5>{*DV=ZJo}+5cZN>>u^PB3uPd%weaCoNiA!`>CB$)VV6vW*Wge*dbwvrmQGAqW<8fNtyw`0X1WgR2UXo4%&ruA#`7j$G zMq7o^Zdk!+TJFS6A=%rS7_^&x6SvbEc$*2!V327cAZB-k)CLJHc zLp!8|#!qdgQi55*I{^$KEAZ20uy$vV+j*5p#Xk|oHp}RXjrS^p@WJCur3oK_zv3YO+AAqN6A&J_1v3u*#Is&y47E znG!ok)4(DWbql)n3YxgpPilyr=ck@HWpz&6a#63?OD1WW_}y_U_XV_@+`}iSTAHmV zJ*f6vh}l>FAhnKEH(YazOZG<0{!FEG=fbUVQr+H)Q;*?xK=<%B6mGxp_nv^8eBNrz zG2t;sH$+-2?qtN(D?j@}ulJw-Vz2j62zKVX;bxv)VEOrKIJ|;jWs6i;$OtZ`o85f3 zT5WdQtGm?{AiW+=cSDopy=338IMgtR#b61wRtGJPt0tuLi(Nd{NL1`FA~n!bB-C=| zKAxn~HFgXzpe1*?q~)YbX0LK~rwbi=QZ7sCopJYA1?+#p6a3}B|77o!W8=vB96@n> zt;v^Sm1@KH4NVg;Xg9@0cbkz=LkYD1-fc4k zi7_}SDn*50Yg`0tsq1x3+U3&VaxvHV4XPR|?w(ZR`V_-pCx;xP1`-EDq(yn4Wl3($ zs#JpY!HhS+*1s|b<=B+f@mucKbP=pj&}hpYa(!$MI8NgBgyPBekSDusNt8AE23&A@ zK#Jfp`_x6{cNR;WQ?O{UD!plFa12`gXX73Sjyp-X#zI~SiQY)n89%0pnu2+O>^o#9&vn4C1tx2q0<#RN0vp1%E&{RR=wAyV~I{DaC z!2z?=iYHxmM$iWEV3PO{r}$yuq`KpCJ-DuW$vJw->E4b{pR0|u9uxO`5O0APM7TO0y9}yFn-)(;3X?%54N@ANW8Adj9kC$w_&&|q&$#Xtul=BWBVzYxv#(b}_zv#Db9mXH=Ptes&bZJ?~uuT7!DJ#aHT=Z7-u9QX@S!JAh?4T;UXX`-_Ervka@|xa^X>g`T z9>ZuIbkl@tu0*OIR`-y_!4KT|@td87t^l?WqLaj*__3=bg<3IC5Jt($YNC(1GtLbd zn9|XWhCLay9wiD@$_9$Kqwtw#Zl6*cqD8w~-YCqW>=<)sk>z}^y_j()T(9~5{bCnt zt8N^%P~#-Pd44aULNo6Ao|o>m>-h=!fT5Op(=~)eWW0<_Zh4!ujhJ`kU@mR2L=<** zFl!*KQtk=~$=gs#LEX-0B=dOkXt~StXYQ*@C{(cweG^Kqo18XK5U!ndw2WsR$FG_SdV8h{o*$elxHq1x zO@q{Lxqn{_@wB4GVq-3JLR8vzQqBfS<2frx(pIZ+WsA$*Xt3K2H`|NNaLij=@KVDG zV8z+vC=ka1c<%>Oh=e9zDqRePv>u7Hur_A!1OlcM%z_rO#Em^R#N)>6#QEV(@vKjs z;;G+qdpGXrTkcOfZv6J*N$g=nnZ#b)P1W91wP6;wFDqzp-PzPD!PPCMrn2<)b|4$D za$l(#U>K8=6GM8qf^AJt^Sqps`biDGSAOcg%~p!kbtp0^_Xd8kBiGVn27^QIm0nG& zmC*pG%=yWt$kRC9*1^HNs&z@BmAa%=LP`*SsW*^QTpD}nfK9%@{D_08sTJh z!)7_0tfc|f`P4E|sKi76cM7!OE8W(#J8-XgVE|u>JxBv|XtF(BG><0cMsS-!q}~*{ z;AU{-FPI!;VycHrI%!tYG|9UP?K)|4RniuhqSA~?PORy~bNy7GI(Lv#KdFvcaq96} z>&>0EZsb)y<@}Jbx|u;62H~@RY~AW7W)Bddx3h@tZc%?A!CEK z-Ksl|{@LTxTk8#gB*S%fs!-#kMu3u^+NvwCXn^*YJ52FtCB(>cX)~F}RG}79%p6sc z2o$T;Bd^4wnMX_OXIPfHuB1>hCVX|}G>4L@VL!m&`^+*;wlra#SYZBPrUMn3q5IN- zXE}6D(6n50uEOvyzx-nW+FN=r^jgjAl=^ROn)4L&%xg3p0+7S;cr_kh-d$YW4KJ?lrhFTr{DN)K=oAuycmoIx?Eo&Mu7IhW z#S{OwT0TvbAf?d2IiH|xfUY&QQUZPgO%KozhLX8kDw!zW6vmo*zO(d}o|mZ>!SdDT z$D)#d3lHzRKX@6Vk{2|DPXWRgXojG6UcEt3oUKY;YF_77ftOV?x9N@j^L_Ja~9Opl%cfD+WZS7*ApqO}9{w zFLOWRoN9|r?uFgmF5nky4K3oEQip&+!(uc2y3IlBs#*gtxPeb(SS|3?M!255C1u&@ z{HR!BVG>~UO`5`y7^9GAkj2b-Xf!VdBCW1vI6!Nr!`LU)R`S=l)2?DOKi(BEXzWHB#ja-+8bq5B?_)Vrdm=DGM}9bG+Fr7b(sN1Z)edXWmh*#a zHA+a_ZYbQ)D^uUbn+9^l41)DGU%&ZlFLx;}yV%9?Qr{0Jeh|e41vZ=GB1T1~db7xxMOn8g{)VVmSG#hmR$V?EE% zT=P=pm%A?7aVN0@$({Hs$T4s!~yZ8|4%*ho zByWm!Cce}AI(++kUxFljINeQmS52PKU^j#hv(5Ej zGX%4|0Cx(nPItJ1tZ!-UWus5c0jg?o>szS7Z<_n2J~%B!CCM5sK8pL0)LL<>cG4}0 znY-SWy5t_qlEOR9I>{*9Q!A$ZWbJ0|47F)r+Yd#8TxN*U4u3d81 z$M7>*HbW@s45FSWJ!V_GMr^*B7FROKg{nr7BB5}0kJ8OjDPwZi4m&=IB6~bkao5qA zR;KO?R4O!-LxIZ4>WB|n;6h`?9j(@ByzVaEk7+yRwpX|1HS?Rw%(mIZAzA;sL= z4A-$qW#;j8O8VIq7}g~=dS&a?Aq58R&_QN!bCWIV!R0Nr%5QGstoLPw+5hN3V2IWF zO(#XJW$Z;9UFzL&ySwY{cbgrUT6T49QI;*J zzO44O>+Tu{9wu!nGl`7H-@?Q$OOs2KLLylwxATa@XR*wlV4ZaQ_rK;5D1EspJ~+2k z-FWFkcoJpt8a@n`KDO!NGes%&)9dRUi@yl~`Mk1n4iY@T@mYMcxT#-wS&W-QG8Ea+ zvfF>qZGD+bgW*vK?<0&-zNI)|?}Q)Qc6-^zR{l@h4SU7voOr3gx!?HId!=t`t7O7_ z1FbKsRL0t@FWa-cxJVhkU(AWOlW@lUoH4GO!0mGiJ^}x~G;WuG+vaj}LAYHj+(y+A zamWh?Bivk!DtoL}$_a$fMrLw_k`{+>1x)ao11f$o<&WA#)c8TZo<-?9&2!lXDa=;d ztx9TAZ+F3YUv{Eyv`&`c(skBuxc1Y|SiV%o@}+YBRUh4;h2?(@!18Y$ph0V0o9rHK z=kIi^a>5ItpGKHGglg=#jYKf!5N7qdE7@X^l$W^b^*IAtnQIeEg5$8LoCALArZb(XA^6wxI@)0Zf4iM(cJ4uSdM1S`$O zq28#%Nnr9|=Y=szJa!#`a&MiqY=mA8Pkz(jU8__9rr>X;mUdD}I8(FY!8$$HCbRlk z4Ox`@+(QY2-d6Co8?2e#-j+>mkv){3LYjL)026M`r{)45{=^}A${{tWze)jQAEO5R$QJ#+q+TlVD``r^ zV1ZGJFPTNl^K?RqZJV@5Nf5t}rm@J_kfj6xifKIs+@2pbwjfyt-JomSY>mTeH(S-- z@_7(F6kdkzfb!69)p`i+GCG$sIxh^dVea^DFMx0V)=P}vZVbWx)oN6<5LUaZ9giAc z-_fY?V5?dPyVY*IqfG@k%TL?|WIa;>OfntT!ci@pV8v<}s}y0%z*`GcwNom=d;Bpd zw`igmzGSrky~>OO^Ngfj^~5*@2TN6dhO~U`Igp~e9yS!b%^1H-^>N+>e00> zLg}wj&CmIM!?yUdMcf=g#?&3dl%i$+N8PrWp&`kE1x&c z6`D(ni)D#)#8S>SW$?@lW6MQTkR{P2)ZZ-$l`mzu9Wxcx)5} z6Q*W{t^lE0&i$*J%veZ^_ zkpmxMjksk9Qne`L%o9ef4%OKmZnn7%Vo!A*Wa_wHRr=}(N=uiqJu5v6$+H+T2`oTh z*%U@rR?=XZEBIYivFsDLLD+0kIf4$}K}7-4mL7BuS-bYUCakFI7eoYP$=fwt;f>!u z`-RUf*jhHU$#kmeBI?iWGm5a-04gYTtjuaKpW^C~n^R|LEK=DY;_Q{W^jo^4aG3h- zxQZ|pCB7w{FTK$-I|OwFWBgPP|Ac5ub%h>m%f`5_GqzA7h^WzTPUb}q@|1zlABL;? zl#WWRGM$M{{u;@4otw$5v)6b%4#!Uo{&Fn-b7cZ^bN+3Mhb$#2^7s!zXEeZ@UVJ`p zp$51>7jtN6APh#kxVH+79Xt2IlS8E5e!&Prd1b;`x@rN@fv%%oj`dZM)`l)ubKy+O zuPYqgWEy+AO*o+nw|Wi_z`|k5`EvmSYyu~dLmW+&t^SZ=@6iNg!!&e$ctb3Vmg(H9 z(_GA$hEFYV`?@TdVMKTaN~ryx8wpg}rIo^Pv973h4z1!N(dA6lvwpNkUTuAAwQF@# zKI|ZymC+qyMKc;6x7eVOOLgvrU~^QkR9_N*Ev-$tO&HwrdGUeX@)CD2!%_{3F3Oo9 zN)QW5F^p>aelGz?S!ibSd1Ck8>fiB*@3XSTbkYY zvAz!Bz_6QcTnK1f@DO^z`X$^{GTrGtB#sK=jk0K4++lT=)WjjZo0rM)VWLqiYpI@i zwIz`xTjs>&7vM2Op0zjj_2wSy4&L!D@c@?zk(zO|{lEc00r5>QUl^Xol7*kEBmdu>HGnNy`50J7k zOsDyc+C6AL{&U{`kaaY9SVmX7;?tz^_zuR?z4I0h+DUv~x_lB9nHYPu&gr*Djuzru zvQpAQHuTG`m2-wfE{Vira^sFeSQU%XO%^CMu3RxX#9WCtzU5@@-%W&Mm_lg+!ir=P zQ}Hr4!O6!aQ1}-rfTo8m({*yAoa)=%5psDAsu zl-InrZprYBS%sV8igd%b=UVh23@=QZqrv8hFaPU9OLdD6`ux4!F#xhKExU+4kv`O- z78bQYRU&P88OC4K<$pzXkUIB_4O@byb7;$s8L6B<+O_n6(g(v;w|jt-DX`T1SL(U{>h?tL zSM>PM!-3{*a3f$j)Mq<)mjAV1lj(t%KLb3Gq`35-6A1Cl{VT8$aKyo_tw$%m`ZD=& z#iSi`Qg~&aWBSnf=HlKw-2k}rl2)&Ion3L#S+%!WY=-j?LlEa{kBROLF@6B|Dr((> z@^+Gu-u^oefE#6b1KcLc43@NRJkGPg4KmSAYE9+Rqr11H^u8H&YasA3`1Z`KV(Ct- zUi5#Y&X%4Td^yfG=Ei=Jq{-|g!t%9w^t@(HzI4~)ceVTEz-_mEy~4EKe7wT+KK8uv zy9sU;uy6m|_0xCze%pV)6D(jvvIEGwpN=1_qf6cSUi&hal^a|6zUQB`XLI4u2|RL4 zE0@oo?it-8YqrLk`yu!qyY44WlIMEFggz=snlK=0yEjG66O)FLL*}RC3*2tYQy~Ir zzlM$yhV0!lX=}v_FXo>dd^&NX+0<8^Os9NL7hWAaUA=jGxAd~~^@RCW<}kUO4zBVo z;6sH0Ejp8SHUI|umot>!u-QkX-)wLRsNPpv@KB+2SdR=gkn>J7qa}YCB&*mYr|Tqb zFKD!@&O%Jv3^j8EZQEM0cvr9vQ-Vt4wJ%5#00DIKj97?6=`u#=A>l&5uyyx+6bX8;NDp6e&n#`M$zPrCl$bgo4G72?en8G zej_b)rgwS1C$aTuzMq&&T#eHuCdFRxm*~Z1S@hvM13bANE-*+Oea3EGG~EPbd=qPq z(Lw8R$fcr!oB(%}=?z94Rh^Fuy%INkxxxD1bU&zBCl2}cFaOqx3~KB(!6%#rr1m6F z_OntqMn%O7$ffB+RqHNb%;T8jW%|6Gqyqt?!1GsW3adzJ^CfBbJr-#Pzp777X7^(k z)c9=dz6ttU&JbB=_QsW$$a#vXda9lO{tbJhKD`y5tNaDUAhhmM($S zipnr+jMe!~=t4w$k-mdR92~U@-Dok5I^VU^af7|mppKfHP1@w8gk!nsb2>hg^q-uY z_^JJ*DN~`JTAIdr7awP~jMxUw^u7z?0@O9}P*qT^1B&JTEAh6y3OnQ*21y7>U|J=? z_lTreBcz~Y^x)_)y)oqk#%7s+X!1xf)dm$^fA^q`prj?^*0Ih?rGALgvswhfOou!?%7y6+_?1m zaoi>#ju%@zW3P2F}QkK!QF2(WW@l&L3Tdsu+BX22rNLO8d52 zs3?y5B8;0#3pjlkhF*h2lD@xN18p<%5`aC*>ZN^Tj=Z_UiKH}04Fb5xE%$ZXCQ&t%RXfN&6v6t418Q;EcJ z=bA_ub4rrufm|=(l_V`96(yYZM zoDwG-#6N_CMM)%FE}|}X(Iv*f!2YUYGBBI8a}tX~#eY+-rSC*%m<$pRu`3c5!@?4d%Ex;0JVSI>3EnDYnnk;^a+1#Lkqdk)plm- z%GHqR;$?GTpnTG3b#TQVN>c?zEzuW$d<#juwU2i4jo7xN4BvrD#GBmNsXCFcQA$9Q z^GSiORosjV!)?_nxW<_(Yff}0SFFeiiFBWa1%2*K4tRRjmfeK3_0Qq$aq6GK>GgH5 zc9zo__;GWt%gIx*!EQLD{y|tSNVYb_fo3E+vv~qAl1FGP^@i+6JR0LLxE17n_(?>@ zSdv|2ysCH{Vrcbws7K&+)Qb%hq+jl9vob9@x))s`=~Y?YL6ugka`>q*MToi_nI-k? zFo6u&k2B!7H8fmdTso*tgana;QUI0q4AZvej4(|_=~L`|*{+s z>7D11WIG*pDJFJJW((Y!Si$cxDOt?_B=0@~N97|L$4ErAm1eNKagOFiKl(rSgE01u#r!tS!7|KWt zh+Ji%4@>&qktbx|1Tp;azvuN+F$rt#!>ft~SHFkUTX8fo=-Lu|FeZ*NtU)>HU1cM! zw8>Vpkw)UA;`O+Y!}kBo(}Ua#k=!et2ss<|7kj)%&#D?~X2m8+4K(k}^0f@c9~(}0 zb@IF_a1;!+Y@oz=mFjQHA|&WCWY7HON1SbGg{YeqadQvl;_D;b4@_k22CB3tR{eor|)VBbMGe zcASfgu76)B_bzKjrics1Tf^t z{L=RLDAmJ)l?eMoGF16t>eAVb^SAcq=6W=>6)-)_(TN#5o)kl&-`lpM_qX*WGx7+t zgB%={RGIEYdSm*YXfNJ`1B`Ssd?FH4)_sC3u@{Px@I~hcheI)lrFk_GC%rMJ4tyGE zuoh!MN!EdWCV^(AkY{GVInMKuB-TUnp|o=0)8Yu)YLxxM`G`NbsNGsh`Mbhr3ogTL zHe;dcS<`cfzAD$^QXJLkU61B&l2LXgGQ!`+%t;h?or-3asY8nJ0QD~gtizp+=v=$P z@F&l#b4=@T9-gIxmi`)7)q2CfC+lNCb)MUNB)d2Sz2_FC_mk>vK;A{>foC>!vw$K? z(TuKrV14*)eSQ54@@>n(lbfGFe}|6!VE*Zi7MyK$-GLk1H1q8X@7;N_()!G<5)KOl zpuls~d3F1TR78j3NLylLL5XHSwls{}1|_0;D8vnfXgGj6Y9FQ@4!XJQ1iHZLn^@`L zAK3F+%+u_d>8zT(-um(6;iV%}N(9cBLsg190jovHP3nT4Q4poy2NybuGhi>da#GRL z3als*R#AWjLCQptIXui5_OEddf+IPu*sjKE4z+9lfG9yj(kpCNrg}A0F66F+5~`H5 zafwrQtrGwraPOj;zEw5kaIUf^j~~DDwnrkuDN1@Y`V13MhH^c%e~$q{F~{&|8+Tm3 z5g|eh-`WW-K6DorEeHPxjG{Jd0&8=Nd8n^|*cQeOviA$Or`DdHR9&O+d~%hL#g~y+ z@u>Psy_HFk5lx4nZbmPJ2VR6#2T+Nc6Ef(gu15-RF~f)E&Sfvg1Lj2mR@p+s!FEruCxIHWoV;Q zUjM?rWu%TUILdk}amFwnkL?l>qaTCKo&El<6NLi@EkB-IZFZviL#KzFkBHCe0KIws zmmYoo8=q*59-{c&R3VVk{%8#o_a9BvVpMPn#tkB4JWw@vJ#)gIY5NHCgkeicyQ9Qs z$dxY-FwcLkviNPUl$~j7&bd>{(4NyVDc>8*Zg-C`zkRNFCk>Ala5ZkWzTWJO>?7E% zcmaIsLmkmsG-?a16O4uh$#Wg1IweuJiMw6EBgk2nDW;S;*&q($5hij}ow}wpF%3q0 z|COCVpAT);u^#_?vz5X0AhuY814HL0Xp-;RaRriUOhM$l5jK|JyCjYfMF;g}byR6gAPbpl z_*lOr;bWEI7A^*o!BhRTSqs3xg+T3j0^u;d>FP8EB|3qrhwprW6RRQ-L&l{G^STHk zYI$&(>pfJJR3x9^#}MV(OEi?S+`r))1zJ6f35W^qP7ko8zZO^_A$0n8$s4MGY`kmTe-waV~T zXe^Q;2+%%SWe@-63fLNgu{@BdovV%6kT~GrB>?997S=E^fFL2ANA?twe;3b~hv(%v zSL{bAf)3&E@EBl8?0^LrHI_K6L6L{84mIN?Wu_kM)!h#EEO?evb*FMk-?+z4$>L1@ zeFi&M9KiDsWP$-PnAdcR#&mEW)7lVT^2hT-1!{C+EHUIXAJ)88S)Uc(EQ= zHh4aJw5sw`g|>3C?B<-&6ty;w6iAg8){Z?|h@k6)ZcQTOW3$W*QlbW3Cy*pVuq*EF&=GW5>;AjlU&*|(Wo`o35JRm@dDd({F*=Ia_qf;}|M@fmZC|Q> zz@TM%IHu)%#~xXkRb#IDC>ZZ={xI0DS4{s z$l0jsADcn;-No}%rZxWUc9_^pb9Fs)2l>lMr~?;FkY@gH>uCLG7*D$oaP;iJ^KURX zn)=X59YhA0w@IO;uc^5aKcK~$Iy!ZsE{4x}CIy%q|8frFK;}%J4sQvbsGBC67j!b; z2w-mIOc2BkY>72YulO^gKtytfYgb{O@W*|=LFd$8tmgUN`4_E;a;)N@SUcFazVT+h$fwZg$1iijBkxritF7gpXwSRil(Oy>pcmh(h|S zQ#unnKb{Qgz{LlBVMSFsdQ6p~QONWmWjT|V%C2u_DX+p|fKdKB(I9Rz2TDconQ9aS z;nE>|EcJH6@)9x&6f6dJYJ+x$qY-QvQ0ew!?BO_+rS5{KRn|(zBZXXz;B2nx5vnB_ z!?7!Mxr2E^ca4FNsQJD4S4q_;6p#l}9h7|<1?F7*(cy36I#^qK;!?J!pb*VRLz0B1 ziFTFzvA9$j$#v_u3W;N?S7csRKb+=Wnl#_NZkEKvRz3`&QKL+MXc#u*nyRd=^wywCvXpqAE>*IND?j#X zyYBmAP%p?I(Eb6lbvLi}^-bMe zsHBU2L`O-Yl43|Yl|_|2CpM#GUXY5(i5APJ!x0|fbJB<;GkPREWT_~c?o!%Q3+ z&bI|5ptxvNhSf~z&v8W9K}b5Znprhk3mlbrIj4>nd*KZc6^uCZ<)|n8UcXhCEC?MFenz?Q=xbvhU zV|`B5<|(H%6FSJj4d$`L9N_^xX~$1pc`%cb`JQHsP}{*rhl8wJ+LGt830Bcs9heM^I$t7Dz%b4fOz$x1OoO&GomsDuDIdYgK&H?5W%#*V)B zw2S!0>CYShy)%p0kIH{4J0)COew1CF-|sIrip(eo$haXkkzO~QKGbV)x@61}0~*pO z7R<4!*I#`P3SRSvR*+9j0{0l$OQdj?z}+EM{4;}I1!jEv&HUqwdX6B26}X`hvDGea z%(j$2L?om7Wp{t7hGfhtBcPUJ{uW1{?jks8*PJYk2JsCEqRtN3A-@Nb`jGzny4SEX zbussD=$W>vyB=Ig9P_bd!rZ$UA6{WtULW{)_O`PxVfIqlN*YG-Q+9<{_beInl^EFq zso||Gu{`>Q5H0d!0x^y++;fCIx{;`cylK_#Al*XE5+5gQy|w3hl8jKuKX61oSl|zJ zo4n0b(CoCkHF~>cYRmr!OQh5iWbC$EfZ{l9-){LBe|?i=ize|-buDE|*;HL%ZjCvI z%!$`=aq@seu>=%Zp`RRkuppo2fvo>2pOSM8AaSUVKCcQJjG3d;u9DdvVz9^B-kXu} z(aJKtDEbvr=5QWw{3VBd^PaB2^C^{4`>=B#DYuel+vXdcv_OSkVsI?d)k4UZhn zHK?8}a=_p(%WXP0CG0U!0h3579fL{>g&BEGer8-oLkm;o?&n|+1sjyY!(Fu`-LHpN zhedS=$P_J+Q-wJ`GT^xb8sPz`Xje$O>H90h`D0>6Ic-NNENrdX9T z>`WiTqT{F-#$a22&#jT;vM%ID+6=ziw`C_jK!k5`^k)p3)?dfW(m-}3%WWy5-PeUU z3*i{%C*pe3_x++%6G^+2)hgp?+*xgi>S)!H90j_kO%p9HiT*{1(=!KSdO;F(dVoZx zG!R&Gs#5SJ<{~+ANsBQnO=xFwyN2o6@GcIQ+>7UX_AE~ov;^ZVAj%BwEwpeQ3#_YR zpPzNBZ5Kb2%4($ayp2uH^4R^ToVOmXBmH)ZN&99+Mj+=z9<#T!%{RXKhxxs}aJ~(P zxV#YB7>!Q^&r5}wqx9MGu0;1upOts%fsm)91)cC-XkeNKgDL&<(Wz*95xZXM=rCj- zPgpQ$fATLSaL*QV&!bUZIC5V$xnba6ZIlE()ufx2YOfW%3Ls-(!oLt` z9T=zTgOWy^aK~|IY?ty&V_G<7RS1Rh@#5T#) zy=6&%+Bi)vNMZwBa^8~V5T`Ow>*Tl$=XlEV?UR)=ql&a_D>GLcOGf;U9G>1As{(pH z1$UqV4UU8LF#VhX>R_J1WhtY}#k7YD!W>0JQSJo3?1D z@s(oUO{e|^&Z1RWmZvRm{T)p=v{J>u+)ULZEHnK-eToB;$E?)(!|diPoS{zzxzY7#as@r_z>GbP4WZ zFpcB)qaXnThWE#XyXT8ORYlxnryY4qKik9!9sD=p%7zsBX>w+Q$nrt*&z2%`mtDew zTE|d`g9v!TD#nvwRAg0{p#9nQD|CKqx}r7HLuYKDkqnNzkcMfrPiTR@XmdD9LfmC` z?hO+1qU3;bv?X$!)#h6M%gHUCbWu!Lgnf#HUwW%4Zw*LrE1A{eb($-2^UdZBN*^9Z zQj=2M_jgGt`g1zu1fo}O*@3paj#5M3udI6V48ZQpqGNx*5_552DHz-S^0b(+j5{dt zxpHpk?dV(2h^ZGIT{*OG#>|X9OtvS)7mmI-&le>D&+l#+ZKxx13eJ zM?c}t)aUvTNbve|dk~BPZ|~eZnatiUXrg*5;}0o`!0I{zDN9X~r8XkE*zxP!QL21Y zvpJfjRiK+e%HR=dD{}YYVq%y4Iz@ES^3Y4Rx`&*qT$b`_nUUi*+2kr;&!d*-?P7q6 zUa{5>-0gi9y^SvzHq4kAppa8)PBx&%dOjVKGzP6qXyf-=%H*6~|0Eg$YSfclX&%jF zRE#*U<3=vHUF-D+DLMY)U0a4?raz7BwqE?cHosbuF{y4=Mvzp48^3~=WI2A@c{FJg zM@AG;jt$2G^OPL1|0k~gY?$0jI#O&XqbN9{gjrBO??|~fJo-V8#jGhS>oH*s+H>@f z83ea;3 zj7lg>_WR=3Htm?*)Nm)vp@nILL)|(fX*Hzj@+n|#<6H1(D8&u+#XelTt;8AIN0ci2S~XUj(Y48v6Yjf~K77iW4VILSp7R`OMYF6bCcQj@nNS zmlN97%{Uk$N*YVXch9|Ys_^H6K-eyXu#t_TKKv1BjtaNeyLz9A=`^E9`{^m%O} z2z=*$CO*~+Ht_XZ(61V+`(bn7#z|3+9B3i1Y-uUP%!e@}0;`bEmWHss!+Y6lGT2MS!D4WrXE! z$Ya-}FE$8XjZm=dEMAQb$pXfSPCG+NE^&90FOjc|f1n>sjkQ9J5F*Ft zlS$GYV+E7cSpC4H(DAq{y9g;?7U6l=f?Xc_4lLD=pVPo4yP8Nz@?|9t1_ryPBj(AziE|2c+mI9MJ`{@BOgU|MlPZP+CZ|WLfurJlvXCW-XB1 z0y-VciKkj|Y=}NwE}`1GW`UZy#7fbkHe=jSr5cMd@=zL)Phsrf@wq&uHNmCkipU57 z{Ouw|Yy5>y|6J#0AjiIulVI8;=I-ORI8!kN)Th6r(A1hd)g#E{_+fLO#$>jxY^BsG zdy_bowotUTxMt$jApM8TCXP39EU&1J5Ggb3XCXYaCA#Gfq5vc|M`GyH*g;kDC z^hYSCN9%)q;R|}8wUucIF<@q8T25`tN^|{u54Z^9Um|}5*DY!3)P-E~H zAUjlC&b%oCS!Fx@)V0LsUbEg_F~Sq2#T9@bF7XVfiW}iX&(_wA{EJ_Gy868Qe!phI z(VtG8blbQC36%?1Pnzp`~8j+S+Q}f z^dko&=@mmQ*J#A5q$p2}!NcftII(2aK|RdRak=;h{2@f>;IEKOQ6p?y5xB-F?>mId z-yqMwZy`r)R+9<&Ov{&`OsxEVw;3)&6)PFZX18`wmj#A^-cToIi?7I+g~lju9D-S5d6}n^1y6?N z(`ig@(Y-KR<-5t6Eb~%jfy7;>e<_@+F1HoETB}~CZ~Or`qGVxA?GvU?3DJia4Xzk= zh>5X<2q%*J$HOzfbakBRRGonIR{oZpr`PQMongL4DmF%k?2sG}ST?zm#M6f3xIZoC zjMb(bnGeJ8RAhh=8m8EAFUI`)#ngTZ+<-hLxJQg-0zqF~RiqmIOq=vJnFV*+Re(2o z#Ts~#lCsc9u1l#_YkZ|i>B3?yg2L~fd(A900SN6Tgf7rTB$4EJ{iY7=CqEACH}QM> zh&*!Z+S1W814zGn#Vf`!^-QI#lTviI6n?;ok`flR6B5yn5<7CJe1=%wv{ty^vJo2) z({2|JlGNJntr8M8ilkRM)kI(-n^QBuVLuH#F?2TTu@zvVfrVB6B$#qhlt~V}@nDC5n1c zLApP*nO*hryd!#;$n%sZwz@vYWl3#Qlw~F4N)tD;o>;C^f@|1HSZ)%gM0+f_{FAM#)7L=%6urnAp z9}`+pM0{K|Ttv@$-*-^~EiA2XIj0G&8Taop(m*wIYvk@zp#HGRskk_auZzxJ70IRC z)+Wu4%B}ISnzYZ+dWYyd1V`?SWB@XYgv-|W2Z#xBH5&@4u&c}W!xjl8g~aG<=PmOV zIxFPT!R}0ku>zE<`?cksax^mTYl9v7)}5RF-r{r`SnvHLp}(SX(0)^8Vtosvkh zpI6>trdZqbO#TO9{--~J`{3UO_X&tQ-ao_r27LH7V*el?6{rUy!Tk&n7)Ozv&^s<-T=T>mF;ISA=C+|j#mZ8D6Bn%IGRrAql3DpE{v3Kna z6-#r8dN9XS<;q7J5$TJYa!gqj*(!2zN|nLH8H(U;+v|(s%!{Et5sXAiiY$hn{z#s_ z^zHafXCHnR8R>_3hQ9&*%ZLVEOLVfl#rdz&4D!JnnDiydWUWy;cUe6&Aogx3$Q!&| zLWaizX_R>fER1YkN_0aclD(+5p_|y`&6;&jK@Y(}@+0=Qr&q z%o`X=k32}mmy$naTq2L+mSV{r-==HuH&dyValBhp6gkp|&YaLgEL5;I+ALvPb?e*M zob9)y7(h4`ox_#3Fsr{B0+J5uW%SlYJvV5 zp~hJW>h?9ZTY_Om`0colNPD!OoalAjl0_OcJilMfx))Oos9R=1o$sRm5xR>dMEb6d z*f3=8qw=@!(j?8${-w=BIFR{=ns>gQIOTH0C3WgIeMt}@{6P*3D9~90EM1_+%G^rk zHf%GPI5DEaKM2bl?Q2nsY39%+f!=N%<%>w0d{-Vv`Z^RGPx(sg3QOiwM?Q;5-dRm@%okCSz8*&P*n&;r21|FjxMDEB|n_G&Y4QG+z{U?3Us=+z6AybUTuNhuGy4hX=wpX4&Z+;&`f zQtW`7hh&ai4l4H>O=r92mO)>wUe}%4RHL`qEQf3MF~L^+ zD1GxnHGcK)^_W)1BZOO4A@3M)c65N`8*wRifM}t(?ZqLc5KEkVP!|Hvi?8&g)^}=u zEb0Y#wM+VJakcPbR{cQ9ibG(LrQNl;tBGbc5~8&6m>Z~JPS?J$W{lIM4Pk0VeNrqN zOxO~Y^lQtdXC|}rV0y%Ww(+1tVn{yz>l(jQ1ZO}CA4+stROnUOSQM6Ph&Y5CQ^a9B z3>0ShYOM)96q$kO<#=S3{sS#y?doqPI!-!zLM-6qA)^k5}wLe=4kde8I9hILfYd8HYr`TWh%DI)` zxxx_tQ>*U@@k#n^vM>L~Fd?J@b;bUinEITk#;K_|rb>gsGCUjmO#pJ^JUBIwcNr-M z`H~ckvWChG+AraOo>!t$J|p2@DFao_(Oeq-xljCw4HJE7G6mTkUfY6Y(d$zAh0W_7 zM{CXj7Ej_xUA9k5zh$r)tK~N_xMN97MXwc~ZCs+lBv@#JIS{a~A>=cV$k&Cjx%c$< zqWM<-aH^Z?X_ED4;3A_jm8Q#WkJif<=`m5Cw%AiMS;aDx&GF%$K;Zskmw3dbXMJ&4 zc3+b#yf8giCh3){RmI&nf)eJ_#FUH__+OAzRcE}IAMuc^Z9kCdgAzxzSxX_u9Ib5F zDKrquA`ib|7iLHCJNMTL77+~;*;x{+JeHJdYZgIttxDGz`ro5OXtLRH&3mZCl9Fwr zsn5w8%78|dx^ZCQSe)Asm%XCBJK#{}PLEjwg4TN|PIMQ0>QxSOcn4}2t+MQ~YHY>v zP7e{%@_#EOhX~akN-XXB?jW6b1q`&5a0GY9;v}0tmCBe5AZ_r3Qgo6?OUsNWQw)#b zE=*44=9PJy&`z8bB~$TKypwDoh+N1&LEasK#UM(0gw$QghK_M+=d|S%8C@;|kj5$m zo+|(3`v^ir+ zl@g+XTDNMY6hEm-R{R<<$f4)OJ|35z`9?gfB4c9{&2DKZFNfIVtS=Qwr7Nm%01lO3 zx9@?*TpncCp1c$tOMa@SqGvR6`(-;?$@zxsQ}^dw9%6y*K~78hqDgDRpJI{<8s<5w ziA9BIcM|7s1jky*UqcVxWKXWvH`i=}u$aXJ!%6M)CjsaO71SzYm|9&;{mcKFfmrNJ z=Af!(=xEr^Eq&Qos*;T18?_Wb#O~?Dr~vzk0*qBj6^?+Q@(x$aY1wIOdn2u=W6A2a zUXCXU(DWYC$EblH6YEb4M!1iax4DDNW@uCB0Y}&pcsM*!J8RX z+QbD)uSBD4#4b={E%SINal7|kM7T7csK4TVn>Ge^WoS`m$@3T%tGhk zA7$~b{h2K|UqIb<0tan41&22RG>&;qoQafk9Hvff~P_uMrX=%=?e z7c2{cEOj&QUtJ_Cw=Fx}$<$xkT5g`Ko83{~5nm}G67_lo5?b>+J;2i{Bs7Z~c3^+NqJFrBxa7YQ1S={vKV(mRcjc!`O<%7u(n2V=B%yd)8C+K-!;lTvrjI!D%9nFVw+pu?{JMy^rQ%W$F zKzS7P5IpG^XswHAHR-m6={tH>v!gX(#Apv?MYCmUlWoY$@HM)!?-;)>^Q}5gHWEyI z+MT;M4#^kztf`;nCi@rGTw&DrmVelOG6YvN99}w9<0bA>NmH0S8*}CS}|LxJY8p(a0RfIp^2Z>OfC{Pr?ee&?Pg#pbm!y6sk-KmME5ChsL&;tsOIk zZ{!hL$g3BP`zXm^)#PUd)nq$O4d5-R7z^$$}W6x^IjH*bbm`id0w?J(xkuj zvk+cM4_Yw1u1tb>df&k0vEBSzZf7kZGmgZ&`GRVlS%9YKPziv z5ZNW*jDqT-i5izbV6|RUkGSNX0uwLaLoMKrdfDu(71&hi;DGn+;S_hD25m|;$ru=Z z0?ZBt{9LKPH3eb9CM0~cCABHP-mh5s*ODg5xwo_~B_`^Ep;xG;jp!^#DT^?`7~rxG zmAxW4%~HeD{3od(`1ypFnp`hJF}4rXXg#$@TV3j1+rQmNH*+LRKnvR$0`N!uqi8gU z;>e0C-gBcz+ecWL;evoEijTP3RG$@!J^9K7uh?h=1AlGvq!qktoQ@K z+J`eI(9x`G>)*une4JBRb|cSEwbGxl5VZUbUQ5_gzY8%6iaCh_PjA?U={_{lB!u)( zrz9WysKxrMF{Y>}7Gr)u`eAIgd-qP5lb*e_5}mtk#d1-Vpi1gStX9W8MMct!6LglT zIi^L4^Nc^~AT-HyHXv-lp3?P-DFFcSIH#auHLPWjq{Kqi@8Bjivkd-+qk7_*_=~rz zTRuWh$Z@QRf`_DZLKmVFB)>PTfMC0$=WF+V3$N{xLisH@$C?a3ZDoN)fIx1!i%apj zz-gRVg)OzbE#8J}j#${vH3k+zdx~qk6XbPThQ)Z2jDHbzL<5X9i_9P(Vigf*YFXi` z&)+7lS&%7E#(RZ33^LTV4L=JyekAnzbK*1zDZ&bv^a+-V>UQMPev2o0005lVcx za!QPTC>{Jjtcom}G06lrYk|WJzwE9mcw!X8b&L zw4#48+ygKgMv@R60WxA*x~OJBcq!uyl81I&Fchx*ZEC`U_#U5OgdLVT?w-FV3mrLZ z=zOo11E_1YxfsBWQwjPON#_rqYt*SeTRfDadY4HR*T3xqTt4B?uG61m-?pD10+V_U z3ox}Gje%0Q$h~jZAE&;^dS_mNFOdcRgL5(f$bGQcKPqq=$$*JWzK-lytUA;g z>Eqv{~AamB}n?IkRd*X`_gEVRwjjSGZlprF561L8@)`N@}nq$DX6CJkSv0TTk zToB7&yp_XoYJ{3nd+w5lzYUor9{4aUo#X*QM;9c;M>9BsKn`S7$&y$gSd!2Uf6nKJG$#uOrC5W!7&x!GHKIffuv z7GCYr#@fQtXneRdj{?uM)3_<|1IcU9Jv47wDLTLkcJ6u=0=-&Y_Xip?mjvV9rQ9F( zKLpvEY(eG!!$e#MnNI4r1PCC8jLK1RLT*43L3dQC(vQ*7txPl+?@d~&Rwt6tHKT=z zB3O7^!h1u&5ZnWu!92o2^R_6QR?n@Twuf#nt~X1F;3{P)k1lxpahcxzImoSLox3~f zf|sW+IY#FRYE9Em%AYCG{17cYXw*Y}-7-&_f!Dw#AQ112pZ#RGAGGZ&AcW-t_5y&s zexvLb!PnWfUknFZc}O?bp@+Mh(P4Zl)<<;@F!ADhN)ppW6d2~vwjZ%2sF6xUsx2~A zE>YO%(r=-o1|lMzq>nKXk*4b0MmYkDYfKg)DwSO_j6-x1=AL}!jlK*1xupp@V8NhI z?5+Ly)DT^H^9E}~tb$zX$n7+Pc{F(IE#_y$mdLu;b{(MW) zmgu9DEAytYlmxUxRw?#`#oZtUWN12ON@hA%4pg}Um-IF$ogX&4W~BnErby%sEpl+# z-{L{#d8!@1n-aMNK6HdW{O5A4VWZz@!anl=u^RN3&#cV!iZa+^l@^-Y{E#sjkwbP0 zKBI&~Y}SJcUAN6h8CvB+B6v~^{f=Hp#03m-;Cg+(Q^G3tOh28acNlIlnZ^038}m(G z^O7itCrZ-CX_P@CWSAc~QTnsnmMb^UC#+XHiP$w6wTqb@JiUkEi1F@Kycf{eYZln# zS=na*gI|UjEql`AS1;}REoV(8QEx~vApFoChGRh2)gh$)lnp3a?zVmnChS~x_Ld*_VP60GZ^sq~ K2o&i5&wl{L!5-NF literal 0 HcmV?d00001 diff --git a/src/sctools/test/test_bam.py b/src/sctools/test/test_bam.py index 84b1502..61f2b15 100644 --- a/src/sctools/test/test_bam.py +++ b/src/sctools/test/test_bam.py @@ -1,3 +1,4 @@ +from copy import deepcopy import glob import os @@ -9,6 +10,43 @@ data_dir = os.path.split(__file__)[0] + '/data/' +TAG_KEYS = ['FOO', 'BAR', 'BAZ'] + +SORTED_VALUES = [(['A', 'A', 'A'], 'A'), + (['A', 'A', 'A'], 'B'), + (['A', 'A', 'B'], 'A'), + (['A', 'A', 'B'], 'B'), + (['A', 'B', 'A'], 'A'), + (['A', 'B', 'A'], 'B'), + (['A', 'B', 'B'], 'A'), + (['A', 'B', 'B'], 'B'), + (['B', 'A', 'A'], 'A'), + (['B', 'A', 'A'], 'B'), + (['B', 'A', 'B'], 'A'), + (['B', 'A', 'B'], 'B'), + (['B', 'B', 'A'], 'A'), + (['B', 'B', 'A'], 'B'), + (['B', 'B', 'B'], 'A'), + (['B', 'B', 'B'], 'B')] + +UNSORTED_VALUES = [(['B', 'B', 'A'], 'A'), + (['B', 'B', 'B'], 'A'), + (['B', 'A', 'A'], 'A'), + (['B', 'B', 'A'], 'B'), + (['B', 'B', 'B'], 'B'), + (['A', 'A', 'B'], 'A'), + (['A', 'B', 'A'], 'B'), + (['A', 'B', 'B'], 'A'), + (['A', 'A', 'A'], 'B'), + (['B', 'A', 'B'], 'A'), + (['A', 'B', 'A'], 'A'), + (['B', 'A', 'A'], 'B'), + (['A', 'A', 'A'], 'A'), + (['A', 'A', 'B'], 'B'), + (['A', 'B', 'B'], 'B'), + (['B', 'A', 'B'], 'B')] + + # TEST SUBSETALIGNMENTS @pytest.fixture(scope='module', params=[data_dir + 'test.sam', data_dir + 'test.bam']) @@ -157,3 +195,171 @@ def test_split_succeeds_with_raise_missing_false_and_no_cr_barcode_passed(tagged os.remove(tagged_bam) # clean up for f in glob.glob('test_output_*'): os.remove(f) + + +# TEST SORTING + +def test_tag_sortable_records_compare_correctly(): + records = make_records_from_values(TAG_KEYS, SORTED_VALUES) + num_records = len(SORTED_VALUES) + for i in range(num_records): + for j in range(num_records): + if i < j: + assert records[i] < records[j] + elif i == j: + assert records[i] == records[j] + else: + assert records[i] > records[j] + + +def test_tag_sortable_records_raises_error_on_different_tag_lists(): + r1 = bam.TagSortableRecord(['FOO', 'BAR'], ['A', 'A'], 'A') + r2 = bam.TagSortableRecord(['BAR', 'BAZ'], ['A', 'A'], 'A') + with pytest.raises(ValueError): + r1 == r2 + + +def test_tag_sortable_records_str(): + record = bam.TagSortableRecord(TAG_KEYS, SORTED_VALUES[0][0], SORTED_VALUES[0][1]) + s = record.__str__() + assert 'TagSortableRecord' in s + assert "['FOO', 'BAR', 'BAZ']" in s + + +def test_verify_sort_on_unsorted_records_raises_error(): + records = make_records_from_values(TAG_KEYS, UNSORTED_VALUES) + with pytest.raises(bam.SortError) as e: + bam.verify_sort(records, TAG_KEYS) + assert 'are not in correct order' in str(e) + + +def test_verify_sort_raises_no_error_on_sorted_records(): + records = make_records_from_values(TAG_KEYS, SORTED_VALUES) + bam.verify_sort(records, TAG_KEYS) + + +def test_sort_by_tags_and_queryname_sorts_correctly_from_file(): + tag_keys = ['UB', 'CB', 'GE'] + with pysam.AlignmentFile(data_dir + 'unsorted.bam', 'rb') as f: + records = f.fetch(until_eof=True) + sorted_records = bam.sort_by_tags_and_queryname(records, tag_keys) + tag_sortable_records = (bam.TagSortableRecord.from_aligned_segment(r, tag_keys) for r in sorted_records) + bam.verify_sort(tag_sortable_records, tag_keys) + + +def test_sort_by_tags_and_queryname_sorts_correctly_from_file_no_tag_keys(): + tag_keys = [] + with pysam.AlignmentFile(data_dir + 'unsorted.bam', 'rb') as f: + records = f.fetch(until_eof=True) + sorted_records = bam.sort_by_tags_and_queryname(records, tag_keys) + tag_sortable_records = (bam.TagSortableRecord.from_aligned_segment(r, tag_keys) for r in sorted_records) + bam.verify_sort(tag_sortable_records, tag_keys) + + +def test_tag_sortable_records_sort_correctly(): + tag_keys = TAG_KEYS + records = make_records_from_values(tag_keys, deepcopy(UNSORTED_VALUES)) + sorted_records = sorted(records) + bam.verify_sort(sorted_records, tag_keys) + + +def test_tag_sortable_records_sort_correctly_when_already_sorted(): + # This is to a bit paranoid, but just make sure sorted stays correct if already sorted + tag_keys = TAG_KEYS + records = make_records_from_values(tag_keys, deepcopy(SORTED_VALUES)) + sorted_records = sorted(records) + bam.verify_sort(sorted_records, tag_keys) + + +def test_sort_by_tags_and_queryname_sorts_correctly_no_tag_keys(): + tag_keys = [] + records = make_records_from_values(tag_keys, deepcopy(UNSORTED_VALUES)) + sorted_records = sorted(records) + bam.verify_sort(sorted_records, tag_keys) + + +def test_tag_sortable_record_missing_tag_value_is_empty_string(): + tags = ['_NOT_REAL_TAG_'] + with pysam.AlignmentFile(data_dir + 'unsorted.bam', 'rb') as f: + records = f.fetch(until_eof=True) + first_record = next(iter(records)) + sortable_record = bam.TagSortableRecord.from_aligned_segment(first_record, tags) + assert sortable_record.tag_values[0] == '' + + +def test_tag_sortable_record_lt_is_false_for_equal_records(): + r1 = bam.TagSortableRecord(tag_keys=TAG_KEYS, tag_values=['A', 'A', 'A'], query_name='A') + r2 = bam.TagSortableRecord(tag_keys=TAG_KEYS, tag_values=['A', 'A', 'A'], query_name='A') + assert not r1 < r2 + + +def test_tag_sortable_record_lt_is_true_for_smaller_query_name(): + r1 = bam.TagSortableRecord(tag_keys=TAG_KEYS, tag_values=['A', 'A', 'A'], query_name='A') + r2 = bam.TagSortableRecord(tag_keys=TAG_KEYS, tag_values=['A', 'A', 'A'], query_name='B') + assert r1 < r2 + + +def test_tag_sortable_record_lt_is_true_for_smaller_tag(): + r1 = bam.TagSortableRecord(tag_keys=TAG_KEYS, tag_values=['A', 'A', 'A'], query_name='A') + r2 = bam.TagSortableRecord(tag_keys=TAG_KEYS, tag_values=['A', 'A', 'B'], query_name='A') + assert r1 < r2 + r1 = bam.TagSortableRecord(tag_keys=TAG_KEYS, tag_values=['A', 'A', 'A'], query_name='A') + r2 = bam.TagSortableRecord(tag_keys=TAG_KEYS, tag_values=['A', 'B', 'A'], query_name='A') + assert r1 < r2 + r1 = bam.TagSortableRecord(tag_keys=TAG_KEYS, tag_values=['A', 'A', 'A'], query_name='A') + r2 = bam.TagSortableRecord(tag_keys=TAG_KEYS, tag_values=['B', 'A', 'A'], query_name='A') + assert r1 < r2 + + +def test_tag_sortable_record_lt_is_true_for_smaller_tag_regardless_of_query_name(): + r1 = bam.TagSortableRecord(tag_keys=TAG_KEYS, tag_values=['A', 'A', 'A'], query_name='B') + r2 = bam.TagSortableRecord(tag_keys=TAG_KEYS, tag_values=['A', 'A', 'B'], query_name='A') + assert r1 < r2 + + +def test_tag_sortable_record_lt_empty_query_name_is_smaller(): + r1 = bam.TagSortableRecord(tag_keys=TAG_KEYS, tag_values=['A', 'A', 'A'], query_name='') + r2 = bam.TagSortableRecord(tag_keys=TAG_KEYS, tag_values=['A', 'A', 'A'], query_name='A') + assert r1 < r2 + + +def test_tag_sortable_record_lt_empty_tag_is_smaller(): + r1 = bam.TagSortableRecord(tag_keys=TAG_KEYS, tag_values=['A', 'A', ''], query_name='A') + r2 = bam.TagSortableRecord(tag_keys=TAG_KEYS, tag_values=['A', 'A', 'A'], query_name='A') + assert r1 < r2 + r1 = bam.TagSortableRecord(tag_keys=TAG_KEYS, tag_values=['A', '', 'A'], query_name='A') + r2 = bam.TagSortableRecord(tag_keys=TAG_KEYS, tag_values=['A', 'A', 'A'], query_name='A') + assert r1 < r2 + r1 = bam.TagSortableRecord(tag_keys=TAG_KEYS, tag_values=['', 'A', 'A'], query_name='A') + r2 = bam.TagSortableRecord(tag_keys=TAG_KEYS, tag_values=['A', 'A', 'A'], query_name='A') + assert r1 < r2 + + +def test_tag_sortable_record_eq_is_true_for_identical_records(): + r1 = bam.TagSortableRecord(tag_keys=TAG_KEYS, tag_values=['A', 'A', 'A'], query_name='A') + r2 = bam.TagSortableRecord(tag_keys=TAG_KEYS, tag_values=['A', 'A', 'A'], query_name='A') + assert r1 == r2 + + +def test_tag_sortable_record_eq_is_false_when_any_difference_exists(): + r1 = bam.TagSortableRecord(tag_keys=TAG_KEYS, tag_values=['A', 'A', 'A'], query_name='A') + r2 = bam.TagSortableRecord(tag_keys=TAG_KEYS, tag_values=['A', 'A', 'A'], query_name='B') + assert not r1 == r2 + r1 = bam.TagSortableRecord(tag_keys=TAG_KEYS, tag_values=['A', 'A', 'A'], query_name='A') + r2 = bam.TagSortableRecord(tag_keys=TAG_KEYS, tag_values=['A', 'A', 'B'], query_name='A') + assert not r1 == r2 + r1 = bam.TagSortableRecord(tag_keys=TAG_KEYS, tag_values=['A', 'A', 'A'], query_name='A') + r2 = bam.TagSortableRecord(tag_keys=TAG_KEYS, tag_values=['A', 'B', 'A'], query_name='A') + assert not r1 == r2 + r1 = bam.TagSortableRecord(tag_keys=TAG_KEYS, tag_values=['A', 'A', 'A'], query_name='A') + r2 = bam.TagSortableRecord(tag_keys=TAG_KEYS, tag_values=['B', 'A', 'A'], query_name='A') + assert not r1 == r2 + + +def make_records_from_values(tag_keys, tags_and_query_name): + records = [] + for i in range(len(tags_and_query_name)): + r = bam.TagSortableRecord(tag_keys=tag_keys, tag_values=tags_and_query_name[i][0], query_name=tags_and_query_name[i][1]) + records.append(r) + return records + diff --git a/src/sctools/test/test_count.py b/src/sctools/test/test_count.py index 05f5ccc..ce58c06 100644 --- a/src/sctools/test/test_count.py +++ b/src/sctools/test/test_count.py @@ -53,7 +53,7 @@ import operator import os import tempfile -from typing import Optional, List, Set, Tuple, Dict, Generator +from typing import Callable, Optional, List, Set, Tuple, Dict, Generator import numpy as np import pysam @@ -100,6 +100,34 @@ def __repr__(self): f'{consts.GENE_NAME_TAG_KEY}: {self.gene_name}' +class CellMoleculeGeneQueryNameSortOrder(bam.AlignmentSortOrder): + """Hierarchical alignment record sort order (cell barcode >= molecule barcode >= gene name >= query name).""" + def __init__( + self, + cell_barcode_tag_key: str = consts.CELL_BARCODE_TAG_KEY, + molecule_barcode_tag_key: str = consts.MOLECULE_BARCODE_TAG_KEY, + gene_name_tag_key: str = consts.GENE_NAME_TAG_KEY) -> None: + assert cell_barcode_tag_key, "Cell barcode tag key can not be None" + assert molecule_barcode_tag_key, "Molecule barcode tag key can not be None" + assert gene_name_tag_key, "Gene name tag key can not be None" + self.cell_barcode_tag_key = cell_barcode_tag_key + self.molecule_barcode_tag_key = molecule_barcode_tag_key + self.gene_name_tag_key = gene_name_tag_key + + def _get_sort_key(self, alignment: pysam.AlignedSegment) -> Tuple[str, str, str, str]: + return (bam.get_tag_or_default(alignment, self.cell_barcode_tag_key, default='N'), + bam.get_tag_or_default(alignment, self.molecule_barcode_tag_key, default='N'), + bam.get_tag_or_default(alignment, self.gene_name_tag_key, default='N'), + alignment.query_name) + + @property + def key_generator(self) -> Callable[[pysam.AlignedSegment], Tuple[str, str, str, str]]: + return self._get_sort_key + + def __repr__(self) -> str: + return 'hierarchical__cell_molecule_gene_query_name' + + class SyntheticTaggedBAMGenerator: """This class generates a synthetic count matrix and an accompanying synthetic tagged BAM file as described in the preamble documentation block. @@ -162,7 +190,7 @@ def __init__( def generate_synthetic_bam_and_counts_matrix( self, output_path: str, num_duplicates: int, num_missing_some_tags: int, num_multiple_gene_alignments: int, max_gene_hits_per_multiple_gene_alignments: int, - alignment_sort_order: bam.AlignmentSortOrder = bam.CellMoleculeGeneQueryNameSortOrder()): + alignment_sort_order: bam.AlignmentSortOrder = CellMoleculeGeneQueryNameSortOrder()): """Generates synthetic count matrix and BAM file and writes them to disk. Parameters @@ -565,7 +593,7 @@ def _get_sorted_count_matrix(count_matrix: np.ndarray, row_index: np.ndarray, co @pytest.mark.parametrize( 'alignment_sort_order', - [bam.QueryNameSortOrder(), bam.CellMoleculeGeneQueryNameSortOrder()], + [bam.QueryNameSortOrder(), CellMoleculeGeneQueryNameSortOrder()], ids=['query_name_sort_order', 'cell_molecule_gene_query_name_sort_order']) def test_count_matrix_from_bam(alignment_sort_order: bam.AlignmentSortOrder, gene_name_to_index): # instantiate a test data generator diff --git a/src/sctools/test/test_entrypoints.py b/src/sctools/test/test_entrypoints.py index e4b82a4..ca49f9b 100644 --- a/src/sctools/test/test_entrypoints.py +++ b/src/sctools/test/test_entrypoints.py @@ -4,9 +4,10 @@ import numpy as np import pysam +import pytest import scipy.sparse as sp -from sctools import platform, count, consts +from sctools import bam, platform, count, consts data_dir = os.path.split(__file__)[0] + '/data/' @@ -81,6 +82,92 @@ def test_split_bam(): os.remove(f) +def test_tag_sort_bam(): + args = [ + '-i', data_dir + 'unsorted.bam', + '-o', 'test_sorted.bam', + '-t', + consts.CELL_BARCODE_TAG_KEY, + consts.GENE_NAME_TAG_KEY, + consts.MOLECULE_BARCODE_TAG_KEY] + + return_call = platform.GenericPlatform.tag_sort_bam(args) + assert return_call == 0 + + tag_keys = [consts.CELL_BARCODE_TAG_KEY, consts.GENE_NAME_TAG_KEY, consts.MOLECULE_BARCODE_TAG_KEY] + with pysam.AlignmentFile('test_sorted.bam', 'rb') as f: + segments = f.fetch(until_eof=True) + tag_sortable_records = (bam.TagSortableRecord.from_aligned_segment(s, tag_keys) for s in segments) + bam.verify_sort(tag_sortable_records, tag_keys) + + for f in glob.glob('test_sorted*'): + os.remove(f) + + +def test_tag_sort_bam_dash_t_specified_multiple_times(): + args = [ + '-i', data_dir + 'unsorted.bam', + '-o', 'test_sorted.bam', + '-t', consts.CELL_BARCODE_TAG_KEY, + '-t', consts.GENE_NAME_TAG_KEY, + '-t', consts.MOLECULE_BARCODE_TAG_KEY] + + return_call = platform.GenericPlatform.tag_sort_bam(args) + assert return_call == 0 + + tag_keys = [consts.CELL_BARCODE_TAG_KEY, consts.GENE_NAME_TAG_KEY, consts.MOLECULE_BARCODE_TAG_KEY] + with pysam.AlignmentFile('test_sorted.bam', 'rb') as f: + segments = f.fetch(until_eof=True) + tag_sortable_record_generator = (bam.TagSortableRecord.from_aligned_segment(s, tag_keys) for s in segments) + bam.verify_sort(tag_sortable_record_generator, tag_keys) + + for f in glob.glob('test_sorted*'): + os.remove(f) + + +def test_tag_sort_bam_no_tags(): + args = [ + '-i', data_dir + 'unsorted.bam', + '-o', 'test_sorted.bam'] + + return_call = platform.GenericPlatform.tag_sort_bam(args) + assert return_call == 0 + + tag_keys = [] + with pysam.AlignmentFile('test_sorted.bam', 'rb') as f: + segments = f.fetch(until_eof=True) + tag_sortable_records = (bam.TagSortableRecord.from_aligned_segment(s, tag_keys) for s in segments) + bam.verify_sort(tag_sortable_records, tag_keys) + + for f in glob.glob('test_sorted*'): + os.remove(f) + + +def test_verify_bam_sort(): + args = [ + '-i', data_dir + 'cell-gene-umi-queryname-sorted.bam', + '-t', + consts.CELL_BARCODE_TAG_KEY, + consts.GENE_NAME_TAG_KEY, + consts.MOLECULE_BARCODE_TAG_KEY] + + return_call = platform.GenericPlatform.verify_bam_sort(args) + assert return_call == 0 + + +def test_verify_bam_sort_raises_error_on_unsorted(): + args = [ + '-i', data_dir + 'unsorted.bam', + '-t', + consts.CELL_BARCODE_TAG_KEY, + consts.GENE_NAME_TAG_KEY, + consts.MOLECULE_BARCODE_TAG_KEY] + + with pytest.raises(bam.SortError) as e: + platform.GenericPlatform.verify_bam_sort(args) + assert 'are not in correct order' in str(e) + + def test_count_merge(): tmp = tempfile.mkdtemp()